OpenVDB  1.1.0
LevelSetFilter.h
Go to the documentation of this file.
1 
2 //
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 //
39 
40 #ifndef OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
41 #define OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
42 
43 #include "LevelSetTracker.h"
44 
45 namespace openvdb {
47 namespace OPENVDB_VERSION_NAME {
48 namespace tools {
49 
54 template<typename GridT,
55  typename InterruptT = util::NullInterrupter>
56 class LevelSetFilter : public LevelSetTracker<GridT, InterruptT>
57 {
58 public:
59  typedef boost::shared_ptr<LevelSetFilter> Ptr;
61  typedef GridT GridType;
62  typedef typename GridType::TreeType TreeType;
63  typedef typename TreeType::LeafNodeType LeafType;
64  typedef typename LeafType::ValueType ValueType;
68 
70  LevelSetFilter(GridType& grid, InterruptT* interrupt = NULL)
71  : BaseType(grid, interrupt), mTask(0)
72  {
73  }
76  : BaseType(other), mTask(other.mTask)
77  {
78  }
79  virtual ~LevelSetFilter() {};
80 
82  void operator()(const RangeType& r) const
83  {
84  if (mTask) mTask(const_cast<LevelSetFilter*>(this), r);
85  else OPENVDB_THROW(ValueError, "task is undefined - call offset(), etc");
86  }
87 
89  void meanCurvature();
90 
92  void laplacian();
93 
98  void gaussian(int width = 1);
99 
101  void offset(ValueType offset);
102 
106  void median(int width = 1);
107 
111  void mean(int width = 1);
112 
113 private:
114  typedef typename boost::function<void (LevelSetFilter*, const RangeType&)> FuncType;
115 
116  FuncType mTask;
117 
118  // Private cook method calling tbb::parallel_for
119  void cook(bool swap)
120  {
121  const int n = BaseType::getGrainSize();
122  if (n>0) {
123  tbb::parallel_for(BaseType::mLeafs->getRange(n), *this);
124  } else {
125  (*this)(BaseType::mLeafs->getRange());
126  }
127  if (swap) BaseType::mLeafs->swapLeafBuffer(1, n==0);
128  }
129 
130  // Private methods called by tbb::parallel_for threads
131  void doBoxX(const RangeType&, Int32);
132  void doBoxY(const RangeType&, Int32);
133  void doBoxZ(const RangeType&, Int32);
134  void doMedian(const RangeType&, int);
135  void doMeanCurvature(const RangeType&);
136  void doLaplacian(const RangeType&);
137  void doOffset(const RangeType&, ValueType);
138 
139 }; // end of LevelSetFilter class
140 
141 
143 
144 template<typename GridT, typename InterruptT>
145 inline void
147 {
148  BaseType::startInterrupter("Median-value flow of level set");
149 
150  BaseType::mLeafs->rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
151 
152  mTask = boost::bind(&LevelSetFilter::doMedian, _1, _2, std::max(1, width));
153  this->cook(true);
154 
155  BaseType::track();
156 
157  BaseType::endInterrupter();
158 }
159 
160 
161 template<typename GridT, typename InterruptT>
162 inline void
164 {
165  BaseType::startInterrupter("Mean-value flow of level set");
166 
167  width = std::max(1, width);
168 
169  BaseType::mLeafs->rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
170 
171  mTask = boost::bind(&LevelSetFilter::doBoxX, _1, _2, width);
172  this->cook(true);
173 
174  mTask = boost::bind(&LevelSetFilter::doBoxY, _1, _2, width);
175  this->cook(true);
176 
177  mTask = boost::bind(&LevelSetFilter::doBoxZ, _1, _2, width);
178  this->cook(true);
179 
180  BaseType::track();
181 
182  BaseType::endInterrupter();
183 }
184 
185 template<typename GridT, typename InterruptT>
186 inline void
188 {
189  BaseType::startInterrupter("Gaussian flow of level set");
190 
191  width = std::max(1, width);
192 
193  BaseType::mLeafs->rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
194 
195  for (int n=0; n<4; ++n) {
196 
197  mTask = boost::bind(&LevelSetFilter::doBoxX, _1, _2, width);
198  this->cook(true);
199 
200  mTask = boost::bind(&LevelSetFilter::doBoxY, _1, _2, width);
201  this->cook(true);
202 
203  mTask = boost::bind(&LevelSetFilter::doBoxZ, _1, _2, width);
204  this->cook(true);
205  }
206 
207  BaseType::track();
208 
209  BaseType::endInterrupter();
210 }
211 
212 
213 
214 template<typename GridT, typename InterruptT>
215 inline void
217 {
218  BaseType::startInterrupter("Mean-curvature flow of level set");
219 
220  BaseType::mLeafs->rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
221 
222  mTask = boost::bind(&LevelSetFilter::doMeanCurvature, _1, _2);
223  this->cook(true);
224 
225  BaseType::track();
226 
227  BaseType::endInterrupter();
228 }
229 
230 template<typename GridT, typename InterruptT>
231 inline void
233 {
234  BaseType::startInterrupter("Laplacian flow of level set");
235 
236  BaseType::mLeafs->rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
237 
238  mTask = boost::bind(&LevelSetFilter::doLaplacian, _1, _2);
239  this->cook(true);
240 
241  BaseType::track();
242 
243  BaseType::endInterrupter();
244 }
245 
246 
247 template<typename GridT, typename InterruptT>
248 inline void
250 {
251  BaseType::startInterrupter("Offsetting level set");
252 
253  BaseType::mLeafs->removeAuxBuffers();// no auxiliary buffers required
254 
255  const ValueType CFL = ValueType(0.5) * BaseType::voxelSize(), offset = openvdb::math::Abs(value);
256  ValueType dist = 0.0;
257  while (offset-dist > ValueType(0.001)*CFL && BaseType::checkInterrupter()) {
258  const ValueType delta = openvdb::math::Min(offset-dist, CFL);
259  dist += delta;
260 
261  mTask = boost::bind(&LevelSetFilter::doOffset, _1, _2, copysign(delta,value));
262  this->cook(false);
263 
264  BaseType::track();
265  }
266 
267  BaseType::endInterrupter();
268 }
269 
270 
272 
274 template<typename GridT, typename InterruptT>
275 inline void
277 {
278  typedef typename LeafType::ValueOnCIter VoxelIterT;
279  BaseType::checkInterrupter();
280  //const float CFL = 0.9f, dt = CFL * mDx * mDx / 6.0f;
281  const ValueType dx = BaseType::voxelSize(),dt = math::Pow2(dx) / ValueType(3.0);
282  math::CurvatureStencil<GridType> stencil(BaseType::mGrid, dx);
283  for (size_t n=range.begin(), e=range.end(); n != e; ++n) {
284  BufferType& buffer = BaseType::mLeafs->getBuffer(n,1);
285  for (VoxelIterT iter = BaseType::mLeafs->leaf(n).cbeginValueOn(); iter; ++iter) {
286  stencil.moveTo(iter);
287  buffer.setValue(iter.pos(), stencil.getValue() + dt * stencil.meanCurvatureNormGrad());
288  }
289  }
290 }
291 
299 template<typename GridT, typename InterruptT>
300 inline void
301 LevelSetFilter<GridT, InterruptT>::doLaplacian(const RangeType& range)
302 {
303  typedef typename LeafType::ValueOnCIter VoxelIterT;
304  BaseType::checkInterrupter();
305  //const float CFL = 0.9f, half_dt = CFL * mDx * mDx / 12.0f;
306  const ValueType dx = BaseType::voxelSize(), dt = math::Pow2(dx) / ValueType(6.0);
307  math::GradStencil<GridType> stencil(BaseType::mGrid, dx);
308  for (size_t n=range.begin(), e=range.end(); n != e; ++n) {
309  BufferType& buffer = BaseType::mLeafs->getBuffer(n,1);
310  for (VoxelIterT iter = BaseType::mLeafs->leaf(n).cbeginValueOn(); iter; ++iter) {
311  stencil.moveTo(iter);
312  buffer.setValue(iter.pos(), stencil.getValue() + dt * stencil.laplacian());
313  }
314  }
315 }
316 
318 template<typename GridT, typename InterruptT>
319 inline void
320 LevelSetFilter<GridT, InterruptT>::doOffset(const RangeType& range, ValueType value)
321 {
322  BaseType::checkInterrupter();
323  for (size_t n=range.begin(), e=range.end(); n != e; ++n)
324  BaseType::mLeafs->leaf(n).addValue(value);
325 }
326 
328 template<typename GridT, typename InterruptT>
329 inline void
330 LevelSetFilter<GridT, InterruptT>::doMedian(const RangeType& range, int width)
331 {
332  typedef typename LeafType::ValueOnCIter VoxelIterT;
333  BaseType::checkInterrupter();
334  typename math::DenseStencil<GridType> stencil(BaseType::mGrid, width);//creates local cache!
335  for (size_t n=range.begin(), e=range.end(); n != e; ++n) {
336  BufferType& buffer = BaseType::mLeafs->getBuffer(n,1);
337  for (VoxelIterT iter=BaseType::mLeafs->leaf(n).cbeginValueOn(); iter; ++iter) {
338  stencil.moveTo(iter);
339  buffer.setValue(iter.pos(), stencil.median());
340  }
341  }
342 }
343 
345 template<typename GridT, typename InterruptT>
346 inline void
347 LevelSetFilter<GridT, InterruptT>::doBoxX(const RangeType& range, Int32 w)
348 {
349  this->checkInterrupter();
350  const ValueType frac = ValueType(1)/ValueType(2*w+1);
351  typename GridT::ConstAccessor acc = BaseType::mGrid.getConstAccessor();
352  for (size_t n=range.begin(), nLast = range.end(); n != nLast; ++n) {
353  const LeafType& leaf = BaseType::mLeafs->leaf(n);
354  BufferType& buffer = BaseType::mLeafs->getBuffer(n, 1);
355  for (typename LeafType::ValueOnCIter iter = leaf.cbeginValueOn(); iter; ++iter) {
356  ValueType sum = zeroVal<ValueType>();
357  math::Coord xyz = iter.getCoord();
358  for (Int32 x = xyz.x()-w, xLast = xyz.x()+w; x <= xLast; ++x) {
359  sum += acc.getValue(xyz.setX(x));
360  }
361  buffer.setValue(iter.pos(), sum*frac);
362  }
363  }
364 }
365 
367 template<typename GridT, typename InterruptT>
368 inline void
369 LevelSetFilter<GridT, InterruptT>::doBoxY(const RangeType& range, Int32 w)
370 {
371  this->checkInterrupter();
372  const ValueType frac = ValueType(1)/ValueType(2*w+1);
373  typename GridT::ConstAccessor acc = BaseType::mGrid.getConstAccessor();
374  for (size_t n=range.begin(), nLast = range.end(); n != nLast; ++n) {
375  const LeafType& leaf = BaseType::mLeafs->leaf(n);
376  BufferType& buffer = BaseType::mLeafs->getBuffer(n, 1);
377  for (typename LeafType::ValueOnCIter iter = leaf.cbeginValueOn(); iter; ++iter) {
378  ValueType sum = zeroVal<ValueType>();
379  math::Coord xyz = iter.getCoord();
380  for (Int32 y = xyz.y()-w, yLast = xyz.y()+w; y <= yLast; ++y) {
381  sum += acc.getValue(xyz.setY(y));
382  }
383  buffer.setValue(iter.pos(), sum*frac);
384  }
385  }
386 }
387 
389 template<typename GridT, typename InterruptT>
390 inline void
391 LevelSetFilter<GridT, InterruptT>::doBoxZ(const RangeType& range, Int32 w)
392 {
393  this->checkInterrupter();
394  const ValueType frac = ValueType(1)/ValueType(2*w+1);
395  typename GridT::ConstAccessor acc = BaseType::mGrid.getConstAccessor();
396  for (size_t n=range.begin(), nLast = range.end(); n != nLast; ++n) {
397  const LeafType& leaf = BaseType::mLeafs->leaf(n);
398  BufferType& buffer = BaseType::mLeafs->getBuffer(n, 1);
399  for (typename LeafType::ValueOnCIter iter = leaf.cbeginValueOn(); iter; ++iter) {
400  ValueType sum = zeroVal<ValueType>();
401  math::Coord xyz = iter.getCoord();
402  for (Int32 z = xyz.z()-w, zLast = xyz.z()+w; z <= zLast; ++z) {
403  sum += acc.getValue(xyz.setZ(z));
404  }
405  buffer.setValue(iter.pos(), sum*frac);
406  }
407  }
408 }
409 
410 } // namespace tools
411 } // namespace OPENVDB_VERSION_NAME
412 } // namespace openvdb
413 
414 #endif // OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
415 
416 // Copyright (c) 2012-2013 DreamWorks Animation LLC
417 // All rights reserved. This software is distributed under the
418 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )