1 #include "subdivide.hpp"
2
3 #include <fstream>
4 #include <queue>
5
6 #include "dedge.hpp"
7 #include "disajoint-tree.hpp"
8 #include "field-math.hpp"
9 #include "parametrizer.hpp"
10
11 namespace qflow {
12
subdivide(MatrixXi & F,MatrixXd & V,VectorXd & rho,VectorXi & V2E,VectorXi & E2E,VectorXi & boundary,VectorXi & nonmanifold,double maxLength)13 void subdivide(MatrixXi &F, MatrixXd &V, VectorXd& rho, VectorXi &V2E, VectorXi &E2E, VectorXi &boundary,
14 VectorXi &nonmanifold, double maxLength) {
15 typedef std::pair<double, int> Edge;
16
17 std::priority_queue<Edge> queue;
18
19 maxLength *= maxLength;
20
21 for (int i = 0; i < E2E.size(); ++i) {
22 int v0 = F(i % 3, i / 3), v1 = F((i + 1) % 3, i / 3);
23 if (nonmanifold[v0] || nonmanifold[v1]) continue;
24 double length = (V.col(v0) - V.col(v1)).squaredNorm();
25 if (length > maxLength || length > std::max(maxLength * 0.75, std::min(rho[v0], rho[v1]) * 1.0)) {
26 int other = E2E[i];
27 if (other == -1 || other > i) queue.push(Edge(length, i));
28 }
29 }
30
31 int nV = V.cols(), nF = F.cols(), nSplit = 0;
32 /*
33 / v0 \
34 v1p 1 | 0 v0p
35 \ v1 /
36
37 / v0 \
38 / 1 | 0 \
39 v1p - vn - v0p
40 \ 2 | 3 /
41 \ v1 /
42
43 f0: vn, v0p, v0
44 f1: vn, v0, v1p
45 f2: vn, v1p, v1
46 f3: vn, v1, v0p
47 */
48 int counter = 0;
49 while (!queue.empty()) {
50 counter += 1;
51 Edge edge = queue.top();
52 queue.pop();
53 int e0 = edge.second, e1 = E2E[e0];
54 bool is_boundary = e1 == -1;
55 int f0 = e0 / 3, f1 = is_boundary ? -1 : (e1 / 3);
56 int v0 = F(e0 % 3, f0), v0p = F((e0 + 2) % 3, f0), v1 = F((e0 + 1) % 3, f0);
57 if ((V.col(v0) - V.col(v1)).squaredNorm() != edge.first) {
58 continue;
59 }
60 int v1p = is_boundary ? -1 : F((e1 + 2) % 3, f1);
61 int vn = nV++;
62 nSplit++;
63 /* Update V */
64 if (nV > V.cols()) {
65 V.conservativeResize(V.rows(), V.cols() * 2);
66 rho.conservativeResize(V.cols() * 2);
67 V2E.conservativeResize(V.cols());
68 boundary.conservativeResize(V.cols());
69 nonmanifold.conservativeResize(V.cols());
70 }
71
72 /* Update V */
73 V.col(vn) = (V.col(v0) + V.col(v1)) * 0.5f;
74 rho[vn] = 0.5f * (rho[v0], rho[v1]);
75 nonmanifold[vn] = false;
76 boundary[vn] = is_boundary;
77
78 /* Update F and E2E */
79 int f2 = is_boundary ? -1 : (nF++);
80 int f3 = nF++;
81 if (nF > F.cols()) {
82 F.conservativeResize(F.rows(), std::max(nF, (int)F.cols() * 2));
83 E2E.conservativeResize(F.cols() * 3);
84 }
85
86 /* Update F */
87 F.col(f0) << vn, v0p, v0;
88 if (!is_boundary) {
89 F.col(f1) << vn, v0, v1p;
90 F.col(f2) << vn, v1p, v1;
91 }
92 F.col(f3) << vn, v1, v0p;
93
94 /* Update E2E */
95 const int e0p = E2E[dedge_prev_3(e0)], e0n = E2E[dedge_next_3(e0)];
96
97 #define sE2E(a, b) \
98 E2E[a] = b; \
99 if (b != -1) E2E[b] = a;
100 sE2E(3 * f0 + 0, 3 * f3 + 2);
101 sE2E(3 * f0 + 1, e0p);
102 sE2E(3 * f3 + 1, e0n);
103 if (is_boundary) {
104 sE2E(3 * f0 + 2, -1);
105 sE2E(3 * f3 + 0, -1);
106 } else {
107 const int e1p = E2E[dedge_prev_3(e1)], e1n = E2E[dedge_next_3(e1)];
108 sE2E(3 * f0 + 2, 3 * f1 + 0);
109 sE2E(3 * f1 + 1, e1n);
110 sE2E(3 * f1 + 2, 3 * f2 + 0);
111 sE2E(3 * f2 + 1, e1p);
112 sE2E(3 * f2 + 2, 3 * f3 + 0);
113 }
114 #undef sE2E
115
116 /* Update V2E */
117 V2E[v0] = 3 * f0 + 2;
118 V2E[vn] = 3 * f0 + 0;
119 V2E[v1] = 3 * f3 + 1;
120 V2E[v0p] = 3 * f0 + 1;
121 if (!is_boundary) V2E[v1p] = 3 * f1 + 2;
122
123 auto schedule = [&](int f) {
124 for (int i = 0; i < 3; ++i) {
125 double length = (V.col(F(i, f)) - V.col(F((i + 1) % 3, f))).squaredNorm();
126 if (length > maxLength
127 || length > std::max(maxLength * 0.75, std::min(rho[F(i, f)], rho[F((i + 1) % 3, f)]) * 1.0))
128 queue.push(Edge(length, f * 3 + i));
129 }
130 };
131
132 schedule(f0);
133 if (!is_boundary) {
134 schedule(f2);
135 schedule(f1);
136 };
137 schedule(f3);
138 }
139 F.conservativeResize(F.rows(), nF);
140 V.conservativeResize(V.rows(), nV);
141 rho.conservativeResize(nV);
142 V2E.conservativeResize(nV);
143 boundary.conservativeResize(nV);
144 nonmanifold.conservativeResize(nV);
145 E2E.conservativeResize(nF * 3);
146 }
147
subdivide_edgeDiff(MatrixXi & F,MatrixXd & V,MatrixXd & N,MatrixXd & Q,MatrixXd & O,MatrixXd * S,VectorXi & V2E,VectorXi & E2E,VectorXi & boundary,VectorXi & nonmanifold,std::vector<Vector2i> & edge_diff,std::vector<DEdge> & edge_values,std::vector<Vector3i> & face_edgeOrients,std::vector<Vector3i> & face_edgeIds,std::vector<int> & sharp_edges,std::map<int,int> & singularities,int max_len)148 void subdivide_edgeDiff(MatrixXi &F, MatrixXd &V, MatrixXd &N, MatrixXd &Q, MatrixXd &O, MatrixXd* S,
149 VectorXi &V2E, VectorXi &E2E, VectorXi &boundary, VectorXi &nonmanifold,
150 std::vector<Vector2i> &edge_diff, std::vector<DEdge> &edge_values,
151 std::vector<Vector3i> &face_edgeOrients, std::vector<Vector3i> &face_edgeIds,
152 std::vector<int>& sharp_edges, std::map<int, int> &singularities, int max_len) {
153 struct EdgeLink {
154 int id;
155 double length;
156 Vector2i diff;
157 int maxlen() const { return std::max(abs(diff[0]), abs(diff[1])); }
158 bool operator<(const EdgeLink &link) const { return maxlen() < link.maxlen(); }
159 };
160
161 struct FaceOrient {
162 int orient;
163 Vector3i d;
164 Vector3d q;
165 Vector3d n;
166 };
167
168 std::vector<FaceOrient> face_spaces(F.cols());
169 std::priority_queue<EdgeLink> queue;
170 std::vector<Vector2i> diffs(E2E.size());
171 for (int i = 0; i < F.cols(); ++i) {
172 for (int j = 0; j < 3; ++j) {
173 int eid = i * 3 + j;
174 diffs[eid] = rshift90(edge_diff[face_edgeIds[i][j]], face_edgeOrients[i][j]);
175 }
176 }
177 for (int i = 0; i < F.cols(); ++i) {
178 FaceOrient orient{};
179 orient.q = Q.col(F(0, i));
180 orient.n = N.col(F(0, i));
181 int orient_diff[3];
182 for (int j = 0; j < 3; ++j) {
183 int final_orient = face_edgeOrients[i][j];
184 int eid = face_edgeIds[i][j];
185 auto value = compat_orientation_extrinsic_index_4(
186 Q.col(edge_values[eid].x), N.col(edge_values[eid].x), orient.q, orient.n);
187 int target_orient = (value.second - value.first + 4) % 4;
188 if (F(j, i) == edge_values[eid].y) target_orient = (target_orient + 2) % 4;
189 orient_diff[j] = (final_orient - target_orient + 4) % 4;
190 }
191 if (orient_diff[0] == orient_diff[1])
192 orient.orient = orient_diff[0];
193 else if (orient_diff[0] == orient_diff[2])
194 orient.orient = orient_diff[2];
195 else if (orient_diff[1] == orient_diff[2])
196 orient.orient = orient_diff[1];
197 orient.d = Vector3i((orient_diff[0] - orient.orient + 4) % 4,
198 (orient_diff[1] - orient.orient + 4) % 4,
199 (orient_diff[2] - orient.orient + 4) % 4);
200 face_spaces[i] = (orient);
201 }
202 for (int i = 0; i < E2E.size(); ++i) {
203 int v0 = F(i % 3, i / 3), v1 = F((i + 1) % 3, i / 3);
204 if (nonmanifold[v0] || nonmanifold[v1]) continue;
205 double length = (V.col(v0) - V.col(v1)).squaredNorm();
206 Vector2i diff = diffs[i];
207 if (abs(diff[0]) > max_len || abs(diff[1]) > max_len) {
208 int other = E2E[i];
209 if (other == -1 || other > i) {
210 EdgeLink e;
211 e.id = i;
212 e.length = length;
213 e.diff = diff;
214 queue.push(e);
215 }
216 }
217 }
218 auto AnalyzeOrient = [&](int f0, const Vector3i &d) {
219 for (int j = 0; j < 3; ++j) {
220 int orient = face_spaces[f0].orient + d[j];
221 int v = std::min(F(j, f0), F((j + 1) % 3, f0));
222 auto value = compat_orientation_extrinsic_index_4(
223 Q.col(v), N.col(v), face_spaces[f0].q, face_spaces[f0].n);
224 if (F(j, f0) != v) orient += 2;
225 face_edgeOrients[f0][j] = (orient + value.second - value.first + 4) % 4;
226 }
227 face_spaces[f0].d = d;
228 for (int j = 0; j < 3; ++j) {
229 int eid = face_edgeIds[f0][j];
230 int orient = face_edgeOrients[f0][j];
231 auto diff = rshift90(diffs[f0 * 3 + j], (4 - orient) % 4);
232 edge_diff[eid] = diff;
233 }
234 };
235 auto FixOrient = [&](int f0) {
236 for (int j = 0; j < 3; ++j) {
237 auto diff = edge_diff[face_edgeIds[f0][j]];
238 if (rshift90(diff, face_edgeOrients[f0][j]) != diffs[f0 * 3 + j]) {
239 int orient = 0;
240 while (orient < 4 && rshift90(diff, orient) != diffs[f0 * 3 + j]) orient += 1;
241 face_spaces[f0].d[j] =
242 (face_spaces[f0].d[j] + orient - face_edgeOrients[f0][j]) % 4;
243 face_edgeOrients[f0][j] = orient;
244 }
245 }
246 };
247 /*
248 auto Length = [&](int f0) {
249 int l = 0;
250 for (int j = 0; j < 3; ++j) {
251 for (int k = 0; k < 2; ++k) {
252 l += abs(diffs[f0*3+j][k]);
253 }
254 printf("<%d %d> ", diffs[f0*3+j][0], diffs[f0*3+j][1]);
255 }
256 printf("\n");
257 return l;
258 };
259 */
260 int nV = V.cols(), nF = F.cols(), nSplit = 0;
261 /*
262 / v0 \
263 v1p 1 | 0 v0p
264 \ v1 /
265
266 / v0 \
267 / 1 | 0 \
268 v1p - vn - v0p
269 \ 2 | 3 /
270 \ v1 /
271
272 f0: vn, v0p, v0
273 f1: vn, v0, v1p
274 f2: vn, v1p, v1
275 f3: vn, v1, v0p
276 */
277 int counter = 0;
278 while (!queue.empty()) {
279 counter += 1;
280 EdgeLink edge = queue.top();
281 queue.pop();
282
283 int e0 = edge.id, e1 = E2E[e0];
284 bool is_boundary = e1 == -1;
285 int f0 = e0 / 3, f1 = is_boundary ? -1 : (e1 / 3);
286 int v0 = F(e0 % 3, f0), v0p = F((e0 + 2) % 3, f0), v1 = F((e0 + 1) % 3, f0);
287 if ((V.col(v0) - V.col(v1)).squaredNorm() != edge.length) {
288 continue;
289 }
290 if (abs(diffs[e0][0]) < 2 && abs(diffs[e0][1]) < 2) continue;
291 if (f1 != -1) {
292 face_edgeOrients.push_back(Vector3i());
293 sharp_edges.push_back(0);
294 sharp_edges.push_back(0);
295 sharp_edges.push_back(0);
296 face_edgeIds.push_back(Vector3i());
297 }
298 int v1p = is_boundary ? -1 : F((e1 + 2) % 3, f1);
299 int vn = nV++;
300 nSplit++;
301 if (nV > V.cols()) {
302 V.conservativeResize(V.rows(), V.cols() * 2);
303 N.conservativeResize(N.rows(), N.cols() * 2);
304 Q.conservativeResize(Q.rows(), Q.cols() * 2);
305 O.conservativeResize(O.rows(), O.cols() * 2);
306 if (S)
307 S->conservativeResize(S->rows(), S->cols() * 2);
308 V2E.conservativeResize(V.cols());
309 boundary.conservativeResize(V.cols());
310 nonmanifold.conservativeResize(V.cols());
311 }
312
313 V.col(vn) = (V.col(v0) + V.col(v1)) * 0.5;
314 N.col(vn) = N.col(v0);
315 Q.col(vn) = Q.col(v0);
316 O.col(vn) = (O.col(v0) + O.col(v1)) * 0.5;
317 if (S)
318 S->col(vn) = S->col(v0);
319
320 nonmanifold[vn] = false;
321 boundary[vn] = is_boundary;
322
323 int eid = face_edgeIds[f0][e0 % 3];
324 int sharp_eid = sharp_edges[e0];
325 int eid01 = face_edgeIds[f0][(e0 + 1) % 3];
326 int sharp_eid01 = sharp_edges[f0 * 3 + (e0 + 1) % 3];
327 int eid02 = face_edgeIds[f0][(e0 + 2) % 3];
328 int sharp_eid02 = sharp_edges[f0 * 3 + (e0 + 2) % 3];
329
330 int eid0, eid1, eid0p, eid1p;
331 int sharp_eid0, sharp_eid1, sharp_eid0p, sharp_eid1p;
332
333 eid0 = eid;
334 sharp_eid0 = sharp_eid;
335 edge_values[eid0] = DEdge(v0, vn);
336
337 eid1 = edge_values.size();
338 sharp_eid1 = sharp_eid;
339 edge_values.push_back(DEdge(vn, v1));
340 edge_diff.push_back(Vector2i());
341
342 eid0p = edge_values.size();
343 sharp_eid0p = 0;
344 edge_values.push_back(DEdge(vn, v0p));
345 edge_diff.push_back(Vector2i());
346
347 int f2 = is_boundary ? -1 : (nF++);
348 int f3 = nF++;
349 sharp_edges.push_back(0);
350 sharp_edges.push_back(0);
351 sharp_edges.push_back(0);
352 face_edgeIds.push_back(Vector3i());
353 face_edgeOrients.push_back(Vector3i());
354
355 if (nF > F.cols()) {
356 F.conservativeResize(F.rows(), std::max(nF, (int)F.cols() * 2));
357 face_spaces.resize(F.cols());
358 E2E.conservativeResize(F.cols() * 3);
359 diffs.resize(F.cols() * 3);
360 }
361
362 auto D01 = diffs[e0];
363 auto D1p = diffs[e0 / 3 * 3 + (e0 + 1) % 3];
364 auto Dp0 = diffs[e0 / 3 * 3 + (e0 + 2) % 3];
365
366 Vector2i D0n = D01 / 2;
367
368 auto orients1 = face_spaces[f0];
369 F.col(f0) << vn, v0p, v0;
370 face_edgeIds[f0] = Vector3i(eid0p, eid02, eid0);
371 sharp_edges[f0 * 3] = sharp_eid0p;
372 sharp_edges[f0 * 3 + 1] = sharp_eid02;
373 sharp_edges[f0 * 3 + 2] = sharp_eid0;
374
375 diffs[f0 * 3] = D01 + D1p - D0n;
376 diffs[f0 * 3 + 1] = Dp0;
377 diffs[f0 * 3 + 2] = D0n;
378 int o1 = e0 % 3, o2 = e1 % 3;
379 AnalyzeOrient(f0, Vector3i(0, orients1.d[(o1 + 2) % 3], orients1.d[o1]));
380 if (!is_boundary) {
381 auto orients2 = face_spaces[f1];
382 int eid11 = face_edgeIds[f1][(e1 + 1) % 3];
383 int sharp_eid11 = sharp_edges[f1 * 3 + (e1 + 1) % 3];
384 int eid12 = face_edgeIds[f1][(e1 + 2) % 3];
385 int sharp_eid12 = sharp_edges[f1 * 3 + (e1 + 2) % 3];
386
387 auto Ds10 = diffs[e1];
388 auto Ds0p = diffs[e1 / 3 * 3 + (e1 + 1) % 3];
389
390 auto Dsp1 = diffs[e1 / 3 * 3 + (e1 + 2) % 3];
391 int orient = 0;
392 while (rshift90(D01, orient) != Ds10) orient += 1;
393 Vector2i Dsn0 = rshift90(D0n, orient);
394
395 F.col(f1) << vn, v0, v1p;
396 eid1p = edge_values.size();
397 sharp_eid1p = 0;
398 edge_values.push_back(DEdge(vn, v1p));
399 edge_diff.push_back(Vector2i());
400
401 sharp_edges[f1 * 3] = sharp_eid0;
402 sharp_edges[f1 * 3 + 1] = sharp_eid11;
403 sharp_edges[f1 * 3 + 2] = sharp_eid1p;
404 face_edgeIds[f1] = (Vector3i(eid0, eid11, eid1p));
405 diffs[f1 * 3] = Dsn0;
406 diffs[f1 * 3 + 1] = Ds0p;
407 diffs[f1 * 3 + 2] = Dsp1 + (Ds10 - Dsn0);
408
409 AnalyzeOrient(f1, Vector3i(orients2.d[o2], orients2.d[(o2 + 1) % 3], 0));
410
411 face_spaces[f2] = face_spaces[f1];
412 sharp_edges[f2 * 3] = sharp_eid1p;
413 sharp_edges[f2 * 3 + 1] = sharp_eid12;
414 sharp_edges[f2 * 3 + 2] = sharp_eid1;
415 face_edgeIds[f2] = (Vector3i(eid1p, eid12, eid1));
416 F.col(f2) << vn, v1p, v1;
417 diffs[f2 * 3] = -Dsp1 - (Ds10 - Dsn0);
418 diffs[f2 * 3 + 1] = Dsp1;
419 diffs[f2 * 3 + 2] = Ds10 - Dsn0;
420
421 AnalyzeOrient(f2, Vector3i(0, orients2.d[(o2 + 2) % 3], orients2.d[o2]));
422 }
423 face_spaces[f3] = face_spaces[f0];
424 sharp_edges[f3 * 3] = sharp_eid1;
425 sharp_edges[f3 * 3 + 1] = sharp_eid01;
426 sharp_edges[f3 * 3 + 2] = sharp_eid0p;
427 face_edgeIds[f3] = (Vector3i(eid1, eid01, eid0p));
428 F.col(f3) << vn, v1, v0p;
429 diffs[f3 * 3] = D01 - D0n;
430 diffs[f3 * 3 + 1] = D1p;
431 diffs[f3 * 3 + 2] = D0n - (D01 + D1p);
432
433 AnalyzeOrient(f3, Vector3i(orients1.d[o1], orients1.d[(o1 + 1) % 3], 0));
434
435 FixOrient(f0);
436 if (!is_boundary) {
437 FixOrient(f1);
438 FixOrient(f2);
439 }
440 FixOrient(f3);
441
442 const int e0p = E2E[dedge_prev_3(e0)], e0n = E2E[dedge_next_3(e0)];
443
444 #define sE2E(a, b) \
445 E2E[a] = b; \
446 if (b != -1) E2E[b] = a;
447 sE2E(3 * f0 + 0, 3 * f3 + 2);
448 sE2E(3 * f0 + 1, e0p);
449 sE2E(3 * f3 + 1, e0n);
450 if (is_boundary) {
451 sE2E(3 * f0 + 2, -1);
452 sE2E(3 * f3 + 0, -1);
453 } else {
454 const int e1p = E2E[dedge_prev_3(e1)], e1n = E2E[dedge_next_3(e1)];
455 sE2E(3 * f0 + 2, 3 * f1 + 0);
456 sE2E(3 * f1 + 1, e1n);
457 sE2E(3 * f1 + 2, 3 * f2 + 0);
458 sE2E(3 * f2 + 1, e1p);
459 sE2E(3 * f2 + 2, 3 * f3 + 0);
460 }
461 #undef sE2E
462
463 V2E[v0] = 3 * f0 + 2;
464 V2E[vn] = 3 * f0 + 0;
465 V2E[v1] = 3 * f3 + 1;
466 V2E[v0p] = 3 * f0 + 1;
467 if (!is_boundary) V2E[v1p] = 3 * f1 + 2;
468
469 auto schedule = [&](int f) {
470 for (int i = 0; i < 3; ++i) {
471 if (abs(diffs[f * 3 + i][0]) > max_len || abs(diffs[f * 3 + i][1]) > max_len) {
472 EdgeLink e;
473 e.id = f * 3 + i;
474 e.length = (V.col(F((i + 1) % 3, f)) - V.col(F(i, f))).squaredNorm();
475 e.diff = diffs[f * 3 + i];
476 queue.push(e);
477 }
478 }
479 };
480
481 schedule(f0);
482 if (!is_boundary) {
483 schedule(f2);
484 schedule(f1);
485 };
486 schedule(f3);
487 }
488 F.conservativeResize(F.rows(), nF);
489 V.conservativeResize(V.rows(), nV);
490 N.conservativeResize(V.rows(), nV);
491 Q.conservativeResize(V.rows(), nV);
492 O.conservativeResize(V.rows(), nV);
493 if (S)
494 S->conservativeResize(S->rows(), nV);
495 V2E.conservativeResize(nV);
496 boundary.conservativeResize(nV);
497 nonmanifold.conservativeResize(nV);
498 E2E.conservativeResize(nF * 3);
499 for (int i = 0; i < F.cols(); ++i) {
500 for (int j = 0; j < 3; ++j) {
501 auto diff = edge_diff[face_edgeIds[i][j]];
502 if (abs(diff[0]) > 1 || abs(diff[1]) > 1) {
503 printf("wrong init %d %d!\n", face_edgeIds[i][j], i * 3 + j);
504 exit(0);
505 }
506 }
507 }
508 for (int i = 0; i < edge_diff.size(); ++i) {
509 if (abs(edge_diff[i][0]) > 1 || abs(edge_diff[i][1]) > 1) {
510 printf("wrong...\n");
511 exit(0);
512 }
513 }
514 }
515
516 } // namespace qflow
517