Skip to content

Commit

Permalink
Moved the old 'sphereScale' rasterizer param to be a uniform scale fo…
Browse files Browse the repository at this point in the history 
…r non-anisotropic neighbourhoods on the PCA methods

Signed-off-by: Nick Avramoussis <[email protected]>
  • Loading branch information
@Idclip
Idclip committed Oct 25, 2023
1 parent de71777 commit faf1b1c
Show file tree
Hide file tree
Showing 6 changed files with 92 additions and 177 deletions.
6 changes: 0 additions & 6 deletions openvdb/openvdb/points/PointRasterizeSDF.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -289,12 +289,6 @@ struct EllipsoidSettings
using BaseT::filter;
using BaseT::interrupter;

/// @brief The sphere scale. Points which are not in the inclusion group
/// specified by the pca attributes have their world space radius scaled
/// by this amount. Typically you'd want this value to be <= 1.0 to
/// produce smaller spheres for isolated points.
Real sphereScale = 1.0;

/// @brief The required principal component analysis attributes which are
/// required to exist on the points being rasterized. These attributes
/// define the rotational and affine transformations which can be used to
Expand Down
6 changes: 5 additions & 1 deletion openvdb/openvdb/points/PrincipalComponentAnalysis.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,6 @@
#include <vector>
#include <memory>
#include <cmath> // std::cbrt
#include <algorithm> // std::accumulate

namespace openvdb {
OPENVDB_USE_VERSION_NAMESPACE
Expand Down Expand Up @@ -88,6 +87,11 @@ struct PcaSettings
/// so a reasonable minimum should be set. Values equal to or less than
/// 0.0, or values greater than 1.0 have undefined results.
float allowedAnisotropyRatio = 0.25f;
/// @param nonAnisotropicStretch The stretch coefficient that should be
/// used for points which have no anisotropic neighbourhood (due to
/// being isolated or not having enough neighbours to reach the
/// specified @sa neighbourThreshold).
float nonAnisotropicStretch = 1.0;
/// @param neighbourThreshold the number of neighbours a point must have
/// to be classified as having an elliptical distribution. Points with
/// less neighbours than this will end up with uniform stretch values of
Expand Down
148 changes: 49 additions & 99 deletions openvdb/openvdb/points/impl/PointRasterizeEllipsoidsSDFImpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,11 +48,9 @@ struct EllipsIndicies
const PcaAttributes& attribs)
: stretch(EllipsIndicies::getAttributeIndex<Vec3f>(desc, attribs.stretch))
, rotation(EllipsIndicies::getAttributeIndex<Mat3s>(desc, attribs.rotation))
, position(EllipsIndicies::getAttributeIndex<Vec3d>(desc, attribs.positionWS))
, ellipsesGroup(desc.groupIndex(attribs.ellipses)) {}
, position(EllipsIndicies::getAttributeIndex<Vec3d>(desc, attribs.positionWS)) {}

const size_t stretch, rotation, position;
const points::AttributeSet::Descriptor::GroupIndex ellipsesGroup;

private:
template<typename ValueT>
Expand Down Expand Up @@ -97,34 +95,28 @@ struct EllipsoidTransfer :
const RadiusType& rt,
const math::Transform& source,
SdfT& surface,
const Real sphereScale,
const EllipsIndicies& indices,
Int64Tree* cpg = nullptr,
const std::unordered_map<const PointDataTree::LeafNodeType*, Index>* ids = nullptr)
: BaseT(pidx, width, rt, source, surface, cpg, ids)
, mSphereScale(sphereScale)
, mIndices(indices)
, mStretchHandle()
, mRotationHandle()
, mPositionWSHandle()
, mEllipsesGroupHandle() {}
, mPositionWSHandle() {}

EllipsoidTransfer(const EllipsoidTransfer& other)
: BaseT(other)
, mSphereScale(other.mSphereScale)
, mIndices(other.mIndices)
, mStretchHandle()
, mRotationHandle()
, mPositionWSHandle()
, mEllipsesGroupHandle() {}
, mPositionWSHandle() {}

inline bool startPointLeaf(const PointDataTree::LeafNodeType& leaf)
{
const bool ret = this->BaseT::startPointLeaf(leaf);
mStretchHandle.reset(new StretchHandleT(leaf.constAttributeArray(mIndices.stretch)));
mRotationHandle.reset(new RotationHandleT(leaf.constAttributeArray(mIndices.rotation)));
mPositionWSHandle.reset(new PwsHandleT(leaf.constAttributeArray(mIndices.position)));
mEllipsesGroupHandle.reset(new points::GroupHandle(leaf.groupHandle(mIndices.ellipsesGroup)));
return ret;
}

