OpenVDB  7.0.0
LevelSetFilter.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
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 {
27 namespace OPENVDB_VERSION_NAME {
28 namespace tools {
29 
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:
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 
54  LevelSetFilter(GridType& grid, InterruptT* interrupt = nullptr)
55  : BaseType(grid, interrupt)
56  , mMinMask(0)
57  , mMaxMask(1)
58  , mInvertMask(false)
59  {
60  }
62  ~LevelSetFilter() override {}
63 
66  AlphaType minMask() const { return mMinMask; }
69  AlphaType maxMask() const { return mMaxMask; }
78  {
79  if (!(min < max)) OPENVDB_THROW(ValueError, "Invalid mask range (expects min < max)");
80  mMinMask = min;
81  mMaxMask = max;
82  }
83 
86  bool isMaskInverted() const { return mInvertMask; }
89  void invertMask(bool invert=true) { mInvertMask = invert; }
90 
93  void meanCurvature(const MaskType* mask = nullptr)
94  {
95  Filter f(this, mask); f.meanCurvature();
96  }
97 
100  void laplacian(const MaskType* mask = nullptr)
101  {
102  Filter f(this, mask); f.laplacian();
103  }
104 
111  void gaussian(int width = 1, const MaskType* mask = nullptr)
112  {
113  Filter f(this, mask); f.gaussian(width);
114  }
115 
119  void offset(ValueType offset, const MaskType* mask = nullptr)
120  {
121  Filter f(this, mask); f.offset(offset);
122  }
123 
130  void median(int width = 1, const MaskType* mask = nullptr)
131  {
132  Filter f(this, mask); f.median(width);
133  }
134 
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;
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;
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 
228 
229 template<typename GridT, typename MaskT, typename InterruptT>
230 inline void
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
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
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
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
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
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
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 
346 
348 template<typename GridT, typename MaskT, typename InterruptT>
349 inline void
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 
388 template<typename GridT, typename MaskT, typename InterruptT>
389 inline void
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 
422 template<typename GridT, typename MaskT, typename InterruptT>
423 inline void
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 
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 
478 template<typename GridT, typename MaskT, typename InterruptT>
479 template <typename AvgT>
480 inline void
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 } // namespace tools
507 } // namespace OPENVDB_VERSION_NAME
508 } // namespace openvdb
509 
510 #endif // OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
LevelSetFilter(GridType &grid, InterruptT *interrupt=nullptr)
Main constructor from a grid.
Definition: LevelSetFilter.h:54
ValueType operator()(Coord xyz)
Definition: LevelSetFilter.h:192
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
GridT::ConstAccessor acc
Definition: LevelSetFilter.h:199
void invertMask(bool invert=true)
Invert the optional mask, i.e. min->max in the original mask maps to 1->0 in the inverted alpha mask...
Definition: LevelSetFilter.h:89
const Int32 width
Definition: LevelSetFilter.h:200
typename TreeType::ValueType ValueType
Definition: LevelSetFilter.h:46
typename MaskType::ValueType AlphaType
Definition: LevelSetFilter.h:47
bool checkInterrupter()
Definition: LevelSetTracker.h:382
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:82
Coord Abs(const Coord &xyz)
Definition: Coord.h:515
AlphaType minMask() const
Return the minimum value of the mask to be used for the derivation of a smooth alpha value...
Definition: LevelSetFilter.h:66
typename CopyConstness< TreeType, NonConstBufferType >::Type BufferType
Definition: LeafManager.h:93
bool isMaskInverted() const
Return true if the mask is inverted, i.e. min->max in the original mask maps to 1->0 in the inverted ...
Definition: LevelSetFilter.h:86
FloatT FloatType
Definition: Interpolation.h:551
Type Pow2(Type x)
Return x2.
Definition: Math.h:495
~LevelSetFilter() override
Default destructor.
Definition: LevelSetFilter.h:62
Axis
Definition: Math.h:849
void meanCurvature(const MaskType *mask=nullptr)
One iteration of mean-curvature flow of the level set.
Definition: LevelSetFilter.h:93
const ValueType frac
Definition: LevelSetFilter.h:201
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:603
Definition: Exceptions.h:65
Definition: Interpolation.h:543
GridType::Ptr meanCurvature(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the mean curvature of the given grid.
Definition: GridOperators.h:1030
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:102
GridT GridType
Definition: LevelSetFilter.h:43
int32_t Int32
Definition: Types.h:33
Dense stencil of a given width.
Definition: Stencils.h:1755
void offset(ValueType offset, const MaskType *mask=nullptr)
Offset the level set by the specified (world) distance.
Definition: LevelSetFilter.h:119
void laplacian(const MaskType *mask=nullptr)
One iteration of Laplacian flow of the level set.
Definition: LevelSetFilter.h:100
Definition: Stencils.h:1506
void startInterrupter(const char *msg)
Definition: LevelSetTracker.h:366
Definition: Exceptions.h:13
Performs multi-threaded interface tracking of narrow band level sets.
Definition: LevelSetTracker.h:55
void mean(int width=1, const MaskType *mask=nullptr)
One iteration of mean-value flow of the level set.
Definition: LevelSetFilter.h:140
void setMaskRange(AlphaType min, AlphaType max)
Define the range for the (optional) scalar mask.
Definition: LevelSetFilter.h:77
void gaussian(int width=1, const MaskType *mask=nullptr)
One iteration of a fast separable Gaussian filter.
Definition: LevelSetFilter.h:111
GridType::Ptr laplacian(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the Laplacian of the given scalar grid.
Definition: GridOperators.h:1013
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:102
AlphaType maxMask() const
Return the maximum value of the mask to be used for the derivation of a smooth alpha value...
Definition: LevelSetFilter.h:69
typename GridType::TreeType TreeType
Definition: LevelSetFilter.h:45
Definition: Stencils.h:1225
typename LeafManagerType::LeafRange LeafRange
Definition: LevelSetTracker.h:65
MaskT MaskType
Definition: LevelSetFilter.h:44
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:154
Filtering (e.g. diffusion) of narrow-band level sets. An optional scalar field can be used to produce...
Definition: LevelSetFilter.h:39
Avg(const GridT &grid, Int32 w)
Definition: LevelSetFilter.h:190
LeafManagerType & leafs()
Definition: LevelSetTracker.h:187
void median(int width=1, const MaskType *mask=nullptr)
One iteration of median-value flow of the level set.
Definition: LevelSetFilter.h:130
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:106