8000 Improve smooth by elalish · Pull Request #828 · elalish/manifold · GitHub
[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Skip to content

Improve smooth #828

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Jun 4, 2024
Merged

Improve smooth #828

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
9f580fb
add max curvature tests
elalish Jun 3, 2024
3d26400
back to radial basis functions
elalish Jun 3, 2024
4c25fca
fix longWay slerp
elalish Jun 4, 2024
5ee4431
cusps eliminated
elalish Jun 4, 2024
f1e4782
cleanup
elalish Jun 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions bindings/wasm/examples/worker.test.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,8 +83,8 @@ suite('Examples', () => {
test('Scallop', async () => {
const result = await runExample('Scallop');
expect(result.genus).to.equal(0, 'Genus');
expect(result.volume).to.be.closeTo(40900, 100, 'Volume');
expect(result.surfaceArea).to.be.closeTo(7740, 10, 'Surface Area');
expect(result.volume).to.be.closeTo(39900, 100, 'Volume');
expect(result.surfaceArea).to.be.closeTo(7930, 10, 'Surface Area');
});

test('Torus Knot', async () => {
Expand Down
132 changes: 73 additions & 59 deletions src/manifold/src/smoothing.cpp
View file 10000
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,13 +41,9 @@ float AngleBetween(glm::vec3 a, glm::vec3 b) {
glm::vec4 CircularTangent(const glm::vec3& tangent, const glm::vec3& edgeVec) {
const glm::vec3 dir = SafeNormalize(tangent);

float weight = glm::abs(glm::dot(dir, SafeNormalize(edgeVec)));
if (weight == 0) {
weight = 1;
}
float weight = glm::max(0.5f, glm::dot(dir, SafeNormalize(edgeVec)));
// Quadratic weighted bezier for circular interpolation
const glm::vec4 bz2 =
weight * glm::vec4(dir * glm::length(edgeVec) / (2 * weight), 1);
const glm::vec4 bz2 = glm::vec4(dir * 0.5f * glm::length(edgeVec), weight);
// Equivalent cubic weighted bezier
const glm::vec4 bz3 = glm::mix(glm::vec4(0, 0, 0, 1), bz2, 2 / 3.0f);
// Convert from homogeneous form to geometric form
Expand Down Expand Up @@ -82,53 +78,72 @@ struct SmoothBezier {
struct InterpTri {
const Manifold::Impl* impl;

glm::vec4 Homogeneous(glm::vec4 v) const {
static glm::vec4 Homogeneous(glm::vec4 v) {
v.x *= v.w;
v.y *= v.w;
v.z *= v.w;
return v;
}

glm::vec4 Homogeneous(glm::vec3 v) const { return glm::vec4(v, 1.0f); }
static glm::vec4 Homogeneous(glm::vec3 v) { return glm::vec4(v, 1.0f); }

glm::vec3 HNormalize(glm::vec4 v) const {
static glm::vec3 HNormalize(glm::vec4 v) {
return v.w == 0 ? v : (glm::vec3(v) / v.w);
}

glm::vec4 Scale(glm::vec4 v, float scale) const {
static glm::vec4 Scale(glm::vec4 v, float scale) {
return glm::vec4(scale * glm::vec3(v), v.w);
}

glm::vec4 Bezier(glm::vec3 point, glm::vec4 tangent) const {
static glm::vec4 Bezier(glm::vec3 point, glm::vec4 tangent) {
return Homogeneous(glm::vec4(point, 0) + tangent);
}

glm::mat2x4 CubicBezier2Linear(glm::vec4 p0, glm::vec4 p1, glm::vec4 p2,
glm::vec4 p3, float x) const {
static glm::mat2x4 CubicBezier2Linear(glm::vec4 p0, glm::vec4 p1,
glm::vec4 p2, glm::vec4 p3, float x) {
glm::mat2x4 out;
glm::vec4 p12 = glm::mix(p1, p2, x);
out[0] = glm::mix(glm::mix(p0, p1, x), p12, x);
out[1] = glm::mix(p12, glm::mix(p2, p3, x), x);
return out;
}

glm::vec3 BezierPoint(glm::mat2x4 points, float x) const {
static glm::vec3 BezierPoint(glm::mat2x4 points, float x) {
return HNormalize(glm::mix(points[0], points[1], x));
}

glm::vec3 BezierTangent(glm::mat2x4 points) const {
static glm::vec3 BezierTangent(glm::mat2x4 points) {
return SafeNormalize(HNormalize(points[1]) - HNormalize(points[0]));
}

glm::vec3 RotateFromTo(glm::vec3 v, glm::quat start, glm::quat end) const {
static glm::vec3 RotateFromTo(glm::vec3 v, glm::quat start, glm::quat end) {
return end * glm::conjugate(start) * v;
}

glm::mat2x4 Bezier2Bezier(const glm::mat2x3& corners,
const glm::mat2x4& tangentsX,
const glm::mat2x4& tangentsY, float x,
const glm::bvec2& pointedEnds,
const glm::vec3& anchor) const {
static glm::quat Slerp(const glm::quat& x, const glm::quat& y, float a,
bool longWay) {
glm::quat z = y;
float cosTheta = glm::dot(x, y);

// Take the long way around the sphere only when requested
if ((cosTheta < 0) != longWay) {
z = -y;
cosTheta = -cosTheta;
}

if (cosTheta > 1.0f - glm::epsilon<float>()) {
return glm::lerp(x, z, a); // for numerical stability
} else {
float angle = glm::acos(cosTheta);
return (glm::sin((1.0f - a) * angle) * x + glm::sin(a * angle) * z) /
glm::sin(angle);
}
}

static glm::mat2x4 Bezier2Bezier(const glm::mat2x3& corners,
const glm::mat2x4& tangentsX,
const glm::mat2x4& tangentsY, float x,
const glm::vec3& anchor) {
const glm::mat2x4 bez = CubicBezier2Linear(
Homogeneous(corners[0]), Bezier(corners[0], tangentsX[0]),
Bezier(corners[1], tangentsX[1]), Homogeneous(corners[1]), x);
Expand All @@ -151,45 +166,39 @@ struct InterpTri {
const glm::quat q1 =
glm::quat_cast(glm::mat3(nTangentsX[1], biTangents[1],
glm::cross(nTangentsX[1], biTangents[1])));
const glm::quat qTmp = glm::slerp(q0, q1, x);
const glm::vec3 edge = corners[1] - corners[0];
const bool longWay =
glm::dot(nTangentsX[0], edge) + glm::dot(nTangentsX[1], edge) < 0;
const glm::quat qTmp = Slerp(q0, q1, x, longWay);
const glm::quat q =
glm::rotation(qTmp * glm::vec3(1, 0, 0), tangent) * qTmp;

const glm::vec3 end0 = pointedEnds[0]
? glm::vec3(0)
: RotateFromTo(glm::vec3(tangentsY[0]), q0, q);
const glm::vec3 end1 = pointedEnds[1]
? glm::vec3(0)
: RotateFromTo(glm::vec3(tangentsY[1]), q1, q);
const glm::vec3 delta = glm::mix(end0, end1, x);
const glm::vec3 delta =
glm::mix(RotateFromTo(glm::vec3(tangentsY[0]), q0, q),
RotateFromTo(glm::vec3(tangentsY[1]), q1, q), x);
const float deltaW = glm::mix(tangentsY[0].w, tangentsY[1].w, x);

return {Homogeneous(end), glm::vec4(delta, deltaW)};
}

glm::vec3 Bezier2D(const glm::mat4x3& corners, const glm::mat4& tangentsX,
const glm::mat4& tangentsY, float x, float y,
const glm::vec3& centroid, bool isTriangle) const {
glm::mat2x4 bez0 = Bezier2Bezier(
{corners[0], corners[1]}, {tangentsX[0], tangentsX[1]},
{tangentsY[0], tangentsY[1]}, x, {isTriangle, false}, centroid);
glm::mat2x4 bez1 = Bezier2Bezier(
{corners[2], corners[3]}, {tangentsX[2], tangentsX[3]},
{tangentsY[2], tangentsY[3]}, 1 - x, {false, isTriangle}, centroid);

const float flatLen =
isTriangle ? x * glm::length(corners[1] - corners[2])
: glm::length(glm::mix(corners[0], corners[1], x) -
glm::mix(corners[3], corners[2], x));
const float scale = glm::length(glm::vec3(bez0[0] - bez1[0])) / flatLen;

const glm::mat2x4 bez = CubicBezier2Linear(
bez0[0], Bezier(glm::vec3(bez0[0]), Scale(bez0[1], scale)),
Bezier(glm::vec3(bez1[0]), Scale(bez1[1], scale)), bez1[0], y);
static glm::vec3 Bezier2D(const glm::mat4x3& corners,
const glm::mat4& tangentsX,
const glm::mat4& tangentsY, float x, float y,
const glm::vec3& centroid) {
glm::mat2x4 bez0 =
Bezier2Bezier({corners[0], corners[1]}, {tangentsX[0], tangentsX[1]},
{tangentsY[0], tangentsY[1]}, x, centroid);
glm::mat2x4 bez1 =
Bezier2Bezier({corners[2], corners[3]}, {tangentsX[2], tangentsX[3]},
{tangentsY[2], tangentsY[3]}, 1 - x, centroid);

const glm::mat2x4 bez =
CubicBezier2Linear(bez0[0], Bezier(glm::vec3(bez0[0]), bez0[1]),
Bezier(glm::vec3(bez1[0]), bez1[1]), bez1[0], y);
return BezierPoint(bez, y);
}

void operator()(thrust::tuple<glm::vec3&, Barycentric> inOut) {
void operator()(thrust::tuple<glm::vec3&, Barycentric> inOut) const {
glm::vec3& pos = thrust::get<0>(inOut);
const int tri = thrust::get<1>(inOut).tri;
const glm::vec4 uvw = thrust::get<1>(inOut).uvw;
Expand Down Expand Up @@ -226,12 +235,17 @@ struct InterpTri {
const int j = Next3(i);
const int k = Prev3(i);
const float x = uvw[k] / (1 - uvw[i]);
const glm::vec3 p =
Bezier2D({corners[i], corners[j], corners[k], corners[i]},
{tangentR[i], tangentL[j], tangentR[k], tangentL[i]},
{tangentL[i], tangentR[j], tangentL[k], tangentR[i]},
1 - uvw[i], x, centroid, true);
posH += Homogeneous(glm::vec4(p, uvw[i]));

const glm::mat2x4 bez =
Bezier2Bezier({corners[j], corners[k]}, {tangentR[j], tangentL[k]},
{tangentL[j], tangentR[k]}, x, centroid);

const glm::mat2x4 bez1 = CubicBezier2Linear(
bez[0], Bezier(glm::vec3(bez[0]), bez[1]),
Bezier(corners[i], glm::mix(tangentR[i], tangentL[i], x)),
Homogeneous(corners[i]), uvw[i]);
const glm::vec3 p = BezierPoint(bez1, uvw[i]);
posH += Homogeneous(glm::vec4(p, uvw[j] * uvw[k]));
}
} else { // quad
const glm::mat4 tangentsX = {
Expand All @@ -248,12 +262,12 @@ struct InterpTri {
const float x = uvw[1] + uvw[2];
const float y = uvw[2] + uvw[3];
const glm::vec3 pX =
Bezier2D(corners, tangentsX, tangentsY, x, y, centroid, false);
Bezier2D(corners, tangentsX, tangentsY, x, y, centroid);
const glm::vec3 pY =
Bezier2D({corners[1], corners[2], corners[3], corners[0]},
{tangentsY[1], tangentsY[2], tangentsY[3], tangentsY[0]},
{tangentsX[1], tangentsX[2], tangentsX[3], tangentsX[0]}, y,
1 - x, centroid, false);
1 - x, centroid);
posH += Homogeneous(glm::vec4(pX, x * (1 - x)));
posH += Homogeneous(glm::vec4(pY, y * (1 - y)));
}
Expand Down Expand Up @@ -885,7 +899,7 @@ void Manifold::Impl::CreateTangents(std::vector<Smoothness> sharpenedEdges) {
&triIsFlatFace](int v) {
auto it = vertTangents.find(v);
if (it == vertTangents.end()) {
fixedHalfedge[vertHalfedge[v]] == true;
fixedHalfedge[vertHalfedge[v]] = true;
return;
}
const std::vector<Pair>& vert = it->second;
Expand Down Expand Up @@ -918,7 +932,7 @@ void Manifold::Impl::CreateTangents(std::vector<Smoothness> sharpenedEdges) {
}
});
} else { // Sharpen vertex uniformly
fixedHalfedge[vertHalfedge[v]] == true;
fixedHalfedge[vertHalfedge[v]] = true;
float smoothness = 0;
float denom = 0;
for (const Pair& pair : vert) {
Expand Down
4 changes: 2 additions & 2 deletions test/samples_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -116,8 +116,8 @@ TEST(Samples, Scallop) {
3, colorCurvature);
CheckNormals(scallop);
auto prop = scallop.GetProperties();
EXPECT_NEAR(prop.volume, 40.9, 0.1);
EXPECT_NEAR(prop.surfaceArea, 77.4, 0.1);
EXPECT_NEAR(prop.volume, 39.9, 0.1);
EXPECT_NEAR(prop.surfaceArea, 79.3, 0.1);
CheckGL(scallop);

#ifdef MANIFOLD_EXPORT
Expand Down
65 changes: 45 additions & 20 deletions test/smooth_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,8 +30,16 @@ TEST(Smooth, Tetrahedron) {
smooth = smooth.Refine(n);
ExpectMeshes(smooth, {{2 * n * n + 2, 4 * n * n}});
auto prop = smooth.GetProperties();
EXPECT_NEAR(prop.volume, 18.7, 0.1);
EXPECT_NEAR(prop.surfaceArea, 34.5, 0.1);
EXPECT_NEAR(prop.volume, 17.0, 0.1);
EXPECT_NEAR(prop.surfaceArea, 32.9, 0.1);

MeshGL out = smooth.CalculateCurvature(-1, 0).GetMeshGL();
float maxMeanCurvature = 0;
for (int i = 3; i < out.vertProperties.size(); i += 4) {
maxMeanCurvature =
glm::max(maxMeanCurvature, glm::abs(out.vertProperties[i]));
}
EXPECT_NEAR(maxMeanCurvature, 4.73, 0.01);

#ifdef MANIFOLD_EXPORT
if (options.exportModels) ExportMesh("smoothTet.glb", smooth.GetMesh(), {});
Expand Down Expand Up @@ -72,8 +80,8 @@ TEST(Smooth, TruncatedCone) {
Manifold cone = Manifold::Cylinder(5, 10, 5, 12).SmoothOut();
Manifold smooth = cone.RefineToLength(0.5).CalculateNormals(0);
auto prop = smooth.GetProperties();
EXPECT_NEAR(prop.volume, 1063.61, 0.01);
EXPECT_NEAR(prop.surfaceArea, 751.78, 0.01);
EXPECT_NEAR(prop.volume, 1062.27, 0.01);
EXPECT_NEAR(prop.surfaceArea, 751.46, 0.01);
MeshGL out = smooth.GetMeshGL();
CheckGL(out);

Expand All @@ -95,8 +103,16 @@ TEST(Smooth, ToLength) {
smooth = smooth.RefineToLength(0.1);
ExpectMeshes(smooth, {{85250, 170496}});
auto prop = smooth.GetProperties();
EXPECT_NEAR(prop.volume, 4842, 1);
EXPECT_NEAR(prop.surfaceArea, 1400, 1);
EXPECT_NEAR(prop.volume, 4604, 1);
EXPECT_NEAR(prop.surfaceArea, 1356, 1);

MeshGL out = smooth.CalculateCurvature(-1, 0).GetMeshGL();
float maxMeanCurvature = 0;
for (int i = 3; i < out.vertProperties.size(); i += 4) {
maxMeanCurvature =
glm::max(maxMeanCurvature, glm::abs(out.vertProperties[i]));
}
EXPECT_NEAR(maxMeanCurvature, 1.67, 0.01);

#ifdef MANIFOLD_EXPORT
if (options.exportModels)
Expand All @@ -106,7 +122,7 @@ TEST(Smooth, ToLength) {

TEST(Smooth, Sphere) {
int n[5] = {4, 8, 16, 32, 64};
float precision[5] = {0.006, 0.003, 0.003, 0.0005, 0.00006};
float precision[5] = {0.04, 0.003, 0.003, 0.0005, 0.00006};
for (int i = 0; i < 5; ++i) {
Manifold sphere = Manifold::Sphere(1, n[i]);
// Refine(odd) puts a center point in the triangle, which is the worst case.
Expand Down Expand Up @@ -156,8 +172,8 @@ TEST(Smooth, Manual) {

ExpectMeshes(interp, {{40002, 80000}});
auto prop = interp.GetProperties();
EXPECT_NEAR(prop.volume, 3.83, 0.01);
EXPECT_NEAR(prop.surfaceArea, 11.95, 0.01);
EXPECT_NEAR(prop.volume, 3.74, 0.01);
EXPECT_NEAR(prop.surfaceArea, 11.78, 0.01);

#ifdef MANIFOLD_EXPORT
if (options.exportModels) {
Expand Down Expand Up @@ -200,8 +216,8 @@ TEST(Smooth, Csaszar) {
csaszar = csaszar.Refine(100);
ExpectMeshes(csaszar, {{70000, 140000}});
auto prop = csaszar.GetProperties();
EXPECT_NEAR(prop.volume, 89873, 10);
EXPECT_NEAR(prop.surfaceArea, 15301, 10);
EXPECT_NEAR(prop.volume, 79890, 10);
EXPECT_NEAR(prop.surfaceArea, 11950, 10);

#ifdef MANIFOLD_EXPORT
if (options.exportModels) {
Expand Down Expand Up @@ -276,19 +292,28 @@ TEST(Smooth, Torus) {
}
}

Manifold smooth = Manifold(torusMesh).RefineToLength(0.1).CalculateNormals(0);
Mesh out = smooth.GetMesh();
for (glm::vec3 v : out.vertPos) {
Manifold smooth = Manifold(torusMesh)
.RefineToLength(0.1)
.CalculateCurvature(-1, 0)
.CalculateNormals(1);
MeshGL out = smooth.GetMeshGL();
float maxMeanCurvature = 0;
for (int i = 0; i < out.vertProperties.size(); i += 7) {
glm::vec3 v(out.vertProperties[i], out.vertProperties[i + 1],
out.vertProperties[i + 2]);
glm::vec3 p(v.x, v.y, 0);
p = glm::normalize(p) * 2.0f;
float r = glm::length(v - p);
ASSERT_NEAR(r, 1, 0.005);
ASSERT_NEAR(r, 1, 0.006);
maxMeanCurvature =
glm::max(maxMeanCurvature, glm::abs(out.vertProperties[i + 3]));
}
EXPECT_NEAR(maxMeanCurvature, 1.63, 0.01);

#ifdef MANIFOLD_EXPORT
ExportOptions options2;
options2.faceted = false;
options2.mat.normalChannels = {3, 4, 5};
options2.mat.normalChannels = {4, 5, 6};
options2.mat.roughness = 0;
if (options.exportModels) ExportMesh("smoothTorus.glb", out, options2);
#endif
Expand All @@ -307,7 +332,7 @@ TEST(Smooth, SineSurface) {
Manifold(surface).CalculateNormals(0, 50).SmoothByNormals(0).Refine(8);
auto prop = smoothed.GetProperties();
EXPECT_NEAR(prop.volume, 7.89, 0.01);
EXPECT_NEAR(prop.surfaceArea, 30.61, 0.01);
EXPECT_NEAR(prop.surfaceArea, 30.60, 0.01);
EXPECT_EQ(smoothed.Genus(), 0);

Manifold smoothed1 = Manifold(surface).SmoothOut(50).Refine(8);
Expand All @@ -318,14 +343,14 @@ TEST(Smooth, SineSurface) {

Manifold smoothed2 = Manifold(surface).SmoothOut(180, 1).Refine(8);
auto prop2 = smoothed2.GetProperties();
EXPECT_NEAR(prop2.volume, 9.06, 0.01);
EXPECT_NEAR(prop2.surfaceArea, 33.81, 0.01);
EXPECT_NEAR(prop2.volume, 9.02, 0.01);
EXPECT_NEAR(prop2.surfaceArea, 33.56, 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.77, 0.01);
EXPECT_NEAR(prop3.surfaceArea, 31.66, 0.01);
EXPECT_EQ(smoothed3.Genus(), 0);

#ifdef MANIFOLD_EXPORT
Expand Down
Loading
0