From faf1b1c56e618242c356817fb0435882d593ed10 Mon Sep 17 00:00:00 2001 From: Nick Avramoussis <4256455+Idclip@users.noreply.github.com> Date: Wed, 25 Oct 2023 13:51:25 +1300 Subject: [PATCH] Moved the old 'sphereScale' rasterizer param to be a uniform scale for non-anisotropic neighbourhoods on the PCA methods Signed-off-by: Nick Avramoussis <4256455+Idclip@users.noreply.github.com> --- openvdb/openvdb/points/PointRasterizeSDF.h | 6 - .../points/PrincipalComponentAnalysis.h | 6 +- .../impl/PointRasterizeEllipsoidsSDFImpl.h | 148 ++++++------------ .../impl/PrincipalComponentAnalysisImpl.h | 2 +- openvdb/openvdb/unittest/TestPCA.cc | 8 +- .../openvdb/unittest/TestPointRasterizeSDF.cc | 99 ++++-------- 6 files changed, 92 insertions(+), 177 deletions(-) diff --git a/openvdb/openvdb/points/PointRasterizeSDF.h b/openvdb/openvdb/points/PointRasterizeSDF.h index 8883ae385c..324a59af91 100644 --- a/openvdb/openvdb/points/PointRasterizeSDF.h +++ b/openvdb/openvdb/points/PointRasterizeSDF.h @@ -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 diff --git a/openvdb/openvdb/points/PrincipalComponentAnalysis.h b/openvdb/openvdb/points/PrincipalComponentAnalysis.h index 95c8c7afd1..9f7b4e2a9d 100644 --- a/openvdb/openvdb/points/PrincipalComponentAnalysis.h +++ b/openvdb/openvdb/points/PrincipalComponentAnalysis.h @@ -32,7 +32,6 @@ #include #include #include // std::cbrt -#include // std::accumulate namespace openvdb { OPENVDB_USE_VERSION_NAMESPACE @@ -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 diff --git a/openvdb/openvdb/points/impl/PointRasterizeEllipsoidsSDFImpl.h b/openvdb/openvdb/points/impl/PointRasterizeEllipsoidsSDFImpl.h index 4fed1ac67c..e97be627ab 100644 --- a/openvdb/openvdb/points/impl/PointRasterizeEllipsoidsSDFImpl.h +++ b/openvdb/openvdb/points/impl/PointRasterizeEllipsoidsSDFImpl.h @@ -48,11 +48,9 @@ struct EllipsIndicies const PcaAttributes& attribs) : stretch(EllipsIndicies::getAttributeIndex(desc, attribs.stretch)) , rotation(EllipsIndicies::getAttributeIndex(desc, attribs.rotation)) - , position(EllipsIndicies::getAttributeIndex(desc, attribs.positionWS)) - , ellipsesGroup(desc.groupIndex(attribs.ellipses)) {} + , position(EllipsIndicies::getAttributeIndex(desc, attribs.positionWS)) {} const size_t stretch, rotation, position; - const points::AttributeSet::Descriptor::GroupIndex ellipsesGroup; private: template @@ -97,26 +95,21 @@ struct EllipsoidTransfer : const RadiusType& rt, const math::Transform& source, SdfT& surface, - const Real sphereScale, const EllipsIndicies& indices, Int64Tree* cpg = nullptr, const std::unordered_map* 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) { @@ -124,7 +117,6 @@ struct EllipsoidTransfer : 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; } @@ -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()); @@ -222,12 +214,10 @@ struct EllipsoidTransfer : } private: - const RealT mSphereScale; const EllipsIndicies& mIndices; std::unique_ptr mStretchHandle; std::unique_ptr mRotationHandle; std::unique_ptr mPositionWSHandle; - std::unique_ptr mEllipsesGroupHandle; }; @@ -245,12 +235,10 @@ 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) {} @@ -258,7 +246,6 @@ struct EllipsSurfaceMaskOp EllipsSurfaceMaskOp(const EllipsSurfaceMaskOp& other, tbb::split) : BaseT(other) , mRadius(other.mRadius) - , mSphereScale(other.mSphereScale) , mHalfband(other.mHalfband) , mIndices(other.mIndices) , mMaxDist(0) {} @@ -280,10 +267,8 @@ struct EllipsSurfaceMaskOp // @todo detect this from PcaSettings and avoid checking if no smoothing points::AttributeHandle pwsHandle(leaf.constAttributeArray(mIndices.position)); points::AttributeHandle 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 @@ -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 @@ -375,17 +325,15 @@ struct EllipsSurfaceMaskOp points::AttributeHandle pwsHandle(leaf.constAttributeArray(mIndices.position)); points::AttributeHandle stretchHandle(leaf.constAttributeArray(mIndices.stretch)); points::AttributeHandle 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::lowest()), minPos(std::numeric_limits::max()); RadiusType rad(mRadius); rad.reset(leaf); - bool hasIsolatedPoints = false; for (Index i = 0; i < stretchHandle.size(); ++i) { @@ -393,57 +341,61 @@ struct EllipsSurfaceMaskOp 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 @@ -473,7 +425,6 @@ struct EllipsSurfaceMaskOp private: const RadiusType& mRadius; - const Real mSphereScale; const Real mHalfband; const EllipsIndicies& mIndices; Vec3i mMaxDist; @@ -501,7 +452,6 @@ rasterizeEllipsoids(const PointDataGridT& points, const std::vector& 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; @@ -547,7 +497,7 @@ rasterizeEllipsoids(const PointDataGridT& points, // pass radius scale as index space EllipsSurfaceMaskOp, 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 @@ -562,7 +512,7 @@ rasterizeEllipsoids(const PointDataGridT& points, grids = doRasterizeSurface (points, attributes, filter, *surface, interrupter, - width, rad, points.transform(), *surface, sphereScale, indices); // args + width, rad, points.transform(), *surface, indices); // args } else { @@ -593,7 +543,7 @@ rasterizeEllipsoids(const PointDataGridT& points, // pass radius scale as index space EllipsSurfaceMaskOp, 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 @@ -608,7 +558,7 @@ rasterizeEllipsoids(const PointDataGridT& points, grids = doRasterizeSurface (points, attributes, filter, *surface, interrupter, - width, rad, points.transform(), *surface, sphereScale, indices); // args + width, rad, points.transform(), *surface, indices); // args } if (interrupter) interrupter->end(); diff --git a/openvdb/openvdb/points/impl/PrincipalComponentAnalysisImpl.h b/openvdb/openvdb/points/impl/PrincipalComponentAnalysisImpl.h index 2695fbbd44..d6ef0b649f 100644 --- a/openvdb/openvdb/points/impl/PrincipalComponentAnalysisImpl.h +++ b/openvdb/openvdb/points/impl/PrincipalComponentAnalysisImpl.h @@ -586,7 +586,7 @@ pca(PointDataGridT& points, // 1) Create persisting attributes const size_t pwsIdx = initAttribute(attrs.positionWS, zeroVal()); const size_t rotIdx = initAttribute(attrs.rotation, zeroVal()); - 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(); diff --git a/openvdb/openvdb/unittest/TestPCA.cc b/openvdb/openvdb/unittest/TestPCA.cc index 0bc24784b8..ab2e0d8aa0 100644 --- a/openvdb/openvdb/unittest/TestPCA.cc +++ b/openvdb/openvdb/unittest/TestPCA.cc @@ -268,6 +268,7 @@ TEST_F(TestPCA, testPCA) s.searchRadius = std::numeric_limits::max(); s.averagePositions = 0.0f; // disable position smoothing s.neighbourThreshold = 3; // more than 2, points should end up as spheres + s.nonAnisotropicStretch = 2.0f; points::pca(*points, s, a); ASSERT_TRUE(points->tree().cbeginLeaf()); @@ -281,7 +282,7 @@ TEST_F(TestPCA, testPCA) CheckPCAAttributeValues( posistions, - {Vec3f(1.0f), Vec3f(1.0f), Vec3f(1.0f)}, + {Vec3f(s.nonAnisotropicStretch), Vec3f(s.nonAnisotropicStretch), Vec3f(s.nonAnisotropicStretch)}, {Mat3s::identity(), Mat3s::identity(), Mat3s::identity()}, posistionsWs, {false, false, false}, @@ -294,6 +295,7 @@ TEST_F(TestPCA, testPCA) s.searchRadius = 0.001f; s.neighbourThreshold = 1; + s.nonAnisotropicStretch = 1.0f; points::pca(*points, s, a); ASSERT_TRUE(points->tree().cbeginLeaf()); @@ -323,6 +325,7 @@ TEST_F(TestPCA, testPCA) s.searchRadius = std::nextafter(0.01f, 1.0f); s.neighbourThreshold = 2; // only center point s.averagePositions = 0.0f; // disable position smoothing + s.nonAnisotropicStretch = 1.3f; points::pca(*points, s, a); ASSERT_TRUE(points->tree().cbeginLeaf()); @@ -336,7 +339,7 @@ TEST_F(TestPCA, testPCA) CheckPCAAttributeValues( posistions, - {Vec3f(1.0), Vec3f(2.51984f, 0.62996f, 0.62996f), Vec3f(1.0)}, + {Vec3f(s.nonAnisotropicStretch), Vec3f(2.51984f, 0.62996f, 0.62996f), Vec3f(s.nonAnisotropicStretch)}, {Mat3s::identity(), Mat3s(0,1,0, 1,0,0, 0,0,1), Mat3s::identity()}, posistionsWs, {false, true, false}, @@ -356,6 +359,7 @@ TEST_F(TestPCA, testPCA) s.neighbourThreshold = 1; // make sure they get classified s.averagePositions = 0.0f; // disable position smoothing s.allowedAnisotropyRatio = 0.25f; + s.nonAnisotropicStretch = 1.0f; points::pca(*points, s, a); ASSERT_TRUE(points->tree().cbeginLeaf()); diff --git a/openvdb/openvdb/unittest/TestPointRasterizeSDF.cc b/openvdb/openvdb/unittest/TestPointRasterizeSDF.cc index 07ff48ab1b..d456901200 100644 --- a/openvdb/openvdb/unittest/TestPointRasterizeSDF.cc +++ b/openvdb/openvdb/unittest/TestPointRasterizeSDF.cc @@ -943,7 +943,6 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) { points::EllipsoidSettings<> s; s.radiusScale = 0.2f; - s.sphereScale = 1.0f; s.halfband = 3; s.transform = nullptr; @@ -959,19 +958,17 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) } // Test single point which is not treated as an ellips. This should give - // identicle results to SphereSettings<> (as long as the sphereScale is 1) + // identicle results to SphereSettings<> { points::EllipsoidSettings<> s; s.radiusScale = 0.2f; - s.sphereScale = 1.0f; s.halfband = 3; s.transform = nullptr; - /// 1) test with a single point but that is not in the ellips ellipses group + /// 1) test with a single point with uniform stretch auto points = PointBuilder({Vec3f(0)}) .voxelsize(0.1) - .group({0}, s.pca.ellipses) // surface as a sphere - .attribute(points::PcaAttributes::StretchT(1.0), s.pca.stretch) + .attribute(points::PcaAttributes::StretchT(1.0), s.pca.stretch) // uniform stretch .attribute(points::PcaAttributes::RotationT::identity(), s.pca.rotation) .attribute(points::PcaAttributes::PosWsT(0.0), s.pca.positionWS) .get(); @@ -999,46 +996,17 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) EXPECT_EQ(sdf->background(), *iter); } - /// 2) test with a single point that is marked as an ellips, but has an identity - /// transformation matrix (and so is essentially a sphere) - points = PointBuilder({Vec3f(0)}) - .voxelsize(0.1) - .group({1}, s.pca.ellipses) // surface as an ellips - .attribute(points::PcaAttributes::StretchT(1.0), s.pca.stretch) - .attribute(points::PcaAttributes::RotationT::identity(), s.pca.rotation) - .attribute(points::PcaAttributes::PosWsT(0.0), s.pca.positionWS) - .get(); - - grids = points::rasterizeSdf(*points, s); - sdf = StaticPtrCast(grids.front()); - EXPECT_TRUE(sdf); - EXPECT_TRUE(sdf->transform() == points->transform()); - EXPECT_EQ(GRID_LEVEL_SET, sdf->getGridClass()); - EXPECT_EQ(float(s.halfband * points->voxelSize()[0]), sdf->background()); - - EXPECT_EQ(Index32(8), sdf->tree().leafCount()); - EXPECT_EQ(Index64(0), sdf->tree().activeTileCount()); - EXPECT_EQ(Index64(485), sdf->tree().activeVoxelCount()); - - for (auto iter = sdf->cbeginValueOn(); iter; ++iter) { - const Vec3d ws = sdf->transform().indexToWorld(iter.getCoord()); - float length = float(ws.length()); // dist to center - length -= float(s.radiusScale); // account for radius - EXPECT_NEAR(length, *iter, 1e-6f); - } - // should only have exterior background for (auto iter = sdf->cbeginValueOff(); iter; ++iter) { EXPECT_EQ(sdf->background(), *iter); } - /// 3) larger radius, larger voxel size, with sphere scale + /// 2) larger radius, larger voxel size s.radiusScale = 1.3f; - s.sphereScale = 1.6f; + float stretch = 1.6f; points = PointBuilder({Vec3f(0)}) .voxelsize(0.5) - .group({0}, s.pca.ellipses) // surface as a sphere - .attribute(points::PcaAttributes::StretchT(1.0), s.pca.stretch) + .attribute(points::PcaAttributes::StretchT(stretch), s.pca.stretch) // uniform stretch .attribute(points::PcaAttributes::RotationT::identity(), s.pca.rotation) .attribute(points::PcaAttributes::PosWsT(0.0), s.pca.positionWS) .get(); @@ -1058,7 +1026,7 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) for (auto iter = sdf->cbeginValueOn(); iter; ++iter) { const Vec3d ws = sdf->transform().indexToWorld(iter.getCoord()); float length = float(ws.length()); // dist to center - length -= float(s.radiusScale * s.sphereScale); // account for radius and sphere scale + length -= float(s.radiusScale * stretch); // account for radius and sphere scale EXPECT_NEAR(length, *iter, 1e-6f); } @@ -1068,7 +1036,7 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) const Vec3d ws = sdf->transform().indexToWorld(iter.getCoord()); float length = float(ws.length()); // dist to center // if length is <= the (rad - halfbandws), voxel is inside the surface - const bool interior = (length <= ((s.radiusScale * s.sphereScale) - (s.halfband * sdf->voxelSize()[0]))); + const bool interior = (length <= ((s.radiusScale * stretch) - (s.halfband * sdf->voxelSize()[0]))); if (interior) EXPECT_EQ(-(sdf->background()), *iter); else EXPECT_EQ(sdf->background(), *iter); interior ? ++interiorOff : ++exteriorOff; @@ -1077,18 +1045,19 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) EXPECT_EQ(size_t(7), interiorOff); EXPECT_EQ(size_t(297441), exteriorOff); - /// 4) offset position, different transform, larger half band - s.sphereScale = 1.0f; + /// 3) offset position, different transform, larger half band + stretch = 1.0f; s.radiusScale = 2.0f; s.halfband = 4; s.transform = math::Transform::createLinearTransform(0.3); + points::PcaAttributes::RotationT rot; + rot.setToRotation({0,1,0}, 45); // arbitrary rotation should have no effect const Vec3f center(-1.2f, 3.4f,-5.6f); points = PointBuilder({Vec3f(center)}) .voxelsize(0.1) - .group({0}, s.pca.ellipses) // surface as a sphere - .attribute(points::PcaAttributes::StretchT(1.0), s.pca.stretch) - .attribute(points::PcaAttributes::RotationT::identity(), s.pca.rotation) + .attribute(points::PcaAttributes::StretchT(stretch), s.pca.stretch) + .attribute(rot, s.pca.rotation) .attribute(points::PcaAttributes::PosWsT(center), s.pca.positionWS) .get(); @@ -1133,7 +1102,6 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) { points::EllipsoidSettings<> s; s.radiusScale = 0.8f; - s.sphereScale = 1.0f; // should have no effect s.halfband = 3; s.transform = nullptr; @@ -1149,7 +1117,6 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) /// 1) test with a single ellips with Y stretch/rotation auto points = PointBuilder({Vec3f(0)}) .voxelsize(0.1) - .group({1}, s.pca.ellipses) // surface as an ellips .attribute(stretch, s.pca.stretch) .attribute(rot, s.pca.rotation) .attribute(points::PcaAttributes::PosWsT(0.0), s.pca.positionWS) @@ -1192,8 +1159,7 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) EXPECT_EQ(size_t(306067), exteriorOff); // Run again with a different Y rotation, should be the same as the above - // as we're only squashing in Y. Also test sphereScale has no impact - s.sphereScale = 10.0f; // should have no effect + // as we're only squashing in Y rot.setToRotation({0,1,0}, -88); points = PointBuilder({Vec3f(0)}) .voxelsize(0.1) @@ -1220,7 +1186,6 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) // along multiple axis { points::EllipsoidSettings<> s; - s.sphereScale = 3.0f; // should have no effect s.halfband = 5; // Design an ellips that is squashed in XYZ and then rotated @@ -1239,7 +1204,6 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) const Vec3f center(-1.2f, 3.4f,-5.6f); auto points = PointBuilder({center}) .voxelsize(0.2) - .group({1}, s.pca.ellipses) // surface as an ellips .attribute(stretch, s.pca.stretch) .attribute(rot, s.pca.rotation) .attribute(points::PcaAttributes::PosWsT(center), s.pca.positionWS) @@ -1286,20 +1250,18 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) { points::EllipsoidSettings<> s; s.radiusScale = 2.0f; - s.sphereScale = 0.8f; s.halfband = 5; s.transform = nullptr; const auto positions = getBoxPoints(); const std::vector positionsVec3d(positions.begin(), positions.end()); - const std::vector ellips {0,1,0,1, 0,1,1,1}; - const std::vector stretches { - {0.0f, 0.0f, 0.0f}, // sphere (should be ignored) + std::vector stretches { + {0.8f, 0.8f, 0.8f}, // sphere, 0.8f {1.0f, 1.0f, 1.0f}, // ellips - {5.0f, 5.0f, 5.0f}, // sphere (should be ignored) + {0.7f, 0.7f, 0.7f}, // sphere, 0.7f {2.0f, 0.3f, 1.5f}, // ellips - {1.0f, 1.0f, 1.0f}, // sphere (should be ignored) + {0.8f, 0.8f, 0.8f}, // sphere, 0.8f {1.1f, 0.8f, 1.0f}, // ellips {0.8f, 4.0f, 1.5f}, // ellips {0.4f, 1.1f, 1.0f} // ellips @@ -1319,7 +1281,6 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) /// 1) test with uniform radius auto points = PointBuilder(positions) .voxelsize(0.3) - .group(ellips, s.pca.ellipses) // surface as an ellips .attribute(stretches, s.pca.stretch) .attribute(rotations, s.pca.rotation) .attribute(positionsVec3d, s.pca.positionWS) @@ -1334,7 +1295,7 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) EXPECT_EQ(Index32(147), sdf->tree().leafCount()); EXPECT_EQ(Index64(0), sdf->tree().activeTileCount()); - EXPECT_EQ(Index64(33252), sdf->tree().activeVoxelCount()); + EXPECT_EQ(Index64(33213), sdf->tree().activeVoxelCount()); // Small lambda that finds the cloest length to a point that could // either be an ellips of sphere, for a given voxel @@ -1347,13 +1308,14 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) math::Mat3s inv = math::Mat3s::identity(); double scale = s.radiusScale / sdf->voxelSize()[0]; if (radii) scale *= double((*radii)[idx]); + auto stretch = stretches[idx]; - if (ellips[idx]) { + if ((stretch.x() != stretch.y()) || (stretch.x() != stretch.z())) { inv = rotations[idx].timesDiagonal(1.0 / stretches[idx]) * rotations[idx].transpose(); } else { - scale *= s.sphereScale; + scale *= stretch.x(); } const Vec3d is = inv * (ijk.asVec3d() - sdf->transform().worldToIndex(pos)); @@ -1384,23 +1346,24 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) } EXPECT_EQ(size_t(147), interiorOff); - EXPECT_EQ(size_t(336622), exteriorOff); + EXPECT_EQ(size_t(336661), exteriorOff); /// 2) test with varying radius - const std::vector radii { 1.0f, 0.0f, 2.0f, 1.1f, 0.2f, 0.5f, 0.8f, 3.0f }; s.radiusScale = 1.2; - s.sphereScale = 0.3; s.radius = "pscale"; s.halfband = 1; + // new sphere radii, note that getClosestLength() reads from this vector + stretches[0] = points::PcaAttributes::StretchT(0.3f); + stretches[4] = points::PcaAttributes::StretchT(0.3f); + points = PointBuilder(positions) .voxelsize(0.3) - .group(ellips, s.pca.ellipses) // surface as an ellips .attribute(radii, s.radius) .attribute(stretches, s.pca.stretch) .attribute(rotations, s.pca.rotation) @@ -1416,7 +1379,7 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) EXPECT_EQ(Index32(39), sdf->tree().leafCount()); EXPECT_EQ(Index64(0), sdf->tree().activeTileCount()); - EXPECT_EQ(Index64(2613), sdf->tree().activeVoxelCount()); + EXPECT_EQ(Index64(2754), sdf->tree().activeVoxelCount()); for (auto iter = sdf->cbeginValueOn(); iter; ++iter) { // get closest dist from all points @@ -1436,8 +1399,8 @@ TEST_F(TestPointRasterizeSDF, testRasterizeEllipsoids) interior ? ++interiorOff : ++exteriorOff; } - EXPECT_EQ(size_t(2137), interiorOff); - EXPECT_EQ(size_t(310083), exteriorOff); + EXPECT_EQ(size_t(2456), interiorOff); + EXPECT_EQ(size_t(309623), exteriorOff); } }