Skip to content

Commit a122bde

Browse files
elalishelalish
authored
Quaternion smoothing (#801)
* revert smoothness change * updated to quaternions - not right yet * fixed triangles, mostly * fixed numerical stability problem * fixed more nans * fixed several bugs * refactor bezier2D * improved circularity * added test * updated tests * fixed mirror test * updated manual smooth model * updated JS test * cleanup
1 parent c9810fd commit a122bde

File tree

9 files changed

+267
-112
lines changed

9 files changed

+267
-112
lines changed

‎.vscode/launch.json

+1-1
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@
1212
"program": "${workspaceFolder}/build/test/manifold_test",
1313
"args": [
1414
"--gtest_break_on_failure",
15-
"--gtest_filter=Manifold.RefineQuads"
15+
"--gtest_filter=SDF.SineSurface"
1616
],
1717
"stopAtEntry": false,
1818
"cwd": "${workspaceFolder}/build/test",

‎bindings/wasm/examples/public/examples.js

+1-1
Original file line numberDiff line numberDiff line change
@@ -221,7 +221,7 @@ export const examples = {
221221
const radius = 30;
222222
const offset = 20;
223223
const wiggles = 12;
224-
const sharpness = 0.6;
224+
const sharpness = 0.8;
225225
const n = 50;
226226

227227
const positions = [];

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

+2-2
Original file line numberDiff line numberDiff line change
@@ -116,8 +116,8 @@ suite('Examples', () => {
116116
test('Scallop', async () => {
117117
const result = await runExample('Scallop');
118118
expect(result.genus).to.equal(0, 'Genus');
119-
expect(result.volume).to.be.closeTo(41500, 100, 'Volume');
120-
expect(result.surfaceArea).to.be.closeTo(7880, 10, 'Surface Area');
119+
expect(result.volume).to.be.closeTo(41400, 100, 'Volume');
120+
expect(result.surfaceArea).to.be.closeTo(7770, 10, 'Surface Area');
121121
});
122122

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

‎samples/src/scallop.cpp

+1-1
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ Manifold Scallop() {
2525
constexpr float radius = 3;
2626
constexpr float offset = 2;
2727
constexpr int wiggles = 12;
28-
constexpr float sharpness = 0.6;
28+
constexpr float sharpness = 0.8;
2929

3030
Mesh scallop;
3131
std::vector<Smoothness> sharpenedEdges;

‎src/manifold/src/impl.h

+1
Original file line numberDiff line numberDiff line change
@@ -180,6 +180,7 @@ struct Manifold::Impl {
180180
Vec<int> VertFlatFace(const Vec<bool>&) const;
181181
std::vector<Smoothness> SharpenEdges(float minSharpAngle,
182182
float minSmoothness) const;
183+
void SharpenTangent(int halfedge, float smoothness);
183184
void SetNormals(int normalIdx, float minSharpAngle);
184185
void LinearizeFlatTangents();
185186
void CreateTangents(int normalIdx);

‎src/manifold/src/smoothing.cpp

+122-64
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,8 @@
1212
// See the License for the specific language governing permissions and
1313
// limitations under the License.
1414

15+
#include <glm/gtx/quaternion.hpp>
16+
1517
#include "impl.h"
1618
#include "par.h"
1719

@@ -81,7 +83,13 @@ struct InterpTri {
8183

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

84-
glm::vec3 HNormalize(glm::vec4 v) const { return glm::vec3(v) / v.w; }
86+
glm::vec3 HNormalize(glm::vec4 v) const {
87+
return v.w == 0 ? v : (glm::vec3(v) / v.w);
88+
}
89+
90+
glm::vec4 Scale(glm::vec4 v, float scale) const {
91+
return glm::vec4(scale * glm::vec3(v), v.w);
92+
}
8593

8694
glm::vec4 Bezier(glm::vec3 point, glm::vec4 tangent) const {
8795
return Homogeneous(glm::vec4(point, 0) + tangent);
@@ -101,43 +109,75 @@ struct InterpTri {
101109
}
102110

103111
glm::vec3 BezierTangent(glm::mat2x4 points) const {
104-
return glm::normalize(HNormalize(points[1]) - HNormalize(points[0]));
112+
return SafeNormalize(HNormalize(points[1]) - HNormalize(points[0]));
113+
}
114+
115+
glm::vec3 RotateFromTo(glm::vec3 v, glm::quat start, glm::quat end) const {
116+
return end * glm::conjugate(start) * v;
105117
}
106118

107119
glm::mat2x4 Bezier2Bezier(const glm::mat2x3& corners,
108120
const glm::mat2x4& tangentsX,
109-
const glm::mat2x4& tangentsY, float x) const {
121+
const glm::mat2x4& tangentsY, float x,
122+
const glm::bvec2& pointedEnds,
123+
const glm::vec3& anchor) const {
110124
const glm::mat2x4 bez = CubicBezier2Linear(
111125
Homogeneous(corners[0]), Bezier(corners[0], tangentsX[0]),
112126
Bezier(corners[1], tangentsX[1]), Homogeneous(corners[1]), x);
113127
const glm::vec3 end = BezierPoint(bez, x);
114128
const glm::vec3 tangent = BezierTangent(bez);
115129

130+
const glm::mat2x3 nTangentsX(SafeNormalize(glm::vec3(tangentsX[0])),
131+
-SafeNormalize(glm::vec3(tangentsX[1])));
116132
const glm::mat2x3 biTangents = {
117-
SafeNormalize(OrthogonalTo(glm::vec3(tangentsY[0]),
118-
SafeNormalize(glm::vec3(tangentsX[0])))),
119-
SafeNormalize(OrthogonalTo(glm::vec3(tangentsY[1]),
120-
-SafeNormalize(glm::vec3(tangentsX[1]))))};
121-
const glm::vec3 normal = SafeNormalize(
122-
glm::cross(glm::mix(biTangents[0], biTangents[1], x), tangent));
123-
const glm::vec3 delta = OrthogonalTo(
124-
glm::mix(glm::vec3(tangentsY[0]), glm::vec3(tangentsY[1]), x), normal);
133+
SafeNormalize(OrthogonalTo(
134+
glm::vec3(tangentsY[0]) + kTolerance * (anchor - corners[0]),
135+
nTangentsX[0])),
136+
SafeNormalize(OrthogonalTo(
137+
glm::vec3(tangentsY[1]) + kTolerance * (anchor - corners[1]),
138+
nTangentsX[1]))};
139+
140+
const glm::quat q0 =
141+
glm::quat_cast(glm::mat3(nTangentsX[0], biTangents[0],
142+
glm::cross(nTangentsX[0], biTangents[0])));
143+
const glm::quat q1 =
144+
glm::quat_cast(glm::mat3(nTangentsX[1], biTangents[1],
145+
glm::cross(nTangentsX[1], biTangents[1])));
146+
const glm::quat qTmp = glm::slerp(q0, q1, x);
147+
const glm::quat q =
148+
glm::rotation(qTmp * glm::vec3(1, 0, 0), tangent) * qTmp;
149+
150+
const glm::vec3 end0 = pointedEnds[0]
151+
? glm::vec3(0)
152+
: RotateFromTo(glm::vec3(tangentsY[0]), q0, q);
153+
const glm::vec3 end1 = pointedEnds[1]
154+
? glm::vec3(0)
155+
: RotateFromTo(glm::vec3(tangentsY[1]), q1, q);
156+
const glm::vec3 delta = glm::mix(end0, end1, x);
125157
const float deltaW = glm::mix(tangentsY[0].w, tangentsY[1].w, x);
126158

127-
return {Homogeneous(end), Homogeneous(glm::vec4(end + delta, deltaW))};
159+
return {Homogeneous(end), glm::vec4(delta, deltaW)};
128160
}
129161

130162
glm::vec3 Bezier2D(const glm::mat4x3& corners, const glm::mat4& tangentsX,
131-
const glm::mat4& tangentsY, float x, float y) const {
132-
glm::mat2x4 bez0 =
133-
Bezier2Bezier({corners[0], corners[1]}, {tangentsX[0], tangentsX[1]},
134-
{tangentsY[0], tangentsY[1]}, x);
135-
glm::mat2x4 bez1 =
136-
Bezier2Bezier({corners[2], corners[3]}, {tangentsX[2], tangentsX[3]},
137-
{tangentsY[2], tangentsY[3]}, 1 - x);
138-
139-
const glm::mat2x4 bez =
140-
CubicBezier2Linear(bez0[0], bez0[1], bez1[1], bez1[0], y);
163+
const glm::mat4& tangentsY, float x, float y,
164+
const glm::vec3& centroid, bool isTriangle) const {
165+
glm::mat2x4 bez0 = Bezier2Bezier(
166+
{corners[0], corners[1]}, {tangentsX[0], tangentsX[1]},
167+
{tangentsY[0], tangentsY[1]}, x, {isTriangle, false}, centroid);
168+
glm::mat2x4 bez1 = Bezier2Bezier(
169+
{corners[2], corners[3]}, {tangentsX[2], tangentsX[3]},
170+
{tangentsY[2], tangentsY[3]}, 1 - x, {false, isTriangle}, centroid);
171+
172+
const float flatLen =
173+
isTriangle ? x * glm::length(corners[1] - corners[2])
174+
: glm::length(glm::mix(corners[0], corners[1], x) -
175+
glm::mix(corners[3], corners[2], x));
176+
const float scale = glm::length(glm::vec3(bez0[0] - bez1[0])) / flatLen;
177+
178+
const glm::mat2x4 bez = CubicBezier2Linear(
179+
bez0[0], Bezier(glm::vec3(bez0[0]), Scale(bez0[1], scale)),
180+
Bezier(glm::vec3(bez1[0]), Scale(bez1[1], scale)), bez1[0], y);
141181
return BezierPoint(bez, y);
142182
}
143183

@@ -172,24 +212,18 @@ struct InterpTri {
172212
impl->halfedgeTangent_[impl->halfedge_[halfedges[2]].pairedHalfedge],
173213
impl->halfedgeTangent_[impl->halfedge_[halfedges[0]].pairedHalfedge],
174214
impl->halfedgeTangent_[impl->halfedge_[halfedges[1]].pairedHalfedge]};
215+
const glm::vec3 centroid = glm::mat3(corners) * glm::vec3(1.0f / 3);
175216

176217
for (const int i : {0, 1, 2}) {
177-
const int j = (i + 1) % 3;
178-
const int k = (i + 2) % 3;
218+
const int j = Next3(i);
219+
const int k = Prev3(i);
179220
const float x = uvw[k] / (1 - uvw[i]);
180-
181-
const glm::mat2x4 bez =
182-
Bezier2Bezier({corners[j], corners[k]}, {tangentR[j], tangentL[k]},
183-
{tangentL[j], tangentR[k]}, x);
184-
185-
const glm::mat2x4 bez1 = CubicBezier2Linear(
186-
bez[0], bez[1],
187-
// Homogeneous(end), Homogeneous(glm::vec4(end + delta, deltaW)),
188-
Bezier(corners[i], glm::mix(tangentR[i], tangentL[i], x)),
189-
Homogeneous(corners[i]), uvw[i]);
190-
const glm::vec3 p = BezierPoint(bez1, uvw[i]);
191-
float w = uvw[j] * uvw[j] * uvw[k] * uvw[k];
192-
posH += Homogeneous(glm::vec4(p, w));
221+
const glm::vec3 p =
222+
Bezier2D({corners[i], corners[j], corners[k], corners[i]},
223+
{tangentR[i], tangentL[j], tangentR[k], tangentL[i]},
224+
{tangentL[i], tangentR[j], tangentL[k], tangentR[i]},
225+
1 - uvw[i], x, centroid, true);
226+
posH += Homogeneous(glm::vec4(p, uvw[i]));
193227
}
194228
} else { // quad
195229
const glm::mat4 tangentsX = {
@@ -202,13 +236,16 @@ struct InterpTri {
202236
impl->halfedgeTangent_[halfedges[1]],
203237
impl->halfedgeTangent_[impl->halfedge_[halfedges[1]].pairedHalfedge],
204238
impl->halfedgeTangent_[halfedges[3]]};
239+
const glm::vec3 centroid = corners * glm::vec4(0.25);
205240
const float x = uvw[1] + uvw[2];
206241
const float y = uvw[2] + uvw[3];
207-
const glm::vec3 pX = Bezier2D(corners, tangentsX, tangentsY, x, y);
208-
const glm::vec3 pY = Bezier2D(
209-
{corners[1], corners[2], corners[3], corners[0]},
210-
{tangentsY[1], tangentsY[2], tangentsY[3], tangentsY[0]},
211-
{tangentsX[1], tangentsX[2], tangentsX[3], tangentsX[0]}, y, 1 - x);
242+
const glm::vec3 pX =
243+
Bezier2D(corners, tangentsX, tangentsY, x, y, centroid, false);
244+
const glm::vec3 pY =
245+
Bezier2D({corners[1], corners[2], corners[3], corners[0]},
246+
{tangentsY[1], tangentsY[2], tangentsY[3], tangentsY[0]},
247+
{tangentsX[1], tangentsX[2], tangentsX[3], tangentsX[0]}, y,
248+
1 - x, centroid, false);
212249
posH += Homogeneous(glm::vec4(pX, x * (1 - x)));
213250
posH += Homogeneous(glm::vec4(pY, y * (1 - y)));
214251
}
@@ -359,6 +396,17 @@ std::vector<Smoothness> Manifold::Impl::SharpenEdges(
359396
return sharpenedEdges;
360397
}
361398

399+
/**
400+
* Sharpen tangents that intersect an edge to sharpen that edge. The weight is
401+
* unchanged, as this has a squared effect on radius of curvature, except
402+
* in the case of zero radius, which is marked with weight = 0.
403+
*/
404+
void Manifold::Impl::SharpenTangent(int halfedge, float smoothness) {
405+
halfedgeTangent_[halfedge] =
406+
glm::vec4(smoothness * glm::vec3(halfedgeTangent_[halfedge]),
407+
smoothness == 0 ? 0 : halfedgeTangent_[halfedge].w);
408+
}
409+
362410
/**
363411
* Instead of calculating the internal shared normals like CalculateNormals
364412
* does, this method fills in vertex properties, unshared across edges that
@@ -552,21 +600,31 @@ void Manifold::Impl::SetNormals(int normalIdx, float minSharpAngle) {
552600
*/
553601
void Manifold::Impl::LinearizeFlatTangents() {
554602
const int n = halfedgeTangent_.size();
555-
for_each_n(autoPolicy(n), zip(halfedgeTangent_.begin(), countAt(0)), n,
556-
[this](thrust::tuple<glm::vec4&, int> inOut) {
557-
glm::vec4& tangent = thrust::get<0>(inOut);
558-
const int halfedge = thrust::get<1>(inOut);
559-
if (tangent.w != 0) {
560-
return;
561-
}
603+
for_each_n(
604+
autoPolicy(n), zip(halfedgeTangent_.begin(), countAt(0)), n,
605+
[this](thrust::tuple<glm::vec4&, int> inOut) {
606+
glm::vec4& tangent = thrust::get<0>(inOut);
607+
const int halfedge = thrust::get<1>(inOut);
608+
glm::vec4& otherTangent =
609+
halfedgeTangent_[halfedge_[halfedge].pairedHalfedge];
610+
611+
const glm::bvec2 flat(tangent.w == 0, otherTangent.w == 0);
612+
if (!halfedge_[halfedge].IsForward() || (!flat[0] && !flat[1])) {
613+
return;
614+
}
562615

563-
const glm::vec4 otherTangent =
564-
halfedgeTangent_[halfedge_[halfedge].pairedHalfedge];
565-
glm::vec3 edge = vertPos_[halfedge_[halfedge].endVert] +
566-
glm::vec3(otherTangent) -
567-
vertPos_[halfedge_[halfedge].startVert];
568-
tangent = glm::vec4(edge / 3.0f, 1);
569-
});
616+
const glm::vec3 edgeVec = vertPos_[halfedge_[halfedge].endVert] -
617+
vertPos_[halfedge_[halfedge].startVert];
618+
619+
if (flat[0] && flat[1]) {
620+
tangent = glm::vec4(edgeVec / 3.0f, 1);
621+
otherTangent = glm::vec4(-edgeVec / 3.0f, 1);
622+
} else if (flat[0]) {
623+
tangent = glm::vec4((edgeVec + glm::vec3(otherTangent)) / 2.0f, 1);
624+
} else {
625+
otherTangent = glm::vec4((-edgeVec + glm::vec3(tangent)) / 2.0f, 1);
626+
}
627+
});
570628
}
571629

572630
/**
@@ -605,10 +663,10 @@ void Manifold::Impl::CreateTangents(int normalIdx) {
605663
const glm::vec3& nextNormal) {
606664
const glm::vec3 diff = nextNormal - normal;
607665
if (glm::dot(diff, diff) > kTolerance * kTolerance) {
608-
if (idx > 1) {
609-
sharpHalfedge[0] = -1;
610-
} else {
666+
if (idx < 2) {
611667
sharpHalfedge[idx++] = halfedge;
668+
} else {
669+
sharpHalfedge[0] = -1; // marks more than 2 sharp edges
612670
}
613671
}
614672
lastNormal = normal;
@@ -638,13 +696,13 @@ void Manifold::Impl::CreateTangents(int normalIdx) {
638696
ForVert(first, [this, first, second](int current) {
639697
if (current != first && current != second &&
640698
!IsMarkedInsideQuad(current)) {
641-
halfedgeTangent_[current] = glm::vec4(0);
699+
SharpenTangent(current, 0);
642700
}
643701
});
644702
} else { // Sharpen vertex uniformly
645-
ForVert(first, [this](int current) {
703+
ForVert(second, [this](int current) {
646704
if (!IsMarkedInsideQuad(current)) {
647-
halfedgeTangent_[current] = glm::vec4(0);
705+
SharpenTangent(current, 0);
648706
}
649707
});
650708
}
@@ -745,7 +803,7 @@ void Manifold::Impl::CreateTangents(std::vector<Smoothness> sharpenedEdges) {
745803
smoothness =
746804
(vert[1].second.smoothness + vert[0].first.smoothness) / 2;
747805
} else if (current != first && !IsMarkedInsideQuad(current)) {
748-
halfedgeTangent_[current] = smoothness * halfedgeTangent_[current];
806+
SharpenTangent(current, smoothness);
749807
}
750808
});
751809
} else { // Sharpen vertex uniformly
@@ -758,7 +816,7 @@ void Manifold::Impl::CreateTangents(std::vector<Smoothness> sharpenedEdges) {
758816

759817
ForVert(vert[0].first.halfedge, [this, smoothness](int current) {
760818
if (!IsMarkedInsideQuad(current)) {
761-
halfedgeTangent_[current] = smoothness * halfedgeTangent_[current];
819+
SharpenTangent(current, smoothness);
762820
}
763821
});
764822
}

‎src/manifold/src/subdivision.cpp

+6-4
Original file line numberDiff line numberDiff line change
@@ -398,8 +398,8 @@ class Partition {
398398
namespace manifold {
399399

400400
/**
401-
* Returns the tri index of the other side of this quad if this tri is part of a
402-
* quad, or -1 otherwise.
401+
* Returns the tri side index (0-2) connected to the other side of this quad if
402+
* this tri is part of a quad, or -1 otherwise.
403403
*/
404404
int Manifold::Impl::GetNeighbor(int tri) const {
405405
int neighbor = -1;
@@ -450,7 +450,7 @@ glm::ivec4 Manifold::Impl::GetHalfedges(int tri) const {
450450
* Returns the BaryIndices, which gives the tri and indices (0-3), such that
451451
* GetHalfedges(val.tri)[val.start4] points back to this halfedge, and val.end4
452452
* will point to the next one. This function handles this for both triangles and
453-
* quads.
453+
* quads. Returns {-1, -1, -1} if the edge is the interior of a quad.
454454
*/
455455
Manifold::Impl::BaryIndices Manifold::Impl::GetIndices(int halfedge) const {
456456
int tri = halfedge / 3;
@@ -468,6 +468,8 @@ Manifold::Impl::BaryIndices Manifold::Impl::GetIndices(int halfedge) const {
468468
tri = pair / 3;
469469
const int j = pair % 3;
470470
idx = Next3(neighbor) == idx ? j : (j + 1) % 4;
471+
} else if (idx > neighbor) {
472+
++idx;
471473
}
472474
return {tri, idx, (idx + 1) % 4};
473475
}
@@ -484,7 +486,7 @@ void Manifold::Impl::FillRetainedVerts(Vec<Barycentric>& vertBary) const {
484486
for (int tri = 0; tri < numTri; ++tri) {
485487
for (const int i : {0, 1, 2}) {
486488
const BaryIndices indices = GetIndices(3 * tri + i);
487-
if (indices.start4 < 0) continue;
489+
if (indices.start4 < 0) continue; // skip quad interiors
488490
glm::vec4 uvw(0);
489491
uvw[indices.start4] = 1;
490492
vertBary[halfedge_[3 * tri + i].startVert] = {indices.tri, uvw};

0 commit comments

Comments
 (0)