diff --git a/openvdb/openvdb/CMakeLists.txt b/openvdb/openvdb/CMakeLists.txt index 1540e06b06..d659dec5a4 100644 --- a/openvdb/openvdb/CMakeLists.txt +++ b/openvdb/openvdb/CMakeLists.txt @@ -485,6 +485,7 @@ set(OPENVDB_LIBRARY_POINTS_IMPL_INCLUDE_FILES set(OPENVDB_LIBRARY_TOOLS_INCLUDE_FILES tools/Activate.h + tools/Blend.h tools/ChangeBackground.h tools/Clip.h tools/Composite.h diff --git a/openvdb/openvdb/tools/Blend.h b/openvdb/openvdb/tools/Blend.h new file mode 100644 index 0000000000..208ec4ef2f --- /dev/null +++ b/openvdb/openvdb/tools/Blend.h @@ -0,0 +1,643 @@ +// Copyright Contributors to the OpenVDB Project +// SPDX-License-Identifier: MPL-2.0 +// +/// @file Blend.h +/// +/// @author Andre Pradhana +/// +/// @brief Define methods to blend two level-sets together. One such approach +/// is by carving an excess fillet so that the resulting blended +/// level-sets can appear to be smoother than a regular union. +/// +/// @details The algorithm used in the function unionFillet is based on +/// a 2007 SIGGRAPH talk titled "Levelsets in production: Spider-man 3" +/// by Allen et al. Paper is available here: +/// https://dl.acm.org/doi/10.1145/1278780.1278815 + +#ifndef OPENVDB_TOOLS_BLEND_HAS_BEEN_INCLUDED +#define OPENVDB_TOOLS_BLEND_HAS_BEEN_INCLUDED + +#include +#include +#include // for isExactlyEqual() +#include +#include + +namespace openvdb { +OPENVDB_USE_VERSION_NAMESPACE +namespace OPENVDB_VERSION_NAME { +namespace tools { + +/// @brief Threaded VDB union with fillet that produces a new grid or tree from +/// immutable inputs. +/// +/// @param lhs Level-set grid to be combined with a second input. +/// +/// @param rhs Level-set grid to be combined with the first input. +/// +/// @param mask Optional float grid that controls the strength of the blending. +/// +/// @param alpha Controls the blend-radius. +/// +/// @param beta Controls the exponent used to control the strength of the blend. +/// +/// @param gamma Controls the strength of the blend by providing a multiplier component. +/// +/// @return The filleted union of the @lhs and @rhs level set inputs. +/// +/// @throw If the transforms of @lsh, @rhs, and @mask do not match. +template::Type> +typename GridT::Ptr +unionFillet(const GridT& lhs, + const GridT& rhs, + typename MaskT::ConstPtr mask, + typename GridT::ValueType alpha, + typename GridT::ValueType beta, + typename GridT::ValueType gamma); + +// +// Main class to handle UnionWithFillet +// +template::Type> +struct UnionWithFillet { + using TreeT = typename GridT::TreeType; + using ValueType = typename TreeT::ValueType; + using TreePtrType = typename TreeT::Ptr; + using LeafNodeType = typename TreeT::LeafNodeType; + using NodeMaskType = typename LeafNodeType::NodeMaskType; + using RootNodeType = typename TreeT::RootNodeType; + using NodeChainType = typename RootNodeType::NodeChainType; + using InternalNodeType = typename NodeChainType::template Get<1>; + using MaskTreeType = typename MaskT::TreeType; + using MaskValueType = typename MaskT::ValueType; + + UnionWithFillet(const GridT& lhsGrid, + const GridT& rhsGrid, + typename MaskT::ConstPtr mask, + const ValueType& bandwidth, + const ValueType& exponent, + const ValueType& multiplier) + : mLhsGrid(&lhsGrid) + , mRhsGrid(&rhsGrid) + , mLhsTree(&(lhsGrid.tree())) + , mRhsTree(&(rhsGrid.tree())) + , mBandwidth(bandwidth) + , mExponent(exponent) + , mMultiplier(multiplier) + { + mMaskTree = mask ? mask->treePtr() : nullptr; + static_assert(std::is_floating_point::value, + "assert in UnionFillet Constructor: " + "level set grids must have scalar/floating-point value types."); + } + + /// @brief Perform blending of two level-sets. + typename GridT::Ptr blend(); + + TreePtrType& segment() { return mSegment; } + +private: + struct FilletParms { + float mAlpha; // Bandwidth Falloff + float mBeta; // Exponent + float mGamma; // Amplitude/Multiplier + }; + + struct BuildPrimarySegment; + struct BuildSecondarySegment; + + TreePtrType mSegment; + GridT const * const mLhsGrid; + GridT const * const mRhsGrid; + TreeT const * const mLhsTree; + TreeT const * const mRhsTree; + typename MaskT::TreeType::ConstPtr mMaskTree = nullptr; + ValueType mBandwidth; + ValueType mExponent; + ValueType mMultiplier; +}; + + +/// @cond OPENVDB_DOCS_INTERNAL +/// @brief Go through the lhs nodes (both internal and leaf nodes). +// TODO: should just be only for float +template +struct UnionWithFillet::BuildPrimarySegment { + using MaskTreeType = typename MaskT::TreeType; + using ValueType = typename TreeT::ValueType; + using TreePtrType = typename TreeT::Ptr; + using LeafNodeType = typename TreeT::LeafNodeType; + using NodeMaskType = typename LeafNodeType::NodeMaskType; + using RootNodeType = typename TreeT::RootNodeType; + using NodeChainType = typename RootNodeType::NodeChainType; + using InternalNodeType = typename NodeChainType::template Get<1>; + + BuildPrimarySegment(TreeT const * const lhs, + TreeT const * const rhs, + typename MaskT::TreeType::ConstPtr mask, + FilletParms parms) + : mSegment(new TreeT(lhs->background())) + , mLhsTree(lhs) + , mRhsTree(rhs) + , mMaskTree(mask) + , mParms(parms) { } + + void operator()() const + { + std::vector leafNodes; + + { + std::vector internalNodes; + mLhsTree->getNodes(internalNodes); + + ProcessInternalNodes op(internalNodes, *mRhsTree, *mSegment, leafNodes, mParms); + tbb::parallel_reduce(tbb::blocked_range(0, internalNodes.size()), op); + } + if (mMaskTree) { + ProcessLeafNodesMask op(leafNodes, *mRhsTree, mMaskTree, *mSegment, mParms); + tbb::parallel_reduce(tbb::blocked_range(0, leafNodes.size()), op); + } else { + ProcessLeafNodes op(leafNodes, *mRhsTree, *mSegment, mParms); + tbb::parallel_reduce(tbb::blocked_range(0, leafNodes.size()), op); + } + } + + TreePtrType& segment() { return mSegment; } + +private: + + struct ProcessInternalNodes { + + ProcessInternalNodes(std::vector& lhsNodes, + const TreeT& rhsTree, TreeT& outputTree, + std::vector& outputLeafNodes, FilletParms parms) + : mLhsNodes(lhsNodes) + , mRhsTree(&rhsTree) + , mLocalTree(mRhsTree->background()) + , mOutputTree(&outputTree) + , mLocalLeafNodes() + , mOutputLeafNodes(&outputLeafNodes) + , mParms(parms) + { } + + ProcessInternalNodes(ProcessInternalNodes& other, tbb::split) + : mLhsNodes(other.mLhsNodes) + , mRhsTree(other.mRhsTree) + , mLocalTree(mRhsTree->background()) + , mOutputTree(&mLocalTree) + , mLocalLeafNodes() + , mOutputLeafNodes(&mLocalLeafNodes) + , mParms(other.mParms) + { } + + void join(ProcessInternalNodes& other) + { + mOutputTree->merge(*other.mOutputTree); + mOutputLeafNodes->insert(mOutputLeafNodes->end(), + other.mOutputLeafNodes->begin(), other.mOutputLeafNodes->end()); + } + + void operator()(const tbb::blocked_range& range) + { + openvdb::tree::ValueAccessor rhsAcc(*mRhsTree); + openvdb::tree::ValueAccessor outputAcc(*mOutputTree); + + std::vector tmpLeafNodes; + + for (size_t n = range.begin(), N = range.end(); n < N; ++n) { + + const InternalNodeType& lhsNode = *mLhsNodes[n]; + const openvdb::math::Coord& ijk = lhsNode.origin(); + const InternalNodeType * rhsNode = + rhsAcc.template probeConstNode(ijk); + + if (rhsNode) { + lhsNode.getNodes(*mOutputLeafNodes); + } else { + if (!(rhsAcc.getValue(ijk) < ValueType(0))) { + tmpLeafNodes.clear(); + lhsNode.getNodes(tmpLeafNodes); + for (const LeafNodeType* leaf : tmpLeafNodes) { + outputAcc.addLeaf(new LeafNodeType(*leaf)); + } + } + } + } // end range loop + } + + const std::vector& mLhsNodes; + TreeT const * const mRhsTree; + TreeT mLocalTree; + TreeT * const mOutputTree; + + std::vector mLocalLeafNodes; + std::vector * const mOutputLeafNodes; + FilletParms mParms; + }; // struct ProcessInternalNodes + + struct ProcessLeafNodes { + + ProcessLeafNodes(std::vector& lhsNodes, + const TreeT& rhsTree, TreeT& output, FilletParms parms) + : mLhsNodes(lhsNodes) + , mRhsTree(&rhsTree) + , mLocalTree(mRhsTree->background()) + , mOutputTree(&output) + , mParms(parms) + { } + + ProcessLeafNodes(ProcessLeafNodes& other, tbb::split) + : mLhsNodes(other.mLhsNodes) + , mRhsTree(other.mRhsTree) + , mLocalTree(mRhsTree->background()) + , mOutputTree(&mLocalTree) + , mParms(other.mParms) + { } + + void join(ProcessLeafNodes& rhs) { mOutputTree->merge(*rhs.mOutputTree); } + + void operator()(const tbb::blocked_range& range) + { + openvdb::tree::ValueAccessor rhsAcc(*mRhsTree); + openvdb::tree::ValueAccessor outputAcc(*mOutputTree); + + for (size_t n = range.begin(), N = range.end(); n < N; ++n) { + + const LeafNodeType& lhsNode = *mLhsNodes[n]; + const openvdb::math::Coord& ijk = lhsNode.origin(); + + const LeafNodeType* rhsNodePt = rhsAcc.probeConstLeaf(ijk); + + if (rhsNodePt) { // combine overlapping nodes + + LeafNodeType* outputNode = outputAcc.touchLeaf(ijk); + ValueType * outputData = outputNode->buffer().data(); + NodeMaskType& outputMask = outputNode->getValueMask(); + + const ValueType * lhsData = lhsNode.buffer().data(); + const NodeMaskType& lhsMask = lhsNode.getValueMask(); + + const ValueType * rhsData = rhsNodePt->buffer().data(); + const NodeMaskType& rhsMask = rhsNodePt->getValueMask(); + + const float alpha = mParms.mAlpha; + const float beta = mParms.mBeta; + const float gamma = mParms.mGamma; + + for (openvdb::Index pos = 0; pos < LeafNodeType::SIZE; ++pos) { + const float A = lhsData[pos]; + const float B = rhsData[pos]; + const float m = openvdb::math::Clamp((alpha - A) / alpha, 0.f, 1.f) * + openvdb::math::Clamp((alpha - B) / alpha, 0.f, 1.f); + const float offset = openvdb::math::Pow(m, beta) * gamma; + const bool isAMin = A < B; + // - sign here. Otherwise, we'll get empty space between geos + outputData[pos] = isAMin ? A - offset : B - offset; + outputMask.set(pos, isAMin ? lhsMask.isOn(pos) : rhsMask.isOn(pos)); + } + } else { + if (!(rhsAcc.getValue(ijk) < ValueType(0.0))) { + outputAcc.addLeaf(new LeafNodeType(lhsNode)); + } + } + } // end range loop + } + + const std::vector& mLhsNodes; + TreeT const * const mRhsTree; + TreeT mLocalTree; + TreeT * const mOutputTree; + FilletParms mParms; + }; // struct ProcessLeafNodes + + struct ProcessLeafNodesMask { + using MaskTreeType = typename MaskT::TreeType; + using MaskLeafNodeType = typename MaskT::TreeType::LeafNodeType; + using MaskNodeMaskType = typename MaskLeafNodeType::NodeMaskType; + using MaskValueType = typename MaskT::ValueType; + + ProcessLeafNodesMask(std::vector& lhsNodes, + const TreeT& rhsTree, typename MaskT::TreeType::ConstPtr maskTree, TreeT& output, FilletParms parms) + : mLhsNodes(lhsNodes) + , mRhsTree(&rhsTree) + , mLocalTree(mRhsTree->background()) + , mMaskTree(maskTree) + , mOutputTree(&output) + , mParms(parms) + { } + + ProcessLeafNodesMask(ProcessLeafNodesMask& other, tbb::split) + : mLhsNodes(other.mLhsNodes) + , mRhsTree(other.mRhsTree) + , mLocalTree(mRhsTree->background()) + , mMaskTree(other.mMaskTree) + , mOutputTree(&mLocalTree) + , mParms(other.mParms) + { } + + void join(ProcessLeafNodesMask& rhs) { mOutputTree->merge(*rhs.mOutputTree); } + + void operator()(const tbb::blocked_range& range) + { + openvdb::tree::ValueAccessor rhsAcc(*mRhsTree); + openvdb::tree::ValueAccessor outputAcc(*mOutputTree); + openvdb::tree::ValueAccessor maskAcc(*mMaskTree); + const float maskBackground = mMaskTree->background(); + + for (size_t n = range.begin(), N = range.end(); n < N; ++n) { + + const LeafNodeType& lhsNode = *mLhsNodes[n]; + const openvdb::math::Coord& ijk = lhsNode.origin(); + + const LeafNodeType* rhsNodePt = rhsAcc.probeConstLeaf(ijk); + const MaskLeafNodeType* maskNodePt = maskAcc.probeConstLeaf(ijk); + + if (rhsNodePt) { // combine overlapping nodes + LeafNodeType* outputNode = outputAcc.touchLeaf(ijk); + ValueType * outputData = outputNode->buffer().data(); + NodeMaskType& outputMask = outputNode->getValueMask(); + + const ValueType* lhsData = lhsNode.buffer().data(); + const NodeMaskType& lhsMask = lhsNode.getValueMask(); + + const ValueType* rhsData = rhsNodePt->buffer().data(); + const NodeMaskType& rhsMask = rhsNodePt->getValueMask(); + + const float alpha = mParms.mAlpha; + const float beta = mParms.mBeta; + const float gamma = mParms.mGamma; + + MaskValueType const * maskData = nullptr; + if (maskNodePt) maskData = maskNodePt->buffer().data(); + + for (openvdb::Index pos = 0; pos < LeafNodeType::SIZE; ++pos) { + const float A = lhsData[pos]; + const float B = rhsData[pos]; + const float m = openvdb::math::Clamp((alpha - A) / alpha, 0.f, 1.f) * + openvdb::math::Clamp((alpha - B) / alpha, 0.f, 1.f); + const float offset = openvdb::math::Pow(m, beta) * gamma; + const bool isAMin = A < B; + // multiply by mask value if it exists + const float multiplier = maskData ? maskData[pos] : maskBackground; + // - sign here. Otherwise, we'll get empty space between geos + outputData[pos] = isAMin ? (A - offset * multiplier) : (B - offset * multiplier); + outputMask.set(pos, isAMin ? lhsMask.isOn(pos) : rhsMask.isOn(pos)); + } + } else { + if (!(rhsAcc.getValue(ijk) < ValueType(0.0))) { + outputAcc.addLeaf(new LeafNodeType(lhsNode)); + } + } + } // end range loop + } + + const std::vector& mLhsNodes; + TreeT const * const mRhsTree; + TreeT mLocalTree; + typename MaskT::TreeType::ConstPtr mMaskTree; + TreeT * const mOutputTree; + FilletParms mParms; + }; // struct ProcessLeafNodesMask + + +private: + TreePtrType mSegment; + TreeT const * const mLhsTree; + TreeT const * const mRhsTree; + typename MaskT::TreeType::ConstPtr mMaskTree; + FilletParms mParms; +}; + + +template +struct UnionWithFillet::BuildSecondarySegment { + using MaskTreeType = typename MaskT::TreeType; + using ValueType = typename TreeT::ValueType; + using TreePtrType = typename TreeT::Ptr; + using LeafNodeType = typename TreeT::LeafNodeType; + using NodeMaskType = typename LeafNodeType::NodeMaskType; + using RootNodeType = typename TreeT::RootNodeType; + using NodeChainType = typename RootNodeType::NodeChainType; + using InternalNodeType = typename NodeChainType::template Get<1>; + using MaskValueType = typename MaskTreeType::ValueType; + + BuildSecondarySegment(TreeT const * const lhs, + TreeT const * const rhs, + FilletParms parms) + : mSegment(new TreeT(lhs->background())) + , mLhsTree(lhs) + , mRhsTree(rhs) + , mParms(parms) { } + + + void operator()() const + { + std::vector leafNodes; + + { + std::vector internalNodes; + mRhsTree->getNodes(internalNodes); + + ProcessInternalNodes op(internalNodes, *mLhsTree, *mSegment, leafNodes, mParms); + tbb::parallel_reduce(tbb::blocked_range(0, internalNodes.size()), op); + } + + ProcessLeafNodes op(leafNodes, *mLhsTree, *mSegment, mParms); + tbb::parallel_reduce(tbb::blocked_range(0, leafNodes.size()), op); + } + + TreePtrType& segment() { return mSegment; } + +private: + + struct ProcessInternalNodes { + + ProcessInternalNodes(std::vector& rhsNodes, + const TreeT& lhsTree, TreeT& outputTree, + std::vector& outputLeafNodes, + FilletParms parms) + : mRhsNodes(rhsNodes) + , mLhsTree(&lhsTree) + , mLocalTree(mLhsTree->background()) + , mOutputTree(&outputTree) + , mLocalLeafNodes() + , mOutputLeafNodes(&outputLeafNodes) + , mParms(parms) + { } + + ProcessInternalNodes(ProcessInternalNodes& other, tbb::split) + : mRhsNodes(other.mRhsNodes) + , mLhsTree(other.mLhsTree) + , mLocalTree(mLhsTree->background()) + , mOutputTree(&mLocalTree) + , mLocalLeafNodes() + , mOutputLeafNodes(&mLocalLeafNodes) + , mParms(other.mParms) + { } + + void join(ProcessInternalNodes& other) + { + mOutputTree->merge(*other.mOutputTree); + mOutputLeafNodes->insert(mOutputLeafNodes->end(), + other.mOutputLeafNodes->begin(), other.mOutputLeafNodes->end()); + } + + void operator()(const tbb::blocked_range& range) + { + openvdb::tree::ValueAccessor lhsAcc(*mLhsTree); + openvdb::tree::ValueAccessor outputAcc(*mOutputTree); + + std::vector tmpLeafNodes; + + for (size_t n = range.begin(), N = range.end(); n < N; ++n) { + + const InternalNodeType& rhsNode = *mRhsNodes[n]; + const openvdb::math::Coord& ijk = rhsNode.origin(); + const InternalNodeType * lhsNode = + lhsAcc.template probeConstNode(ijk); + + if (lhsNode) { + rhsNode.getNodes(*mOutputLeafNodes); + } else { + if (!(lhsAcc.getValue(ijk) < ValueType(0))) { + tmpLeafNodes.clear(); + rhsNode.getNodes(tmpLeafNodes); + for (const LeafNodeType* leaf : tmpLeafNodes) { + outputAcc.addLeaf(new LeafNodeType(*leaf)); + } + } + } + } // end range loop + } + + const std::vector& mRhsNodes; + TreeT const * const mLhsTree; + TreeT mLocalTree; + TreeT * const mOutputTree; + + std::vector mLocalLeafNodes; + std::vector * const mOutputLeafNodes; + FilletParms mParms; + }; // struct ProcessInternalNodes + + struct ProcessLeafNodes { + + ProcessLeafNodes(std::vector& rhsNodes, + const TreeT& lhsTree, TreeT& output, + FilletParms parms) + : mRhsNodes(rhsNodes) + , mLhsTree(&lhsTree) + , mLocalTree(mLhsTree->background()) + , mOutputTree(&output) + , mParms(parms) + { } + + ProcessLeafNodes(ProcessLeafNodes& rhs, tbb::split) + : mRhsNodes(rhs.mRhsNodes) + , mLhsTree(rhs.mLhsTree) + , mLocalTree(mLhsTree->background()) + , mOutputTree(&mLocalTree) + , mParms(rhs.mParms) + { } + + void join(ProcessLeafNodes& rhs) { mOutputTree->merge(*rhs.mOutputTree); } + + void operator()(const tbb::blocked_range& range) + { + openvdb::tree::ValueAccessor lhsAcc(*mLhsTree); + openvdb::tree::ValueAccessor outputAcc(*mOutputTree); + + for (size_t n = range.begin(), N = range.end(); n < N; ++n) { + + const LeafNodeType& rhsNode = *mRhsNodes[n]; + const openvdb::math::Coord& ijk = rhsNode.origin(); + + const LeafNodeType* lhsNode = lhsAcc.probeConstLeaf(ijk); + + if (!lhsNode) { + if (!(lhsAcc.getValue(ijk) < ValueType(0))) { + outputAcc.addLeaf(new LeafNodeType(rhsNode)); + } + } + } // end range loop + } + + const std::vector& mRhsNodes; + TreeT const * const mLhsTree; + TreeT mLocalTree; + TreeT * const mOutputTree; + FilletParms mParms; + }; // struct ProcessLeafNodes + + +private: + TreePtrType mSegment; + TreeT const * const mLhsTree; + TreeT const * const mRhsTree; + FilletParms mParms; +}; + + +template +typename GridT::Ptr UnionWithFillet::blend() { + FilletParms parms; + parms.mAlpha = mBandwidth; + parms.mBeta = mExponent; + parms.mGamma = mMultiplier; + BuildPrimarySegment primary(mLhsTree, mRhsTree, mMaskTree, parms); + BuildSecondarySegment secondary(mLhsTree, mRhsTree, parms); + + // Exploiting nested parallelism + tbb::task_group tasks; + tasks.run(primary); + tasks.run(secondary); + tasks.wait(); + + primary.segment()->merge(*secondary.segment()); + + // The leafnode (level = 0) sign is set in the segment construction. + openvdb::tools::signedFloodFill(*primary.segment(), /*threaded=*/true, /*grainSize=*/1, /*minLevel=*/1); + + mSegment = primary.segment(); + + typename GridT::Ptr ret = GridT::create(mSegment); + ret->setTransform((mRhsGrid->transform()).copy()); + ret->setGridClass(openvdb::GRID_LEVEL_SET); + return ret; +} + + +template::Type> +typename GridT::Ptr +unionFillet(const GridT& lhs, + const GridT& rhs, + typename MaskT::ConstPtr mask, + typename GridT::ValueType alpha, + typename GridT::ValueType beta, + typename GridT::ValueType gamma) +{ + static_assert(std::is_floating_point::value, + "assert in unionFillet: " + "level set grids must have scalar/floating-point value types."); + + // sanitizer + const openvdb::math::Transform& lhsXform = lhs.constTransform(); + const openvdb::math::Transform& rhsXform = rhs.constTransform(); + if (lhsXform != rhsXform) throw std::runtime_error("The two grids need to have the same transforms."); + if (mask) { + const openvdb::math::Transform& maskXform = mask->constTransform(); + if (lhsXform != maskXform) throw std::runtime_error("The grids and the mask need to have the same transforms."); + } + + UnionWithFillet uf(lhs, rhs, mask, alpha, beta, gamma); + return uf.blend(); +} +} // tools +} // OPENVDB_VERSION_NAME +} // openvdb + +#endif // OPENVDB_TOOLS_BLEND_HAS_BEEN_INCLUDED + diff --git a/openvdb_houdini/openvdb_houdini/CMakeLists.txt b/openvdb_houdini/openvdb_houdini/CMakeLists.txt index e12e84d083..eec9321ec0 100644 --- a/openvdb_houdini/openvdb_houdini/CMakeLists.txt +++ b/openvdb_houdini/openvdb_houdini/CMakeLists.txt @@ -242,6 +242,7 @@ if(NOT OPENVDB_DSO_NAMES) SOP_OpenVDB_Diagnostics SOP_OpenVDB_Extrapolate SOP_OpenVDB_Fill + SOP_OpenVDB_Fillet SOP_OpenVDB_Filter SOP_OpenVDB_Filter_Level_Set SOP_OpenVDB_Fracture diff --git a/openvdb_houdini/openvdb_houdini/SOP_OpenVDB_Fillet.cc b/openvdb_houdini/openvdb_houdini/SOP_OpenVDB_Fillet.cc new file mode 100644 index 0000000000..67e53b3f2f --- /dev/null +++ b/openvdb_houdini/openvdb_houdini/SOP_OpenVDB_Fillet.cc @@ -0,0 +1,649 @@ +// Copyright Contributors to the OpenVDB Project +// SPDX-License-Identifier: MPL-2.0 +// +// @file SOP_OpenVDB_Fillet.cc +// +// @author FX R&D OpenVDB team +// +// @brief OpenVDB Fillet to combine two level-sets together. + +#include +#include +#include + +#include +#include +#include // for isExactlyEqual() +#include +#include +#include // for resampleToMatch() +#include // for levelSetRebuild() +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace hvdb = openvdb_houdini; +namespace hutil = houdini_utils; + +namespace { + +// +// Resampling options +// +enum ResampleMode { + RESAMPLE_OFF, // don't auto-resample grids + RESAMPLE_B, // resample B to match A + RESAMPLE_A, // resample A to match B + RESAMPLE_HI_RES, // resample higher-res grid to match lower-res + RESAMPLE_LO_RES // resample lower-res grid to match higher-res +}; + + +enum { RESAMPLE_MODE_FIRST = RESAMPLE_OFF, RESAMPLE_MODE_LAST = RESAMPLE_LO_RES }; + + +// +// VDBFilletParms struct to be used in SOP +// +struct VDBFilletParms { + float mAlpha = 10.f; // Falloff + float mBeta = 100.f; // Exponent + float mGamma = 10.f; // Amplitude/Multiplier + ResampleMode mResampleMode = ResampleMode::RESAMPLE_OFF; // Resample mode + int mSamplingOrder = 0; + hvdb::GridCPtr mABaseGrid = nullptr; + hvdb::GridCPtr mBBaseGrid = nullptr; + openvdb::FloatGrid::ConstPtr mMaskPtr = nullptr; +}; + +// Helper function to compare both scalar and vector values +template +static bool approxEq(const ValueT& a, const ValueT& b) { + return openvdb::math::isRelOrApproxEqual(a, b, ValueT(1e-6f), ValueT(1e-8f)); +} + +const char* const sResampleModeMenuItems[] = { + "off", "Off", + "btoa", "B to Match A", + "atob", "A to Match B", + "hitolo", "Higher-res to Match Lower-res", + "lotohi", "Lower-res to Match Higher-res", + nullptr +}; + +inline ResampleMode +asResampleMode(exint i, ResampleMode defaultMode = RESAMPLE_B) +{ + return (i >= RESAMPLE_MODE_FIRST && i <= RESAMPLE_MODE_LAST) + ? static_cast(i) : defaultMode; +} + +using StringVec = std::vector; + +// Split a string into group patterns separated by whitespace. +// For example, given '@name=d* @id="1 2" {grp1 grp2}', return +// ['@name=d*', '@id="1 2"', '{grp1 grp2}']. +// (This is nonstandard. Normally, multiple patterns are unioned +// to define a single group.) +// Nesting of quotes and braces is not supported. +inline StringVec +splitPatterns(const std::string& str) +{ + StringVec patterns; + bool quoted = false, braced = false; + std::string pattern; + for (const auto c: str) { + if (isspace(c)) { + if (pattern.empty()) continue; // skip whitespace between patterns + if (quoted || braced) { + pattern.push_back(c); // keep whitespace within quotes or braces + } else { + // At the end of a pattern. Start a new pattern. + patterns.push_back(pattern); + pattern.clear(); + quoted = braced = false; + } + } else { + switch (c) { + case '"': quoted = !quoted; break; + case '{': braced = true; break; + case '}': braced = false; break; + default: break; + } + pattern.push_back(c); + } + } + if (!pattern.empty()) { patterns.push_back(pattern); } // add the final pattern + + // If no patterns were found, add an empty pattern, which matches everything. + if (patterns.empty()) { patterns.push_back(""); } + + return patterns; +} + +inline UT_String +getGridName(const GU_PrimVDB* vdb, const UT_String& defaultName = "") +{ + UT_String name{UT_String::ALWAYS_DEEP}; + if (vdb != nullptr) { + name = vdb->getGridName(); + if (!name.isstring()) name = defaultName; + } + return name; +} + +template +struct ResampleOp { + ResampleOp(typename GridT::ConstPtr &aGrid, typename GridT::ConstPtr &bGrid, const VDBFilletParms &parms) + : mAGrid(aGrid) + , mBGrid(bGrid) + , mParms(parms) + { } + + // Function to return a scalar grid's background value as a double + static double getScalarBackgroundValue(const hvdb::Grid& baseGrid) + { + double bg = 0.0; + baseGrid.apply([&bg](const auto &grid) { + bg = static_cast(grid.background()); + }); + return bg; + } + + typename GridT::Ptr resampleToMatch(const GridT& src, const hvdb::Grid& ref, int order) + { + using ValueT = typename GridT::ValueType; + const ValueT ZERO = openvdb::zeroVal(); + + const openvdb::math::Transform& refXform = ref.constTransform(); + + typename GridT::Ptr dest; + if (src.getGridClass() == openvdb::GRID_LEVEL_SET) { + // For level set grids, use the level set rebuild tool to both resample the + // source grid to match the reference grid and to rebuild the resulting level set. + const bool refIsLevelSet = ref.getGridClass() == openvdb::GRID_LEVEL_SET; + OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN + const ValueT halfWidth = refIsLevelSet + ? ValueT(ZERO + getScalarBackgroundValue(ref) * (1.0 / ref.voxelSize()[0])) + : ValueT(src.background() * (1.0 / src.voxelSize()[0])); + OPENVDB_NO_TYPE_CONVERSION_WARNING_END + + if (!openvdb::math::isFinite(halfWidth)) { + std::stringstream msg; + msg << "Resample to match: Illegal narrow band width = " << halfWidth + << ", caused by grid '" << ref.getName() << "' with background " + << getScalarBackgroundValue(ref); + throw std::invalid_argument(msg.str()); + } + + try { + dest = openvdb::tools::levelSetRebuild(src, /*iso=*/ZERO, + /*exWidth=*/halfWidth, /*inWidth=*/halfWidth, &refXform/*, &mInterrupt.interrupter*/); + } catch (openvdb::TypeError&) { + dest.reset(); + std::stringstream msg; + msg << "skipped rebuild of level set grid " << src.getName() << " of type " << src.type(); + throw std::runtime_error(msg.str().c_str()); + } + } + if (!dest && src.constTransform() != refXform) { + // If level set rebuild failed due to an unsupported grid type, + // use the grid transformer tool to resample the source grid to match + // the reference grid. + dest = src.copyWithNewTree(); + dest->setTransform(refXform.copy()); + using namespace openvdb; + switch (order) { + case 0: tools::resampleToMatch(src, *dest/*, mInterrupt.interrupter()*/); break; + case 1: tools::resampleToMatch(src, *dest/*, mInterrupt.interrupter()*/); break; + case 2: tools::resampleToMatch(src, *dest/*, mInterrupt.interrupter()*/); break; + // note: no default case because sampling order is guaranteed to be 0, 1, or 2 in evalParms. + } + } + return dest; + } + + void resampleMask(typename GridT::ConstPtr& mask) + { + if (!mask) return; + const openvdb::math::Transform& maskXform = mask->constTransform(); + const openvdb::math::Transform& aGrdXform = mAGrid->constTransform(); + if (maskXform != aGrdXform) { + typename GridT::ConstPtr maskRsmpl = this->resampleToMatch(*mask /* src */, *mAGrid /* ref */, mParms.mSamplingOrder); + mask = maskRsmpl; + } + } + + void resampleGrids() + { + if (!mAGrid || !mBGrid) return; + + // One of RESAMPLE_A, RESAMPLE_B or RESAMPLE_OFF, specifying whether + // grid A, grid B or neither grid was resampled + int resampleWhich = RESAMPLE_OFF; + + // Determine which of the two grids should be resampled. + if (mParms.mResampleMode == RESAMPLE_HI_RES || mParms.mResampleMode == RESAMPLE_LO_RES) { + const openvdb::Vec3d + aVoxSize = mAGrid->voxelSize(), + bVoxSize = mBGrid->voxelSize(); + const double + aVoxVol = aVoxSize[0] * aVoxSize[1] * aVoxSize[2], + bVoxVol = bVoxSize[0] * bVoxSize[1] * bVoxSize[2]; + resampleWhich = ((aVoxVol > bVoxVol && mParms.mResampleMode == RESAMPLE_LO_RES) + || (aVoxVol < bVoxVol && mParms.mResampleMode == RESAMPLE_HI_RES)) + ? RESAMPLE_A : RESAMPLE_B; + } else { + resampleWhich = mParms.mResampleMode; + } + + if (mAGrid->constTransform() != mBGrid->constTransform()) { + // If the A and B grid transforms don't match, one of the grids + // should be resampled into the other's index space. + if (mParms.mResampleMode == RESAMPLE_OFF) { + // Resampling is disabled. Just log a warning. + std::ostringstream msg; + msg << mAGrid->getName() << " and " << mBGrid->getName() << " transforms don't match"; + throw std::runtime_error(msg.str().c_str()); + } else { + if (resampleWhich == RESAMPLE_A) { + // Resample grid A into grid B's index space. + mAGrid = this->resampleToMatch(*mAGrid, *mBGrid, mParms.mSamplingOrder); + } else if (resampleWhich == RESAMPLE_B) { + // Resample grid B into grid A's index space. + mBGrid = this->resampleToMatch(*mBGrid, *mAGrid, mParms.mSamplingOrder); + } + } + } + + // If both grids are level sets, ensure that their background values match. + // (If one of the grids was resampled, then the background values should + // already match.) + if (mAGrid->getGridClass() == openvdb::GRID_LEVEL_SET && + mBGrid->getGridClass() == openvdb::GRID_LEVEL_SET) + { + const double + a = this->getScalarBackgroundValue(*mAGrid), + b = this->getScalarBackgroundValue(*mBGrid); + if (!approxEq(a, b)) { + if (mParms.mResampleMode == RESAMPLE_OFF) { + // Resampling/rebuilding is disabled. Just log a warning. + std::ostringstream msg; + msg << mAGrid->getName() << " and " << mBGrid->getName() + << " background values don't match (" + << std::setprecision(3) << a << " vs. " << b << ");\n" + << " the output grid will not be a valid level set"; + throw std::runtime_error(msg.str().c_str()); + } else { + // One of the two grids needs a level set rebuild. + if (resampleWhich == RESAMPLE_A) { + // Rebuild A to match B's background value. + mAGrid = this->resampleToMatch(*mAGrid, *mBGrid, mParms.mSamplingOrder); + } else if (resampleWhich == RESAMPLE_B) { + // Rebuild B to match A's background value. + mBGrid = this->resampleToMatch(*mBGrid, *mAGrid, mParms.mSamplingOrder); + } + } + } + } + } + + typename GridT::ConstPtr &mAGrid; + typename GridT::ConstPtr &mBGrid; + const VDBFilletParms mParms; + ResampleMode mResampleMode; + int mSamplingOrder; +}; +} + +class SOP_OpenVDB_Fillet: public hvdb::SOP_NodeVDB +{ +public: + SOP_OpenVDB_Fillet(OP_Network*, const char* name, OP_Operator*); + ~SOP_OpenVDB_Fillet() override {} + + static OP_Node* factory(OP_Network*, const char* name, OP_Operator*); + + class Cache: public SOP_VDBCacheOptions + { + protected: + OP_ERROR cookVDBSop(OP_Context&) override; + OP_ERROR evalParms(OP_Context&, VDBFilletParms&); + private: + hvdb::GridPtr blendLevelSets( + hvdb::GridCPtr aGrid, + hvdb::GridCPtr bGrid, + VDBFilletParms parm); + }; // class Cache + + // Return true for a given input if the connector to the input + // should be drawn dashed rather than solid. + int isRefInput(unsigned idx) const override { return (idx == 1); } + +protected: + unsigned disableParms() override; +}; + + +//////////////////////////////////////// + + +OP_ERROR +SOP_OpenVDB_Fillet::Cache::evalParms(OP_Context& context, VDBFilletParms& parms) +{ + const fpreal time = context.getTime(); + + parms.mAlpha = static_cast(evalFloat("alpha", 0, time)); + parms.mBeta = static_cast(evalFloat("beta", 0, time)); + parms.mGamma = static_cast(evalFloat("gamma", 0, time)); + parms.mResampleMode = asResampleMode(evalInt("resample", 0, time)); + parms.mSamplingOrder = static_cast(evalInt("resampleinterp", 0, time)); + + if (parms.mSamplingOrder < 0 || parms.mSamplingOrder > 2) { + addWarning(SOP_MESSAGE, "Sampling order should be 0, 1, or 2."); + return UT_ERROR_ABORT; + } + + // mask + int useMask = static_cast(evalInt("mask", 0, time)); + + openvdb::FloatGrid::ConstPtr maskGrid; + + if (useMask && this->nInputs() == 3) { + const GU_Detail* maskGeo = inputGeo(2); + const auto maskName = evalStdString("maskname", time); + if (!maskGeo) { + addWarning(SOP_MESSAGE, "The mask input is empty."); + return UT_ERROR_ABORT; + } + + const GA_PrimitiveGroup* maskGroup = parsePrimitiveGroups(maskName.c_str(), GroupCreator(maskGeo)); + if (!maskGroup && !maskName.empty()) { + addWarning(SOP_MESSAGE, "Mask not found."); + return UT_ERROR_ABORT; + } + + hvdb::VdbPrimCIterator maskIt(maskGeo, maskGroup); + if (!maskIt) { + addWarning(SOP_MESSAGE, "The mask input is empty."); + return UT_ERROR_ABORT; + } else if (maskIt->getStorageType() != UT_VDB_FLOAT) { + addWarning(SOP_MESSAGE, "The mask grid has to be a Float Grid."); + return UT_ERROR_ABORT; + } + + maskGrid = openvdb::gridConstPtrCast(maskIt->getGridPtr()); + if (!maskGrid) { + addWarning(SOP_MESSAGE, "The mask grid has to be a Float Grid."); + return UT_ERROR_ABORT; + } + + parms.mMaskPtr = maskGrid; + + } else if (useMask && this->nInputs() < 3) { + addWarning(SOP_MESSAGE, "Need a mask grid as a third input to SOP."); + return UT_ERROR_ABORT; + } // if not using mask, do nothing. + + return error(); +} + + +//////////////////////////////////////// + + +// Build UI and register this operator. +void +newSopOperator(OP_OperatorTable* table) +{ + if (table == nullptr) return; + + hutil::ParmList parms; + + // Group A + parms.add(hutil::ParmFactory(PRM_STRING, "agroup", "Group A") + .setChoiceList(&hutil::PrimGroupMenuInput1) + .setTooltip("Use a subset of the first input as the A VDB(s).") + .setDocumentation( + "The VDBs to be used from the first input" + " (see [specifying volumes|/model/volumes#group])")); + + // Group B + parms.add(hutil::ParmFactory(PRM_STRING, "bgroup", "Group B") + .setChoiceList(&hutil::PrimGroupMenuInput2) + .setTooltip("Use a subset of the second input as the B VDB(s).") + .setDocumentation( + "The VDBs to be used from the second input" + " (see [specifying volumes|/model/volumes#group])")); + + // Mask multiplier + parms.add(hutil::ParmFactory(PRM_TOGGLE, "mask", "") + .setDefault(PRMzeroDefaults) + .setTypeExtended(PRM_TYPE_TOGGLE_JOIN) + .setTooltip("Enable / disable the mask.")); + + parms.add(hutil::ParmFactory(PRM_STRING, "maskname", "Alpha Mask") + .setChoiceList(&hutil::PrimGroupMenuInput3) + .setTooltip("Optional scalar VDB used for alpha masking\n\n" + "Values are assumed to be between 0 and 1.")); + + // Band radius of influence / falloff width, i.e. alpha + parms.add(hutil::ParmFactory(PRM_FLT_J, "alpha", "Band Radius") + .setDefault(10.f) + .setRange(PRM_RANGE_UI, 0.f, PRM_RANGE_UI, 1000.f) + .setTooltip( + "Band radius of influence measures the distance from the zero\n" + "iso-contour of the intersection that is going to be modified.\n" + "This is measured in world-space.")); + + // Exponent, i.e. beta + parms.add(hutil::ParmFactory(PRM_FLT_J, "beta", "Exponent") + .setDefault(100.f) + .setRange(PRM_RANGE_UI, 0.f, PRM_RANGE_UI, 1000.f) + .setTooltip( + "Blending curve exponential used in the model.")); + + // Amplitude + parms.add(hutil::ParmFactory(PRM_FLT_J, "gamma", "Multiplier") + .setDefault(10.f) + .setRange(PRM_RANGE_UI, 0.f, PRM_RANGE_UI, 1000.f) + .setTooltip( + "Amplitude provides a multiplier to make the blended\n" + "influence weaker or stronger.")); + + // Menu of resampling options + parms.add(hutil::ParmFactory(PRM_ORD, "resample", "Resample") + .setDefault(PRMoneDefaults) + .setChoiceListItems(PRM_CHOICELIST_SINGLE, sResampleModeMenuItems) + .setTooltip( + "If the A and B VDBs have different transforms, one VDB should\n" + "be resampled to match the other before the two are combined.\n" + "Also, level set VDBs should have matching background values\n" + "(i.e., matching narrow band widths).")); + + // Menu of resampling interpolation order options + parms.add(hutil::ParmFactory(PRM_ORD, "resampleinterp", "Interpolation") + .setDefault(PRMoneDefaults) + .setChoiceListItems(PRM_CHOICELIST_SINGLE, { + "point", "Nearest", + "linear", "Linear", + "quadratic", "Quadratic" + }) + .setTooltip( + "Specify the type of interpolation to be used when\n" + "resampling one VDB to match the other's transform.") + .setDocumentation( + "The type of interpolation to be used when resampling one VDB" + " to match the other's transform\n\n" + "Nearest neighbor interpolation is fast but can introduce noticeable" + " sampling artifacts. Quadratic interpolation is slow but high-quality." + " Linear interpolation is intermediate in speed and quality.")); + + + // Register this operator. + // (See houdini_utils/Utils.h for OpFactory details.) + hvdb::OpenVDBOpFactory("VDB Fillet", SOP_OpenVDB_Fillet::factory, parms, *table) + .addInput("A VDBs") + .addOptionalInput("B VDBs") + .addOptionalInput("Mask VDB") + .setVerb(SOP_NodeVerb::COOK_INPLACE, []() { return new SOP_OpenVDB_Fillet::Cache; }); +} + + +//////////////////////////////////////// + + +OP_Node* +SOP_OpenVDB_Fillet::factory(OP_Network* net, + const char* name, OP_Operator* op) +{ + return new SOP_OpenVDB_Fillet(net, name, op); +} + + +SOP_OpenVDB_Fillet::SOP_OpenVDB_Fillet(OP_Network* net, + const char* name, OP_Operator* op): + hvdb::SOP_NodeVDB(net, name, op) +{ +} + + +//////////////////////////////////////// + + +// Enable/disable or show/hide parameters in the UI. +unsigned +SOP_OpenVDB_Fillet::disableParms() +{ + unsigned changed = 0; + + // Disable parms. + //changed += enableParm("dummy", + // evalInt("dummy", 0, /*time=*/0)); + + + //setVisibleState("dummy", getEnableState("dummy")); + + return changed; +} + + +//////////////////////////////////////// + +hvdb::GridPtr +SOP_OpenVDB_Fillet::Cache::blendLevelSets( + hvdb::GridCPtr aGrid, + hvdb::GridCPtr bGrid, + VDBFilletParms parms) +{ + openvdb::FloatGrid::ConstPtr a = openvdb::gridConstPtrCast(aGrid); + openvdb::FloatGrid::ConstPtr b = openvdb::gridConstPtrCast(bGrid); + + // sanitizers + if (!a) throw std::runtime_error("Missing the first grid"); + if (!b) throw std::runtime_error("Missing the second grid"); + if (a->getGridClass() != openvdb::GRID_LEVEL_SET) + throw std::runtime_error("First grid needs to be a level-set."); + if (b->getGridClass() != openvdb::GRID_LEVEL_SET) + throw std::runtime_error("Second grid needs to be a level-set."); + + ResampleOp rsmpl(a, b, parms); + rsmpl.resampleGrids(); + if (parms.mMaskPtr) rsmpl.resampleMask(parms.mMaskPtr); + + hvdb::GridPtr ret = openvdb::tools::unionFillet( + *a, *b, parms.mMaskPtr, parms.mAlpha, parms.mBeta, parms.mGamma); + + return ret; +} + +OP_ERROR +SOP_OpenVDB_Fillet::Cache::cookVDBSop(OP_Context& context) +{ + try { + UT_AutoInterrupt progress{"Processing VDB grids"}; + + const fpreal time = context.getTime(); + + // Get gdps + GU_Detail* aGdp = gdp; + const GU_Detail* bGdp = inputGeo(1, context); + + // Get the group of grids to process. + const auto aGroupStr = evalStdString("agroup", time); + const auto bGroupStr = evalStdString("bgroup", time); + const auto* bGroup = (!bGdp ? nullptr : matchGroup(*bGdp, bGroupStr)); + + // Fill-in VDBFilletParms + VDBFilletParms parms; + if (evalParms(context, parms) >= UT_ERROR_ABORT) return error(); + + // The collation pattern that we follow is 'pairs' as in + // SOP_OpenVDB_Combine. This means to combine pairs of A and B VDBs + // in the order in which they appear in their respective groups. + std::vector aGroupVec; + for (const auto& pattern: splitPatterns(aGroupStr)) { + aGroupVec.push_back(matchGroup(*aGdp, pattern)); + } + + // Iterate over one or more A groups. + for (const auto* aGroup: aGroupVec) { + hvdb::VdbPrimIterator aIt{aGdp, GA_Range::safedeletions{}, aGroup}; + hvdb::VdbPrimCIterator bIt{bGdp, bGroup}; + + // Populate two vectors of primitives, one comprising the A grids + // and the other the B grids. (In the case of flattening operations, + // these grids might be taken from the same input.) + // Note: the following relies on exhausted iterators returning nullptr + // and on incrementing an exhausted iterator being a no-op. + std::vector aVdbVec; + std::vector bVdbVec; + for ( ; aIt && bIt; ++aIt, ++bIt) { + aVdbVec.push_back(*aIt); + bVdbVec.push_back(*bIt); + } + + std::set vdbsToRemove; + // Iterate over A and, optionally, B grids. + for (size_t i = 0, N = std::min(aVdbVec.size(), bVdbVec.size()); i < N; ++i) { + if (progress.wasInterrupted()) { throw std::runtime_error{"interrupted"}; } + + // Note: even if needA is false, we still need to delete A grids. + GU_PrimVDB* aVdb = aVdbVec[i]; + const GU_PrimVDB* bVdb = bVdbVec[i]; + + hvdb::GridPtr aGrid; + hvdb::GridCPtr bGrid; + if (aVdb) aGrid = aVdb->getGridPtr(); + if (bVdb) bGrid = bVdb->getConstGridPtr(); + + parms.mABaseGrid = aGrid; + parms.mBBaseGrid = bGrid; + + if (hvdb::GridPtr outGrid = blendLevelSets(aGrid, bGrid, parms)) + { + // Name the output grid after the A grid if the A grid is used, + // or after the B grid otherwise. + UT_String outGridName = getGridName(aVdb); + // Add a new VDB primitive for the output grid to the output gdp. + GU_PrimVDB::buildFromGrid(*gdp, outGrid, /*copyAttrsFrom=*/ aVdb, outGridName); + vdbsToRemove.insert(aVdb); + } + } // iterate over aVdbVec and bVdbVec + + // Remove primitives that were copied from input 0. + for (GU_PrimVDB* vdb: vdbsToRemove) { + if (vdb) gdp->destroyPrimitive(*vdb, /*andPoints=*/true); + } + } // aGroup : aGroupVec + } catch (std::exception& e) { + addError(SOP_MESSAGE, e.what()); + } + return error(); +}