Expand All @@ -135,18 +127,18 @@ struct EllipsoidTransfer :
// Position may have been smoothed so we need to use the ws handle
const Vec3d PWS = this->mPositionWSHandle->get(id);
const Vec3d P = this->targetTransform().worldToIndex(PWS);
// group membership denotes non-isolated points
const bool isolated = !mEllipsesGroupHandle->get(id);
const Vec3f stretch = mStretchHandle->get(id);

if (isolated) {
// If we have a uniform stretch, treat as a sphere
// @todo worth using a tolerance here?
if ((stretch.x() == stretch.y()) && (stretch.x() == stretch.z())) {
this->BaseT::rasterizePoint(P, id, bounds,
this->mRadius.eval(id, typename RadiusType::ValueType(mSphereScale)));
this->mRadius.eval(id, typename RadiusType::ValueType(stretch.x())));
return;
}

const auto& r = this->mRadius.eval(id);
const RealT radius = r.get(); // index space radius
const Vec3f stretch = mStretchHandle->get(id);
const math::Mat3s rotation = mRotationHandle->get(id);
const math::Mat3s ellipsoidTransform = rotation.timesDiagonal(stretch);
const Vec3d max = calcEllipsoidBoundMax(ellipsoidTransform, r.max());
Expand Down Expand Up @@ -222,12 +214,10 @@ struct EllipsoidTransfer :
}

private:
const RealT mSphereScale;
const EllipsIndicies& mIndices;
std::unique_ptr<StretchHandleT> mStretchHandle;
std::unique_ptr<RotationHandleT> mRotationHandle;
std::unique_ptr<PwsHandleT> mPositionWSHandle;
std::unique_ptr<points::GroupHandle> mEllipsesGroupHandle;
};


Expand All @@ -245,20 +235,17 @@ struct EllipsSurfaceMaskOp
const math::Transform& src,
const math::Transform& trg,
const RadiusType& rad,
const Real sphereScale,
const Real halfband,
const EllipsIndicies& indices)
: BaseT(src, trg, nullptr)
, mRadius(rad)
, mSphereScale(sphereScale)
, mHalfband(halfband)
, mIndices(indices)
, mMaxDist(0) {}

EllipsSurfaceMaskOp(const EllipsSurfaceMaskOp& other, tbb::split)
: BaseT(other)
, mRadius(other.mRadius)
, mSphereScale(other.mSphereScale)
, mHalfband(other.mHalfband)
, mIndices(other.mIndices)
, mMaxDist(0) {}
Expand All @@ -280,10 +267,8 @@ struct EllipsSurfaceMaskOp
// @todo detect this from PcaSettings and avoid checking if no smoothing
points::AttributeHandle<Vec3d> pwsHandle(leaf.constAttributeArray(mIndices.position));
points::AttributeHandle<Vec3f> stretchHandle(leaf.constAttributeArray(mIndices.stretch));
points::GroupHandle ellipses(leaf.groupHandle(mIndices.ellipsesGroup));
if (stretchHandle.size() == 0) return;

RadiusT maxr = 0;
// The max stretch coefficient. We can't analyze each xyz component
// individually as we don't take into account the ellips rotation, so
// have to expand the worst case uniformly
Expand All @@ -293,57 +278,22 @@ struct EllipsSurfaceMaskOp

RadiusType rad(mRadius);
rad.reset(leaf);
bool hasIsolatedPoints = false;

for (Index i = 0; i < Index(stretchHandle.size()); ++i)
{
const Vec3d Pws = pwsHandle.get(i);
maxPos = math::maxComponent(maxPos, Pws);
minPos = math::minComponent(minPos, Pws);
const auto stretch = stretchHandle.get(i);

if constexpr(!RadiusType::Fixed) {
// This is doing an unnecessary multi by the scale for every
// point as we could defer the radius scale multi till after
// all min/max operations
const auto r = rad.eval(i).get(); // index space radius
maxr = std::max(maxr, r);

if (!ellipses.get(i)) {
hasIsolatedPoints = true;
}
else {
const auto stretch = stretchHandle.get(i);
const float maxcoef = std::max(stretch[0], std::max(stretch[1], stretch[2]));
// scaling by r puts the stretch values in index space
maxs = std::max(maxs, maxcoef * r);
}
}
else {
if (!ellipses.get(i)) {
hasIsolatedPoints = true;
}
else {
const auto stretch = stretchHandle.get(i);
const float maxcoef = std::max(stretch[0], std::max(stretch[1], stretch[2]));
maxs = std::max(maxs, maxcoef);
}
}
float maxcoef = std::max(stretch[0], std::max(stretch[1], stretch[2]));
if constexpr(!RadiusType::Fixed) maxcoef *= rad.eval(i).get(); // index space radius
maxs = std::max(maxs, maxcoef);
}

if constexpr(RadiusType::Fixed) {
maxr = rad.get();
// scaling by r puts the stretch values in index space
maxs *= float(maxr);
}

// If there are isolated points that we will surface as spheres, use
// the max between the stretch values and scaled sphere radii
if (hasIsolatedPoints) {
// scale radius by sphere scale
maxr *= RadiusT(mSphereScale);
// choose max component scale - compare the ellips to isolated
// points and choose the largest radius of the two
maxs = std::max(maxs, float(maxr));
maxs *= float(rad.get());
}

// @note This addition of the halfband here doesn't take into account
Expand Down Expand Up @@ -375,75 +325,77 @@ struct EllipsSurfaceMaskOp
points::AttributeHandle<Vec3d> pwsHandle(leaf.constAttributeArray(mIndices.position));
points::AttributeHandle<Vec3f> stretchHandle(leaf.constAttributeArray(mIndices.stretch));
points::AttributeHandle<math::Mat3s> rotHandle(leaf.constAttributeArray(mIndices.rotation));
points::GroupHandle ellipses(leaf.groupHandle(mIndices.ellipsesGroup));
if (stretchHandle.size() == 0) return;

RadiusT maxr = 0;
RadiusT maxr = 0, maxs = 0;
Vec3d maxBounds(0),
maxPos(std::numeric_limits<Real>::lowest()),
minPos(std::numeric_limits<Real>::max());

RadiusType rad(mRadius);
rad.reset(leaf);
bool hasIsolatedPoints = false;

for (Index i = 0; i < stretchHandle.size(); ++i)
{
const Vec3d Pws = pwsHandle.get(i);
maxPos = math::maxComponent(maxPos, Pws);
minPos = math::minComponent(minPos, Pws);

if constexpr(!RadiusType::Fixed) {
// This is doing an unnecessary multi by the scale for every
// point as we could defer the radius scale multi till after
// all min/max operations
const auto r = rad.eval(i).get(); // index space radius
maxr = std::max(maxr, r);
const Vec3f stretch = stretchHandle.get(i);
const bool isEllips = (stretch.x() != stretch.y()) || (stretch.x() != stretch.z());

if (!ellipses.get(i)) {
hasIsolatedPoints = true;
if constexpr(RadiusType::Fixed)
{
if (!isEllips) {
maxs = std::max(maxs, RadiusT(stretch.x()));
}
else {
const Vec3f stretch = stretchHandle.get(i);
const math::Mat3s rotation = rotHandle.get(i);
const math::Mat3s ellipsoidTransform = rotation.timesDiagonal(stretch);
// scaling by r puts the bounds values in index space
const Vec3d bounds = calcEllipsoidBoundMax(ellipsoidTransform, r);
// For fixed radii, compared the squared distances - we sqrt at the end
const Vec3d bounds = calcUnitEllipsoidBoundMaxSq(ellipsoidTransform);
maxBounds = math::maxComponent(maxBounds, bounds);
}
}
else {
if (!ellipses.get(i)) {
hasIsolatedPoints = true;
else
{
// This is doing an unnecessary multi by the scale for every
// point as we could defer the radius scale multi till after
// all min/max operations
const auto r = rad.eval(i).get(); // index space radius
maxr = std::max(maxr, r);

if (!isEllips) {
// for variable radii, scale each stretch by r
maxs = std::max(maxs, r * RadiusT(stretch.x()));
}
else {
const Vec3f stretch = stretchHandle.get(i);
const math::Mat3s rotation = rotHandle.get(i);
const math::Mat3s ellipsoidTransform = rotation.timesDiagonal(stretch);
// For fixed radii, compared the squared distances - we sqrt at the end
const Vec3d bounds = calcUnitEllipsoidBoundMaxSq(ellipsoidTransform);
// scaling by r puts the bounds values in index space
const Vec3d bounds = calcEllipsoidBoundMax(ellipsoidTransform, r);
maxBounds = math::maxComponent(maxBounds, bounds);
}
}
}

if constexpr(RadiusType::Fixed) {
// We don't do the sqrt per point for fixed radii - resolve the
// actual maxBounds now
if constexpr(RadiusType::Fixed)
{
maxr = rad.get(); // index space radius
maxs *= maxr; // Also scale maxs by the radius
// scaling by r puts the bounds values in index space
maxBounds[0] = std::sqrt(maxBounds[0]) * double(maxr);
maxBounds[1] = std::sqrt(maxBounds[1]) * double(maxr);
maxBounds[2] = std::sqrt(maxBounds[2]) * double(maxr);
}

if (hasIsolatedPoints) {
// scale radius by sphere scale
maxr *= RadiusT(mSphereScale);
// choose max component scale - compare the ellips to isolated
// points and choose the largest radius of the two
maxBounds[0] = std::max(double(maxr), maxBounds[0]);
maxBounds[1] = std::max(double(maxr), maxBounds[1]);
maxBounds[2] = std::max(double(maxr), maxBounds[2]);
}
// Account for uniform stretch values - compare the ellips to isolated
// points and choose the largest radius of the two
maxBounds[0] = std::max(double(maxs), maxBounds[0]);
maxBounds[1] = std::max(double(maxs), maxBounds[1]);
maxBounds[2] = std::max(double(maxs), maxBounds[2]);

// @note This addition of the halfband here doesn't take into account
// the squash on the halfband itself. The subsequent rasterizer
Expand Down Expand Up @@ -473,7 +425,6 @@ struct EllipsSurfaceMaskOp

private:
const RadiusType& mRadius;
const Real mSphereScale;
const Real mHalfband;
const EllipsIndicies& mIndices;
Vec3i mMaxDist;
Expand Down Expand Up @@ -501,7 +452,6 @@ rasterizeEllipsoids(const PointDataGridT& points,
const std::vector<std::string>& attributes = settings.attributes;
const Real halfband = settings.halfband;
const Real radiusScale = settings.radiusScale;
const Real sphereScale = settings.sphereScale;
auto* interrupter = settings.interrupter;

math::Transform::Ptr transform = settings.transform;
Expand Down Expand Up @@ -547,7 +497,7 @@ rasterizeEllipsoids(const PointDataGridT& points,
// pass radius scale as index space

EllipsSurfaceMaskOp<FixedBandRadius<Real>, MaskTreeT, InterrupterType>
op(points.transform(), *transform, rad, sphereScale, halfband, indices);
op(points.transform(), *transform, rad, halfband, indices);
tbb::parallel_reduce(manager.leafRange(), op);

surface = rasterize_sdf_internal::initSdfFromMasks<SdfT, MaskTreeT>
Expand All @@ -562,7 +512,7 @@ rasterizeEllipsoids(const PointDataGridT& points,

grids = doRasterizeSurface<SdfT, EllipsoidTransfer, AttributeTypes, InterrupterType>
(points, attributes, filter, *surface, interrupter,
width, rad, points.transform(), *surface, sphereScale, indices); // args
width, rad, points.transform(), *surface, indices); // args
}
else
{
Expand Down Expand Up @@ -593,7 +543,7 @@ rasterizeEllipsoids(const PointDataGridT& points,

// pass radius scale as index space
EllipsSurfaceMaskOp<VaryingBandRadius<RadiusT>, MaskTreeT, InterrupterType>
op(points.transform(), *transform, rad, sphereScale, halfband, indices);
op(points.transform(), *transform, rad, halfband, indices);
tbb::parallel_reduce(manager.leafRange(), op);

surface = rasterize_sdf_internal::initSdfFromMasks<SdfT, MaskTreeT>
Expand All @@ -608,7 +558,7 @@ rasterizeEllipsoids(const PointDataGridT& points,

grids = doRasterizeSurface<SdfT, EllipsoidTransfer, AttributeTypes, InterrupterType>
(points, attributes, filter, *surface, interrupter,
width, rad, points.transform(), *surface, sphereScale, indices); // args
width, rad, points.transform(), *surface, indices); // args
}

if (interrupter) interrupter->end();
Expand Down
2 changes: 1 addition & 1 deletion openvdb/openvdb/points/impl/PrincipalComponentAnalysisImpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -586,7 +586,7 @@ pca(PointDataGridT& points,
// 1) Create persisting attributes
const size_t pwsIdx = initAttribute(attrs.positionWS, zeroVal<PcaAttributes::PosWsT>());
const size_t rotIdx = initAttribute(attrs.rotation, zeroVal<PcaAttributes::RotationT>());
const size_t strIdx = initAttribute(attrs.stretch, PcaAttributes::StretchT(1));
const size_t strIdx = initAttribute(attrs.stretch, PcaAttributes::StretchT(settings.nonAnisotropicStretch));

// 2) Create temporary attributes
const auto& descriptor = leaf->attributeSet().descriptor();
Expand Down
Loading

0 comments on commit faf1b1c

Please sign in to comment.