diff options
Diffstat (limited to 'extern/openvdb/internal/openvdb/tools/LevelSetMeasure.h')
| -rw-r--r-- | extern/openvdb/internal/openvdb/tools/LevelSetMeasure.h | 560 |
1 files changed, 560 insertions, 0 deletions
diff --git a/extern/openvdb/internal/openvdb/tools/LevelSetMeasure.h b/extern/openvdb/internal/openvdb/tools/LevelSetMeasure.h new file mode 100644 index 00000000000..600900e0166 --- /dev/null +++ b/extern/openvdb/internal/openvdb/tools/LevelSetMeasure.h @@ -0,0 +1,560 @@ +/////////////////////////////////////////////////////////////////////////// +// +// Copyright (c) 2012-2013 DreamWorks Animation LLC +// +// All rights reserved. This software is distributed under the +// Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ ) +// +// Redistributions of source code must retain the above copyright +// and license notice and the following restrictions and disclaimer. +// +// * Neither the name of DreamWorks Animation nor the names of +// its contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL, +// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE +// LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00. +// +/////////////////////////////////////////////////////////////////////////// +// +/// @author Ken Museth +/// +/// @file LevelSetMeasure.h + +#ifndef OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED +#define OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED + +#include <tbb/parallel_for.h> +#include <tbb/parallel_sort.h> +#include <boost/bind.hpp> +#include <boost/function.hpp> +#include <boost/type_traits/is_floating_point.hpp> +#include <boost/utility/enable_if.hpp> +#include <boost/math/constants/constants.hpp>//for Pi +#include <openvdb/math/Math.h> +#include <openvdb/Types.h> +#include <openvdb/Grid.h> +#include <openvdb/tree/LeafManager.h> +#include <openvdb/tree/ValueAccessor.h> +#include <openvdb/math/FiniteDifference.h> +#include <openvdb/math/Operators.h> +#include <openvdb/util/NullInterrupter.h> + +namespace openvdb { +OPENVDB_USE_VERSION_NAMESPACE +namespace OPENVDB_VERSION_NAME { +namespace tools { + +/// @brief Return the surface area of a narrow-band level set. +/// +/// @param grid a scalar, floating-point grid with one or more disjoint, +/// closed isosurfaces at the given @a isovalue +/// @param useWorldSpace if true the area is computed in +/// world space units, else in voxel units. +/// +/// @throw TypeError if @a grid is not scalar or not floating-point or not a level set. +template<class GridType> +inline Real +levelSetArea(const GridType& grid, bool useWorldSpace = true); + +/// @brief Return the volume of a narrow-band level set surface. +/// +/// @param grid a scalar, floating-point grid with one or more disjoint, +/// closed isosurfaces at the given @a isovalue +/// @param useWorldSpace if true the volume is computed in +/// world space units, else in voxel units. +/// +/// @throw TypeError if @a grid is not scalar or not floating-point or not a level set. +template<class GridType> +inline Real +levelSetVolume(const GridType& grid, bool useWorldSpace = true); + +/// @brief Compute the surface area and volume of a narrow-band level set. +/// +/// @param grid a scalar, floating-point grid with one or more disjoint, +/// closed isosurfaces at the given @a isovalue +/// @param area surface area of the level set +/// @param volume volume of the level set surface +/// @param useWorldSpace if true the area and volume are computed in +/// world space units, else in voxel units. +/// +/// @throw TypeError if @a grid is not scalar or not floating-point or not a level set. +template<class GridType> +inline void +levelSetMeasure(const GridType& grid, Real& area, Real& volume, bool useWorldSpace = true); + +/// @brief Compute the surface area and volume of a narrow-band level set. +/// +/// @param grid a scalar, floating-point grid with one or more disjoint, +/// closed isosurfaces at the given @a isovalue +/// @param area surface area of the level set +/// @param volume volume of the level set surface +/// @param avgCurvature average mean curvature of the level set surface +/// @param useWorldSpace if true the area, volume and curvature are computed in +/// world space units, else in voxel units. +/// +/// @throw TypeError if @a grid is not scalar or not floating-point or not a level set. +template<class GridType> +inline void +levelSetMeasure(const GridType& grid, Real& area, Real& volume, Real& avgCurvature, + bool useWorldSpace = true); + +/// @brief Smeared-out and continuous Dirac Delta function. +template<typename RealT> +class DiracDelta +{ +public: + DiracDelta(RealT eps) : mC(0.5/eps), mD(2*boost::math::constants::pi<RealT>()*mC), mE(eps) {} + inline RealT operator()(RealT phi) const { return math::Abs(phi) > mE ? 0 : mC*(1+cos(mD*phi)); } +private: + const RealT mC, mD, mE; +}; + + +/// @brief Multi-threaded computation of surface area, volume and +/// average mean-curvature for narrow band level sets. +/// +/// @details To reduce the risk of round-off errors (primarily due to +/// catastrophic cancellation) and guarantee determinism during +/// multi-threading this class is implemented using parallel_for, and +/// delayed reduction of a sorted list. +template<typename GridT, + typename InterruptT = util::NullInterrupter> +class LevelSetMeasure +{ +public: + typedef GridT GridType; + typedef typename GridType::TreeType TreeType; + typedef typename TreeType::ValueType ValueType; + typedef typename tree::LeafManager<const TreeType> ManagerType; + typedef typename ManagerType::LeafRange RangeType; + BOOST_STATIC_ASSERT(boost::is_floating_point<ValueType>::value); + + /// @brief Main constructor from a grid + /// @param grid The level set to be measured. + /// @param interrupt Optional interrupter. + /// @throw RuntimeError if the grid is not a level set. + LevelSetMeasure(const GridType& grid, InterruptT* interrupt = NULL); + + LevelSetMeasure(ManagerType& leafs, Real Dx, InterruptT* interrupt); + + /// @brief Re-initialize using the specified grid. + void reinit(const GridType& grid); + + /// @brief Re-initialize using the specified LeafManager and voxelSize. + void reinit(ManagerType& leafs, Real dx); + + /// @brief Destructor + ~LevelSetMeasure() {} + + /// @return the grain-size used for multi-threading + int getGrainSize() const { return mGrainSize; } + + /// @brief Set the grain-size used for multi-threading. + /// @note A grainsize of 0 or less disables multi-threading! + void setGrainSize(int grainsize) { mGrainSize = grainsize; } + + /// @brief Compute the surface area and volume of the level + /// set. Use the last argument to specify the result in world or + /// voxel units. + /// @note This method is faster (about 3x) then the measure method + /// below that also computes the average mean-curvature. + void measure(Real& area, Real& volume, bool useWorldUnits = true); + + /// @brief Compute the surface area, volume, and average + /// mean-curvatue of the level set. Use the last argument to + /// specify the result in world or voxel units. + /// @note This method is slower (about 3x) then the measure method + /// above that only computes the area and volume. + void measure(Real& area, Real& volume, Real& avgMeanCurvature, bool useWorldUnits = true); + + /// @brief Used internally by tbb::parallel_reduce(). + /// @param range The range over which to perform multi-threading. + /// @warning Never call this method directly! + void operator()(const RangeType& range) const + { + if (mTask) mTask(const_cast<LevelSetMeasure*>(this), range); + else OPENVDB_THROW(ValueError, "task is undefined"); + } + +private: + typedef typename GridT::ConstAccessor AccT; + typedef typename TreeType::LeafNodeType LeafT; + typedef typename LeafT::ValueOnCIter VoxelCIterT; + typedef typename ManagerType::BufferType BufferT; + typedef typename RangeType::Iterator LeafIterT; + + AccT mAcc; + ManagerType* mLeafs; + InterruptT* mInterrupter; + double mDx; + double* mArray; + typename boost::function<void (LevelSetMeasure*, const RangeType&)> mTask; + int mGrainSize; + + // @brief Return false if the process was interrupted + bool checkInterrupter(); + + // Private methods called by tbb::parallel_reduce threads + void measure2( const RangeType& ); + + // Private methods called by tbb::parallel_reduce threads + void measure3( const RangeType& ); + + inline double reduce(double* first, double scale) + { + double* last = first + mLeafs->leafCount(); + tbb::parallel_sort(first, last);//reduces catastrophic cancellation + Real sum = 0.0; + while(first != last) sum += *first++; + return scale * sum; + } + +}; // end of LevelSetMeasure class + + +template<typename GridT, typename InterruptT> +inline +LevelSetMeasure<GridT, InterruptT>::LevelSetMeasure(const GridType& grid, InterruptT* interrupt) + : mAcc(grid.tree()) + , mLeafs(NULL) + , mInterrupter(interrupt) + , mDx(grid.voxelSize()[0]) + , mArray(NULL) + , mTask(0) + , mGrainSize(1) +{ + if (!grid.hasUniformVoxels()) { + OPENVDB_THROW(RuntimeError, + "The transform must have uniform scale for the LevelSetMeasure to function"); + } + if (grid.getGridClass() != GRID_LEVEL_SET) { + OPENVDB_THROW(RuntimeError, + "LevelSetMeasure only supports level sets;" + " try setting the grid class to \"level set\""); + } +} + + +template<typename GridT, typename InterruptT> +inline +LevelSetMeasure<GridT, InterruptT>::LevelSetMeasure( + ManagerType& leafs, Real dx, InterruptT* interrupt) + : mAcc(leafs.tree()) + , mLeafs(&leafs) + , mInterrupter(interrupt) + , mDx(dx) + , mArray(NULL) + , mTask(0) + , mGrainSize(1) +{ +} + +template<typename GridT, typename InterruptT> +inline void +LevelSetMeasure<GridT, InterruptT>::reinit(const GridType& grid) +{ + if (!grid.hasUniformVoxels()) { + OPENVDB_THROW(RuntimeError, + "The transform must have uniform scale for the LevelSetMeasure to function"); + } + if (grid.getGridClass() != GRID_LEVEL_SET) { + OPENVDB_THROW(RuntimeError, + "LevelSetMeasure only supports level sets;" + " try setting the grid class to \"level set\""); + } + mLeafs = NULL; + mAcc = grid.getConstAccessor(); + mDx = grid.voxelSize()[0]; +} + + +template<typename GridT, typename InterruptT> +inline void +LevelSetMeasure<GridT, InterruptT>::reinit(ManagerType& leafs, Real dx) +{ + mLeafs = &leafs; + mAcc = AccT(leafs.tree()); + mDx = dx; +} + +//////////////////////////////////////// + + +template<typename GridT, typename InterruptT> +inline void +LevelSetMeasure<GridT, InterruptT>::measure(Real& area, Real& volume, bool useWorldUnits) +{ + if (mInterrupter) mInterrupter->start("Measuring level set"); + mTask = boost::bind(&LevelSetMeasure::measure2, _1, _2); + + const bool newLeafs = mLeafs == NULL; + if (newLeafs) mLeafs = new ManagerType(mAcc.tree()); + const size_t leafCount = mLeafs->leafCount(); + if (leafCount == 0) { + area = volume = 0; + return; + } + mArray = new double[2*leafCount]; + + if (mGrainSize>0) { + tbb::parallel_for(mLeafs->leafRange(mGrainSize), *this); + } else { + (*this)(mLeafs->leafRange()); + } + + const double dx = useWorldUnits ? mDx : 1.0; + area = this->reduce(mArray, math::Pow2(dx)); + volume = this->reduce(mArray + leafCount, math::Pow3(dx) / 3.0); + + if (newLeafs) { + delete mLeafs; + mLeafs = NULL; + } + delete [] mArray; + + if (mInterrupter) mInterrupter->end(); +} + + +template<typename GridT, typename InterruptT> +inline void +LevelSetMeasure<GridT, InterruptT>::measure(Real& area, Real& volume, Real& avgMeanCurvature, + bool useWorldUnits) +{ + if (mInterrupter) mInterrupter->start("Measuring level set"); + mTask = boost::bind(&LevelSetMeasure::measure3, _1, _2); + + const bool newLeafs = mLeafs == NULL; + if (newLeafs) mLeafs = new ManagerType(mAcc.tree()); + const size_t leafCount = mLeafs->leafCount(); + if (leafCount == 0) { + area = volume = avgMeanCurvature = 0; + return; + } + mArray = new double[3*leafCount]; + + if (mGrainSize>0) { + tbb::parallel_for(mLeafs->leafRange(mGrainSize), *this); + } else { + (*this)(mLeafs->leafRange()); + } + + const double dx = useWorldUnits ? mDx : 1.0; + area = this->reduce(mArray, math::Pow2(dx)); + volume = this->reduce(mArray + leafCount, math::Pow3(dx) / 3.0); + avgMeanCurvature = this->reduce(mArray + 2*leafCount, dx/area); + + if (newLeafs) { + delete mLeafs; + mLeafs = NULL; + } + delete [] mArray; + + if (mInterrupter) mInterrupter->end(); +} + + +///////////////////////// PRIVATE METHODS ////////////////////// + + +template<typename GridT, typename InterruptT> +inline bool +LevelSetMeasure<GridT, InterruptT>::checkInterrupter() +{ + if (util::wasInterrupted(mInterrupter)) { + tbb::task::self().cancel_group_execution(); + return false; + } + return true; +} + +template<typename GridT, typename InterruptT> +inline void +LevelSetMeasure<GridT, InterruptT>::measure2(const RangeType& range) +{ + typedef math::Vec3<ValueType> Vec3T; + typedef math::ISGradient<math::CD_2ND> Grad; + this->checkInterrupter(); + const Real invDx = 1.0/mDx; + const DiracDelta<Real> DD(1.5); + const size_t leafCount = mLeafs->leafCount(); + for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) { + Real sumA = 0, sumV = 0;//reduce risk of catastrophic cancellation + for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) { + const Real dd = DD(invDx * (*voxelIter)); + if (dd > 0.0) { + const Coord p = voxelIter.getCoord(); + const Vec3T g = invDx*Grad::result(mAcc, p);//voxel units + sumA += dd * g.dot(g); + sumV += dd * (g[0]*p[0]+g[1]*p[1]+g[2]*p[2]); + } + } + double* v = mArray + leafIter.pos(); + *v = sumA; + v += leafCount; + *v = sumV; + } +} + + +template<typename GridT, typename InterruptT> +inline void +LevelSetMeasure<GridT, InterruptT>::measure3(const RangeType& range) +{ + typedef math::Vec3<ValueType> Vec3T; + typedef math::ISGradient<math::CD_2ND> Grad; + typedef math::ISMeanCurvature<math::CD_SECOND, math::CD_2ND> Curv; + this->checkInterrupter(); + const Real invDx = 1.0/mDx; + const DiracDelta<Real> DD(1.5); + ValueType alpha, beta; + const size_t leafCount = mLeafs->leafCount(); + for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) { + Real sumA = 0, sumV = 0, sumC = 0;//reduce risk of catastrophic cancellation + for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) { + const Real dd = DD(invDx * (*voxelIter)); + if (dd > 0.0) { + const Coord p = voxelIter.getCoord(); + const Vec3T g = invDx*Grad::result(mAcc, p);//voxel units + const Real dA = dd * g.dot(g); + sumA += dA; + sumV += dd * (g[0]*p[0]+g[1]*p[1]+g[2]*p[2]); + Curv::result(mAcc, p, alpha, beta); + sumC += dA * alpha/(2*math::Pow2(beta))*invDx; + } + } + double* v = mArray + leafIter.pos(); + *v = sumA; + v += leafCount; + *v = sumV; + v += leafCount; + *v = sumC; + } +} + +//////////////////////////////////////// + +template<class GridT> +inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType>, Real>::type +doLevelSetArea(const GridT& grid, bool useWorldSpace) +{ + Real area, volume; + LevelSetMeasure<GridT> m(grid); + m.measure(area, volume, useWorldSpace); + return area; +} + +template<class GridT> +inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType>, Real>::type +doLevelSetArea(const GridT&, bool) +{ + OPENVDB_THROW(TypeError, + "level set area is supported only for scalar, floating-point grids"); +} + +template<class GridT> +inline Real +levelSetArea(const GridT& grid, bool useWorldSpace) +{ + return doLevelSetArea<GridT>(grid, useWorldSpace); +} + +//////////////////////////////////////// + +template<class GridT> +inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType>, Real>::type +doLevelSetVolume(const GridT& grid, bool useWorldSpace) +{ + Real area, volume; + LevelSetMeasure<GridT> m(grid); + m.measure(area, volume, useWorldSpace); + return volume; +} + +template<class GridT> +inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType>, Real>::type +doLevelSetVolume(const GridT&, bool) +{ + OPENVDB_THROW(TypeError, + "level set volume is supported only for scalar, floating-point grids"); +} + +template<class GridT> +inline Real +levelSetVolume(const GridT& grid, bool useWorldSpace) +{ + return doLevelSetVolume<GridT>(grid, useWorldSpace); +} + +//////////////////////////////////////// + +template<class GridT> +inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType> >::type +doLevelSetMeasure(const GridT& grid, Real& area, Real& volume, bool useWorldSpace) +{ + LevelSetMeasure<GridT> m(grid); + m.measure(area, volume, useWorldSpace); +} + +template<class GridT> +inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType> >::type +doLevelSetMeasure(const GridT&, Real&, Real&, bool) +{ + OPENVDB_THROW(TypeError, + "level set measure is supported only for scalar, floating-point grids"); +} + +template<class GridT> +inline void +levelSetMeasure(const GridT& grid, Real& area, Real& volume, bool useWorldSpace) +{ + doLevelSetMeasure<GridT>(grid, area, volume, useWorldSpace); +} + +//////////////////////////////////////// + +template<class GridT> +inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType> >::type +doLevelSetMeasure(const GridT& grid, Real& area, Real& volume, Real& avgCurvature, + bool useWorldSpace) +{ + LevelSetMeasure<GridT> m(grid); + m.measure(area, volume, avgCurvature, useWorldSpace); +} + +template<class GridT> +inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType> >::type +doLevelSetMeasure(const GridT&, Real&, Real&, Real&, bool) +{ + OPENVDB_THROW(TypeError, + "level set measure is supported only for scalar, floating-point grids"); +} + +template<class GridT> +inline void +levelSetMeasure(const GridT& grid, Real& area, Real& volume, Real& avgCurvature, bool useWorldSpace) +{ + doLevelSetMeasure<GridT>(grid, area, volume, avgCurvature, useWorldSpace); +} + +} // namespace tools +} // namespace OPENVDB_VERSION_NAME +} // namespace openvdb + +#endif // OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED + +// Copyright (c) 2012-2013 DreamWorks Animation LLC +// All rights reserved. This software is distributed under the +// Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ ) |