OpenVDB  2.0.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 "LevelSetTracker.h"
46 #include "Interpolation.h"
47 
48 namespace openvdb {
50 namespace OPENVDB_VERSION_NAME {
51 namespace tools {
52 
57 template<typename GridT,
58  typename InterruptT = util::NullInterrupter>
59 class LevelSetFilter : public LevelSetTracker<GridT, InterruptT>
60 {
61 public:
63  typedef GridT GridType;
64  typedef typename GridType::TreeType TreeType;
65  typedef typename TreeType::ValueType ValueType;
67 
71  LevelSetFilter(GridType& grid, InterruptT* interrupt = NULL)
72  : BaseType(grid, interrupt), mTask(0), mMask(NULL)
73  {
74  }
79  : BaseType(other), mTask(other.mTask), mMask(other.mMask)
80  {
81  }
83  virtual ~LevelSetFilter() {};
84 
88  void operator()(const RangeType& range) const
89  {
90  if (mTask) mTask(const_cast<LevelSetFilter*>(this), range);
91  else OPENVDB_THROW(ValueError, "task is undefined - call offset(), etc");
92  }
93 
96  void meanCurvature(const GridT* mask = NULL);
97 
100  void laplacian(const GridT* mask = NULL);
101 
108  void gaussian(int width = 1, const GridT* mask = NULL);
109 
113  void offset(ValueType offset, const GridT* mask = NULL);
114 
121  void median(int width = 1, const GridT* mask = NULL);
122 
128  void mean(int width = 1, const GridT* mask = NULL);
129 
130 private:
131  typedef typename GridT::ConstAccessor AccT;
132  typedef typename TreeType::LeafNodeType LeafT;
133  typedef typename LeafT::ValueOnIter VoxelIterT;
134  typedef typename LeafT::ValueOnCIter VoxelCIterT;
135  typedef typename tree::LeafManager<TreeType>::BufferType BufferT;
136  typedef typename RangeType::Iterator LeafIterT;
137 
138  // Only two private member data
139  typename boost::function<void (LevelSetFilter*, const RangeType&)> mTask;
140  const GridT* mMask;
141 
142  // Private cook method calling tbb::parallel_for
143  void cook(bool swap)
144  {
145  const int n = BaseType::getGrainSize();
146  if (n>0) {
147  tbb::parallel_for(BaseType::leafs().leafRange(n), *this);
148  } else {
149  (*this)(BaseType::leafs().leafRange());
150  }
151  if (swap) BaseType::leafs().swapLeafBuffer(1, n==0);
152  }
153 
154  // Private class to perform interpolated alpha-blending
155  struct AffineCombine
156  {
157  AffineCombine(const GridT& grid, const GridT* mask)
158  : mGrid(&grid), mMask(mask), mAcc(mask ? new AccT(mask->tree()) : NULL),
159  mIdentical(mask ? grid.transform() == mask->transform() : false)
160  {
161  }
162  ~AffineCombine() { delete mAcc; }
163  inline ValueType alpha(const Coord& xyz) const
164  {
165  assert(mMask);
166  if (mIdentical) return mAcc->getValue(xyz);
167  const Vec3R world = mGrid->indexToWorld(xyz);
168  const Vec3R voxel = mMask->worldToIndex(world);
169  return tools::BoxSampler::sample<AccT>(*mAcc, voxel);
170  }
171  inline ValueType operator()(const Coord& xyz, ValueType unfiltered, ValueType filtered) const
172  {
173  if (mMask==NULL) return filtered;
174  const ValueType alpha = this->alpha(xyz);
175  return alpha*filtered + (1-alpha)*unfiltered;
176  }
177  const GridT* mGrid;
178  const GridT* mMask;
179  const AccT* mAcc;
180  const bool mIdentical;
181  };
182 
183  // Private driver method for mean and gaussian filtering
184  void box(int width);
185 
186  template <size_t Axis>
187  struct Avg {
188  Avg(const GridT& grid, Int32 w) : acc(grid.tree()), width(w), frac(1/ValueType(2*w+1)) {}
189  ValueType operator()(Coord xyz) {
190  ValueType sum = zeroVal<ValueType>();
191  Int32& i = xyz[Axis], j = i + width;
192  for (i -= width; i <= j; ++i) sum += acc.getValue(xyz);
193  return sum*frac;
194  }
195  typename GridT::ConstAccessor acc;
196  const Int32 width;
197  const ValueType frac;
198  };
199 
200  // Private methods called by tbb::parallel_for threads
201  template <typename AvgT>
202  void doBox( const RangeType& r, Int32 w);
203  void doBoxX(const RangeType& r, Int32 w) { this->doBox<Avg<0> >(r,w); }
204  void doBoxZ(const RangeType& r, Int32 w) { this->doBox<Avg<1> >(r,w); }
205  void doBoxY(const RangeType& r, Int32 w) { this->doBox<Avg<2> >(r,w); }
206  void doMedian(const RangeType&, int);
207  void doMeanCurvature(const RangeType&);
208  void doLaplacian(const RangeType&);
209  void doOffset(const RangeType&, ValueType);
210 
211 }; // end of LevelSetFilter class
212 
213 
215 
216 template<typename GridT, typename InterruptT>
217 inline void
218 LevelSetFilter<GridT, InterruptT>::median(int width, const GridT* mask)
219 {
220  mMask = mask;
221 
222  BaseType::startInterrupter("Median-value flow of level set");
223 
224  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
225 
226  mTask = boost::bind(&LevelSetFilter::doMedian, _1, _2, std::max(1, width));
227  this->cook(true);
228 
229  BaseType::track();
230 
231  BaseType::endInterrupter();
232 }
233 
234 
235 template<typename GridT, typename InterruptT>
236 inline void
237 LevelSetFilter<GridT, InterruptT>::mean(int width, const GridT* mask)
238 {
239  mMask = mask;
240 
241  BaseType::startInterrupter("Mean-value flow of level set");
242 
243  this->box(width);
244 
245  BaseType::endInterrupter();
246 }
247 
248 template<typename GridT, typename InterruptT>
249 inline void
250 LevelSetFilter<GridT, InterruptT>::gaussian(int width, const GridT* mask)
251 {
252  mMask = mask;
253 
254  BaseType::startInterrupter("Gaussian flow of level set");
255 
256  for (int n=0; n<4; ++n) this->box(width);
257 
258  BaseType::endInterrupter();
259 }
260 
261 template<typename GridT, typename InterruptT>
262 inline void
264 {
265  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
266 
267  width = std::max(1, width);
268 
269  mTask = boost::bind(&LevelSetFilter::doBoxX, _1, _2, width);
270  this->cook(true);
271 
272  mTask = boost::bind(&LevelSetFilter::doBoxY, _1, _2, width);
273  this->cook(true);
274 
275  mTask = boost::bind(&LevelSetFilter::doBoxZ, _1, _2, width);
276  this->cook(true);
277 
278  BaseType::track();
279 }
280 
281 
282 template<typename GridT, typename InterruptT>
283 inline void
285 {
286  mMask = mask;
287 
288  BaseType::startInterrupter("Mean-curvature flow of level set");
289 
290  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
291 
292  mTask = boost::bind(&LevelSetFilter::doMeanCurvature, _1, _2);
293  this->cook(true);
294 
295  BaseType::track();
296 
297  BaseType::endInterrupter();
298 }
299 
300 template<typename GridT, typename InterruptT>
301 inline void
303 {
304  mMask = mask;
305 
306  BaseType::startInterrupter("Laplacian flow of level set");
307 
308  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
309 
310  mTask = boost::bind(&LevelSetFilter::doLaplacian, _1, _2);
311  this->cook(true);
312 
313  BaseType::track();
314 
315  BaseType::endInterrupter();
316 }
317 
318 
319 template<typename GridT, typename InterruptT>
320 inline void
322 {
323  mMask = mask;
324 
325  BaseType::startInterrupter("Offsetting level set");
326 
327  BaseType::leafs().removeAuxBuffers();// no auxiliary buffers required
328 
329  const ValueType CFL = ValueType(0.5) * BaseType::voxelSize(), offset = openvdb::math::Abs(value);
330  ValueType dist = 0.0;
331  while (offset-dist > ValueType(0.001)*CFL && BaseType::checkInterrupter()) {
332  const ValueType delta = openvdb::math::Min(offset-dist, CFL);
333  dist += delta;
334 
335  mTask = boost::bind(&LevelSetFilter::doOffset, _1, _2, copysign(delta, value));
336  this->cook(false);
337 
338  BaseType::track();
339  }
340 
341  BaseType::endInterrupter();
342 }
343 
344 
346 
348 template<typename GridT, typename InterruptT>
349 inline void
351 {
352  BaseType::checkInterrupter();
353  //const float CFL = 0.9f, dt = CFL * mDx * mDx / 6.0f;
354  const ValueType dx = BaseType::voxelSize(),dt = math::Pow2(dx) / ValueType(3.0);
355  math::CurvatureStencil<GridType> stencil(BaseType::grid(), dx);
356  AffineCombine comb(BaseType::grid(), mMask);
357  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
358  BufferT& buffer = leafIter.buffer(1);
359  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
360  stencil.moveTo(iter);
361  ValueType A = stencil.getCenterValue(), B = A + dt * stencil.meanCurvatureNormGrad();
362  buffer.setValue(iter.pos(), comb(stencil.getCenterCoord(), A, B));
363  }
364  }
365 }
366 
374 template<typename GridT, typename InterruptT>
375 inline void
376 LevelSetFilter<GridT, InterruptT>::doLaplacian(const RangeType& range)
377 {
378  BaseType::checkInterrupter();
379  //const float CFL = 0.9f, half_dt = CFL * mDx * mDx / 12.0f;
380  const ValueType dx = BaseType::voxelSize(), dt = math::Pow2(dx) / ValueType(6.0);
381  math::GradStencil<GridType> stencil(BaseType::grid(), dx);
382  AffineCombine comb(BaseType::grid(), mMask);
383  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
384  BufferT& buffer = leafIter.buffer(1);
385  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
386  stencil.moveTo(iter);
387  const ValueType A = stencil.getCenterValue(), B = A + dt * stencil.laplacian();
388  buffer.setValue(iter.pos(), comb(stencil.getCenterCoord(), A, B));
389  }
390  }
391 }
392 
394 template<typename GridT, typename InterruptT>
395 inline void
396 LevelSetFilter<GridT, InterruptT>::doOffset(const RangeType& range, ValueType offset)
397 {
398  BaseType::checkInterrupter();
399  if (mMask) {
400  AffineCombine comb(BaseType::grid(), mMask);
401  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
402  for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
403  iter.setValue(*iter + comb.alpha(iter.getCoord())*offset);
404  }
405  }
406  } else {
407  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
408  for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
409  iter.setValue(*iter + offset);
410  }
411  }
412  }
413 }
414 
416 template<typename GridT, typename InterruptT>
417 inline void
418 LevelSetFilter<GridT, InterruptT>::doMedian(const RangeType& range, int width)
419 {
420  BaseType::checkInterrupter();
421  typename math::DenseStencil<GridType> stencil(BaseType::grid(), width);//creates local cache!
422  AffineCombine comb(BaseType::grid(), mMask);
423  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
424  BufferT& buffer = leafIter.buffer(1);
425  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
426  stencil.moveTo(iter);
427  const ValueType A = stencil.getCenterValue(), B = stencil.median();
428  buffer.setValue(iter.pos(), comb(stencil.getCenterCoord(), A, B));
429  }
430  }
431 }
432 
434 template<typename GridT, typename InterruptT>
435 template <typename AvgT>
436 inline void
437 LevelSetFilter<GridT, InterruptT>::doBox(const RangeType& range, Int32 w)
438 {
439  BaseType::checkInterrupter();
440  AvgT avg(BaseType::grid(), w);
441  AffineCombine comb(BaseType::grid(), mMask);
442  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
443  BufferT& buffer = leafIter.buffer(1);
444  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
445  const Coord xyz = iter.getCoord();
446  buffer.setValue(iter.pos(), comb(xyz, *iter, avg(xyz)));
447  }
448  }
449 }
450 
451 } // namespace tools
452 } // namespace OPENVDB_VERSION_NAME
453 } // namespace openvdb
454 
455 #endif // OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
456 
457 // Copyright (c) 2012-2013 DreamWorks Animation LLC
458 // All rights reserved. This software is distributed under the
459 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
OPENVDB_API Hermite max(const Hermite &, const Hermite &)
min and max operations done directly on the compressed data.
Definition: Stencils.h:1452
GridType::TreeType TreeType
Definition: LevelSetFilter.h:64
tree::LeafManager< TreeType >::LeafRange RangeType
Definition: LevelSetFilter.h:66
Performs multi-threaded interface tracking of narrow band level sets. This is the building-block for ...
LevelSetFilter(GridType &grid, InterruptT *interrupt=NULL)
Main constructor from a grid.
Definition: LevelSetFilter.h:71
int32_t Int32
Definition: Types.h:58
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Type Pow2(Type x)
Return .
Definition: Math.h:460
LevelSetFilter(const LevelSetFilter &other)
Shallow copy constructor called by tbb::parallel_for() threads during filtering.
Definition: LevelSetFilter.h:78
TreeType::ValueType ValueType
Definition: LevelSetFilter.h:65
LevelSetTracker< GridT, InterruptT > BaseType
Definition: LevelSetFilter.h:62
AccessorT mAcc
Definition: GridOperators.h:259
GridT GridType
Definition: LevelSetFilter.h:63
Performs multi-threaded interface tracking of narrow band level sets.
Definition: LevelSetTracker.h:64
#define OPENVDB_VERSION_NAME
Definition: version.h:45
math::Vec3< Real > Vec3R
Definition: Types.h:74
Filtering (i.e. diffusion) of narrow-band level sets.
Definition: LevelSetFilter.h:59
CopyConstness< TreeType, NonConstBufferType >::Type BufferType
Definition: LeafManager.h:117
GridType::Ptr laplacian(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the Laplacian of the given scalar grid.
Definition: GridOperators.h:700
virtual ~LevelSetFilter()
Destructor.
Definition: LevelSetFilter.h:83
Definition: Exceptions.h:88
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:246
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:56
Axis
Definition: Math.h:762
const Type & Min(const Type &a, const Type &b)
Return the minimum of two values.
Definition: Math.h:568
void operator()(const RangeType &range) const
Used internally by tbb::parallel_for().
Definition: LevelSetFilter.h:88
GridType::Ptr meanCurvature(const GridType &grid, bool threaded, InterruptT *interrupt)
Compute the mean curvature of the given grid.
Definition: GridOperators.h:708