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