OpenVDB  2.1.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  : mSampler(mask, grid), mMin(min), mInvNorm(1/(max-min)), mInvert(invert)
206  {
207  assert(min < max);
208  }
209  inline bool operator()(const Coord& xyz, AlphaType& a, AlphaType& b) const
210  {
211  a = mSampler(xyz);
212  const AlphaType t = (a-mMin)*mInvNorm;
213  a = t > 0 ? t < 1 ? (3-2*t)*t*t : 1 : 0;//smooth mapping to 0->1
214  b = 1 - a;
215  if (mInvert) std::swap(a,b);
216  return a>0;
217  }
218  tools::DualGridSampler<MaskType, GridType, tools::BoxSampler> mSampler;
219  const AlphaType mMin, mInvNorm;
220  const bool mInvert;
221  };
222 
223  // Private driver method for mean and gaussian filtering
224  void box(int width);
225 
226  template <size_t Axis>
227  struct Avg {
228  Avg(const GridT& grid, Int32 w) :
229  acc(grid.tree()), width(w), frac(1/ValueType(2*w+1)) {}
230  ValueType operator()(Coord xyz) {
231  ValueType sum = zeroVal<ValueType>();
232  Int32& i = xyz[Axis], j = i + width;
233  for (i -= width; i <= j; ++i) sum += acc.getValue(xyz);
234  return sum*frac;
235  }
236  typename GridT::ConstAccessor acc;
237  const Int32 width;
238  const ValueType frac;
239  };
240 
241  // Private methods called by tbb::parallel_for threads
242  template <typename AvgT>
243  void doBox( const RangeType& r, Int32 w);
244  void doBoxX(const RangeType& r, Int32 w) { this->doBox<Avg<0> >(r,w); }
245  void doBoxZ(const RangeType& r, Int32 w) { this->doBox<Avg<1> >(r,w); }
246  void doBoxY(const RangeType& r, Int32 w) { this->doBox<Avg<2> >(r,w); }
247  void doMedian(const RangeType&, int);
248  void doMeanCurvature(const RangeType&);
249  void doLaplacian(const RangeType&);
250  void doOffset(const RangeType&, ValueType);
251 
252 }; // end of LevelSetFilter class
253 
254 
256 
257 template<typename GridT, typename MaskT, typename InterruptT>
258 inline void
260 {
261  mMask = mask;
262 
263  BaseType::startInterrupter("Median-value flow of level set");
264 
265  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
266 
267  mTask = boost::bind(&LevelSetFilter::doMedian, _1, _2, std::max(1, width));
268  this->cook(true);
269 
270  BaseType::track();
271 
272  BaseType::endInterrupter();
273 }
274 
275 template<typename GridT, typename MaskT, typename InterruptT>
276 inline void
278 {
279  mMask = mask;
280 
281  BaseType::startInterrupter("Mean-value flow of level set");
282 
283  this->box(width);
284 
285  BaseType::endInterrupter();
286 }
287 
288 template<typename GridT, typename MaskT, typename InterruptT>
289 inline void
291 {
292  mMask = mask;
293 
294  BaseType::startInterrupter("Gaussian flow of level set");
295 
296  for (int n=0; n<4; ++n) this->box(width);
297 
298  BaseType::endInterrupter();
299 }
300 
301 template<typename GridT, typename MaskT, typename InterruptT>
302 inline void
304 {
305  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
306 
307  width = std::max(1, width);
308 
309  mTask = boost::bind(&LevelSetFilter::doBoxX, _1, _2, width);
310  this->cook(true);
311 
312  mTask = boost::bind(&LevelSetFilter::doBoxY, _1, _2, width);
313  this->cook(true);
314 
315  mTask = boost::bind(&LevelSetFilter::doBoxZ, _1, _2, width);
316  this->cook(true);
317 
318  BaseType::track();
319 }
320 
321 template<typename GridT, typename MaskT, typename InterruptT>
322 inline void
324 {
325  mMask = mask;
326 
327  BaseType::startInterrupter("Mean-curvature flow of level set");
328 
329  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
330 
331  mTask = boost::bind(&LevelSetFilter::doMeanCurvature, _1, _2);
332  this->cook(true);
333 
334  BaseType::track();
335 
336  BaseType::endInterrupter();
337 }
338 
339 template<typename GridT, typename MaskT, typename InterruptT>
340 inline void
342 {
343  mMask = mask;
344 
345  BaseType::startInterrupter("Laplacian flow of level set");
346 
347  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
348 
349  mTask = boost::bind(&LevelSetFilter::doLaplacian, _1, _2);
350  this->cook(true);
351 
352  BaseType::track();
353 
354  BaseType::endInterrupter();
355 }
356 
357 template<typename GridT, typename MaskT, typename InterruptT>
358 inline void
360 {
361  mMask = mask;
362 
363  BaseType::startInterrupter("Offsetting level set");
364 
365  BaseType::leafs().removeAuxBuffers();// no auxiliary buffers required
366 
367  const ValueType CFL = ValueType(0.5) * BaseType::voxelSize(), offset = openvdb::math::Abs(value);
368  ValueType dist = 0.0;
369  while (offset-dist > ValueType(0.001)*CFL && BaseType::checkInterrupter()) {
370  const ValueType delta = openvdb::math::Min(offset-dist, CFL);
371  dist += delta;
372 
373  mTask = boost::bind(&LevelSetFilter::doOffset, _1, _2, copysign(delta, value));
374  this->cook(false);
375 
376  BaseType::track();
377  }
378 
379  BaseType::endInterrupter();
380 }
381 
382 
384 
386 template<typename GridT, typename MaskT, typename InterruptT>
387 inline void
389 {
390  BaseType::checkInterrupter();
391  //const float CFL = 0.9f, dt = CFL * mDx * mDx / 6.0f;
392  const ValueType dx = BaseType::voxelSize(), dt = math::Pow2(dx) / ValueType(3.0);
393  math::CurvatureStencil<GridType> stencil(BaseType::grid(), dx);
394  if (mMask) {
395  AlphaType a, b;
396  AlphaMask alpha(BaseType::grid(), *mMask, mMinMask, mMaxMask, mInvertMask);
397  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
398  BufferT& buffer = leafIter.buffer(1);
399  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
400  if (alpha(iter.getCoord(), a, b)) {
401  stencil.moveTo(iter);
402  const ValueType phi0 = *iter, phi1 = phi0 + dt*stencil.meanCurvatureNormGrad();
403  buffer.setValue(iter.pos(), b*phi0 + a*phi1);
404  }
405  }
406  }
407  } else {
408  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
409  BufferT& buffer = leafIter.buffer(1);
410  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
411  stencil.moveTo(iter);
412  buffer.setValue(iter.pos(), *iter + dt*stencil.meanCurvatureNormGrad());
413  }
414  }
415  }
416 }
417 
425 template<typename GridT, typename MaskT, typename InterruptT>
426 inline void
427 LevelSetFilter<GridT, MaskT, InterruptT>::doLaplacian(const RangeType& range)
428 {
429  BaseType::checkInterrupter();
430  //const float CFL = 0.9f, half_dt = CFL * mDx * mDx / 12.0f;
431  const ValueType dx = BaseType::voxelSize(), dt = math::Pow2(dx) / ValueType(6.0);
432  math::GradStencil<GridType> stencil(BaseType::grid(), dx);
433  if (mMask) {
434  AlphaType a, b;
435  AlphaMask alpha(BaseType::grid(), *mMask, mMinMask, mMaxMask, mInvertMask);
436  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
437  BufferT& buffer = leafIter.buffer(1);
438  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
439  if (alpha(iter.getCoord(), a, b)) {
440  stencil.moveTo(iter);
441  const ValueType phi0 = *iter, phi1 = phi0 + dt*stencil.laplacian();
442  buffer.setValue(iter.pos(), b*phi0 + a*phi1);
443  }
444  }
445  }
446  } else {
447  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
448  BufferT& buffer = leafIter.buffer(1);
449  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
450  stencil.moveTo(iter);
451  buffer.setValue(iter.pos(), *iter + dt*stencil.laplacian());
452  }
453  }
454  }
455 }
456 
458 template<typename GridT, typename MaskT, typename InterruptT>
459 inline void
460 LevelSetFilter<GridT, MaskT, InterruptT>::doOffset(const RangeType& range, ValueType offset)
461 {
462  BaseType::checkInterrupter();
463  if (mMask) {
464  AlphaType a, b;
465  AlphaMask alpha(BaseType::grid(), *mMask, mMinMask, mMaxMask, mInvertMask);
466  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
467  for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
468  if (alpha(iter.getCoord(), a, b)) iter.setValue(*iter + a*offset);
469  }
470  }
471  } else {
472  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
473  for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
474  iter.setValue(*iter + offset);
475  }
476  }
477  }
478 }
479 
481 template<typename GridT, typename MaskT, typename InterruptT>
482 inline void
483 LevelSetFilter<GridT, MaskT, InterruptT>::doMedian(const RangeType& range, int width)
484 {
485  BaseType::checkInterrupter();
486  typename math::DenseStencil<GridType> stencil(BaseType::grid(), width);//creates local cache!
487  if (mMask) {
488  AlphaType a, b;
489  AlphaMask alpha(BaseType::grid(), *mMask, mMinMask, mMaxMask, mInvertMask);
490  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
491  BufferT& buffer = leafIter.buffer(1);
492  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
493  if (alpha(iter.getCoord(), a, b)) {
494  stencil.moveTo(iter);
495  buffer.setValue(iter.pos(), b*(*iter) + a*stencil.median());
496  }
497  }
498  }
499  } else {
500  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
501  BufferT& buffer = leafIter.buffer(1);
502  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
503  stencil.moveTo(iter);
504  buffer.setValue(iter.pos(), stencil.median());
505  }
506  }
507  }
508 }
509 
511 template<typename GridT, typename MaskT, typename InterruptT>
512 template <typename AvgT>
513 inline void
514 LevelSetFilter<GridT, MaskT, InterruptT>::doBox(const RangeType& range, Int32 w)
515 {
516  BaseType::checkInterrupter();
517  AvgT avg(BaseType::grid(), w);
518  if (mMask) {
519  AlphaType a, b;
520  AlphaMask alpha(BaseType::grid(), *mMask, mMinMask, mMaxMask, mInvertMask);
521  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
522  BufferT& buffer = leafIter.buffer(1);
523  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
524  const Coord xyz = iter.getCoord();
525  if (alpha(xyz, a, b)) buffer.setValue(iter.pos(), b*(*iter)+ a*avg(xyz));
526  }
527  }
528  } else {
529  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
530  BufferT& buffer = leafIter.buffer(1);
531  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
532  buffer.setValue(iter.pos(), avg(iter.getCoord()));
533  }
534  }
535  }
536 }
537 
538 } // namespace tools
539 } // namespace OPENVDB_VERSION_NAME
540 } // namespace openvdb
541 
542 #endif // OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
543 
544 // Copyright (c) 2012-2013 DreamWorks Animation LLC
545 // All rights reserved. This software is distributed under the
546 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
void setMaskRange(AlphaType min, AlphaType max)
Define the range for the (optional) scalar mask.
Definition: LevelSetFilter.h:125
bool isMaskInverted() const
Return true if the mask is inverted, i.e. min-&gt;max in the original mask maps to 1-&gt;0 in the inverted ...
Definition: LevelSetFilter.h:134
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
OPENVDB_API Hermite min(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
virtual ~LevelSetFilter()
Destructor.
Definition: LevelSetFilter.h:101
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Definition: Exceptions.h:88
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:253
CopyConstness< TreeType, NonConstBufferType >::Type BufferType
Definition: LeafManager.h:117
Type Pow2(Type x)
Return .
Definition: Math.h:467
void invertMask(bool invert=true)
Invert the optional mask, i.e. min-&gt;max in the original mask maps to 1-&gt;0 in the inverted alpha mask...
Definition: LevelSetFilter.h:137
Filtering (e.g. diffusion) of narrow-band level sets. An optional scalar field can be used to produce...
Definition: LevelSetFilter.h:64
MaskT MaskType
Definition: LevelSetFilter.h:69
void operator()(const RangeType &range) const
Used internally by tbb::parallel_for().
Definition: LevelSetFilter.h:106
#define OPENVDB_VERSION_NAME
Definition: version.h:45
tree::LeafManager< TreeType >::LeafRange RangeType
Definition: LevelSetFilter.h:73
Definition: Stencils.h:1452
LevelSetFilter(GridType &grid, InterruptT *interrupt=NULL)
Main constructor from a grid.
Definition: LevelSetFilter.h:79
const MaskGridType * mMask
Definition: GridOperators.h:386
int32_t Int32
Definition: Types.h:58
GridType::Ptr meanCurvature(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the mean curvature of the given grid.
Definition: GridOperators.h:932
OPENVDB_API Hermite max(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
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
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
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:575
TreeType::ValueType ValueType
Definition: LevelSetFilter.h:71
Performs multi-threaded interface tracking of narrow band level sets.
Definition: LevelSetTracker.h:66
GridT GridType
Definition: LevelSetFilter.h:68
Axis
Definition: Math.h:769
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67
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
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:916