OpenVDB  1.2.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  ~AffineCombine() { delete mAcc; }
160  inline ValueType alpha(const Coord& xyz) const
161  {
162  assert(mMask);
163  const Vec3R world = mGrid->indexToWorld(xyz);
164  const Vec3R voxel = mMask->worldToIndex(world);
165  return tools::BoxSampler::sample<AccT>(*mAcc, voxel);
166  }
167  inline ValueType operator()(const Coord& xyz, ValueType unfiltered, ValueType filtered) const
168  {
169  if (mMask==NULL) return filtered;
170  const ValueType alpha = this->alpha(xyz);
171  return alpha*filtered + (1-alpha)*unfiltered;
172  }
173  const GridT* mGrid;
174  const GridT* mMask;
175  const AccT* mAcc;
176  };
177 
178  // Private driver method for mean and gaussian filtering
179  void box(int width);
180 
181  template <size_t Axis>
182  struct Avg {
183  Avg(const GridT& grid, Int32 w) : acc(grid.tree()), width(w), frac(1/ValueType(2*w+1)) {}
184  ValueType operator()(Coord xyz) {
185  ValueType sum = zeroVal<ValueType>();
186  Int32& i = xyz[Axis], j = i + width;
187  for (i -= width; i <= j; ++i) sum += acc.getValue(xyz);
188  return sum*frac;
189  }
190  typename GridT::ConstAccessor acc;
191  const Int32 width;
192  const ValueType frac;
193  };
194 
195  // Private methods called by tbb::parallel_for threads
196  template <typename AvgT>
197  void doBox( const RangeType& r, Int32 w);
198  void doBoxX(const RangeType& r, Int32 w) { this->doBox<Avg<0> >(r,w); }
199  void doBoxZ(const RangeType& r, Int32 w) { this->doBox<Avg<1> >(r,w); }
200  void doBoxY(const RangeType& r, Int32 w) { this->doBox<Avg<2> >(r,w); }
201  void doMedian(const RangeType&, int);
202  void doMeanCurvature(const RangeType&);
203  void doLaplacian(const RangeType&);
204  void doOffset(const RangeType&, ValueType);
205 
206 }; // end of LevelSetFilter class
207 
208 
210 
211 template<typename GridT, typename InterruptT>
212 inline void
213 LevelSetFilter<GridT, InterruptT>::median(int width, const GridT* mask)
214 {
215  mMask = mask;
216 
217  BaseType::startInterrupter("Median-value flow of level set");
218 
219  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
220 
221  mTask = boost::bind(&LevelSetFilter::doMedian, _1, _2, std::max(1, width));
222  this->cook(true);
223 
224  BaseType::track();
225 
226  BaseType::endInterrupter();
227 }
228 
229 
230 template<typename GridT, typename InterruptT>
231 inline void
232 LevelSetFilter<GridT, InterruptT>::mean(int width, const GridT* mask)
233 {
234  mMask = mask;
235 
236  BaseType::startInterrupter("Mean-value flow of level set");
237 
238  this->box(width);
239 
240  BaseType::endInterrupter();
241 }
242 
243 template<typename GridT, typename InterruptT>
244 inline void
245 LevelSetFilter<GridT, InterruptT>::gaussian(int width, const GridT* mask)
246 {
247  mMask = mask;
248 
249  BaseType::startInterrupter("Gaussian flow of level set");
250 
251  for (int n=0; n<4; ++n) this->box(width);
252 
253  BaseType::endInterrupter();
254 }
255 
256 template<typename GridT, typename InterruptT>
257 inline void
259 {
260  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
261 
262  width = std::max(1, width);
263 
264  mTask = boost::bind(&LevelSetFilter::doBoxX, _1, _2, width);
265  this->cook(true);
266 
267  mTask = boost::bind(&LevelSetFilter::doBoxY, _1, _2, width);
268  this->cook(true);
269 
270  mTask = boost::bind(&LevelSetFilter::doBoxZ, _1, _2, width);
271  this->cook(true);
272 
273  BaseType::track();
274 }
275 
276 
277 template<typename GridT, typename InterruptT>
278 inline void
280 {
281  mMask = mask;
282 
283  BaseType::startInterrupter("Mean-curvature flow of level set");
284 
285  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
286 
287  mTask = boost::bind(&LevelSetFilter::doMeanCurvature, _1, _2);
288  this->cook(true);
289 
290  BaseType::track();
291 
292  BaseType::endInterrupter();
293 }
294 
295 template<typename GridT, typename InterruptT>
296 inline void
298 {
299  mMask = mask;
300 
301  BaseType::startInterrupter("Laplacian flow of level set");
302 
303  BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
304 
305  mTask = boost::bind(&LevelSetFilter::doLaplacian, _1, _2);
306  this->cook(true);
307 
308  BaseType::track();
309 
310  BaseType::endInterrupter();
311 }
312 
313 
314 template<typename GridT, typename InterruptT>
315 inline void
317 {
318  mMask = mask;
319 
320  BaseType::startInterrupter("Offsetting level set");
321 
322  BaseType::leafs().removeAuxBuffers();// no auxiliary buffers required
323 
324  const ValueType CFL = ValueType(0.5) * BaseType::voxelSize(), offset = openvdb::math::Abs(value);
325  ValueType dist = 0.0;
326  while (offset-dist > ValueType(0.001)*CFL && BaseType::checkInterrupter()) {
327  const ValueType delta = openvdb::math::Min(offset-dist, CFL);
328  dist += delta;
329 
330  mTask = boost::bind(&LevelSetFilter::doOffset, _1, _2, copysign(delta, value));
331  this->cook(false);
332 
333  BaseType::track();
334  }
335 
336  BaseType::endInterrupter();
337 }
338 
339 
341 
343 template<typename GridT, typename InterruptT>
344 inline void
346 {
347  BaseType::checkInterrupter();
348  //const float CFL = 0.9f, dt = CFL * mDx * mDx / 6.0f;
349  const ValueType dx = BaseType::voxelSize(),dt = math::Pow2(dx) / ValueType(3.0);
350  math::CurvatureStencil<GridType> stencil(BaseType::grid(), dx);
351  AffineCombine comb(BaseType::grid(), mMask);
352  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
353  BufferT& buffer = leafIter.buffer(1);
354  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
355  stencil.moveTo(iter);
356  ValueType A = stencil.getCenterValue(), B = A + dt * stencil.meanCurvatureNormGrad();
357  buffer.setValue(iter.pos(), comb(stencil.getCenterCoord(), A, B));
358  }
359  }
360 }
361 
369 template<typename GridT, typename InterruptT>
370 inline void
371 LevelSetFilter<GridT, InterruptT>::doLaplacian(const RangeType& range)
372 {
373  BaseType::checkInterrupter();
374  //const float CFL = 0.9f, half_dt = CFL * mDx * mDx / 12.0f;
375  const ValueType dx = BaseType::voxelSize(), dt = math::Pow2(dx) / ValueType(6.0);
376  math::GradStencil<GridType> stencil(BaseType::grid(), dx);
377  AffineCombine comb(BaseType::grid(), mMask);
378  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
379  BufferT& buffer = leafIter.buffer(1);
380  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
381  stencil.moveTo(iter);
382  const ValueType A = stencil.getCenterValue(), B = A + dt * stencil.laplacian();
383  buffer.setValue(iter.pos(), comb(stencil.getCenterCoord(), A, B));
384  }
385  }
386 }
387 
389 template<typename GridT, typename InterruptT>
390 inline void
391 LevelSetFilter<GridT, InterruptT>::doOffset(const RangeType& range, ValueType offset)
392 {
393  BaseType::checkInterrupter();
394  if (mMask) {
395  AffineCombine comb(BaseType::grid(), mMask);
396  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
397  for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
398  iter.setValue(*iter + comb.alpha(iter.getCoord())*offset);
399  }
400  }
401  } else {
402  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
403  for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
404  iter.setValue(*iter + offset);
405  }
406  }
407  }
408 }
409 
411 template<typename GridT, typename InterruptT>
412 inline void
413 LevelSetFilter<GridT, InterruptT>::doMedian(const RangeType& range, int width)
414 {
415  BaseType::checkInterrupter();
416  typename math::DenseStencil<GridType> stencil(BaseType::grid(), width);//creates local cache!
417  AffineCombine comb(BaseType::grid(), mMask);
418  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
419  BufferT& buffer = leafIter.buffer(1);
420  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
421  stencil.moveTo(iter);
422  const ValueType A = stencil.getCenterValue(), B = stencil.median();
423  buffer.setValue(iter.pos(), comb(stencil.getCenterCoord(), A, B));
424  }
425  }
426 }
427 
429 template<typename GridT, typename InterruptT>
430 template <typename AvgT>
431 inline void
432 LevelSetFilter<GridT, InterruptT>::doBox(const RangeType& range, Int32 w)
433 {
434  BaseType::checkInterrupter();
435  AvgT avg(BaseType::grid(), w);
436  AffineCombine comb(BaseType::grid(), mMask);
437  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
438  BufferT& buffer = leafIter.buffer(1);
439  for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
440  const Coord xyz = iter.getCoord();
441  buffer.setValue(iter.pos(), comb(xyz, *iter, avg(xyz)));
442  }
443  }
444 }
445 
446 } // namespace tools
447 } // namespace OPENVDB_VERSION_NAME
448 } // namespace openvdb
449 
450 #endif // OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
451 
452 // Copyright (c) 2012-2013 DreamWorks Animation LLC
453 // All rights reserved. This software is distributed under the
454 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )