From 980b83007f2dd75b1acb79ee00a9c34a355085a8 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Thu, 13 Jun 2024 14:32:37 -0700 Subject: [PATCH 1/2] fixed several bugs --- src/manifold/src/smoothing.cpp | 61 ++++++++++++++++++++++------------ test/smooth_test.cpp | 12 +++---- 2 files changed, 45 insertions(+), 28 deletions(-) diff --git a/src/manifold/src/smoothing.cpp b/src/manifold/src/smoothing.cpp index 86f4f942a..bba0863bf 100644 --- a/src/manifold/src/smoothing.cpp +++ b/src/manifold/src/smoothing.cpp @@ -21,16 +21,21 @@ namespace { using namespace manifold; -glm::vec3 OrthogonalTo(glm::vec3 in, glm::vec3 ref) { - in -= glm::dot(in, ref) * ref; - return in; +// Returns a normalized vector orthogonal to ref, in the plane of ref and in, +// unless in and ref are colinear, in which case it falls back to the plane of +// ref and altIn. +glm::vec3 OrthogonalTo(glm::vec3 in, glm::vec3 altIn, glm::vec3 ref) { + glm::vec3 out = in - glm::dot(in, ref) * ref; + if (glm::dot(out, out) < kTolerance * glm::dot(in, in)) { + out = altIn - glm::dot(altIn, ref) * ref; + } + return SafeNormalize(out); } // Get the angle between two unit-vectors. float AngleBetween(glm::vec3 a, glm::vec3 b) { const float dot = glm::dot(a, b); - return dot >= 1 ? kTolerance - : (dot <= -1 ? glm::pi() : glm::acos(dot)); + return dot >= 1 ? 0 : (dot <= -1 ? glm::pi() : glm::acos(dot)); } // Calculate a tangent vector in the form of a weighted cubic Bezier taking as @@ -146,12 +151,10 @@ struct InterpTri { const glm::mat2x3 nTangentsX(SafeNormalize(glm::vec3(tangentsX[0])), -SafeNormalize(glm::vec3(tangentsX[1]))); const glm::mat2x3 biTangents = { - SafeNormalize(OrthogonalTo( - glm::vec3(tangentsY[0]) + kTolerance * (anchor - corners[0]), - nTangentsX[0])), - SafeNormalize(OrthogonalTo( - glm::vec3(tangentsY[1]) + kTolerance * (anchor - corners[1]), - nTangentsX[1]))}; + OrthogonalTo(glm::vec3(tangentsY[0]), (anchor - corners[0]), + nTangentsX[0]), + OrthogonalTo(glm::vec3(tangentsY[1]), (anchor - corners[1]), + nTangentsX[1])}; const glm::quat q0 = glm::quat_cast(glm::mat3(nTangentsX[0], biTangents[0], @@ -580,8 +583,8 @@ void Manifold::Impl::SetNormals(int normalIdx, float minSharpAngle) { // more constant curvature to meet the opposite normal. Achieve // this by pointing the tangent toward the opposite bezier // control point instead of the vert itself. - pos += glm::vec3( - CircularTangent(OrthogonalTo(edgeVec, normal), edgeVec)); + pos += glm::vec3(TangentFromNormal( + normal, halfedge_[current].pairedHalfedge)); } return FaceEdge( {halfedge_[current].face, SafeNormalize(pos - centerPos)}); @@ -720,8 +723,9 @@ void Manifold::Impl::DistributeTangents(const Vec& fixedHalfedges) { const glm::vec3 center = vertPos_[halfedge_[halfedge].startVert]; glm::vec3 lastEdgeVec = SafeNormalize(vertPos_[halfedge_[halfedge].endVert] - center); - glm::vec3 lastTangent = + const glm::vec3 firstTangent = SafeNormalize(glm::vec3(halfedgeTangent_[halfedge])); + glm::vec3 lastTangent = firstTangent; int current = halfedge; do { current = NextHalfedge(halfedge_[current].pairedHalfedge); @@ -730,16 +734,20 @@ void Manifold::Impl::DistributeTangents(const Vec& fixedHalfedges) { SafeNormalize(vertPos_[halfedge_[current].endVert] - center); const glm::vec3 thisTangent = SafeNormalize(glm::vec3(halfedgeTangent_[current])); - const glm::vec3 cp = glm::cross(thisTangent, lastTangent); - normal += cp; + normal += glm::cross(thisTangent, lastTangent); // cumulative sum desiredAngle.push_back( AngleBetween(thisEdgeVec, lastEdgeVec) + (desiredAngle.size() > 0 ? desiredAngle.back() : 0)); - currentAngle.push_back( - AngleBetween(thisTangent, lastTangent) * - glm::sign(glm::dot(cp, approxNormal)) + - (currentAngle.size() > 0 ? currentAngle.back() : 0)); + if (current == halfedge) { + currentAngle.push_back(glm::two_pi()); + } else { + currentAngle.push_back(AngleBetween(thisTangent, firstTangent)); + if (glm::dot(approxNormal, glm::cross(thisTangent, firstTangent)) < + 0) { + currentAngle.back() = glm::two_pi() - currentAngle.back(); + } + } lastEdgeVec = thisEdgeVec; lastTangent = thisTangent; } while (!fixedHalfedges[current]); @@ -760,8 +768,17 @@ void Manifold::Impl::DistributeTangents(const Vec& fixedHalfedges) { do { current = NextHalfedge(halfedge_[current].pairedHalfedge); if (IsMarkedInsideQuad(current)) continue; - const float angle = - currentAngle[i] - scale * desiredAngle[i] - offset; + desiredAngle[i] *= scale; + const float lastAngle = i > 0 ? desiredAngle[i - 1] : 0; + // shrink obtuse angles + if (desiredAngle[i] - lastAngle > glm::pi()) { + desiredAngle[i] = lastAngle + glm::pi(); + } else if (i + 1 < desiredAngle.size() && + scale * desiredAngle[i + 1] - desiredAngle[i] > + glm::pi()) { + desiredAngle[i] = scale * desiredAngle[i + 1] - glm::pi(); + } + const float angle = currentAngle[i] - desiredAngle[i] - offset; glm::vec3 tangent(halfedgeTangent_[current]); halfedgeTangent_[current] = glm::vec4( glm::rotate(tangent, angle, normal), halfedgeTangent_[current].w); diff --git a/test/smooth_test.cpp b/test/smooth_test.cpp index c4477208d..3a9287f15 100644 --- a/test/smooth_test.cpp +++ b/test/smooth_test.cpp @@ -339,8 +339,8 @@ TEST(Smooth, SineSurface) { Manifold smoothed = Manifold(surface).CalculateNormals(0, 50).SmoothByNormals(0).Refine(8); auto prop = smoothed.GetProperties(); - EXPECT_NEAR(prop.volume, 8.07, 0.01); - EXPECT_NEAR(prop.surfaceArea, 30.88, 0.01); + EXPECT_NEAR(prop.volume, 8.09, 0.01); + EXPECT_NEAR(prop.surfaceArea, 30.93, 0.01); EXPECT_EQ(smoothed.Genus(), 0); Manifold smoothed1 = Manifold(surface).SmoothOut(50).Refine(8); @@ -351,14 +351,14 @@ TEST(Smooth, SineSurface) { Manifold smoothed2 = Manifold(surface).SmoothOut(180, 1).Refine(8); auto prop2 = smoothed2.GetProperties(); - EXPECT_NEAR(prop2.volume, 9.02, 0.01); - EXPECT_NEAR(prop2.surfaceArea, 33.56, 0.01); + EXPECT_NEAR(prop2.volume, 8.99, 0.01); + EXPECT_NEAR(prop2.surfaceArea, 33.61, 0.01); EXPECT_EQ(smoothed2.Genus(), 0); Manifold smoothed3 = Manifold(surface).SmoothOut(50, 0.5).Refine(8); auto prop3 = smoothed3.GetProperties(); - EXPECT_NEAR(prop3.volume, 8.46, 0.01); - EXPECT_NEAR(prop3.surfaceArea, 31.66, 0.01); + EXPECT_NEAR(prop3.volume, 8.43, 0.01); + EXPECT_NEAR(prop3.surfaceArea, 31.72, 0.01); EXPECT_EQ(smoothed3.Genus(), 0); #ifdef MANIFOLD_EXPORT From 10d84fecd840117145a568abed9b75500e583938 Mon Sep 17 00:00:00 2001 From: Emmett Lalish Date: Thu, 13 Jun 2024 15:03:47 -0700 Subject: [PATCH 2/2] fixed average rotation --- src/manifold/src/smoothing.cpp | 8 +++++++- test/smooth_test.cpp | 2 +- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/src/manifold/src/smoothing.cpp b/src/manifold/src/smoothing.cpp index bba0863bf..3cc0d9091 100644 --- a/src/manifold/src/smoothing.cpp +++ b/src/manifold/src/smoothing.cpp @@ -32,6 +32,12 @@ glm::vec3 OrthogonalTo(glm::vec3 in, glm::vec3 altIn, glm::vec3 ref) { return SafeNormalize(out); } +float Wrap(float radians) { + return radians < -glm::pi() ? radians + glm::two_pi() + : radians > glm::pi() ? radians - glm::two_pi() + : radians; +} + // Get the angle between two unit-vectors. float AngleBetween(glm::vec3 a, glm::vec3 b) { const float dot = glm::dot(a, b); @@ -758,7 +764,7 @@ void Manifold::Impl::DistributeTangents(const Vec& fixedHalfedges) { float offset = 0; if (current == halfedge) { // only one - find average offset for (int i = 0; i < currentAngle.size(); ++i) { - offset += currentAngle[i] - scale * desiredAngle[i]; + offset += Wrap(currentAngle[i] - scale * desiredAngle[i]); } offset /= currentAngle.size(); } diff --git a/test/smooth_test.cpp b/test/smooth_test.cpp index 3a9287f15..ea8c08b45 100644 --- a/test/smooth_test.cpp +++ b/test/smooth_test.cpp @@ -351,7 +351,7 @@ TEST(Smooth, SineSurface) { Manifold smoothed2 = Manifold(surface).SmoothOut(180, 1).Refine(8); auto prop2 = smoothed2.GetProperties(); - EXPECT_NEAR(prop2.volume, 8.99, 0.01); + EXPECT_NEAR(prop2.volume, 9.00, 0.01); EXPECT_NEAR(prop2.surfaceArea, 33.61, 0.01); EXPECT_EQ(smoothed2.Genus(), 0);