Skip to content

Commit bdae9dd

Browse files
elalishelalish
authored
Improve smooth (#828)
* add max curvature tests * back to radial basis functions * fix longWay slerp * cusps eliminated * cleanup
1 parent 5f4059a commit bdae9dd

File tree

4 files changed

+122
-83
lines changed

4 files changed

+122
-83
lines changed

‎bindings/wasm/examples/worker.test.ts

+2-2
Original file line numberDiff line numberDiff line change
@@ -83,8 +83,8 @@ suite('Examples', () => {
8383
test('Scallop', async () => {
8484
const result = await runExample('Scallop');
8585
expect(result.genus).to.equal(0, 'Genus');
86-
expect(result.volume).to.be.closeTo(40900, 100, 'Volume');
87-
expect(result.surfaceArea).to.be.closeTo(7740, 10, 'Surface Area');
86+
expect(result.volume).to.be.closeTo(39900, 100, 'Volume');
87+
expect(result.surfaceArea).to.be.closeTo(7930, 10, 'Surface Area');
8888
});
8989

9090
test('Torus Knot', async () => {

‎src/manifold/src/smoothing.cpp

+73-59
Original file line numberDiff line numberDiff line change
@@ -41,13 +41,9 @@ float AngleBetween(glm::vec3 a, glm::vec3 b) {
4141
glm::vec4 CircularTangent(const glm::vec3& tangent, const glm::vec3& edgeVec) {
4242
const glm::vec3 dir = SafeNormalize(tangent);
4343

44-
float weight = glm::abs(glm::dot(dir, SafeNormalize(edgeVec)));
45-
if (weight == 0) {
46-
weight = 1;
47-
}
44+
float weight = glm::max(0.5f, glm::dot(dir, SafeNormalize(edgeVec)));
4845
// Quadratic weighted bezier for circular interpolation
49-
const glm::vec4 bz2 =
50-
weight * glm::vec4(dir * glm::length(edgeVec) / (2 * weight), 1);
46+
const glm::vec4 bz2 = glm::vec4(dir * 0.5f * glm::length(edgeVec), weight);
5147
// Equivalent cubic weighted bezier
5248
const glm::vec4 bz3 = glm::mix(glm::vec4(0, 0, 0, 1), bz2, 2 / 3.0f);
5349
// Convert from homogeneous form to geometric form
@@ -82,53 +78,72 @@ struct SmoothBezier {
8278
struct InterpTri {
8379
const Manifold::Impl* impl;
8480

85-
glm::vec4 Homogeneous(glm::vec4 v) const {
81+
static glm::vec4 Homogeneous(glm::vec4 v) {
8682
v.x *= v.w;
8783
v.y *= v.w;
8884
v.z *= v.w;
8985
return v;
9086
}
9187

92-
glm::vec4 Homogeneous(glm::vec3 v) const { return glm::vec4(v, 1.0f); }
88+
static glm::vec4 Homogeneous(glm::vec3 v) { return glm::vec4(v, 1.0f); }
9389

94-
glm::vec3 HNormalize(glm::vec4 v) const {
90+
static glm::vec3 HNormalize(glm::vec4 v) {
9591
return v.w == 0 ? v : (glm::vec3(v) / v.w);
9692
}
9793

98-
glm::vec4 Scale(glm::vec4 v, float scale) const {
94+
static glm::vec4 Scale(glm::vec4 v, float scale) {
9995
return glm::vec4(scale * glm::vec3(v), v.w);
10096
}
10197

102-
glm::vec4 Bezier(glm::vec3 point, glm::vec4 tangent) const {
98+
static glm::vec4 Bezier(glm::vec3 point, glm::vec4 tangent) {
10399
return Homogeneous(glm::vec4(point, 0) + tangent);
104100
}
105101

106-
glm::mat2x4 CubicBezier2Linear(glm::vec4 p0, glm::vec4 p1, glm::vec4 p2,
107-
glm::vec4 p3, float x) const {
102+
static glm::mat2x4 CubicBezier2Linear(glm::vec4 p0, glm::vec4 p1,
103+
glm::vec4 p2, glm::vec4 p3, float x) {
108104
glm::mat2x4 out;
109105
glm::vec4 p12 = glm::mix(p1, p2, x);
110106
out[0] = glm::mix(glm::mix(p0, p1, x), p12, x);
111107
out[1] = glm::mix(p12, glm::mix(p2, p3, x), x);
112108
return out;
113109
}
114110

115-
glm::vec3 BezierPoint(glm::mat2x4 points, float x) const {
111+
static glm::vec3 BezierPoint(glm::mat2x4 points, float x) {
116112
return HNormalize(glm::mix(points[0], points[1], x));
117113
}
118114

119-
glm::vec3 BezierTangent(glm::mat2x4 points) const {
115+
static glm::vec3 BezierTangent(glm::mat2x4 points) {
120116
return SafeNormalize(HNormalize(points[1]) - HNormalize(points[0]));
121117
}
122118

123-
glm::vec3 RotateFromTo(glm::vec3 v, glm::quat start, glm::quat end) const {
119+
static glm::vec3 RotateFromTo(glm::vec3 v, glm::quat start, glm::quat end) {
124120
return end * glm::conjugate(start) * v;
125121
}
126122

127-
glm::mat2x4 Bezier2Bezier(const glm::mat2x3& corners,
128-
const glm::mat2x4& tangentsX,
129-
const glm::mat2x4& tangentsY, float x,
130-
const glm::bvec2& pointedEnds,
131-
const glm::vec3& anchor) const {
123+
static glm::quat Slerp(const glm::quat& x, const glm::quat& y, float a,
124+
bool longWay) {
125+
glm::quat z = y;
126+
float cosTheta = glm::dot(x, y);
127+
128+
// Take the long way around the sphere only when requested
129+
if ((cosTheta < 0) != longWay) {
130+
z = -y;
131+
cosTheta = -cosTheta;
132+
}
133+
134+
if (cosTheta > 1.0f - glm::epsilon<float>()) {
135+
return glm::lerp(x, z, a); // for numerical stability
136+
} else {
137+
float angle = glm::acos(cosTheta);
138+
return (glm::sin((1.0f - a) * angle) * x + glm::sin(a * angle) * z) /
139+
glm::sin(angle);
140+
}
141+
}
142+
143+
static glm::mat2x4 Bezier2Bezier(const glm::mat2x3& corners,
144+
const glm::mat2x4& tangentsX,
145+
const glm::mat2x4& tangentsY, float x,
146+
const glm::vec3& anchor) {
132147
const glm::mat2x4 bez = CubicBezier2Linear(
133148
Homogeneous(corners[0]), Bezier(corners[0], tangentsX[0]),
134149
Bezier(corners[1], tangentsX[1]), Homogeneous(corners[1]), x);
@@ -151,45 +166,39 @@ struct InterpTri {
151166
const glm::quat q1 =
152167
glm::quat_cast(glm::mat3(nTangentsX[1], biTangents[1],
153168
glm::cross(nTangentsX[1], biTangents[1])));
154-
const glm::quat qTmp = glm::slerp(q0, q1, x);
169+
const glm::vec3 edge = corners[1] - corners[0];
170+
const bool longWay =
171+
glm::dot(nTangentsX[0], edge) + glm::dot(nTangentsX[1], edge) < 0;
172+
const glm::quat qTmp = Slerp(q0, q1, x, longWay);
155173
const glm::quat q =
156174
glm::rotation(qTmp * glm::vec3(1, 0, 0), tangent) * qTmp;
157175

158-
const glm::vec3 end0 = pointedEnds[0]
159-
? glm::vec3(0)
160-
: RotateFromTo(glm::vec3(tangentsY[0]), q0, q);
161-
const glm::vec3 end1 = pointedEnds[1]
162-
? glm::vec3(0)
163-
: RotateFromTo(glm::vec3(tangentsY[1]), q1, q);
164-
const glm::vec3 delta = glm::mix(end0, end1, x);
176+
const glm::vec3 delta =
177+
glm::mix(RotateFromTo(glm::vec3(tangentsY[0]), q0, q),
178+
RotateFromTo(glm::vec3(tangentsY[1]), q1, q), x);
165179
const float deltaW = glm::mix(tangentsY[0].w, tangentsY[1].w, x);
166180

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

170-
glm::vec3 Bezier2D(const glm::mat4x3& corners, const glm::mat4& tangentsX,
171-
const glm::mat4& tangentsY, float x, float y,
172-
const glm::vec3& centroid, bool isTriangle) const {
173-
glm::mat2x4 bez0 = Bezier2Bezier(
174-
{corners[0], corners[1]}, {tangentsX[0], tangentsX[1]},
175-
{tangentsY[0], tangentsY[1]}, x, {isTriangle, false}, centroid);
176-
glm::mat2x4 bez1 = Bezier2Bezier(
177-
{corners[2], corners[3]}, {tangentsX[2], tangentsX[3]},
178-
{tangentsY[2], tangentsY[3]}, 1 - x, {false, isTriangle}, centroid);
179-
180-
const float flatLen =
181-
isTriangle ? x * glm::length(corners[1] - corners[2])
182-
: glm::length(glm::mix(corners[0], corners[1], x) -
183-
glm::mix(corners[3], corners[2], x));
184-
const float scale = glm::length(glm::vec3(bez0[0] - bez1[0])) / flatLen;
185-
186-
const glm::mat2x4 bez = CubicBezier2Linear(
187-
bez0[0], Bezier(glm::vec3(bez0[0]), Scale(bez0[1], scale)),
188-
Bezier(glm::vec3(bez1[0]), Scale(bez1[1], scale)), bez1[0], y);
184+
static glm::vec3 Bezier2D(const glm::mat4x3& corners,
185+
const glm::mat4& tangentsX,
186+
const glm::mat4& tangentsY, float x, float y,
187+
const glm::vec3& centroid) {
188+
glm::mat2x4 bez0 =
189+
Bezier2Bezier({corners[0], corners[1]}, {tangentsX[0], tangentsX[1]},
190+
{tangentsY[0], tangentsY[1]}, x, centroid);
191+
glm::mat2x4 bez1 =
192+
Bezier2Bezier({corners[2], corners[3]}, {tangentsX[2], tangentsX[3]},
193+
{tangentsY[2], tangentsY[3]}, 1 - x, centroid);
194+
195+
const glm::mat2x4 bez =
196+
CubicBezier2Linear(bez0[0], Bezier(glm::vec3(bez0[0]), bez0[1]),
197+
Bezier(glm::vec3(bez1[0]), bez1[1]), bez1[0], y);
189198
return BezierPoint(bez, y);
190199
}
191200

192-
void operator()(thrust::tuple<glm::vec3&, Barycentric> inOut) {
201+
void operator()(thrust::tuple<glm::vec3&, Barycentric> inOut) const {
193202
glm::vec3& pos = thrust::get<0>(inOut);
194203
const int tri = thrust::get<1>(inOut).tri;
195204
const glm::vec4 uvw = thrust::get<1>(inOut).uvw;
@@ -226,12 +235,17 @@ struct InterpTri {
226235
const int j = Next3(i);
227236
const int k = Prev3(i);
228237
const float x = uvw[k] / (1 - uvw[i]);
229-
const glm::vec3 p =
230-
Bezier2D({corners[i], corners[j], corners[k], corners[i]},
231-
{tangentR[i], tangentL[j], tangentR[k], tangentL[i]},
232-
{tangentL[i], tangentR[j], tangentL[k], tangentR[i]},
233-
1 - uvw[i], x, centroid, true);
234-
posH += Homogeneous(glm::vec4(p, uvw[i]));
238+
239+
const glm::mat2x4 bez =
240+
Bezier2Bezier({corners[j], corners[k]}, {tangentR[j], tangentL[k]},
241+
{tangentL[j], tangentR[k]}, x, centroid);
242+
243+
const glm::mat2x4 bez1 = CubicBezier2Linear(
244+
bez[0], Bezier(glm::vec3(bez[0]), bez[1]),
245+
Bezier(corners[i], glm::mix(tangentR[i], tangentL[i], x)),
246+
Homogeneous(corners[i]), uvw[i]);
247+
const glm::vec3 p = BezierPoint(bez1, uvw[i]);
248+
posH += Homogeneous(glm::vec4(p, uvw[j] * uvw[k]));
235249
}
236250
} else { // quad
237251
const glm::mat4 tangentsX = {
@@ -248,12 +262,12 @@ struct InterpTri {
248262
const float x = uvw[1] + uvw[2];
249263
const float y = uvw[2] + uvw[3];
250264
const glm::vec3 pX =
251-
Bezier2D(corners, tangentsX, tangentsY, x, y, centroid, false);
265+
Bezier2D(corners, tangentsX, tangentsY, x, y, centroid);
252266
const glm::vec3 pY =
253267
Bezier2D({corners[1], corners[2], corners[3], corners[0]},
254268
{tangentsY[1], tangentsY[2], tangentsY[3], tangentsY[0]},
255269
{tangentsX[1], tangentsX[2], tangentsX[3], tangentsX[0]}, y,
256-
1 - x, centroid, false);
270+
1 - x, centroid);
257271
posH += Homogeneous(glm::vec4(pX, x * (1 - x)));
258272
posH += Homogeneous(glm::vec4(pY, y * (1 - y)));
259273
}
@@ -885,7 +899,7 @@ void Manifold::Impl::CreateTangents(std::vector<Smoothness> sharpenedEdges) {
885899
&triIsFlatFace](int v) {
886900
auto it = vertTangents.find(v);
887901
if (it == vertTangents.end()) {
888-
fixedHalfedge[vertHalfedge[v]] == true;
902+
fixedHalfedge[vertHalfedge[v]] = true;
889903
return;
890904
}
891905
const std::vector<Pair>& vert = it->second;
@@ -918,7 +932,7 @@ void Manifold::Impl::CreateTangents(std::vector<Smoothness> sharpenedEdges) {
918932
}
919933
});
920934
} else { // Sharpen vertex uniformly
921-
fixedHalfedge[vertHalfedge[v]] == true;
935+
fixedHalfedge[vertHalfedge[v]] = true;
922936
float smoothness = 0;
923937
float denom = 0;
924938
for (const Pair& pair : vert) {

‎test/samples_test.cpp

+2-2
Original file line numberDiff line numberDiff line change
@@ -116,8 +116,8 @@ TEST(Samples, Scallop) {
116116
3, colorCurvature);
117117
CheckNormals(scallop);
118118
auto prop = scallop.GetProperties();
119-
EXPECT_NEAR(prop.volume, 40.9, 0.1);
120-
EXPECT_NEAR(prop.surfaceArea, 77.4, 0.1);
119+
EXPECT_NEAR(prop.volume, 39.9, 0.1);
120+
EXPECT_NEAR(prop.surfaceArea, 79.3, 0.1);
121121
CheckGL(scallop);
122122

123123
#ifdef MANIFOLD_EXPORT

‎test/smooth_test.cpp

+45-20
Original file line numberDiff line numberDiff line change
@@ -30,8 +30,16 @@ TEST(Smooth, Tetrahedron) {
3030
smooth = smooth.Refine(n);
3131
ExpectMeshes(smooth, {{2 * n * n + 2, 4 * n * n}});
3232
auto prop = smooth.GetProperties();
33-
EXPECT_NEAR(prop.volume, 18.7, 0.1);
34-
EXPECT_NEAR(prop.surfaceArea, 34.5, 0.1);
33+
EXPECT_NEAR(prop.volume, 17.0, 0.1);
34+
EXPECT_NEAR(prop.surfaceArea, 32.9, 0.1);
35+
36+
MeshGL out = smooth.CalculateCurvature(-1, 0).GetMeshGL();
37+
float maxMeanCurvature = 0;
38+
for (int i = 3; i < out.vertProperties.size(); i += 4) {
39+
maxMeanCurvature =
40+
glm::max(maxMeanCurvature, glm::abs(out.vertProperties[i]));
41+
}
42+
EXPECT_NEAR(maxMeanCurvature, 4.73, 0.01);
3543

3644
#ifdef MANIFOLD_EXPORT
3745
if (options.exportModels) ExportMesh("smoothTet.glb", smooth.GetMesh(), {});
@@ -72,8 +80,8 @@ TEST(Smooth, TruncatedCone) {
7280
Manifold cone = Manifold::Cylinder(5, 10, 5, 12).SmoothOut();
7381
Manifold smooth = cone.RefineToLength(0.5).CalculateNormals(0);
7482
auto prop = smooth.GetProperties();
75-
EXPECT_NEAR(prop.volume, 1063.61, 0.01);
76-
EXPECT_NEAR(prop.surfaceArea, 751.78, 0.01);
83+
EXPECT_NEAR(prop.volume, 1062.27, 0.01);
84+
EXPECT_NEAR(prop.surfaceArea, 751.46, 0.01);
7785
MeshGL out = smooth.GetMeshGL();
7886
CheckGL(out);
7987

@@ -95,8 +103,16 @@ TEST(Smooth, ToLength) {
95103
smooth = smooth.RefineToLength(0.1);
96104
ExpectMeshes(smooth, {{85250, 170496}});
97105
auto prop = smooth.GetProperties();
98-
EXPECT_NEAR(prop.volume, 4842, 1);
99-
EXPECT_NEAR(prop.surfaceArea, 1400, 1);
106+
EXPECT_NEAR(prop.volume, 4604, 1);
107+
EXPECT_NEAR(prop.surfaceArea, 1356, 1);
108+
109+
MeshGL out = smooth.CalculateCurvature(-1, 0).GetMeshGL();
110+
float maxMeanCurvature = 0;
111+
for (int i = 3; i < out.vertProperties.size(); i += 4) {
112+
maxMeanCurvature =
113+
glm::max(maxMeanCurvature, glm::abs(out.vertProperties[i]));
114+
}
115+
EXPECT_NEAR(maxMeanCurvature, 1.67, 0.01);
100116

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

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

157173
ExpectMeshes(interp, {{40002, 80000}});
158174
auto prop = interp.GetProperties();
159-
EXPECT_NEAR(prop.volume, 3.83, 0.01);
160-
EXPECT_NEAR(prop.surfaceArea, 11.95, 0.01);
175+
EXPECT_NEAR(prop.volume, 3.74, 0.01);
176+
EXPECT_NEAR(prop.surfaceArea, 11.78, 0.01);
161177

162178
#ifdef MANIFOLD_EXPORT
163179
if (options.exportModels) {
@@ -200,8 +216,8 @@ TEST(Smooth, Csaszar) {
200216
csaszar = csaszar.Refine(100);
201217
ExpectMeshes(csaszar, {{70000, 140000}});
202218
auto prop = csaszar.GetProperties();
203-
EXPECT_NEAR(prop.volume, 89873, 10);
204-
EXPECT_NEAR(prop.surfaceArea, 15301, 10);
219+
EXPECT_NEAR(prop.volume, 79890, 10);
220+
EXPECT_NEAR(prop.surfaceArea, 11950, 10);
205221

206222
#ifdef MANIFOLD_EXPORT
207223
if (options.exportModels) {
@@ -276,19 +292,28 @@ TEST(Smooth, Torus) {
276292
}
277293
}
278294

279-
Manifold smooth = Manifold(torusMesh).RefineToLength(0.1).CalculateNormals(0);
280-
Mesh out = smooth.GetMesh();
281-
for (glm::vec3 v : out.vertPos) {
295+
Manifold smooth = Manifold(torusMesh)
296+
.RefineToLength(0.1)
297+
.CalculateCurvature(-1, 0)
298+
.CalculateNormals(1);
299+
MeshGL out = smooth.GetMeshGL();
300+
float maxMeanCurvature = 0;
301+
for (int i = 0; i < out.vertProperties.size(); i += 7) {
302+
glm::vec3 v(out.vertProperties[i], out.vertProperties[i + 1],
303+
out.vertProperties[i + 2]);
282304
glm::vec3 p(v.x, v.y, 0);
283305
p = glm::normalize(p) * 2.0f;
284306
float r = glm::length(v - p);
285-
ASSERT_NEAR(r, 1, 0.005);
307+
ASSERT_NEAR(r, 1, 0.006);
308+
maxMeanCurvature =
309+
glm::max(maxMeanCurvature, glm::abs(out.vertProperties[i + 3]));
286310
}
311+
EXPECT_NEAR(maxMeanCurvature, 1.63, 0.01);
287312

288313
#ifdef MANIFOLD_EXPORT
289314
ExportOptions options2;
290315
options2.faceted = false;
291-
options2.mat.normalChannels = {3, 4, 5};
316+
options2.mat.normalChannels = {4, 5, 6};
292317
options2.mat.roughness = 0;
293318
if (options.exportModels) ExportMesh("smoothTorus.glb", out, options2);
294319
#endif
@@ -307,7 +332,7 @@ TEST(Smooth, SineSurface) {
307332
Manifold(surface).CalculateNormals(0, 50).SmoothByNormals(0).Refine(8);
308333
auto prop = smoothed.GetProperties();
309334
EXPECT_NEAR(prop.volume, 7.89, 0.01);
310-
EXPECT_NEAR(prop.surfaceArea, 30.61, 0.01);
335+
EXPECT_NEAR(prop.surfaceArea, 30.60, 0.01);
311336
EXPECT_EQ(smoothed.Genus(), 0);
312337

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

319344
Manifold smoothed2 = Manifold(surface).SmoothOut(180, 1).Refine(8);
320345
auto prop2 = smoothed2.GetProperties();
321-
EXPECT_NEAR(prop2.volume, 9.06, 0.01);
322-
EXPECT_NEAR(prop2.surfaceArea, 33.81, 0.01);
346+
EXPECT_NEAR(prop2.volume, 9.02, 0.01);
347+
EXPECT_NEAR(prop2.surfaceArea, 33.56, 0.01);
323348
EXPECT_EQ(smoothed2.Genus(), 0);
324349

325350
Manifold smoothed3 = Manifold(surface).SmoothOut(50, 0.5).Refine(8);
326351
auto prop3 = smoothed3.GetProperties();
327352
EXPECT_NEAR(prop3.volume, 8.46, 0.01);
328-
EXPECT_NEAR(prop3.surfaceArea, 31.77, 0.01);
353+
EXPECT_NEAR(prop3.surfaceArea, 31.66, 0.01);
329354
EXPECT_EQ(smoothed3.Genus(), 0);
330355

331356
#ifdef MANIFOLD_EXPORT

0 commit comments

Comments
 (0)