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