GCC Code Coverage Report

Directory: ./
File: openvdb/openvdb/tools/LevelSetFilter.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 3 175 1.7%
Functions: 0 74 0.0%
Branches: 1 290 0.3%
Line Branch Exec Source
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
4 /// @author Ken Museth
5 ///
6 /// @file tools/LevelSetFilter.h
7 ///
8 /// @brief Performs various types of level set deformations with
9 /// interface tracking. These unrestricted deformations include
10 /// surface smoothing (e.g., Laplacian flow), filtering (e.g., mean
11 /// value) and morphological operations (e.g., morphological opening).
12 /// All these operations can optionally be masked with another grid that
13 /// acts as an alpha-mask.
14
15 #ifndef OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
16 #define OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
17
18 #include "LevelSetTracker.h"
19 #include "Interpolation.h"
20 #include <algorithm> // for std::max()
21 #include <functional>
22 #include <type_traits>
23
24
25 namespace openvdb {
26 OPENVDB_USE_VERSION_NAMESPACE
27 namespace OPENVDB_VERSION_NAME {
28 namespace tools {
29
30 /// @brief Filtering (e.g. diffusion) of narrow-band level sets. An
31 /// optional scalar field can be used to produce a (smooth) alpha mask
32 /// for the filtering.
33 ///
34 /// @note This class performs proper interface tracking which allows
35 /// for unrestricted surface deformations
36 template<typename GridT,
37 typename MaskT = typename GridT::template ValueConverter<float>::Type,
38 typename InterruptT = util::NullInterrupter>
39 class LevelSetFilter : public LevelSetTracker<GridT, InterruptT>
40 {
41 public:
42 using BaseType = LevelSetTracker<GridT, InterruptT>;
43 using GridType = GridT;
44 using MaskType = MaskT;
45 using TreeType = typename GridType::TreeType;
46 using ValueType = typename TreeType::ValueType;
47 using AlphaType = typename MaskType::ValueType;
48 static_assert(std::is_floating_point<AlphaType>::value,
49 "LevelSetFilter requires a mask grid with floating-point values");
50
51 /// @brief Main constructor from a grid
52 /// @param grid The level set to be filtered.
53 /// @param interrupt Optional interrupter.
54 2 LevelSetFilter(GridType& grid, InterruptT* interrupt = nullptr)
55 : BaseType(grid, interrupt)
56 , mMinMask(0)
57 , mMaxMask(1)
58
1/4
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
 2 , mInvertMask(false)
59 {
60 }
61 /// @brief Default destructor
62 2 ~LevelSetFilter() override {}
63
64 /// @brief Return the minimum value of the mask to be used for the
65 /// derivation of a smooth alpha value.
66 ✗ AlphaType minMask() const { return mMinMask; }
67 /// @brief Return the maximum value of the mask to be used for the
68 /// derivation of a smooth alpha value.
69 ✗ AlphaType maxMask() const { return mMaxMask; }
70 /// @brief Define the range for the (optional) scalar mask.
71 /// @param min Minimum value of the range.
72 /// @param max Maximum value of the range.
73 /// @details Mask values outside the range maps to alpha values of
74 /// respectfully zero and one, and values inside the range maps
75 /// smoothly to 0->1 (unless of course the mask is inverted).
76 /// @throw ValueError if @a min is not smaller than @a max.
77 ✗ void setMaskRange(AlphaType min, AlphaType max)
78 {
79 ✗ if (!(min < max)) OPENVDB_THROW(ValueError, "Invalid mask range (expects min < max)");
80 ✗ mMinMask = min;
81 ✗ mMaxMask = max;
82 }
83
84 /// @brief Return true if the mask is inverted, i.e. min->max in the
85 /// original mask maps to 1->0 in the inverted alpha mask.
86 ✗ bool isMaskInverted() const { return mInvertMask; }
87 /// @brief Invert the optional mask, i.e. min->max in the original
88 /// mask maps to 1->0 in the inverted alpha mask.
89 ✗ void invertMask(bool invert=true) { mInvertMask = invert; }
90
91 /// @brief One iteration of mean-curvature flow of the level set.
92 /// @param mask Optional alpha mask.
93 ✗ void meanCurvature(const MaskType* mask = nullptr)
94 {
95 ✗ Filter f(this, mask); f.meanCurvature();
96 }
97
98 /// @brief One iteration of Laplacian flow of the level set.
99 /// @param mask Optional alpha mask.
100 ✗ void laplacian(const MaskType* mask = nullptr)
101 {
102 ✗ Filter f(this, mask); f.laplacian();
103 }
104
105 /// @brief One iteration of a fast separable Gaussian filter.
106 /// @param width Width of the Gaussian kernel in voxel units.
107 /// @param mask Optional alpha mask.
108 ///
109 /// @note This is approximated as 4 iterations of a separable mean filter
110 /// which typically leads an approximation that's better than 95%!
111 ✗ void gaussian(int width = 1, const MaskType* mask = nullptr)
112 {
113 ✗ Filter f(this, mask); f.gaussian(width);
114 }
115
116 /// @brief Offset the level set by the specified (world) distance.
117 /// @param offset Value of the offset.
118 /// @param mask Optional alpha mask.
119 ✗ void offset(ValueType offset, const MaskType* mask = nullptr)
120 {
121 ✗ Filter f(this, mask); f.offset(offset);
122 }
123
124 /// @brief One iteration of median-value flow of the level set.
125 /// @param width Width of the median-value kernel in voxel units.
126 /// @param mask Optional alpha mask.
127 ///
128 /// @warning This filter is not separable and is hence relatively
129 /// slow!
130 ✗ void median(int width = 1, const MaskType* mask = nullptr)
131 {
132 ✗ Filter f(this, mask); f.median(width);
133 }
134
135 /// @brief One iteration of mean-value flow of the level set.
136 /// @param width Width of the mean-value kernel in voxel units.
137 /// @param mask Optional alpha mask.
138 ///
139 /// @note This filter is separable so it's fast!
140 ✗ void mean(int width = 1, const MaskType* mask = nullptr)
141 {
142 ✗ Filter f(this, mask); f.mean(width);
143 }
144
145 private:
146 // disallow copy construction and copy by assignment!
147 LevelSetFilter(const LevelSetFilter&);// not implemented
148 LevelSetFilter& operator=(const LevelSetFilter&);// not implemented
149
150 // Private struct that implements all the filtering.
151 struct Filter
152 {
153 using LeafT = typename TreeType::LeafNodeType;
154 using VoxelIterT = typename LeafT::ValueOnIter;
155 using VoxelCIterT = typename LeafT::ValueOnCIter;
156 using BufferT = typename tree::LeafManager<TreeType>::BufferType;
157 using LeafRange = typename tree::LeafManager<TreeType>::LeafRange;
158 using LeafIterT = typename LeafRange::Iterator;
159 using AlphaMaskT = tools::AlphaMask<GridT, MaskT>;
160
161 ✗ Filter(LevelSetFilter* parent, const MaskType* mask) : mParent(parent), mMask(mask) {}
162 ✗ Filter(const Filter&) = default;
163 ✗ virtual ~Filter() {}
164
165 void box(int width);
166 void median(int width);
167 void mean(int width);
168 void gaussian(int width);
169 void laplacian();
170 void meanCurvature();
171 void offset(ValueType value);
172 ✗ void operator()(const LeafRange& r) const
173 {
174 ✗ if (mTask) mTask(const_cast<Filter*>(this), r);
175 ✗ else OPENVDB_THROW(ValueError, "task is undefined - don\'t call this method directly");
176 }
177 ✗ void cook(bool swap)
178 {
179 ✗ const int n = mParent->getGrainSize();
180 ✗ if (n>0) {
181 ✗ tbb::parallel_for(mParent->leafs().leafRange(n), *this);
182 } else {
183 ✗ (*this)(mParent->leafs().leafRange());
184 }
185 ✗ if (swap) mParent->leafs().swapLeafBuffer(1, n==0);
186 }
187
188 template <size_t Axis>
189 struct Avg {
190 ✗ Avg(const GridT& grid, Int32 w) :
191 ✗ acc(grid.tree()), width(w), frac(1/ValueType(2*w+1)) {}
192 ✗ inline ValueType operator()(Coord xyz)
193 {
194 ValueType sum = zeroVal<ValueType>();
195 ✗ Int32& i = xyz[Axis], j = i + width;
196 ✗ for (i -= width; i <= j; ++i) sum += acc.getValue(xyz);
197 ✗ return sum*frac;
198 }
199 typename GridT::ConstAccessor acc;
200 const Int32 width;
201 const ValueType frac;
202 };
203
204 template<typename AvgT>
205 void boxImpl(const LeafRange& r, Int32 w);
206
207 ✗ void boxXImpl(const LeafRange& r, Int32 w) { this->boxImpl<Avg<0> >(r,w); }
208 ✗ void boxZImpl(const LeafRange& r, Int32 w) { this->boxImpl<Avg<1> >(r,w); }
209 ✗ void boxYImpl(const LeafRange& r, Int32 w) { this->boxImpl<Avg<2> >(r,w); }
210
211 void medianImpl(const LeafRange&, int);
212 void meanCurvatureImpl(const LeafRange&);
213 void laplacianImpl(const LeafRange&);
214 void offsetImpl(const LeafRange&, ValueType);
215
216 LevelSetFilter* mParent;
217 const MaskType* mMask;
218 typename std::function<void (Filter*, const LeafRange&)> mTask;
219 }; // end of private Filter struct
220
221 AlphaType mMinMask, mMaxMask;
222 bool mInvertMask;
223
224 }; // end of LevelSetFilter class
225
226
227 ////////////////////////////////////////
228
229 template<typename GridT, typename MaskT, typename InterruptT>
230 inline void
231 ✗ LevelSetFilter<GridT, MaskT, InterruptT>::Filter::median(int width)
232 {
233 ✗ mParent->startInterrupter("Median-value flow of level set");
234
235 ✗ mParent->leafs().rebuildAuxBuffers(1, mParent->getGrainSize()==0);
236
237 ✗ mTask = std::bind(&Filter::medianImpl,
238 ✗ std::placeholders::_1, std::placeholders::_2, std::max(1, width));
239 ✗ this->cook(true);
240
241 ✗ mParent->track();
242
243 ✗ mParent->endInterrupter();
244 }
245
246 template<typename GridT, typename MaskT, typename InterruptT>
247 inline void
248 ✗ LevelSetFilter<GridT, MaskT, InterruptT>::Filter::mean(int width)
249 {
250 ✗ mParent->startInterrupter("Mean-value flow of level set");
251
252 ✗ this->box(width);
253
254 ✗ mParent->endInterrupter();
255 }
256
257 template<typename GridT, typename MaskT, typename InterruptT>
258 inline void
259 ✗ LevelSetFilter<GridT, MaskT, InterruptT>::Filter::gaussian(int width)
260 {
261 ✗ mParent->startInterrupter("Gaussian flow of level set");
262
263 ✗ for (int n=0; n<4; ++n) this->box(width);
264
265 ✗ mParent->endInterrupter();
266 }
267
268 template<typename GridT, typename MaskT, typename InterruptT>
269 inline void
270 ✗ LevelSetFilter<GridT, MaskT, InterruptT>::Filter::box(int width)
271 {
272 ✗ mParent->leafs().rebuildAuxBuffers(1, mParent->getGrainSize()==0);
273
274 ✗ width = std::max(1, width);
275
276 ✗ mTask = std::bind(&Filter::boxXImpl, std::placeholders::_1, std::placeholders::_2, width);
277 ✗ this->cook(true);
278
279 ✗ mTask = std::bind(&Filter::boxYImpl, std::placeholders::_1, std::placeholders::_2, width);
280 ✗ this->cook(true);
281
282 ✗ mTask = std::bind(&Filter::boxZImpl, std::placeholders::_1, std::placeholders::_2, width);
283 ✗ this->cook(true);
284
285 ✗ mParent->track();
286 }
287
288 template<typename GridT, typename MaskT, typename InterruptT>
289 inline void
290 ✗ LevelSetFilter<GridT, MaskT, InterruptT>::Filter::meanCurvature()
291 {
292 ✗ mParent->startInterrupter("Mean-curvature flow of level set");
293
294 ✗ mParent->leafs().rebuildAuxBuffers(1, mParent->getGrainSize()==0);
295
296 ✗ mTask = std::bind(&Filter::meanCurvatureImpl, std::placeholders::_1, std::placeholders::_2);
297 ✗ this->cook(true);
298
299 ✗ mParent->track();
300
301 ✗ mParent->endInterrupter();
302 }
303
304 template<typename GridT, typename MaskT, typename InterruptT>
305 inline void
306 ✗ LevelSetFilter<GridT, MaskT, InterruptT>::Filter::laplacian()
307 {
308 ✗ mParent->startInterrupter("Laplacian flow of level set");
309
310 ✗ mParent->leafs().rebuildAuxBuffers(1, mParent->getGrainSize()==0);
311
312 ✗ mTask = std::bind(&Filter::laplacianImpl, std::placeholders::_1, std::placeholders::_2);
313 ✗ this->cook(true);
314
315 ✗ mParent->track();
316
317 ✗ mParent->endInterrupter();
318 }
319
320 template<typename GridT, typename MaskT, typename InterruptT>
321 inline void
322 ✗ LevelSetFilter<GridT, MaskT, InterruptT>::Filter::offset(ValueType value)
323 {
324 ✗ mParent->startInterrupter("Offsetting level set");
325
326 ✗ mParent->leafs().removeAuxBuffers();// no auxiliary buffers required
327
328 ✗ const ValueType CFL = ValueType(0.5) * mParent->voxelSize(), offset = openvdb::math::Abs(value);
329 ValueType dist = 0.0;
330 ✗ while (offset-dist > ValueType(0.001)*CFL && mParent->checkInterrupter()) {
331 ✗ const ValueType delta = openvdb::math::Min(offset-dist, CFL);
332 ✗ dist += delta;
333
334 ✗ mTask = std::bind(&Filter::offsetImpl,
335 ✗ std::placeholders::_1, std::placeholders::_2, copysign(delta, value));
336 ✗ this->cook(false);
337
338 ✗ mParent->track();
339 }
340
341 ✗ mParent->endInterrupter();
342 }
343
344
345 ///////////////////////// PRIVATE METHODS //////////////////////
346
347 /// Performs parabolic mean-curvature diffusion
348 template<typename GridT, typename MaskT, typename InterruptT>
349 inline void
350 ✗ LevelSetFilter<GridT, MaskT, InterruptT>::Filter::meanCurvatureImpl(const LeafRange& range)
351 {
352 ✗ mParent->checkInterrupter();
353 //const float CFL = 0.9f, dt = CFL * mDx * mDx / 6.0f;
354 ✗ const ValueType dx = mParent->voxelSize(), dt = math::Pow2(dx) / ValueType(3.0);
355 ✗ math::CurvatureStencil<GridType> stencil(mParent->grid(), dx);
356 ✗ if (mMask) {
357 typename AlphaMaskT::FloatType a, b;
358 ✗ AlphaMaskT alpha(mParent->grid(), *mMask, mParent->minMask(),
359 ✗ mParent->maxMask(), mParent->isMaskInverted());
360 ✗ for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
361 ✗ ValueType* buffer = leafIter.buffer(1).data();
362 ✗ for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
363 ✗ if (alpha(iter.getCoord(), a, b)) {
364 ✗ stencil.moveTo(iter);
365 ✗ const ValueType phi0 = *iter, phi1 = phi0 + dt*stencil.meanCurvatureNormGrad();
366 ✗ buffer[iter.pos()] = b * phi0 + a * phi1;
367 }
368 }
369 }
370 } else {
371 ✗ for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
372 ✗ ValueType* buffer = leafIter.buffer(1).data();
373 ✗ for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
374 ✗ stencil.moveTo(iter);
375 ✗ buffer[iter.pos()] = *iter + dt*stencil.meanCurvatureNormGrad();
376 }
377 }
378 }
379 }
380
381 /// Performs Laplacian diffusion. Note if the grids contains a true
382 /// signed distance field (e.g. a solution to the Eikonal equation)
383 /// Laplacian diffusions (e.g. geometric heat equation) is actually
384 /// identical to mean curvature diffusion, yet less computationally
385 /// expensive! In other words if you're performing renormalization
386 /// anyway (e.g. rebuilding the narrow-band) you should consider
387 /// performing Laplacian diffusion over mean curvature flow!
388 template<typename GridT, typename MaskT, typename InterruptT>
389 inline void
390 ✗ LevelSetFilter<GridT, MaskT, InterruptT>::Filter::laplacianImpl(const LeafRange& range)
391 {
392 ✗ mParent->checkInterrupter();
393 //const float CFL = 0.9f, half_dt = CFL * mDx * mDx / 12.0f;
394 ✗ const ValueType dx = mParent->voxelSize(), dt = math::Pow2(dx) / ValueType(6.0);
395 ✗ math::GradStencil<GridType> stencil(mParent->grid(), dx);
396 ✗ if (mMask) {
397 typename AlphaMaskT::FloatType a, b;
398 ✗ AlphaMaskT alpha(mParent->grid(), *mMask, mParent->minMask(),
399 ✗ mParent->maxMask(), mParent->isMaskInverted());
400 ✗ for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
401 ✗ ValueType* buffer = leafIter.buffer(1).data();
402 ✗ for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
403 ✗ if (alpha(iter.getCoord(), a, b)) {
404 ✗ stencil.moveTo(iter);
405 ✗ const ValueType phi0 = *iter, phi1 = phi0 + dt*stencil.laplacian();
406 ✗ buffer[iter.pos()] = b * phi0 + a * phi1;
407 }
408 }
409 }
410 } else {
411 ✗ for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
412 ✗ ValueType* buffer = leafIter.buffer(1).data();
413 ✗ for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
414 ✗ stencil.moveTo(iter);
415 ✗ buffer[iter.pos()] = *iter + dt*stencil.laplacian();
416 }
417 }
418 }
419 }
420
421 /// Offsets the values by a constant
422 template<typename GridT, typename MaskT, typename InterruptT>
423 inline void
424 ✗ LevelSetFilter<GridT, MaskT, InterruptT>::Filter::offsetImpl(
425 const LeafRange& range, ValueType offset)
426 {
427 ✗ mParent->checkInterrupter();
428 ✗ if (mMask) {
429 typename AlphaMaskT::FloatType a, b;
430 ✗ AlphaMaskT alpha(mParent->grid(), *mMask, mParent->minMask(),
431 ✗ mParent->maxMask(), mParent->isMaskInverted());
432 ✗ for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
433 ✗ for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
434 ✗ if (alpha(iter.getCoord(), a, b)) iter.setValue(*iter + a*offset);
435 }
436 }
437 } else {
438 ✗ for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
439 ✗ for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
440 ✗ iter.setValue(*iter + offset);
441 }
442 }
443 }
444 }
445
446 /// Performs simple but slow median-value diffusion
447 template<typename GridT, typename MaskT, typename InterruptT>
448 inline void
449 ✗ LevelSetFilter<GridT, MaskT, InterruptT>::Filter::medianImpl(const LeafRange& range, int width)
450 {
451 ✗ mParent->checkInterrupter();
452 ✗ typename math::DenseStencil<GridType> stencil(mParent->grid(), width);//creates local cache!
453 ✗ if (mMask) {
454 typename AlphaMaskT::FloatType a, b;
455 ✗ AlphaMaskT alpha(mParent->grid(), *mMask, mParent->minMask(),
456 ✗ mParent->maxMask(), mParent->isMaskInverted());
457 ✗ for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
458 ✗ ValueType* buffer = leafIter.buffer(1).data();
459 ✗ for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
460 ✗ if (alpha(iter.getCoord(), a, b)) {
461 stencil.moveTo(iter);
462 ✗ buffer[iter.pos()] = b * (*iter) + a * stencil.median();
463 }
464 }
465 }
466 } else {
467 ✗ for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
468 ✗ ValueType* buffer = leafIter.buffer(1).data();
469 ✗ for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
470 stencil.moveTo(iter);
471 ✗ buffer[iter.pos()] = stencil.median();
472 }
473 }
474 }
475 }
476
477 /// One dimensional convolution of a separable box filter
478 template<typename GridT, typename MaskT, typename InterruptT>
479 template <typename AvgT>
480 inline void
481 ✗ LevelSetFilter<GridT, MaskT, InterruptT>::Filter::boxImpl(const LeafRange& range, Int32 w)
482 {
483 ✗ mParent->checkInterrupter();
484 ✗ AvgT avg(mParent->grid(), w);
485 ✗ if (mMask) {
486 typename AlphaMaskT::FloatType a, b;
487 ✗ AlphaMaskT alpha(mParent->grid(), *mMask, mParent->minMask(),
488 ✗ mParent->maxMask(), mParent->isMaskInverted());
489 ✗ for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
490 ✗ ValueType* buffer = leafIter.buffer(1).data();
491 ✗ for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
492 ✗ const Coord xyz = iter.getCoord();
493 ✗ if (alpha(xyz, a, b)) buffer[iter.pos()] = b * (*iter)+ a * avg(xyz);
494 }
495 }
496 } else {
497 ✗ for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
498 ✗ ValueType* buffer = leafIter.buffer(1).data();
499 ✗ for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
500 ✗ buffer[iter.pos()] = avg(iter.getCoord());
501 }
502 }
503 }
504 }
505
506
507 ////////////////////////////////////////
508
509
510 // Explicit Template Instantiation
511
512 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
513
514 #ifdef OPENVDB_INSTANTIATE_LEVELSETFILTER
515 #include <openvdb/util/ExplicitInstantiation.h>
516 #endif
517
518 OPENVDB_INSTANTIATE_CLASS LevelSetFilter<FloatGrid, FloatGrid, util::NullInterrupter>;
519 OPENVDB_INSTANTIATE_CLASS LevelSetFilter<DoubleGrid, FloatGrid, util::NullInterrupter>;
520
521 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
522
523
524 } // namespace tools
525 } // namespace OPENVDB_VERSION_NAME
526 } // namespace openvdb
527
528 #endif // OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
529