GCC Code Coverage Report

Directory: ./
File: openvdb/openvdb/tools/TopologyToLevelSet.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 32 54 59.3%
Functions: 1 3 33.3%
Branches: 1 32 3.1%
Line Branch Exec Source
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
4 /// @file TopologyToLevelSet.h
5 ///
6 /// @brief This tool generates a narrow-band signed distance field / level set
7 /// from the interface between active and inactive voxels in a vdb grid.
8 ///
9 /// @par Example:
10 /// Combine with @c tools::PointsToVolume for fast point cloud to level set conversion.
11
12 #ifndef OPENVDB_TOOLS_TOPOLOGY_TO_LEVELSET_HAS_BEEN_INCLUDED
13 #define OPENVDB_TOOLS_TOPOLOGY_TO_LEVELSET_HAS_BEEN_INCLUDED
14
15 #include "LevelSetFilter.h"
16 #include "Morphology.h" // for erodeActiveValues and dilateActiveValues
17
18 #include <openvdb/Grid.h>
19 #include <openvdb/Types.h>
20 #include <openvdb/math/FiniteDifference.h> // for math::BiasedGradientScheme
21 #include <openvdb/util/NullInterrupter.h>
22 #include <openvdb/openvdb.h>
23 #include <openvdb/points/PointDataGrid.h>
24
25 #include <tbb/task_group.h>
26
27 #include <algorithm> // for std::min(), std::max()
28 #include <vector>
29
30
31 namespace openvdb {
32 OPENVDB_USE_VERSION_NAMESPACE
33 namespace OPENVDB_VERSION_NAME {
34 namespace tools {
35
36
37 /// @brief Compute the narrow-band signed distance to the interface between
38 /// active and inactive voxels in the input grid.
39 ///
40 /// @return A shared pointer to a new sdf / level set grid of type @c float
41 ///
42 /// @param grid Input grid of arbitrary type whose active voxels are used
43 /// in constructing the level set.
44 /// @param halfWidth Half the width of the narrow band in voxel units.
45 /// @param closingSteps Number of morphological closing steps used to fill gaps
46 /// in the active voxel region.
47 /// @param dilation Number of voxels to expand the active voxel region.
48 /// @param smoothingSteps Number of smoothing interations.
49 template<typename GridT>
50 typename GridT::template ValueConverter<float>::Type::Ptr
51 topologyToLevelSet(const GridT& grid,
52 int halfWidth = 3,
53 int closingSteps = 1,
54 int dilation = 0,
55 int smoothingSteps = 0);
56
57
58 /// @brief Compute the narrow-band signed distance to the interface between
59 /// active and inactive voxels in the input grid.
60 ///
61 /// @return A shared pointer to a new sdf / level set grid of type @c float
62 ///
63 /// @param grid Input grid of arbitrary type whose active voxels are used
64 /// in constructing the level set.
65 /// @param halfWidth Half the width of the narrow band in voxel units.
66 /// @param closingSteps Number of morphological closing steps used to fill gaps
67 /// in the active voxel region.
68 /// @param dilation Number of voxels to expand the active voxel region.
69 /// @param smoothingSteps Number of smoothing interations.
70 /// @param interrupt Optional object adhering to the util::NullInterrupter interface.
71 template<typename GridT, typename InterrupterT>
72 typename GridT::template ValueConverter<float>::Type::Ptr
73 topologyToLevelSet(const GridT& grid,
74 int halfWidth = 3,
75 int closingSteps = 1,
76 int dilation = 0,
77 int smoothingSteps = 0,
78 InterrupterT* interrupt = nullptr);
79
80
81 ////////////////////////////////////////
82
83 /// @cond OPENVDB_DOCS_INTERNAL
84
85 namespace ttls_internal {
86
87
88 template<typename TreeType>
89 struct OffsetAndMinComp
90 {
91 using LeafNodeType = typename TreeType::LeafNodeType;
92 using ValueType = typename TreeType::ValueType;
93
94 ✗ OffsetAndMinComp(std::vector<LeafNodeType*>& lhsNodes,
95 const TreeType& rhsTree,
96 ValueType offset)
97 : mLhsNodes(lhsNodes)
98 , mRhsTree(rhsTree)
99 ✗ , mOffset(offset) {}
100
101 ✗ void operator()(const tbb::blocked_range<size_t>& range) const
102 {
103 using Iterator = typename LeafNodeType::ValueOnIter;
104
105 ✗ tree::ValueAccessor<const TreeType> rhsAcc(mRhsTree);
106 ✗ const ValueType offset = mOffset;
107
108 ✗ for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
109
110 ✗ LeafNodeType& lhsNode = *mLhsNodes[n];
111 const LeafNodeType* rhsNodePt = rhsAcc.probeConstLeaf(lhsNode.origin());
112 ✗ if (!rhsNodePt) continue;
113
114 ✗ auto* const data = lhsNode.buffer().data();
115 ✗ for (Iterator it = lhsNode.beginValueOn(); it; ++it) {
116 ✗ ValueType& val = data[it.pos()];
117 ✗ val = std::min(val, offset + rhsNodePt->getValue(it.pos()));
118 }
119 }
120 }
121
122 private:
123 std::vector<LeafNodeType*>& mLhsNodes;
124 const TreeType& mRhsTree;
125 const ValueType mOffset;
126 }; // struct OffsetAndMinComp
127
128
129 template<typename GridType, typename InterrupterType>
130 inline void
131 2 normalizeLevelSet(GridType& grid, const int halfWidthInVoxels, InterrupterType* interrupt = nullptr)
132 {
133 LevelSetFilter<GridType, GridType, InterrupterType> filter(grid, interrupt);
134 filter.setSpatialScheme(math::FIRST_BIAS);
135 filter.setNormCount(halfWidthInVoxels);
136 filter.normalize();
137
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 filter.prune();
138 2 }
139
140
141 template<typename GridType, typename InterrupterType>
142 inline void
143 ✗ smoothLevelSet(GridType& grid, int iterations, int halfBandWidthInVoxels,
144 InterrupterType* interrupt = nullptr)
145 {
146 using ValueType = typename GridType::ValueType;
147 using TreeType = typename GridType::TreeType;
148 using LeafNodeType = typename TreeType::LeafNodeType;
149
150 ✗ GridType filterGrid(grid);
151
152 LevelSetFilter<GridType, GridType, InterrupterType> filter(filterGrid, interrupt);
153 filter.setSpatialScheme(math::FIRST_BIAS);
154
155 ✗ for (int n = 0; n < iterations; ++n) {
156 ✗ if (interrupt && interrupt->wasInterrupted()) break;
157 ✗ filter.mean(1);
158 }
159
160 std::vector<LeafNodeType*> nodes;
161 grid.tree().getNodes(nodes);
162
163 ✗ const ValueType offset = ValueType(double(0.5) * grid.transform().voxelSize()[0]);
164
165 ✗ tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
166 OffsetAndMinComp<TreeType>(nodes, filterGrid.tree(), -offset));
167
168 // Clean up any damanage that was done by the min operation
169 ✗ normalizeLevelSet(grid, halfBandWidthInVoxels, interrupt);
170 }
171
172 } // namespace ttls_internal
173
174 /// @endcond
175
176 template<typename GridT, typename InterrupterT>
177 typename GridT::template ValueConverter<float>::Type::Ptr
178 2 topologyToLevelSet(const GridT& grid, int halfWidth, int closingSteps, int dilation,
179 int smoothingSteps, InterrupterT* interrupt)
180 {
181 using MaskTreeT = typename GridT::TreeType::template ValueConverter<ValueMask>::Type;
182 using FloatTreeT = typename GridT::TreeType::template ValueConverter<float>::Type;
183 using FloatGridT = Grid<FloatTreeT>;
184
185 // Check inputs
186
187 2 halfWidth = std::max(halfWidth, 1);
188 2 closingSteps = std::max(closingSteps, 0);
189 2 dilation = std::max(dilation, 0);
190
191 2 if (!grid.hasUniformVoxels()) {
192 ✗ OPENVDB_THROW(ValueError, "Non-uniform voxels are not supported!");
193 }
194
195 // Copy the topology into a MaskGrid.
196 4 MaskTreeT maskTree(grid.tree(), false/*background*/, openvdb::TopologyCopy());
197
198 // Morphological closing operation.
199 2 tools::dilateActiveValues(maskTree, closingSteps + dilation, tools::NN_FACE, tools::PRESERVE_TILES);
200 2 tools::erodeActiveValues(maskTree, /*iterations=*/closingSteps, tools::NN_FACE, tools::PRESERVE_TILES);
201 2 tools::pruneInactive(maskTree);
202
203 // Generate a volume with an implicit zero crossing at the boundary
204 // between active and inactive values in the input grid.
205 2 const float background = float(grid.voxelSize()[0]) * float(halfWidth);
206 2 typename FloatTreeT::Ptr lsTree(
207 2 new FloatTreeT(maskTree, /*out=*/background, /*in=*/-background, openvdb::TopologyCopy()));
208
209 tbb::task_group pool;
210 8 pool.run([&]() {
211 2 tools::erodeActiveValues(maskTree, /*iterations=*/halfWidth, tools::NN_FACE, tools::PRESERVE_TILES);
212 2 tools::pruneInactive(maskTree);
213 });
214 4 pool.run([&]() {
215 2 tools::dilateActiveValues(*lsTree, /*iterations=*/halfWidth, tools::NN_FACE, tools::PRESERVE_TILES);
216 });
217 2 pool.wait();// wait for both tasks to complete
218
219 lsTree->topologyDifference(maskTree);
220 2 maskTree.clear();
221
222 lsTree->voxelizeActiveTiles(); // has to come before ::pruneLevelSet
223 2 tools::pruneLevelSet(*lsTree, /*threading=*/true);
224
225 // Create a level set grid from the tree
226 4 typename FloatGridT::Ptr lsGrid = FloatGridT::create(lsTree);
227 4 lsGrid->setTransform(grid.transform().copy());
228 2 lsGrid->setGridClass(openvdb::GRID_LEVEL_SET);
229
230 // Use a PDE based scheme to propagate distance values from the
231 // implicit zero crossing.
232 2 ttls_internal::normalizeLevelSet(*lsGrid, 3*halfWidth, interrupt);
233
234 // Additional filtering
235 2 if (smoothingSteps > 0) {
236 ✗ ttls_internal::smoothLevelSet(*lsGrid, smoothingSteps, halfWidth, interrupt);
237 }
238
239 2 return lsGrid;
240 }
241
242
243 template<typename GridT>
244 typename GridT::template ValueConverter<float>::Type::Ptr
245 2 topologyToLevelSet(const GridT& grid, int halfWidth, int closingSteps, int dilation, int smoothingSteps)
246 {
247 2 util::NullInterrupter interrupt;
248 2 return topologyToLevelSet(grid, halfWidth, closingSteps, dilation, smoothingSteps, &interrupt);
249 }
250
251
252 ////////////////////////////////////////
253
254
255 // Explicit Template Instantiation
256
257 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
258
259 #ifdef OPENVDB_INSTANTIATE_TOPOLOGYTOLEVELSET
260 #include <openvdb/util/ExplicitInstantiation.h>
261 #endif
262
263 #define _FUNCTION(TreeT) \
264 Grid<TreeT>::ValueConverter<float>::Type::Ptr topologyToLevelSet(const Grid<TreeT>&, int, int, int, int, \
265 util::NullInterrupter*)
266 OPENVDB_ALL_TREE_INSTANTIATE(_FUNCTION)
267 #undef _FUNCTION
268
269 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
270
271
272 } // namespace tools
273 } // namespace OPENVDB_VERSION_NAME
274 } // namespace openvdb
275
276 #endif // OPENVDB_TOOLS_TOPOLOGY_TO_LEVELSET_HAS_BEEN_INCLUDED
277
278