OpenVDB  2.3.0
LevelSetMeasure.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 //
34 
35 #ifndef OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
36 #define OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
37 
38 #include <tbb/parallel_for.h>
39 #include <tbb/parallel_sort.h>
40 #include <boost/bind.hpp>
41 #include <boost/function.hpp>
42 #include <boost/type_traits/is_floating_point.hpp>
43 #include <boost/utility/enable_if.hpp>
44 #include <boost/math/constants/constants.hpp>//for Pi
45 #include <openvdb/math/Math.h>
46 #include <openvdb/Types.h>
47 #include <openvdb/Grid.h>
48 #include <openvdb/tree/LeafManager.h>
49 #include <openvdb/tree/ValueAccessor.h>
50 #include <openvdb/math/FiniteDifference.h>
51 #include <openvdb/math/Operators.h>
52 #include <openvdb/util/NullInterrupter.h>
53 
54 namespace openvdb {
56 namespace OPENVDB_VERSION_NAME {
57 namespace tools {
58 
67 template<class GridType>
68 inline Real
69 levelSetArea(const GridType& grid, bool useWorldSpace = true);
70 
79 template<class GridType>
80 inline Real
81 levelSetVolume(const GridType& grid, bool useWorldSpace = true);
82 
93 template<class GridType>
94 inline void
95 levelSetMeasure(const GridType& grid, Real& area, Real& volume, bool useWorldSpace = true);
96 
108 template<class GridType>
109 inline void
110 levelSetMeasure(const GridType& grid, Real& area, Real& volume, Real& avgCurvature,
111  bool useWorldSpace = true);
112 
114 template<typename RealT>
116 {
117 public:
118  DiracDelta(RealT eps) : mC(0.5/eps), mD(2*boost::math::constants::pi<RealT>()*mC), mE(eps) {}
119  inline RealT operator()(RealT phi) const { return math::Abs(phi) > mE ? 0 : mC*(1+cos(mD*phi)); }
120 private:
121  const RealT mC, mD, mE;
122 };
123 
124 
132 template<typename GridT,
133  typename InterruptT = util::NullInterrupter>
135 {
136 public:
137  typedef GridT GridType;
138  typedef typename GridType::TreeType TreeType;
139  typedef typename TreeType::ValueType ValueType;
142  BOOST_STATIC_ASSERT(boost::is_floating_point<ValueType>::value);
143 
148  LevelSetMeasure(const GridType& grid, InterruptT* interrupt = NULL);
149 
150  LevelSetMeasure(ManagerType& leafs, Real Dx, InterruptT* interrupt);
151 
153  void reinit(const GridType& grid);
154 
156  void reinit(ManagerType& leafs, Real dx);
157 
160 
162  int getGrainSize() const { return mGrainSize; }
163 
166  void setGrainSize(int grainsize) { mGrainSize = grainsize; }
167 
173  void measure(Real& area, Real& volume, bool useWorldUnits = true);
174 
180  void measure(Real& area, Real& volume, Real& avgMeanCurvature, bool useWorldUnits = true);
181 
185  void operator()(const RangeType& range) const
186  {
187  if (mTask) mTask(const_cast<LevelSetMeasure*>(this), range);
188  else OPENVDB_THROW(ValueError, "task is undefined");
189  }
190 
191 private:
192  typedef typename GridT::ConstAccessor AccT;
193  typedef typename TreeType::LeafNodeType LeafT;
194  typedef typename LeafT::ValueOnCIter VoxelCIterT;
195  typedef typename ManagerType::BufferType BufferT;
196  typedef typename RangeType::Iterator LeafIterT;
197 
198  AccT mAcc;
199  ManagerType* mLeafs;
200  InterruptT* mInterrupter;
201  double mDx;
202  double* mArray;
203  typename boost::function<void (LevelSetMeasure*, const RangeType&)> mTask;
204  int mGrainSize;
205 
206  // @brief Return false if the process was interrupted
207  bool checkInterrupter();
208 
209  // Private methods called by tbb::parallel_reduce threads
210  void measure2( const RangeType& );
211 
212  // Private methods called by tbb::parallel_reduce threads
213  void measure3( const RangeType& );
214 
215  inline double reduce(double* first, double scale)
216  {
217  double* last = first + mLeafs->leafCount();
218  tbb::parallel_sort(first, last);//reduces catastrophic cancellation
219  Real sum = 0.0;
220  while(first != last) sum += *first++;
221  return scale * sum;
222  }
223 
224 }; // end of LevelSetMeasure class
225 
226 
227 template<typename GridT, typename InterruptT>
228 inline
230  : mAcc(grid.tree())
231  , mLeafs(NULL)
232  , mInterrupter(interrupt)
233  , mDx(grid.voxelSize()[0])
234  , mArray(NULL)
235  , mTask(0)
236  , mGrainSize(1)
237 {
238  if (!grid.hasUniformVoxels()) {
240  "The transform must have uniform scale for the LevelSetMeasure to function");
241  }
242  if (grid.getGridClass() != GRID_LEVEL_SET) {
244  "LevelSetMeasure only supports level sets;"
245  " try setting the grid class to \"level set\"");
246  }
247 }
248 
249 
250 template<typename GridT, typename InterruptT>
251 inline
253  ManagerType& leafs, Real dx, InterruptT* interrupt)
254  : mAcc(leafs.tree())
255  , mLeafs(&leafs)
256  , mInterrupter(interrupt)
257  , mDx(dx)
258  , mArray(NULL)
259  , mTask(0)
260  , mGrainSize(1)
261 {
262 }
263 
264 template<typename GridT, typename InterruptT>
265 inline void
267 {
268  if (!grid.hasUniformVoxels()) {
270  "The transform must have uniform scale for the LevelSetMeasure to function");
271  }
272  if (grid.getGridClass() != GRID_LEVEL_SET) {
274  "LevelSetMeasure only supports level sets;"
275  " try setting the grid class to \"level set\"");
276  }
277  mLeafs = NULL;
278  mAcc = grid.getConstAccessor();
279  mDx = grid.voxelSize()[0];
280 }
281 
282 
283 template<typename GridT, typename InterruptT>
284 inline void
286 {
287  mLeafs = &leafs;
288  mAcc = AccT(leafs.tree());
289  mDx = dx;
290 }
291 
293 
294 
295 template<typename GridT, typename InterruptT>
296 inline void
297 LevelSetMeasure<GridT, InterruptT>::measure(Real& area, Real& volume, bool useWorldUnits)
298 {
299  if (mInterrupter) mInterrupter->start("Measuring level set");
300  mTask = boost::bind(&LevelSetMeasure::measure2, _1, _2);
301 
302  const bool newLeafs = mLeafs == NULL;
303  if (newLeafs) mLeafs = new ManagerType(mAcc.tree());
304  const size_t leafCount = mLeafs->leafCount();
305  if (leafCount == 0) {
306  area = volume = 0;
307  return;
308  }
309  mArray = new double[2*leafCount];
310 
311  if (mGrainSize>0) {
312  tbb::parallel_for(mLeafs->leafRange(mGrainSize), *this);
313  } else {
314  (*this)(mLeafs->leafRange());
315  }
316 
317  const double dx = useWorldUnits ? mDx : 1.0;
318  area = this->reduce(mArray, math::Pow2(dx));
319  volume = this->reduce(mArray + leafCount, math::Pow3(dx) / 3.0);
320 
321  if (newLeafs) {
322  delete mLeafs;
323  mLeafs = NULL;
324  }
325  delete [] mArray;
326 
327  if (mInterrupter) mInterrupter->end();
328 }
329 
330 
331 template<typename GridT, typename InterruptT>
332 inline void
333 LevelSetMeasure<GridT, InterruptT>::measure(Real& area, Real& volume, Real& avgMeanCurvature,
334  bool useWorldUnits)
335 {
336  if (mInterrupter) mInterrupter->start("Measuring level set");
337  mTask = boost::bind(&LevelSetMeasure::measure3, _1, _2);
338 
339  const bool newLeafs = mLeafs == NULL;
340  if (newLeafs) mLeafs = new ManagerType(mAcc.tree());
341  const size_t leafCount = mLeafs->leafCount();
342  if (leafCount == 0) {
343  area = volume = avgMeanCurvature = 0;
344  return;
345  }
346  mArray = new double[3*leafCount];
347 
348  if (mGrainSize>0) {
349  tbb::parallel_for(mLeafs->leafRange(mGrainSize), *this);
350  } else {
351  (*this)(mLeafs->leafRange());
352  }
353 
354  const double dx = useWorldUnits ? mDx : 1.0;
355  area = this->reduce(mArray, math::Pow2(dx));
356  volume = this->reduce(mArray + leafCount, math::Pow3(dx) / 3.0);
357  avgMeanCurvature = this->reduce(mArray + 2*leafCount, dx/area);
358 
359  if (newLeafs) {
360  delete mLeafs;
361  mLeafs = NULL;
362  }
363  delete [] mArray;
364 
365  if (mInterrupter) mInterrupter->end();
366 }
367 
368 
370 
371 
372 template<typename GridT, typename InterruptT>
373 inline bool
375 {
376  if (util::wasInterrupted(mInterrupter)) {
377  tbb::task::self().cancel_group_execution();
378  return false;
379  }
380  return true;
381 }
382 
383 template<typename GridT, typename InterruptT>
384 inline void
385 LevelSetMeasure<GridT, InterruptT>::measure2(const RangeType& range)
386 {
387  typedef math::Vec3<ValueType> Vec3T;
388  typedef math::ISGradient<math::CD_2ND> Grad;
389  this->checkInterrupter();
390  const Real invDx = 1.0/mDx;
391  const DiracDelta<Real> DD(1.5);
392  const size_t leafCount = mLeafs->leafCount();
393  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
394  Real sumA = 0, sumV = 0;//reduce risk of catastrophic cancellation
395  for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
396  const Real dd = DD(invDx * (*voxelIter));
397  if (dd > 0.0) {
398  const Coord p = voxelIter.getCoord();
399  const Vec3T g = invDx*Grad::result(mAcc, p);//voxel units
400  sumA += dd * g.dot(g);
401  sumV += dd * (g[0]*p[0]+g[1]*p[1]+g[2]*p[2]);
402  }
403  }
404  double* v = mArray + leafIter.pos();
405  *v = sumA;
406  v += leafCount;
407  *v = sumV;
408  }
409 }
410 
411 
412 template<typename GridT, typename InterruptT>
413 inline void
414 LevelSetMeasure<GridT, InterruptT>::measure3(const RangeType& range)
415 {
416  typedef math::Vec3<ValueType> Vec3T;
417  typedef math::ISGradient<math::CD_2ND> Grad;
418  typedef math::ISMeanCurvature<math::CD_SECOND, math::CD_2ND> Curv;
419  this->checkInterrupter();
420  const Real invDx = 1.0/mDx;
421  const DiracDelta<Real> DD(1.5);
422  ValueType alpha, beta;
423  const size_t leafCount = mLeafs->leafCount();
424  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
425  Real sumA = 0, sumV = 0, sumC = 0;//reduce risk of catastrophic cancellation
426  for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
427  const Real dd = DD(invDx * (*voxelIter));
428  if (dd > 0.0) {
429  const Coord p = voxelIter.getCoord();
430  const Vec3T g = invDx*Grad::result(mAcc, p);//voxel units
431  const Real dA = dd * g.dot(g);
432  sumA += dA;
433  sumV += dd * (g[0]*p[0]+g[1]*p[1]+g[2]*p[2]);
434  Curv::result(mAcc, p, alpha, beta);
435  sumC += dA * alpha/(2*math::Pow2(beta))*invDx;
436  }
437  }
438  double* v = mArray + leafIter.pos();
439  *v = sumA;
440  v += leafCount;
441  *v = sumV;
442  v += leafCount;
443  *v = sumC;
444  }
445 }
446 
448 
449 template<class GridT>
450 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType>, Real>::type
451 doLevelSetArea(const GridT& grid, bool useWorldSpace)
452 {
453  Real area, volume;
454  LevelSetMeasure<GridT> m(grid);
455  m.measure(area, volume, useWorldSpace);
456  return area;
457 }
458 
459 template<class GridT>
460 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType>, Real>::type
461 doLevelSetArea(const GridT&, bool)
462 {
464  "level set area is supported only for scalar, floating-point grids");
465 }
466 
467 template<class GridT>
468 inline Real
469 levelSetArea(const GridT& grid, bool useWorldSpace)
470 {
471  return doLevelSetArea<GridT>(grid, useWorldSpace);
472 }
473 
475 
476 template<class GridT>
477 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType>, Real>::type
478 doLevelSetVolume(const GridT& grid, bool useWorldSpace)
479 {
480  Real area, volume;
481  LevelSetMeasure<GridT> m(grid);
482  m.measure(area, volume, useWorldSpace);
483  return volume;
484 }
485 
486 template<class GridT>
487 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType>, Real>::type
488 doLevelSetVolume(const GridT&, bool)
489 {
491  "level set volume is supported only for scalar, floating-point grids");
492 }
493 
494 template<class GridT>
495 inline Real
496 levelSetVolume(const GridT& grid, bool useWorldSpace)
497 {
498  return doLevelSetVolume<GridT>(grid, useWorldSpace);
499 }
500 
502 
503 template<class GridT>
504 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType> >::type
505 doLevelSetMeasure(const GridT& grid, Real& area, Real& volume, bool useWorldSpace)
506 {
507  LevelSetMeasure<GridT> m(grid);
508  m.measure(area, volume, useWorldSpace);
509 }
510 
511 template<class GridT>
512 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType> >::type
513 doLevelSetMeasure(const GridT&, Real&, Real&, bool)
514 {
516  "level set measure is supported only for scalar, floating-point grids");
517 }
518 
519 template<class GridT>
520 inline void
521 levelSetMeasure(const GridT& grid, Real& area, Real& volume, bool useWorldSpace)
522 {
523  doLevelSetMeasure<GridT>(grid, area, volume, useWorldSpace);
524 }
525 
527 
528 template<class GridT>
529 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType> >::type
530 doLevelSetMeasure(const GridT& grid, Real& area, Real& volume, Real& avgCurvature,
531  bool useWorldSpace)
532 {
533  LevelSetMeasure<GridT> m(grid);
534  m.measure(area, volume, avgCurvature, useWorldSpace);
535 }
536 
537 template<class GridT>
538 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType> >::type
539 doLevelSetMeasure(const GridT&, Real&, Real&, Real&, bool)
540 {
542  "level set measure is supported only for scalar, floating-point grids");
543 }
544 
545 template<class GridT>
546 inline void
547 levelSetMeasure(const GridT& grid, Real& area, Real& volume, Real& avgCurvature, bool useWorldSpace)
548 {
549  doLevelSetMeasure<GridT>(grid, area, volume, avgCurvature, useWorldSpace);
550 }
551 
552 } // namespace tools
553 } // namespace OPENVDB_VERSION_NAME
554 } // namespace openvdb
555 
556 #endif // OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
557 
558 // Copyright (c) 2012-2013 DreamWorks Animation LLC
559 // All rights reserved. This software is distributed under the
560 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
void setGrainSize(int grainsize)
Set the grain-size used for multi-threading.
Definition: LevelSetMeasure.h:166
Real levelSetArea(const GridT &grid, bool useWorldSpace)
Definition: LevelSetMeasure.h:469
GridType::TreeType TreeType
Definition: LevelSetMeasure.h:138
boost::enable_if< boost::is_floating_point< typename GridT::ValueType >, Real >::type doLevelSetArea(const GridT &grid, bool useWorldSpace)
Definition: LevelSetMeasure.h:451
boost::enable_if< boost::is_floating_point< typename GridT::ValueType >, Real >::type doLevelSetVolume(const GridT &grid, bool useWorldSpace)
Definition: LevelSetMeasure.h:478
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Definition: Mat.h:146
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:244
Real levelSetArea(const GridType &grid, bool useWorldSpace=true)
Return the surface area of a narrow-band level set.
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
Type Pow2(Type x)
Return .
Definition: Math.h:458
DiracDelta(RealT eps)
Definition: LevelSetMeasure.h:118
int getGrainSize() const
Definition: LevelSetMeasure.h:162
Definition: Exceptions.h:88
Smeared-out and continuous Dirac Delta function.
Definition: LevelSetMeasure.h:115
Definition: Types.h:170
TreeType::ValueType ValueType
Definition: LevelSetMeasure.h:139
LevelSetMeasure(const GridType &grid, InterruptT *interrupt=NULL)
Main constructor from a grid.
Definition: LevelSetMeasure.h:229
~LevelSetMeasure()
Destructor.
Definition: LevelSetMeasure.h:159
boost::enable_if< boost::is_floating_point< typename GridT::ValueType > >::type doLevelSetMeasure(const GridT &grid, Real &area, Real &volume, bool useWorldSpace)
Definition: LevelSetMeasure.h:505
Definition: Exceptions.h:87
void measure(Real &area, Real &volume, bool useWorldUnits=true)
Compute the surface area and volume of the level set. Use the last argument to specify the result in ...
Definition: LevelSetMeasure.h:297
#define OPENVDB_VERSION_NAME
Definition: version.h:45
tree::LeafManager< const TreeType > ManagerType
Definition: LevelSetMeasure.h:140
Gradient operators defined in index space of various orders.
Definition: Operators.h:122
void operator()(const RangeType &range) const
Used internally by tbb::parallel_reduce().
Definition: LevelSetMeasure.h:185
size_t leafCount() const
Return the number of leaf nodes.
Definition: LeafManager.h:295
MatType scale(const Vec3< typename MatType::value_type > &s)
Return a matrix that scales by s.
Definition: Mat.h:593
GridT GridType
Definition: LevelSetMeasure.h:137
void levelSetMeasure(const GridT &grid, Real &area, Real &volume, Real &avgCurvature, bool useWorldSpace)
Definition: LevelSetMeasure.h:547
ManagerType::LeafRange RangeType
Definition: LevelSetMeasure.h:141
Type Pow3(Type x)
Return .
Definition: Math.h:462
RealT operator()(RealT phi) const
Definition: LevelSetMeasure.h:119
TreeType & tree()
Return the tree associated with this manager.
Definition: LeafManager.h:298
Definition: Exceptions.h:86
void reinit(const GridType &grid)
Re-initialize using the specified grid.
Definition: LevelSetMeasure.h:266
Real levelSetVolume(const GridType &grid, bool useWorldSpace=true)
Return the volume of a narrow-band level set surface.
double Real
Definition: Types.h:63
void levelSetMeasure(const GridType &grid, Real &area, Real &volume, bool useWorldSpace=true)
Compute the surface area and volume of a narrow-band level set.
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:67
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:109
Multi-threaded computation of surface area, volume and average mean-curvature for narrow band level s...
Definition: LevelSetMeasure.h:134
Real levelSetVolume(const GridT &grid, bool useWorldSpace)
Definition: LevelSetMeasure.h:496