OpenVDB  2.3.0
LevelSetFilter.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2013 DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
30 //
41 
42 #ifndef OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
43 #define OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
44 
45 #include <assert.h>
46 #include <boost/type_traits/is_floating_point.hpp>
47 #include "LevelSetTracker.h"
48 #include "Interpolation.h"
49 
50 namespace openvdb {
52 namespace OPENVDB_VERSION_NAME {
53 namespace tools {
54 
61 template<typename GridT,
62  typename MaskT = typename GridT::template ValueConverter<float>::Type,
63  typename InterruptT = util::NullInterrupter>
64 class LevelSetFilter : public LevelSetTracker<GridT, InterruptT>
65 {
66 public:
68  typedef GridT GridType;
69  typedef MaskT MaskType;
70  typedef typename GridType::TreeType TreeType;
71  typedef typename TreeType::ValueType ValueType;
72  typedef typename MaskType::ValueType AlphaType;
74  BOOST_STATIC_ASSERT(boost::is_floating_point<AlphaType>::value);
75 
79  LevelSetFilter(GridType& grid, InterruptT* interrupt = NULL)
80  : BaseType(grid, interrupt)
81  , mTask(0)
82  , mMask(NULL)
83  , mMinMask(0)
84  , mMaxMask(1)
85  , mInvertMask(false)
86  {
87  }
92  : BaseType(other)
93  , mTask(other.mTask)
94  , mMask(other.mMask)
95  , mMinMask(other.mMinMask)
96  , mMaxMask(other.mMaxMask)
97  , mInvertMask(other.mInvertMask)
98  {
99  }
101  virtual ~LevelSetFilter() {};
102 
106  void operator()(const RangeType& range) const
107  {
108  if (mTask) mTask(const_cast<LevelSetFilter*>(this), range);
109  else OPENVDB_THROW(ValueError, "task is undefined - call offset(), etc");
110  }
111 
114  AlphaType minMask() const { return mMinMask; }
117  AlphaType maxMask() const { return mMaxMask; }
126  {
127  if (!(min < max)) OPENVDB_THROW(ValueError, "Invalid mask range (expects min < max)");
128  mMinMask = min;
129  mMaxMask = max;
130  }
131 
134  bool isMaskInverted() const { return mInvertMask; }
137  void invertMask(bool invert=true) { mInvertMask = invert; }
138 
141  void meanCurvature(const MaskType* mask = NULL);
142 
145  void laplacian(const MaskType* mask = NULL);
146 
153  void gaussian(int width = 1, const MaskType* mask = NULL);
154 
158  void offset(ValueType offset, const MaskType* mask = NULL);
159 
166  void median(int width = 1, const MaskType* mask = NULL);
167 
173  void mean(int width = 1, const MaskType* mask = NULL);
174 
175 private:
176  typedef typename TreeType::LeafNodeType LeafT;
177  typedef typename LeafT::ValueOnIter VoxelIterT;
178  typedef typename LeafT::ValueOnCIter VoxelCIterT;
179  typedef typename tree::LeafManager<TreeType>::BufferType BufferT;
180  typedef typename RangeType::Iterator LeafIterT;
181 
182  // Only two private member data
183  typename boost::function<void (LevelSetFilter*, const RangeType&)> mTask;
184  const MaskType* mMask;
185  AlphaType mMinMask, mMaxMask;
186  bool mInvertMask;
187 
188  // Private cook method calling tbb::parallel_for
189  void cook(bool swap)
190  {
191  const int n = BaseType::getGrainSize();
192  if (n>0) {
193  tbb::parallel_for(BaseType::leafs().leafRange(n), *this);
194  } else {
195  (*this)(BaseType::leafs().leafRange());
196  }
197  if (swap) BaseType::leafs().swapLeafBuffer(1, n==0);
198  }
199 
200  // Private class to derive the normalized alpha mask
201  struct AlphaMask
202  {
203  AlphaMask(const GridType& grid, const MaskType& mask,
204  AlphaType min, AlphaType max, bool invert)
205  : mAcc(mask.tree()), mSampler(mAcc, mask.transform(), grid.transform()),
206  mMin(min), mInvNorm(1/(max-min)), mInvert(invert)
207  {
208  assert(min < max);
209  }
210  inline bool operator()(const Coord& xyz, AlphaType& a, AlphaType& b) const
211  {
212  a = mSampler(xyz);
213  const AlphaType t = (a-mMin)*mInvNorm;
214  a = t > 0 ? t < 1 ? (3-2*t)*t*t : 1 : 0;//smooth mapping to 0->1
215  b = 1 - a;
216  if (mInvert) std::swap(a,b);
217  return a>0;
218  }
219  typedef typename MaskType::ConstAccessor AccType;
220  AccType mAcc;
221  tools::DualGridSampler<AccType, tools::BoxSampler> mSampler;
222  const AlphaType mMin, mInvNorm;
223  const bool mInvert;
224  };
225 
226  // Private driver method for mean and gaussian filtering
227  void box(int width);
228 
229  template <size_t Axis>
230  struct Avg {
231  Avg(const GridT& grid, Int32 w) :
232  acc(grid.tree()), width(w), frac(1/ValueType(2*w+1)) {}
233  ValueType operator()(Coord xyz) {
234  ValueType sum = zeroVal<ValueType>();
235  Int32& i = xyz[Axis], j = i + width;
236  for (i -= width; i <= j; ++i) sum += acc.getValue(xyz);
237  return sum*frac;
238  }
239  typename GridT::ConstAccessor acc;
240  const Int32 width;
241  const ValueType frac;
242  };
243 
244  // Private methods called by tbb::parallel_for threads
245  template <typename AvgT>
246  void doBox( const RangeType& r, Int32 w);
247  void doBoxX(const RangeType& r, Int32 w) { this->doBox<Avg<0> >(r,w); }
248  void doBoxZ(const RangeType& r, Int32 w) { this->doBox<Avg<1> >(r,w); }
249  void doBoxY(const RangeType& r, Int32 w) { this->doBox<Avg<2> >(r,w); }
250  void doMedian(const RangeType&, int);
251  void doMeanCurvature(const RangeType&);
252  void doLaplacian(const RangeType&);
253  void doOffset(const RangeType&, ValueType);
254 
255 }; // end of LevelSetFilter class
256 
257 
259 
260 template<typename GridT, typename MaskT, typename InterruptT>
261 inline void
263 {
264  mMask = mask;
265 
266  BaseType::startInterrupter("Median-value flow of level set");
267 
268  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
269 
270  mTask = boost::bind(&LevelSetFilter::doMedian, _1, _2, std::max(1, width));
271  this->cook(true);
272 
273  BaseType::track();
274 
275  BaseType::endInterrupter();
276 }
277 
278 template<typename GridT, typename MaskT, typename InterruptT>
279 inline void
281 {
282  mMask = mask;
283 
284  BaseType::startInterrupter("Mean-value flow of level set");
285 
286  this->box(width);
287 
288  BaseType::endInterrupter();
289 }
290 
291 template<typename GridT, typename MaskT, typename InterruptT>
292 inline void
294 {
295  mMask = mask;
296 
297  BaseType::startInterrupter("Gaussian flow of level set");
298 
299  for (int n=0; n<4; ++n) this->box(width);
300 
301  BaseType::endInterrupter();
302 }
303 
304 template<typename GridT, typename MaskT, typename InterruptT>
305 inline void
307 {
308  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
309 
310  width = std::max(1, width);
311 
312  mTask = boost::bind(&LevelSetFilter::doBoxX, _1, _2, width);
313  this->cook(true);
314 
315  mTask = boost::bind(&LevelSetFilter::doBoxY, _1, _2, width);
316  this->cook(true);
317 
318  mTask = boost::bind(&LevelSetFilter::doBoxZ, _1, _2, width);
319  this->cook(true);
320 
321  BaseType::track();
322 }
323 
324 template<typename GridT, typename MaskT, typename InterruptT>
325 inline void
327 {
328  mMask = mask;
329 
330  BaseType::startInterrupter("Mean-curvature flow of level set");
331 
332  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
333 
334  mTask = boost::bind(&LevelSetFilter::doMeanCurvature, _1, _2);
335  this->cook(true);
336 
337  BaseType::track();
338 
339  BaseType::endInterrupter();
340 }
341 
342 template<typename GridT, typename MaskT, typename InterruptT>
343 inline void
345 {
346  mMask = mask;
347 
348  BaseType::startInterrupter("Laplacian flow of level set");
349 
350  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
351 
352  mTask = boost::bind(&LevelSetFilter::doLaplacian, _1, _2);
353  this->cook(true);
354 
355  BaseType::track();
356 
357  BaseType::endInterrupter();
358 }
359 
360 template<typename GridT, typename MaskT, typename InterruptT>
361 inline void
363 {
364  mMask = mask;
365 
366  BaseType::startInterrupter("Offsetting level set");
367 
368  BaseType::leafs().removeAuxBuffers();// no auxiliary buffers required
369 
370  const ValueType CFL = ValueType(0.5) * BaseType::voxelSize(), offset = openvdb::math::Abs(value);
371  ValueType dist = 0.0;
372  while (offset-dist > ValueType(0.001)*CFL && BaseType::checkInterrupter()) {
373  const ValueType delta = openvdb::math::Min(offset-dist, CFL);
374  dist += delta;
375 
376  mTask = boost::bind(&LevelSetFilter::doOffset, _1, _2, copysign(delta, value));
377  this->cook(false);
378 
379  BaseType::track();
380  }
381 
382  BaseType::endInterrupter();
383 }
384 
385 
387 
389 template<typename GridT, typename MaskT, typename InterruptT>
390 inline void
392 {
393  BaseType::checkInterrupter();
394  //const float CFL = 0.9f, dt = CFL * mDx * mDx / 6.0f;
395  const ValueType dx = BaseType::voxelSize(), dt = math::Pow2(dx) / ValueType(3.0);
396  math::CurvatureStencil<GridType> stencil(BaseType::grid(), dx);
397  if (mMask) {
398  AlphaType a, b;
399  AlphaMask alpha(BaseType::grid(), *mMask, mMinMask, mMaxMask, mInvertMask);
400  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
401  BufferT& buffer = leafIter.buffer(1);
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.meanCurvatureNormGrad();
406  buffer.setValue(iter.pos(), b*phi0 + a*phi1);
407  }
408  }
409  }
410  } else {
411  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
412  BufferT& buffer = leafIter.buffer(1);
413  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
414  stencil.moveTo(iter);
415  buffer.setValue(iter.pos(), *iter + dt*stencil.meanCurvatureNormGrad());
416  }
417  }
418  }
419 }
420 
428 template<typename GridT, typename MaskT, typename InterruptT>
429 inline void
430 LevelSetFilter<GridT, MaskT, InterruptT>::doLaplacian(const RangeType& range)
431 {
432  BaseType::checkInterrupter();
433  //const float CFL = 0.9f, half_dt = CFL * mDx * mDx / 12.0f;
434  const ValueType dx = BaseType::voxelSize(), dt = math::Pow2(dx) / ValueType(6.0);
435  math::GradStencil<GridType> stencil(BaseType::grid(), dx);
436  if (mMask) {
437  AlphaType a, b;
438  AlphaMask alpha(BaseType::grid(), *mMask, mMinMask, mMaxMask, mInvertMask);
439  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
440  BufferT& buffer = leafIter.buffer(1);
441  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
442  if (alpha(iter.getCoord(), a, b)) {
443  stencil.moveTo(iter);
444  const ValueType phi0 = *iter, phi1 = phi0 + dt*stencil.laplacian();
445  buffer.setValue(iter.pos(), b*phi0 + a*phi1);
446  }
447  }
448  }
449  } else {
450  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
451  BufferT& buffer = leafIter.buffer(1);
452  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
453  stencil.moveTo(iter);
454  buffer.setValue(iter.pos(), *iter + dt*stencil.laplacian());
455  }
456  }
457  }
458 }
459 
461 template<typename GridT, typename MaskT, typename InterruptT>
462 inline void
463 LevelSetFilter<GridT, MaskT, InterruptT>::doOffset(const RangeType& range, ValueType offset)
464 {
465  BaseType::checkInterrupter();
466  if (mMask) {
467  AlphaType a, b;
468  AlphaMask alpha(BaseType::grid(), *mMask, mMinMask, mMaxMask, mInvertMask);
469  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
470  for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
471  if (alpha(iter.getCoord(), a, b)) iter.setValue(*iter + a*offset);
472  }
473  }
474  } else {
475  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
476  for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
477  iter.setValue(*iter + offset);
478  }
479  }
480  }
481 }
482 
484 template<typename GridT, typename MaskT, typename InterruptT>
485 inline void
486 LevelSetFilter<GridT, MaskT, InterruptT>::doMedian(const RangeType& range, int width)
487 {
488  BaseType::checkInterrupter();
489  typename math::DenseStencil<GridType> stencil(BaseType::grid(), width);//creates local cache!
490  if (mMask) {
491  AlphaType a, b;
492  AlphaMask alpha(BaseType::grid(), *mMask, mMinMask, mMaxMask, mInvertMask);
493  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
494  BufferT& buffer = leafIter.buffer(1);
495  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
496  if (alpha(iter.getCoord(), a, b)) {
497  stencil.moveTo(iter);
498  buffer.setValue(iter.pos(), b*(*iter) + a*stencil.median());
499  }
500  }
501  }
502  } else {
503  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
504  BufferT& buffer = leafIter.buffer(1);
505  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
506  stencil.moveTo(iter);
507  buffer.setValue(iter.pos(), stencil.median());
508  }
509  }
510  }
511 }
512 
514 template<typename GridT, typename MaskT, typename InterruptT>
515 template <typename AvgT>
516 inline void
517 LevelSetFilter<GridT, MaskT, InterruptT>::doBox(const RangeType& range, Int32 w)
518 {
519  BaseType::checkInterrupter();
520  AvgT avg(BaseType::grid(), w);
521  if (mMask) {
522  AlphaType a, b;
523  AlphaMask alpha(BaseType::grid(), *mMask, mMinMask, mMaxMask, mInvertMask);
524  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
525  BufferT& buffer = leafIter.buffer(1);
526  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
527  const Coord xyz = iter.getCoord();
528  if (alpha(xyz, a, b)) buffer.setValue(iter.pos(), b*(*iter)+ a*avg(xyz));
529  }
530  }
531  } else {
532  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
533  BufferT& buffer = leafIter.buffer(1);
534  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
535  buffer.setValue(iter.pos(), avg(iter.getCoord()));
536  }
537  }
538  }
539 }
540 
541 } // namespace tools
542 } // namespace OPENVDB_VERSION_NAME
543 } // namespace openvdb
544 
545 #endif // OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
546 
547 // Copyright (c) 2012-2013 DreamWorks Animation LLC
548 // All rights reserved. This software is distributed under the
549 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
GridT GridType
Definition: LevelSetFilter.h:68
OPENVDB_API Hermite min(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
MaskType::ValueType AlphaType
Definition: LevelSetFilter.h:72
LevelSetTracker< GridT, InterruptT > BaseType
Definition: LevelSetFilter.h:67
GridType::Ptr laplacian(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the Laplacian of the given scalar grid.
Definition: GridOperators.h:1016
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
MaskT MaskType
Definition: LevelSetFilter.h:69
void setMaskRange(AlphaType min, AlphaType max)
Define the range for the (optional) scalar mask.
Definition: LevelSetFilter.h:125
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:244
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:134
Type Pow2(Type x)
Return .
Definition: Math.h:458
Definition: Exceptions.h:88
virtual ~LevelSetFilter()
Destructor.
Definition: LevelSetFilter.h:101
Definition: Stencils.h:1452
GridType::Ptr meanCurvature(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the mean curvature of the given grid.
Definition: GridOperators.h:1033
LevelSetFilter(const LevelSetFilter &other)
Shallow copy constructor called by tbb::parallel_for() threads during filtering.
Definition: LevelSetFilter.h:91
GridType::TreeType TreeType
Definition: LevelSetFilter.h:70
#define OPENVDB_VERSION_NAME
Definition: version.h:45
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:137
OPENVDB_API Hermite max(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
TreeType::ValueType ValueType
Definition: LevelSetFilter.h:71
void operator()(const RangeType &range) const
Used internally by tbb::parallel_for().
Definition: LevelSetFilter.h:106
int32_t Int32
Definition: Types.h:59
tree::LeafManager< TreeType >::LeafRange RangeType
Definition: LevelSetFilter.h:73
LevelSetFilter(GridType &grid, InterruptT *interrupt=NULL)
Main constructor from a grid.
Definition: LevelSetFilter.h:79
AlphaType maxMask() const
Return the maximum value of the mask to be used for the derivation of a smooth alpha value...
Definition: LevelSetFilter.h:117
Filtering (e.g. diffusion) of narrow-band level sets. An optional scalar field can be used to produce...
Definition: LevelSetFilter.h:64
Axis
Definition: Math.h:767
AlphaType minMask() const
Return the minimum value of the mask to be used for the derivation of a smooth alpha value...
Definition: LevelSetFilter.h:114
CopyConstness< TreeType, NonConstBufferType >::Type BufferType
Definition: LeafManager.h:117
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67
Performs multi-threaded interface tracking of narrow band level sets.
Definition: LevelSetTracker.h:66
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:566