From b8406e9ce94e46dc0d3280949649f73c18602215 Mon Sep 17 00:00:00 2001 From: Nick Avramoussis <4256455+Idclip@users.noreply.github.com> Date: Mon, 28 Aug 2023 15:27:51 +1200 Subject: [PATCH] Fixed an issue with smooth position rasterization which wouldn't always activate the correct topology if the search radius was larger than the point radius Signed-off-by: Nick Avramoussis <4256455+Idclip@users.noreply.github.com> --- openvdb/openvdb/points/PointRasterizeSDF.h | 4 + .../points/impl/PointRasterizeSDFImpl.h | 187 +++++++++++++++--- .../openvdb/unittest/TestPointRasterizeSDF.cc | 48 ++++- pendingchanges/surfacing.txt | 2 + 4 files changed, 207 insertions(+), 34 deletions(-) diff --git a/openvdb/openvdb/points/PointRasterizeSDF.h b/openvdb/openvdb/points/PointRasterizeSDF.h index 6a5130a17d..8883ae385c 100644 --- a/openvdb/openvdb/points/PointRasterizeSDF.h +++ b/openvdb/openvdb/points/PointRasterizeSDF.h @@ -257,6 +257,10 @@ struct SmoothSphereSettings /// larger than the point radius. Both this and the `radiusScale` /// parameters are given in world space units and are applied to every /// point to generate a surface mask. + /// @warning If this value is less than the sum of the maximum particle + /// radius and the half band width, the exterior half band width may be + /// smaller than desired. In these cases, consider running a levelset + /// renormalize or a levelset rebuild. Real searchRadius = 1.0; }; diff --git a/openvdb/openvdb/points/impl/PointRasterizeSDFImpl.h b/openvdb/openvdb/points/impl/PointRasterizeSDFImpl.h index a5b819ead2..b9dc26a742 100644 --- a/openvdb/openvdb/points/impl/PointRasterizeSDFImpl.h +++ b/openvdb/openvdb/points/impl/PointRasterizeSDFImpl.h @@ -643,11 +643,15 @@ struct SurfaceMaskOp protected: SurfaceMaskOp(const math::Transform& points, const math::Transform& surface, + // Clip the surface to this bounds. Only used for the smooth + // raster workflow to limit to search radii topology init + const CoordBBox* const maxClipBounds = nullptr, InterrupterT* interrupter = nullptr) : mMask(new MaskTreeT) , mMaskOff(new MaskTreeT) , mPointsTransform(points) , mSurfaceTransform(surface) + , mMaxClipBounds(maxClipBounds ? this->toSurfaceBounds(*maxClipBounds) : CoordBBox::inf()) , mInterrupter(interrupter) {} SurfaceMaskOp(const SurfaceMaskOp& other) @@ -655,6 +659,7 @@ struct SurfaceMaskOp , mMaskOff(new MaskTreeT) , mPointsTransform(other.mPointsTransform) , mSurfaceTransform(other.mSurfaceTransform) + , mMaxClipBounds(other.mMaxClipBounds) , mInterrupter(other.mInterrupter) {} // @brief Sparse fill a tree with activated bounding boxes expanded from @@ -696,8 +701,16 @@ struct SurfaceMaskOp this->deactivate(bounds); } - inline void activate(const CoordBBox& bounds) { mMask->sparseFill(bounds, true, true); } - inline void deactivate(const CoordBBox& bounds) { mMaskOff->sparseFill(bounds, true, true); } + inline void activate(CoordBBox& bounds) + { + bounds.intersect(mMaxClipBounds); + mMask->sparseFill(bounds, true, true); + } + + inline void deactivate(const CoordBBox& bounds) + { + mMaskOff->sparseFill(bounds, true, true); + } inline bool interrupted() { @@ -752,6 +765,7 @@ struct SurfaceMaskOp std::unique_ptr mMaskOff; const math::Transform& mPointsTransform; const math::Transform& mSurfaceTransform; + const CoordBBox mMaxClipBounds; InterrupterT* mInterrupter; }; @@ -769,8 +783,9 @@ struct FixedSurfaceMaskOp const math::Transform& surface, const double minBandRadius, // sdf index space const double maxBandRadius, // sdf index space + const CoordBBox* const maxClipBounds, InterrupterT* interrupter = nullptr) - : BaseT(points, surface, interrupter) + : BaseT(points, surface, maxClipBounds, interrupter) , mMin(), mMax() { // calculate the min interior cube area of activity. this is the side @@ -827,13 +842,15 @@ struct VariableSurfaceMaskOp VariableSurfaceMaskOp(const math::Transform& pointsTransform, const math::Transform& surfaceTransform, - const RadiusTreeT& min, + const RadiusTreeT* const min, const RadiusTreeT& max, - const Real scale, + const Real minScale, + const Real maxScale, const Real halfband, + const CoordBBox* const maxClipBounds, InterrupterT* interrupter = nullptr) - : BaseT(pointsTransform, surfaceTransform, interrupter) - , mMin(min), mMax(max), mScale(scale) + : BaseT(pointsTransform, surfaceTransform, maxClipBounds, interrupter) + , mMin(min), mMax(max), mMinScale(minScale), mMaxScale(maxScale) , mHalfband(halfband) {} VariableSurfaceMaskOp(const VariableSurfaceMaskOp&) = default; @@ -848,7 +865,9 @@ struct VariableSurfaceMaskOp const int32_t max = this->maxDist(maxacc.getValue(leafIter->origin())); this->activate(*leafIter, max); } - const tree::ValueAccessor minacc(mMin); + + if (!mMin) return; + const tree::ValueAccessor minacc(*mMin); for (auto leafIter = range.begin(); leafIter; ++leafIter) { const int32_t min = this->minDist(minacc.getValue(leafIter->origin())); if (min < 0) continue; @@ -860,14 +879,14 @@ struct VariableSurfaceMaskOp inline int32_t maxDist(const typename RadiusTreeT::ValueType& maxRadiusWs) const { // max radius in index space - const Real maxBandRadius = (Real(maxRadiusWs) * mScale) + mHalfband; + const Real maxBandRadius = (Real(maxRadiusWs) * mMaxScale) + mHalfband; return static_cast(math::Round(maxBandRadius)); // furthest voxel } inline int32_t minDist(const typename RadiusTreeT::ValueType& minRadiusWs) const { // min radius in index space - Real minBandRadius = math::Max(0.0, (Real(minRadiusWs) * mScale) - mHalfband); + Real minBandRadius = math::Max(0.0, (Real(minRadiusWs) * mMinScale) - mHalfband); // calculate the min interior cube area of activity. this is the side // of the largest possible cube that fits into the radius "min": // d = 2r -> 2r = 3x^2 -> x = 2r / sqrt(3) @@ -889,9 +908,9 @@ struct VariableSurfaceMaskOp } private: - const RadiusTreeT& mMin; + const RadiusTreeT* const mMin; const RadiusTreeT& mMax; - const Real mScale, mHalfband; + const Real mMinScale, mMaxScale, mHalfband; }; template @@ -945,7 +964,7 @@ initFixedSdf(const PointDataGridT& points, if (interrupter) interrupter->start("Generating uniform surface topology"); FixedSurfaceMaskOp op(points.transform(), - *transform, minBandRadius, maxBandRadius, interrupter); + *transform, minBandRadius, maxBandRadius, /*clipbounds=*/nullptr, interrupter); LeafManagerT manager(points.tree()); tbb::parallel_reduce(manager.leafRange(), op); @@ -977,7 +996,72 @@ initVariableSdf(const PointDataGridT& points, if (interrupter) interrupter->start("Generating variable surface topology"); VariableSurfaceMaskOp - op(points.transform(), *transform, min, max, scale, halfband, interrupter); + op(points.transform(), *transform, &min, max, scale, scale, halfband, + /*clipbounds=*/nullptr, interrupter); + + LeafManagerT manager(points.tree()); + tbb::parallel_reduce(manager.leafRange(), op); + + typename SdfT::Ptr surface = + initSdfFromMasks(transform, bg, op.mask(), op.maskoff()); + + if (interrupter) interrupter->end(); + return surface; +} + +template +inline typename SdfT::Ptr +initFixedSmoothSdf(const PointDataGridT& points, + math::Transform::Ptr transform, + const typename SdfT::ValueType bg, + const Real maxBandRadius, + const CoordBBox& bounds, + InterrupterT* interrupter) +{ + using LeafManagerT = tree::LeafManager; + using MaskTreeT = typename SdfT::TreeType::template ValueConverter::Type; + + if (interrupter) interrupter->start("Generating uniform surface topology"); + + // @note Currently don't use min radii to deactivate, can't compute + // this with the ZB kernel + FixedSurfaceMaskOp op(points.transform(), + *transform, /*minBandRadius=*/0.0, maxBandRadius, &bounds, interrupter); + + LeafManagerT manager(points.tree()); + tbb::parallel_reduce(manager.leafRange(), op); + + typename SdfT::Ptr surface = + initSdfFromMasks(transform, bg, op.mask(), op.maskoff()); + + if (interrupter) interrupter->end(); + return surface; +} + +template +inline typename SdfT::Ptr +initVariableSmoothSdf(const PointDataGridT& points, + math::Transform::Ptr transform, + const typename SdfT::ValueType bg, + const RadiusTreeT& maxTree, + const Real scale, + const Real halfband, + const CoordBBox& bounds, + InterrupterT* interrupter) +{ + using LeafManagerT = tree::LeafManager; + using MaskTreeT = typename SdfT::TreeType::template ValueConverter::Type; + + if (interrupter) interrupter->start("Generating variable surface topology"); + + // @note Currently don't use min radii/tree to deactivate, can't + // compute this with the ZB kernel + VariableSurfaceMaskOp + op(points.transform(), *transform, nullptr, maxTree, /*minscale=*/1.0, + scale, halfband, &bounds, interrupter); LeafManagerT manager(points.tree()); tbb::parallel_reduce(manager.leafRange(), op); @@ -1272,24 +1356,44 @@ rasterizeSmoothSpheres(const PointDataGridT& points, static_cast(vs * halfband); const Real indexSpaceSearch = settings.searchRadius / vs; + // max possible index space search radius + const size_t width = static_cast(math::RoundUp(indexSpaceSearch)); + + // The topology we need to activate is at a distance based on the maximum + // radii of each point and the uniform search radius. Even though we're + // guaranteed to be generating new positions within the distribution of + // point neighbours, these positions may end up outside of active topology + // were we _only_ to use the radius of the particles for topology + // activation. + const Real maxActivationRadius = std::max(settings.searchRadius, settings.radiusScale) / vs; const auto leaf = points.constTree().cbeginLeaf(); + // Compute estimated max bounds for clipping. This is used if the search + // radius is larger than the max particle radius (as we don't need to + // activate topology further outside the bounds of the point data grid). + // This bounds is expanded by the halfband + max radii + CoordBBox bounds; + points.tree().evalLeafBoundingBox(bounds); + typename SdfT::Ptr surface; GridPtrVec grids; if (settings.radius.empty()) { - const FixedBandRadius bands(settings.radiusScale / vs, halfband); + // This is the max possible distance we need to activate, but we'll + // clip this at the edges of the point bounds (as the ZB kernel will + // only create positions in between points). + const FixedBandRadius bands(maxActivationRadius, halfband); const Real max = bands.max(); - surface = initFixedSdf - (points, transform, background, /*min*/0.0, max, interrupter); + // Compute max radius in index space and expand bounds + bounds.expand(halfband + math::Round(settings.radiusScale / vs)); + // init surface + surface = initFixedSmoothSdf + (points, transform, background, max, bounds, interrupter); if (!leaf) return GridPtrVec(1, surface); - // max possible index space search radius - const size_t width = static_cast(math::RoundUp(indexSpaceSearch)); - const FixedRadius rad(settings.radiusScale / vs); if (interrupter) interrupter->start("Rasterizing particles to level set using constant Zhu-Bridson"); @@ -1297,29 +1401,50 @@ rasterizeSmoothSpheres(const PointDataGridT& points, (points, attributes, filter, *surface, interrupter, width, rad, indexSpaceSearch, points.transform(), *surface); // args } - else { + else + { using RadiusT = typename SettingsT::RadiusAttributeType; using PointDataTreeT = typename PointDataGridT::TreeType; using RadTreeT = typename PointDataTreeT::template ValueConverter::Type; + // We currently don't use the min values for the ZB kernel topology + // activation. + // @todo We should be able to deactivate on some metric RadiusT min(0), max(0); - typename RadTreeT::Ptr mintree(new RadTreeT), maxtree(new RadTreeT); + typename RadTreeT::Ptr maxtree(new RadTreeT); points::evalMinMax - (points.tree(), settings.radius, min, max, filter, mintree.get(), maxtree.get()); + (points.tree(), settings.radius, min, max, filter, nullptr, maxtree.get()); + + if ((settings.searchRadius > settings.radiusScale) && (min < settings.searchRadius)) { + // set radius tree values to search distances if they are less, + // just for the surface topology initialization. This is the max + // possible distance we need to activate, but we'll clip this at + // the edges of the point bounds (as the ZB kernel will only + // create positions in-between points). + tools::foreach(maxtree->beginValueOn(), + [&](auto& iter) { + iter.modifyValue([&](auto& r) { + if (r < settings.searchRadius) { + // initVariableSmoothSdf scales radii by (settings.radiusScale/vs). + // We don't want to scale the search radii by the radius scale, so + // cancel it out here. + r = (settings.searchRadius / settings.radiusScale); + } + }); + }, /*thread=*/true, /*shared=*/true); + } - // search distance at the SDF transform const RadiusT indexSpaceScale = RadiusT(settings.radiusScale / vs); - surface = initVariableSdf - (points, transform, background, *mintree, *maxtree, - indexSpaceScale, halfband, interrupter); - mintree.reset(); + // Compute max radius in index space and expand bounds + bounds.expand(halfband + math::Round(max * indexSpaceScale)); + // init surface + surface = initVariableSmoothSdf + (points, transform, background, *maxtree, + indexSpaceScale, halfband, bounds, interrupter); maxtree.reset(); if (!leaf) return GridPtrVec(1, surface); - // max possible index space search radius - const size_t width = static_cast(math::RoundUp(indexSpaceSearch)); - const size_t ridx = leaf->attributeSet().find(settings.radius); if (ridx == AttributeSet::INVALID_POS) { OPENVDB_THROW(RuntimeError, "Failed to find radius attribute"); diff --git a/openvdb/openvdb/unittest/TestPointRasterizeSDF.cc b/openvdb/openvdb/unittest/TestPointRasterizeSDF.cc index f63a95ce83..acb5250f6a 100644 --- a/openvdb/openvdb/unittest/TestPointRasterizeSDF.cc +++ b/openvdb/openvdb/unittest/TestPointRasterizeSDF.cc @@ -704,6 +704,48 @@ TEST_F(TestPointRasterizeSDF, testRasterizeSmoothSpheres) } } + // Test two points qhich create a ghost particle outside of their + // radii (to test that the surface topology correctly accounts for the + // search distance) + { + FixedSurfacing s; + s.halfband = 2; + s.points = PointBuilder({ Vec3f(0), Vec3f(0,5,0) }).voxelsize(0.1).get(); + s.transform = nullptr; + double radius = 0.5, search = 5; // large search + + FloatGrid::Ptr sdf = s.surfaceSmooth(radius, search); + EXPECT_TRUE(sdf && sdf->getName() == "fixed_avg"); + EXPECT_TRUE(sdf->transform() == s.points->transform()); + EXPECT_EQ(GRID_LEVEL_SET, sdf->getGridClass()); + EXPECT_EQ(float(s.halfband * s.points->voxelSize()[0]), sdf->background()); + + // @todo regression test. find a way to better test these values + size_t interiorOn = 0, exteriorOn = 0; + EXPECT_EQ(Index32(36), sdf->tree().leafCount()); + EXPECT_EQ(Index64(0), sdf->tree().activeTileCount()); + EXPECT_EQ(Index64(4188), sdf->tree().activeVoxelCount()); + + for (auto iter = sdf->cbeginValueOn(); iter; ++iter) { + EXPECT_TRUE(*iter > -sdf->background()); + EXPECT_TRUE(*iter < sdf->background()); + if (*iter > 0) ++exteriorOn; + else ++interiorOn; + } + EXPECT_EQ(size_t(1244), interiorOn); + EXPECT_EQ(size_t(2944), exteriorOn); + + size_t interiorOff = 0, exteriorOff = 0; + for (auto iter = sdf->cbeginValueOff(); iter; ++iter) { + EXPECT_TRUE((*iter <= -sdf->background()) || (*iter >= sdf->background())) + << *iter << " " << sdf->background();; + if (*iter > 0) ++exteriorOff; + else ++interiorOff; + } + EXPECT_EQ(size_t(323), interiorOff); + EXPECT_EQ(size_t(308789), exteriorOff); + } + // Test multiple points - 8 points at cube corner positions { FixedSurfacing s; @@ -726,7 +768,7 @@ TEST_F(TestPointRasterizeSDF, testRasterizeSmoothSpheres) // test points from a box with coords at 0.5 and a search radius // large enough to create a smoothed box - // @todo find a way to better test these values + // @todo regression test. find a way to better test these values positions = getBoxPoints(/*scale*/0.5f); radius = 0.2, search = 1.2; s.halfband = 3; @@ -881,9 +923,9 @@ TEST_F(TestPointRasterizeSDF, testRasterizeVariableSmoothSpheres) EXPECT_TRUE(sdf->transform() == s.points->transform()); EXPECT_EQ(GRID_LEVEL_SET, sdf->getGridClass()); EXPECT_EQ(float(s.halfband * s.points->voxelSize()[0]), sdf->background()); - EXPECT_EQ(Index32(64), sdf->tree().leafCount()); + EXPECT_EQ(Index32(60), sdf->tree().leafCount()); EXPECT_EQ(Index64(0), sdf->tree().activeTileCount()); - EXPECT_EQ(Index64(15011), sdf->tree().activeVoxelCount()); + EXPECT_EQ(Index64(14215), sdf->tree().activeVoxelCount()); for (auto iter = sdf->cbeginValueOn(); iter; ++iter) { EXPECT_TRUE(*iter > -sdf->background()); EXPECT_TRUE(*iter < sdf->background()); diff --git a/pendingchanges/surfacing.txt b/pendingchanges/surfacing.txt index 3f97747675..3ce1c30d20 100644 --- a/pendingchanges/surfacing.txt +++ b/pendingchanges/surfacing.txt @@ -22,4 +22,6 @@ OpenVDB: - Fixed a bug in points::rasterizeSpheres and points::rasterizeSmoothSpheres which could cause sections of the generated surface to incorrectly be marked as negative interior values. + - Fixed a bug in points::rasterizeSmoothSpheres which could result in + incomplete surface reconstruction with significantly larger search distances. - Fixed some precision issues in various Matrix methods