Skip to content

Commit

Permalink
Fixed an issue with smooth position rasterization which wouldn't alwa…
Browse files Browse the repository at this point in the history 
…ys activate the correct topology if the search radius was larger than the point radius

Signed-off-by: Nick Avramoussis <[email protected]>
  • Loading branch information
@Idclip
Idclip committed Sep 26, 2023
1 parent e04f8d8 commit b8406e9
Show file tree
Hide file tree
Showing 4 changed files with 207 additions and 34 deletions.
4 changes: 4 additions & 0 deletions openvdb/openvdb/points/PointRasterizeSDF.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
};

Expand Down
187 changes: 156 additions & 31 deletions openvdb/openvdb/points/impl/PointRasterizeSDFImpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -643,18 +643,23 @@ 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)
: mMask(new MaskTreeT)
, 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
Expand Down Expand Up @@ -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()
{
Expand Down Expand Up @@ -752,6 +765,7 @@ struct SurfaceMaskOp
std::unique_ptr<MaskTreeT> mMaskOff;
const math::Transform& mPointsTransform;
const math::Transform& mSurfaceTransform;
const CoordBBox mMaxClipBounds;
InterrupterT* mInterrupter;
};

Expand All @@ -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
Expand Down Expand Up @@ -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;
Expand All @@ -848,7 +865,9 @@ struct VariableSurfaceMaskOp
const int32_t max = this->maxDist(maxacc.getValue(leafIter->origin()));
this->activate(*leafIter, max);
}
const tree::ValueAccessor<const RadiusTreeT> minacc(mMin);

if (!mMin) return;
const tree::ValueAccessor<const RadiusTreeT> minacc(*mMin);
for (auto leafIter = range.begin(); leafIter; ++leafIter) {
const int32_t min = this->minDist(minacc.getValue(leafIter->origin()));
if (min < 0) continue;
Expand All @@ -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<int32_t>(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)
Expand All @@ -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 <typename SdfT, typename MaskTreeT>
Expand Down Expand Up @@ -945,7 +964,7 @@ initFixedSdf(const PointDataGridT& points,
if (interrupter) interrupter->start("Generating uniform surface topology");

FixedSurfaceMaskOp<MaskTreeT, InterrupterT> op(points.transform(),
*transform, minBandRadius, maxBandRadius, interrupter);
*transform, minBandRadius, maxBandRadius, /*clipbounds=*/nullptr, interrupter);

LeafManagerT manager(points.tree());
tbb::parallel_reduce(manager.leafRange(), op);
Expand Down Expand Up @@ -977,7 +996,72 @@ initVariableSdf(const PointDataGridT& points,
if (interrupter) interrupter->start("Generating variable surface topology");

VariableSurfaceMaskOp<RadiusTreeT, MaskTreeT, InterrupterT>
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<SdfT, MaskTreeT>(transform, bg, op.mask(), op.maskoff());

if (interrupter) interrupter->end();
return surface;
}

template <typename SdfT, typename InterrupterT, typename PointDataGridT>
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<const typename PointDataGridT::TreeType>;
using MaskTreeT = typename SdfT::TreeType::template ValueConverter<ValueMask>::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<MaskTreeT, InterrupterT> 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<SdfT, MaskTreeT>(transform, bg, op.mask(), op.maskoff());

if (interrupter) interrupter->end();
return surface;
}

template <typename SdfT,
typename InterrupterT,
typename PointDataGridT,
typename RadiusTreeT>
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<const typename PointDataGridT::TreeType>;
using MaskTreeT = typename SdfT::TreeType::template ValueConverter<ValueMask>::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<RadiusTreeT, MaskTreeT, InterrupterT>
op(points.transform(), *transform, nullptr, maxTree, /*minscale=*/1.0,
scale, halfband, &bounds, interrupter);

LeafManagerT manager(points.tree());
tbb::parallel_reduce(manager.leafRange(), op);
Expand Down Expand Up @@ -1272,54 +1356,95 @@ rasterizeSmoothSpheres(const PointDataGridT& points,
static_cast<typename SdfT::ValueType>(vs * halfband);

const Real indexSpaceSearch = settings.searchRadius / vs;
// max possible index space search radius
const size_t width = static_cast<size_t>(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<Real> 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<Real> bands(maxActivationRadius, halfband);
const Real max = bands.max();

surface = initFixedSdf<SdfT, InterrupterType>
(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<SdfT, InterrupterType>
(points, transform, background, max, bounds, interrupter);

if (!leaf) return GridPtrVec(1, surface);

// max possible index space search radius
const size_t width = static_cast<size_t>(math::RoundUp(indexSpaceSearch));

const FixedRadius<Real> rad(settings.radiusScale / vs);
if (interrupter) interrupter->start("Rasterizing particles to level set using constant Zhu-Bridson");

grids = doRasterizeSurface<SdfT, AveragePositionTransfer, AttributeTypes, InterrupterType>
(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<RadiusT>::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<RadiusT, UnknownCodec>
(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<SdfT, InterrupterType>
(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<SdfT, InterrupterType>
(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<size_t>(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");
Expand Down
48 changes: 45 additions & 3 deletions openvdb/openvdb/unittest/TestPointRasterizeSDF.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<points::NullFilter> 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<points::NullFilter> s;
Expand All @@ -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;
Expand Down Expand Up @@ -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());
Expand Down
Loading

0 comments on commit b8406e9

Please sign in to comment.