1 // This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
2 #include "meshoptimizer.h"
3
4 #include <assert.h>
5 #include <float.h>
6 #include <math.h>
7 #include <string.h>
8
9 #ifndef TRACE
10 #define TRACE 0
11 #endif
12
13 #if TRACE
14 #include <stdio.h>
15 #endif
16
17 // This work is based on:
18 // Michael Garland and Paul S. Heckbert. Surface simplification using quadric error metrics. 1997
19 // Michael Garland. Quadric-based polygonal surface simplification. 1999
20 // Peter Lindstrom. Out-of-Core Simplification of Large Polygonal Models. 2000
21 // Matthias Teschner, Bruno Heidelberger, Matthias Mueller, Danat Pomeranets, Markus Gross. Optimized Spatial Hashing for Collision Detection of Deformable Objects. 2003
22 // Peter Van Sandt, Yannis Chronis, Jignesh M. Patel. Efficiently Searching In-Memory Sorted Arrays: Revenge of the Interpolation Search? 2019
23 namespace meshopt
24 {
25
26 struct EdgeAdjacency
27 {
28 unsigned int* counts;
29 unsigned int* offsets;
30 unsigned int* data;
31 };
32
buildEdgeAdjacency(EdgeAdjacency & adjacency,const unsigned int * indices,size_t index_count,size_t vertex_count,meshopt_Allocator & allocator)33 static void buildEdgeAdjacency(EdgeAdjacency& adjacency, const unsigned int* indices, size_t index_count, size_t vertex_count, meshopt_Allocator& allocator)
34 {
35 size_t face_count = index_count / 3;
36
37 // allocate arrays
38 adjacency.counts = allocator.allocate<unsigned int>(vertex_count);
39 adjacency.offsets = allocator.allocate<unsigned int>(vertex_count);
40 adjacency.data = allocator.allocate<unsigned int>(index_count);
41
42 // fill edge counts
43 memset(adjacency.counts, 0, vertex_count * sizeof(unsigned int));
44
45 for (size_t i = 0; i < index_count; ++i)
46 {
47 assert(indices[i] < vertex_count);
48
49 adjacency.counts[indices[i]]++;
50 }
51
52 // fill offset table
53 unsigned int offset = 0;
54
55 for (size_t i = 0; i < vertex_count; ++i)
56 {
57 adjacency.offsets[i] = offset;
58 offset += adjacency.counts[i];
59 }
60
61 assert(offset == index_count);
62
63 // fill edge data
64 for (size_t i = 0; i < face_count; ++i)
65 {
66 unsigned int a = indices[i * 3 + 0], b = indices[i * 3 + 1], c = indices[i * 3 + 2];
67
68 adjacency.data[adjacency.offsets[a]++] = b;
69 adjacency.data[adjacency.offsets[b]++] = c;
70 adjacency.data[adjacency.offsets[c]++] = a;
71 }
72
73 // fix offsets that have been disturbed by the previous pass
74 for (size_t i = 0; i < vertex_count; ++i)
75 {
76 assert(adjacency.offsets[i] >= adjacency.counts[i]);
77
78 adjacency.offsets[i] -= adjacency.counts[i];
79 }
80 }
81
82 struct PositionHasher
83 {
84 const float* vertex_positions;
85 size_t vertex_stride_float;
86
hashmeshopt::PositionHasher87 size_t hash(unsigned int index) const
88 {
89 // MurmurHash2
90 const unsigned int m = 0x5bd1e995;
91 const int r = 24;
92
93 unsigned int h = 0;
94 const unsigned int* key = reinterpret_cast<const unsigned int*>(vertex_positions + index * vertex_stride_float);
95
96 for (size_t i = 0; i < 3; ++i)
97 {
98 unsigned int k = key[i];
99
100 k *= m;
101 k ^= k >> r;
102 k *= m;
103
104 h *= m;
105 h ^= k;
106 }
107
108 return h;
109 }
110
equalmeshopt::PositionHasher111 bool equal(unsigned int lhs, unsigned int rhs) const
112 {
113 return memcmp(vertex_positions + lhs * vertex_stride_float, vertex_positions + rhs * vertex_stride_float, sizeof(float) * 3) == 0;
114 }
115 };
116
hashBuckets2(size_t count)117 static size_t hashBuckets2(size_t count)
118 {
119 size_t buckets = 1;
120 while (buckets < count)
121 buckets *= 2;
122
123 return buckets;
124 }
125
126 template <typename T, typename Hash>
hashLookup2(T * table,size_t buckets,const Hash & hash,const T & key,const T & empty)127 static T* hashLookup2(T* table, size_t buckets, const Hash& hash, const T& key, const T& empty)
128 {
129 assert(buckets > 0);
130 assert((buckets & (buckets - 1)) == 0);
131
132 size_t hashmod = buckets - 1;
133 size_t bucket = hash.hash(key) & hashmod;
134
135 for (size_t probe = 0; probe <= hashmod; ++probe)
136 {
137 T& item = table[bucket];
138
139 if (item == empty)
140 return &item;
141
142 if (hash.equal(item, key))
143 return &item;
144
145 // hash collision, quadratic probing
146 bucket = (bucket + probe + 1) & hashmod;
147 }
148
149 assert(false && "Hash table is full"); // unreachable
150 return 0;
151 }
152
buildPositionRemap(unsigned int * remap,unsigned int * wedge,const float * vertex_positions_data,size_t vertex_count,size_t vertex_positions_stride,meshopt_Allocator & allocator)153 static void buildPositionRemap(unsigned int* remap, unsigned int* wedge, const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride, meshopt_Allocator& allocator)
154 {
155 PositionHasher hasher = {vertex_positions_data, vertex_positions_stride / sizeof(float)};
156
157 size_t table_size = hashBuckets2(vertex_count);
158 unsigned int* table = allocator.allocate<unsigned int>(table_size);
159 memset(table, -1, table_size * sizeof(unsigned int));
160
161 // build forward remap: for each vertex, which other (canonical) vertex does it map to?
162 // we use position equivalence for this, and remap vertices to other existing vertices
163 for (size_t i = 0; i < vertex_count; ++i)
164 {
165 unsigned int index = unsigned(i);
166 unsigned int* entry = hashLookup2(table, table_size, hasher, index, ~0u);
167
168 if (*entry == ~0u)
169 *entry = index;
170
171 remap[index] = *entry;
172 }
173
174 // build wedge table: for each vertex, which other vertex is the next wedge that also maps to the same vertex?
175 // entries in table form a (cyclic) wedge loop per vertex; for manifold vertices, wedge[i] == remap[i] == i
176 for (size_t i = 0; i < vertex_count; ++i)
177 wedge[i] = unsigned(i);
178
179 for (size_t i = 0; i < vertex_count; ++i)
180 if (remap[i] != i)
181 {
182 unsigned int r = remap[i];
183
184 wedge[i] = wedge[r];
185 wedge[r] = unsigned(i);
186 }
187 }
188
189 enum VertexKind
190 {
191 Kind_Manifold, // not on an attribute seam, not on any boundary
192 Kind_Border, // not on an attribute seam, has exactly two open edges
193 Kind_Seam, // on an attribute seam with exactly two attribute seam edges
194 Kind_Complex, // none of the above; these vertices can move as long as all wedges move to the target vertex
195 Kind_Locked, // none of the above; these vertices can't move
196
197 Kind_Count
198 };
199
200 // manifold vertices can collapse onto anything
201 // border/seam vertices can only be collapsed onto border/seam respectively
202 // complex vertices can collapse onto complex/locked
203 // a rule of thumb is that collapsing kind A into kind B preserves the kind B in the target vertex
204 // for example, while we could collapse Complex into Manifold, this would mean the target vertex isn't Manifold anymore
205 const unsigned char kCanCollapse[Kind_Count][Kind_Count] = {
206 {1, 1, 1, 1, 1},
207 {0, 1, 0, 0, 0},
208 {0, 0, 1, 0, 0},
209 {0, 0, 0, 1, 1},
210 {0, 0, 0, 0, 0},
211 };
212
213 // if a vertex is manifold or seam, adjoining edges are guaranteed to have an opposite edge
214 // note that for seam edges, the opposite edge isn't present in the attribute-based topology
215 // but is present if you consider a position-only mesh variant
216 const unsigned char kHasOpposite[Kind_Count][Kind_Count] = {
217 {1, 1, 1, 0, 1},
218 {1, 0, 1, 0, 0},
219 {1, 1, 1, 0, 1},
220 {0, 0, 0, 0, 0},
221 {1, 0, 1, 0, 0},
222 };
223
hasEdge(const EdgeAdjacency & adjacency,unsigned int a,unsigned int b)224 static bool hasEdge(const EdgeAdjacency& adjacency, unsigned int a, unsigned int b)
225 {
226 unsigned int count = adjacency.counts[a];
227 const unsigned int* data = adjacency.data + adjacency.offsets[a];
228
229 for (size_t i = 0; i < count; ++i)
230 if (data[i] == b)
231 return true;
232
233 return false;
234 }
235
findWedgeEdge(const EdgeAdjacency & adjacency,const unsigned int * wedge,unsigned int a,unsigned int b)236 static unsigned int findWedgeEdge(const EdgeAdjacency& adjacency, const unsigned int* wedge, unsigned int a, unsigned int b)
237 {
238 unsigned int v = a;
239
240 do
241 {
242 if (hasEdge(adjacency, v, b))
243 return v;
244
245 v = wedge[v];
246 } while (v != a);
247
248 return ~0u;
249 }
250
countOpenEdges(const EdgeAdjacency & adjacency,unsigned int vertex,unsigned int * last=0)251 static size_t countOpenEdges(const EdgeAdjacency& adjacency, unsigned int vertex, unsigned int* last = 0)
252 {
253 size_t result = 0;
254
255 unsigned int count = adjacency.counts[vertex];
256 const unsigned int* data = adjacency.data + adjacency.offsets[vertex];
257
258 for (size_t i = 0; i < count; ++i)
259 if (!hasEdge(adjacency, data[i], vertex))
260 {
261 result++;
262
263 if (last)
264 *last = data[i];
265 }
266
267 return result;
268 }
269
classifyVertices(unsigned char * result,unsigned int * loop,size_t vertex_count,const EdgeAdjacency & adjacency,const unsigned int * remap,const unsigned int * wedge)270 static void classifyVertices(unsigned char* result, unsigned int* loop, size_t vertex_count, const EdgeAdjacency& adjacency, const unsigned int* remap, const unsigned int* wedge)
271 {
272 for (size_t i = 0; i < vertex_count; ++i)
273 loop[i] = ~0u;
274
275 #if TRACE
276 size_t lockedstats[4] = {};
277 #define TRACELOCKED(i) lockedstats[i]++;
278 #else
279 #define TRACELOCKED(i) (void)0
280 #endif
281
282 for (size_t i = 0; i < vertex_count; ++i)
283 {
284 if (remap[i] == i)
285 {
286 if (wedge[i] == i)
287 {
288 // no attribute seam, need to check if it's manifold
289 unsigned int v = 0;
290 size_t edges = countOpenEdges(adjacency, unsigned(i), &v);
291
292 // note: we classify any vertices with no open edges as manifold
293 // this is technically incorrect - if 4 triangles share an edge, we'll classify vertices as manifold
294 // it's unclear if this is a problem in practice
295 // also note that we classify vertices as border if they have *one* open edge, not two
296 // this is because we only have half-edges - so a border vertex would have one incoming and one outgoing edge
297 if (edges == 0)
298 {
299 result[i] = Kind_Manifold;
300 }
301 else if (edges == 1)
302 {
303 result[i] = Kind_Border;
304 loop[i] = v;
305 }
306 else
307 {
308 result[i] = Kind_Locked;
309 TRACELOCKED(0);
310 }
311 }
312 else if (wedge[wedge[i]] == i)
313 {
314 // attribute seam; need to distinguish between Seam and Locked
315 unsigned int a = 0;
316 size_t a_count = countOpenEdges(adjacency, unsigned(i), &a);
317 unsigned int b = 0;
318 size_t b_count = countOpenEdges(adjacency, wedge[i], &b);
319
320 // seam should have one open half-edge for each vertex, and the edges need to "connect" - point to the same vertex post-remap
321 if (a_count == 1 && b_count == 1)
322 {
323 unsigned int ao = findWedgeEdge(adjacency, wedge, a, wedge[i]);
324 unsigned int bo = findWedgeEdge(adjacency, wedge, b, unsigned(i));
325
326 if (ao != ~0u && bo != ~0u)
327 {
328 result[i] = Kind_Seam;
329
330 loop[i] = a;
331 loop[wedge[i]] = b;
332 }
333 else
334 {
335 result[i] = Kind_Locked;
336 TRACELOCKED(1);
337 }
338 }
339 else
340 {
341 result[i] = Kind_Locked;
342 TRACELOCKED(2);
343 }
344 }
345 else
346 {
347 // more than one vertex maps to this one; we don't have classification available
348 result[i] = Kind_Locked;
349 TRACELOCKED(3);
350 }
351 }
352 else
353 {
354 assert(remap[i] < i);
355
356 result[i] = result[remap[i]];
357 }
358 }
359
360 #if TRACE
361 printf("locked: many open edges %d, disconnected seam %d, many seam edges %d, many wedges %d\n",
362 int(lockedstats[0]), int(lockedstats[1]), int(lockedstats[2]), int(lockedstats[3]));
363 #endif
364 }
365
366 struct Vector3
367 {
368 float x, y, z;
369 };
370
rescalePositions(Vector3 * result,const float * vertex_positions_data,size_t vertex_count,size_t vertex_positions_stride)371 static void rescalePositions(Vector3* result, const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride)
372 {
373 size_t vertex_stride_float = vertex_positions_stride / sizeof(float);
374
375 float minv[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
376 float maxv[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX};
377
378 for (size_t i = 0; i < vertex_count; ++i)
379 {
380 const float* v = vertex_positions_data + i * vertex_stride_float;
381
382 result[i].x = v[0];
383 result[i].y = v[1];
384 result[i].z = v[2];
385
386 for (int j = 0; j < 3; ++j)
387 {
388 float vj = v[j];
389
390 minv[j] = minv[j] > vj ? vj : minv[j];
391 maxv[j] = maxv[j] < vj ? vj : maxv[j];
392 }
393 }
394
395 float extent = 0.f;
396
397 extent = (maxv[0] - minv[0]) < extent ? extent : (maxv[0] - minv[0]);
398 extent = (maxv[1] - minv[1]) < extent ? extent : (maxv[1] - minv[1]);
399 extent = (maxv[2] - minv[2]) < extent ? extent : (maxv[2] - minv[2]);
400
401 float scale = extent == 0 ? 0.f : 1.f / extent;
402
403 for (size_t i = 0; i < vertex_count; ++i)
404 {
405 result[i].x = (result[i].x - minv[0]) * scale;
406 result[i].y = (result[i].y - minv[1]) * scale;
407 result[i].z = (result[i].z - minv[2]) * scale;
408 }
409 }
410
411 struct Quadric
412 {
413 float a00, a11, a22;
414 float a10, a20, a21;
415 float b0, b1, b2, c;
416 float w;
417 };
418
419 struct Collapse
420 {
421 unsigned int v0;
422 unsigned int v1;
423
424 union {
425 unsigned int bidi;
426 float error;
427 unsigned int errorui;
428 };
429 };
430
normalize(Vector3 & v)431 static float normalize(Vector3& v)
432 {
433 float length = sqrtf(v.x * v.x + v.y * v.y + v.z * v.z);
434
435 if (length > 0)
436 {
437 v.x /= length;
438 v.y /= length;
439 v.z /= length;
440 }
441
442 return length;
443 }
444
quadricAdd(Quadric & Q,const Quadric & R)445 static void quadricAdd(Quadric& Q, const Quadric& R)
446 {
447 Q.a00 += R.a00;
448 Q.a11 += R.a11;
449 Q.a22 += R.a22;
450 Q.a10 += R.a10;
451 Q.a20 += R.a20;
452 Q.a21 += R.a21;
453 Q.b0 += R.b0;
454 Q.b1 += R.b1;
455 Q.b2 += R.b2;
456 Q.c += R.c;
457 Q.w += R.w;
458 }
459
quadricError(const Quadric & Q,const Vector3 & v)460 static float quadricError(const Quadric& Q, const Vector3& v)
461 {
462 float rx = Q.b0;
463 float ry = Q.b1;
464 float rz = Q.b2;
465
466 rx += Q.a10 * v.y;
467 ry += Q.a21 * v.z;
468 rz += Q.a20 * v.x;
469
470 rx *= 2;
471 ry *= 2;
472 rz *= 2;
473
474 rx += Q.a00 * v.x;
475 ry += Q.a11 * v.y;
476 rz += Q.a22 * v.z;
477
478 float r = Q.c;
479 r += rx * v.x;
480 r += ry * v.y;
481 r += rz * v.z;
482
483 float s = Q.w == 0.f ? 0.f : 1.f / Q.w;
484
485 return fabsf(r) * s;
486 }
487
quadricFromPlane(Quadric & Q,float a,float b,float c,float d,float w)488 static void quadricFromPlane(Quadric& Q, float a, float b, float c, float d, float w)
489 {
490 float aw = a * w;
491 float bw = b * w;
492 float cw = c * w;
493 float dw = d * w;
494
495 Q.a00 = a * aw;
496 Q.a11 = b * bw;
497 Q.a22 = c * cw;
498 Q.a10 = a * bw;
499 Q.a20 = a * cw;
500 Q.a21 = b * cw;
501 Q.b0 = a * dw;
502 Q.b1 = b * dw;
503 Q.b2 = c * dw;
504 Q.c = d * dw;
505 Q.w = w;
506 }
507
quadricFromPoint(Quadric & Q,float x,float y,float z,float w)508 static void quadricFromPoint(Quadric& Q, float x, float y, float z, float w)
509 {
510 // we need to encode (x - X) ^ 2 + (y - Y)^2 + (z - Z)^2 into the quadric
511 Q.a00 = w;
512 Q.a11 = w;
513 Q.a22 = w;
514 Q.a10 = 0.f;
515 Q.a20 = 0.f;
516 Q.a21 = 0.f;
517 Q.b0 = -2.f * x * w;
518 Q.b1 = -2.f * y * w;
519 Q.b2 = -2.f * z * w;
520 Q.c = (x * x + y * y + z * z) * w;
521 Q.w = w;
522 }
523
quadricFromTriangle(Quadric & Q,const Vector3 & p0,const Vector3 & p1,const Vector3 & p2,float weight)524 static void quadricFromTriangle(Quadric& Q, const Vector3& p0, const Vector3& p1, const Vector3& p2, float weight)
525 {
526 Vector3 p10 = {p1.x - p0.x, p1.y - p0.y, p1.z - p0.z};
527 Vector3 p20 = {p2.x - p0.x, p2.y - p0.y, p2.z - p0.z};
528
529 // normal = cross(p1 - p0, p2 - p0)
530 Vector3 normal = {p10.y * p20.z - p10.z * p20.y, p10.z * p20.x - p10.x * p20.z, p10.x * p20.y - p10.y * p20.x};
531 float area = normalize(normal);
532
533 float distance = normal.x * p0.x + normal.y * p0.y + normal.z * p0.z;
534
535 // we use sqrtf(area) so that the error is scaled linearly; this tends to improve silhouettes
536 quadricFromPlane(Q, normal.x, normal.y, normal.z, -distance, sqrtf(area) * weight);
537 }
538
quadricFromTriangleEdge(Quadric & Q,const Vector3 & p0,const Vector3 & p1,const Vector3 & p2,float weight)539 static void quadricFromTriangleEdge(Quadric& Q, const Vector3& p0, const Vector3& p1, const Vector3& p2, float weight)
540 {
541 Vector3 p10 = {p1.x - p0.x, p1.y - p0.y, p1.z - p0.z};
542 float length = normalize(p10);
543
544 // p20p = length of projection of p2-p0 onto normalize(p1 - p0)
545 Vector3 p20 = {p2.x - p0.x, p2.y - p0.y, p2.z - p0.z};
546 float p20p = p20.x * p10.x + p20.y * p10.y + p20.z * p10.z;
547
548 // normal = altitude of triangle from point p2 onto edge p1-p0
549 Vector3 normal = {p20.x - p10.x * p20p, p20.y - p10.y * p20p, p20.z - p10.z * p20p};
550 normalize(normal);
551
552 float distance = normal.x * p0.x + normal.y * p0.y + normal.z * p0.z;
553
554 // note: the weight is scaled linearly with edge length; this has to match the triangle weight
555 quadricFromPlane(Q, normal.x, normal.y, normal.z, -distance, length * weight);
556 }
557
fillFaceQuadrics(Quadric * vertex_quadrics,const unsigned int * indices,size_t index_count,const Vector3 * vertex_positions,const unsigned int * remap)558 static void fillFaceQuadrics(Quadric* vertex_quadrics, const unsigned int* indices, size_t index_count, const Vector3* vertex_positions, const unsigned int* remap)
559 {
560 for (size_t i = 0; i < index_count; i += 3)
561 {
562 unsigned int i0 = indices[i + 0];
563 unsigned int i1 = indices[i + 1];
564 unsigned int i2 = indices[i + 2];
565
566 Quadric Q;
567 quadricFromTriangle(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], 1.f);
568
569 quadricAdd(vertex_quadrics[remap[i0]], Q);
570 quadricAdd(vertex_quadrics[remap[i1]], Q);
571 quadricAdd(vertex_quadrics[remap[i2]], Q);
572 }
573 }
574
fillEdgeQuadrics(Quadric * vertex_quadrics,const unsigned int * indices,size_t index_count,const Vector3 * vertex_positions,const unsigned int * remap,const unsigned char * vertex_kind,const unsigned int * loop)575 static void fillEdgeQuadrics(Quadric* vertex_quadrics, const unsigned int* indices, size_t index_count, const Vector3* vertex_positions, const unsigned int* remap, const unsigned char* vertex_kind, const unsigned int* loop)
576 {
577 for (size_t i = 0; i < index_count; i += 3)
578 {
579 static const int next[3] = {1, 2, 0};
580
581 for (int e = 0; e < 3; ++e)
582 {
583 unsigned int i0 = indices[i + e];
584 unsigned int i1 = indices[i + next[e]];
585
586 unsigned char k0 = vertex_kind[i0];
587 unsigned char k1 = vertex_kind[i1];
588
589 // check that i0 and i1 are border/seam and are on the same edge loop
590 // loop[] tracks half edges so we only need to check i0->i1
591 if (k0 != k1 || (k0 != Kind_Border && k0 != Kind_Seam) || loop[i0] != i1)
592 continue;
593
594 unsigned int i2 = indices[i + next[next[e]]];
595
596 // we try hard to maintain border edge geometry; seam edges can move more freely
597 // due to topological restrictions on collapses, seam quadrics slightly improves collapse structure but aren't critical
598 const float kEdgeWeightSeam = 1.f;
599 const float kEdgeWeightBorder = 10.f;
600
601 float edgeWeight = (k0 == Kind_Seam) ? kEdgeWeightSeam : kEdgeWeightBorder;
602
603 Quadric Q;
604 quadricFromTriangleEdge(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], edgeWeight);
605
606 quadricAdd(vertex_quadrics[remap[i0]], Q);
607 quadricAdd(vertex_quadrics[remap[i1]], Q);
608 }
609 }
610 }
611
pickEdgeCollapses(Collapse * collapses,const unsigned int * indices,size_t index_count,const unsigned int * remap,const unsigned char * vertex_kind,const unsigned int * loop)612 static size_t pickEdgeCollapses(Collapse* collapses, const unsigned int* indices, size_t index_count, const unsigned int* remap, const unsigned char* vertex_kind, const unsigned int* loop)
613 {
614 size_t collapse_count = 0;
615
616 for (size_t i = 0; i < index_count; i += 3)
617 {
618 static const int next[3] = {1, 2, 0};
619
620 for (int e = 0; e < 3; ++e)
621 {
622 unsigned int i0 = indices[i + e];
623 unsigned int i1 = indices[i + next[e]];
624
625 // this can happen either when input has a zero-length edge, or when we perform collapses for complex
626 // topology w/seams and collapse a manifold vertex that connects to both wedges onto one of them
627 // we leave edges like this alone since they may be important for preserving mesh integrity
628 if (remap[i0] == remap[i1])
629 continue;
630
631 unsigned char k0 = vertex_kind[i0];
632 unsigned char k1 = vertex_kind[i1];
633
634 // the edge has to be collapsible in at least one direction
635 if (!(kCanCollapse[k0][k1] | kCanCollapse[k1][k0]))
636 continue;
637
638 // manifold and seam edges should occur twice (i0->i1 and i1->i0) - skip redundant edges
639 if (kHasOpposite[k0][k1] && remap[i1] > remap[i0])
640 continue;
641
642 // two vertices are on a border or a seam, but there's no direct edge between them
643 // this indicates that they belong to two different edge loops and we should not collapse this edge
644 // loop[] tracks half edges so we only need to check i0->i1
645 if (k0 == k1 && (k0 == Kind_Border || k0 == Kind_Seam) && loop[i0] != i1)
646 continue;
647
648 // edge can be collapsed in either direction - we will pick the one with minimum error
649 // note: we evaluate error later during collapse ranking, here we just tag the edge as bidirectional
650 if (kCanCollapse[k0][k1] & kCanCollapse[k1][k0])
651 {
652 Collapse c = {i0, i1, {/* bidi= */ 1}};
653 collapses[collapse_count++] = c;
654 }
655 else
656 {
657 // edge can only be collapsed in one direction
658 unsigned int e0 = kCanCollapse[k0][k1] ? i0 : i1;
659 unsigned int e1 = kCanCollapse[k0][k1] ? i1 : i0;
660
661 Collapse c = {e0, e1, {/* bidi= */ 0}};
662 collapses[collapse_count++] = c;
663 }
664 }
665 }
666
667 return collapse_count;
668 }
669
rankEdgeCollapses(Collapse * collapses,size_t collapse_count,const Vector3 * vertex_positions,const Quadric * vertex_quadrics,const unsigned int * remap)670 static void rankEdgeCollapses(Collapse* collapses, size_t collapse_count, const Vector3* vertex_positions, const Quadric* vertex_quadrics, const unsigned int* remap)
671 {
672 for (size_t i = 0; i < collapse_count; ++i)
673 {
674 Collapse& c = collapses[i];
675
676 unsigned int i0 = c.v0;
677 unsigned int i1 = c.v1;
678
679 // most edges are bidirectional which means we need to evaluate errors for two collapses
680 // to keep this code branchless we just use the same edge for unidirectional edges
681 unsigned int j0 = c.bidi ? i1 : i0;
682 unsigned int j1 = c.bidi ? i0 : i1;
683
684 const Quadric& qi = vertex_quadrics[remap[i0]];
685 const Quadric& qj = vertex_quadrics[remap[j0]];
686
687 float ei = quadricError(qi, vertex_positions[i1]);
688 float ej = quadricError(qj, vertex_positions[j1]);
689
690 // pick edge direction with minimal error
691 c.v0 = ei <= ej ? i0 : j0;
692 c.v1 = ei <= ej ? i1 : j1;
693 c.error = ei <= ej ? ei : ej;
694 }
695 }
696
697 #if TRACE > 1
dumpEdgeCollapses(const Collapse * collapses,size_t collapse_count,const unsigned char * vertex_kind)698 static void dumpEdgeCollapses(const Collapse* collapses, size_t collapse_count, const unsigned char* vertex_kind)
699 {
700 size_t ckinds[Kind_Count][Kind_Count] = {};
701 float cerrors[Kind_Count][Kind_Count] = {};
702
703 for (int k0 = 0; k0 < Kind_Count; ++k0)
704 for (int k1 = 0; k1 < Kind_Count; ++k1)
705 cerrors[k0][k1] = FLT_MAX;
706
707 for (size_t i = 0; i < collapse_count; ++i)
708 {
709 unsigned int i0 = collapses[i].v0;
710 unsigned int i1 = collapses[i].v1;
711
712 unsigned char k0 = vertex_kind[i0];
713 unsigned char k1 = vertex_kind[i1];
714
715 ckinds[k0][k1]++;
716 cerrors[k0][k1] = (collapses[i].error < cerrors[k0][k1]) ? collapses[i].error : cerrors[k0][k1];
717 }
718
719 for (int k0 = 0; k0 < Kind_Count; ++k0)
720 for (int k1 = 0; k1 < Kind_Count; ++k1)
721 if (ckinds[k0][k1])
722 printf("collapses %d -> %d: %d, min error %e\n", k0, k1, int(ckinds[k0][k1]), cerrors[k0][k1]);
723 }
724
dumpLockedCollapses(const unsigned int * indices,size_t index_count,const unsigned char * vertex_kind)725 static void dumpLockedCollapses(const unsigned int* indices, size_t index_count, const unsigned char* vertex_kind)
726 {
727 size_t locked_collapses[Kind_Count][Kind_Count] = {};
728
729 for (size_t i = 0; i < index_count; i += 3)
730 {
731 static const int next[3] = {1, 2, 0};
732
733 for (int e = 0; e < 3; ++e)
734 {
735 unsigned int i0 = indices[i + e];
736 unsigned int i1 = indices[i + next[e]];
737
738 unsigned char k0 = vertex_kind[i0];
739 unsigned char k1 = vertex_kind[i1];
740
741 locked_collapses[k0][k1] += !kCanCollapse[k0][k1] && !kCanCollapse[k1][k0];
742 }
743 }
744
745 for (int k0 = 0; k0 < Kind_Count; ++k0)
746 for (int k1 = 0; k1 < Kind_Count; ++k1)
747 if (locked_collapses[k0][k1])
748 printf("locked collapses %d -> %d: %d\n", k0, k1, int(locked_collapses[k0][k1]));
749 }
750 #endif
751
sortEdgeCollapses(unsigned int * sort_order,const Collapse * collapses,size_t collapse_count)752 static void sortEdgeCollapses(unsigned int* sort_order, const Collapse* collapses, size_t collapse_count)
753 {
754 const int sort_bits = 11;
755
756 // fill histogram for counting sort
757 unsigned int histogram[1 << sort_bits];
758 memset(histogram, 0, sizeof(histogram));
759
760 for (size_t i = 0; i < collapse_count; ++i)
761 {
762 // skip sign bit since error is non-negative
763 unsigned int key = (collapses[i].errorui << 1) >> (32 - sort_bits);
764
765 histogram[key]++;
766 }
767
768 // compute offsets based on histogram data
769 size_t histogram_sum = 0;
770
771 for (size_t i = 0; i < 1 << sort_bits; ++i)
772 {
773 size_t count = histogram[i];
774 histogram[i] = unsigned(histogram_sum);
775 histogram_sum += count;
776 }
777
778 assert(histogram_sum == collapse_count);
779
780 // compute sort order based on offsets
781 for (size_t i = 0; i < collapse_count; ++i)
782 {
783 // skip sign bit since error is non-negative
784 unsigned int key = (collapses[i].errorui << 1) >> (32 - sort_bits);
785
786 sort_order[histogram[key]++] = unsigned(i);
787 }
788 }
789
performEdgeCollapses(unsigned int * collapse_remap,unsigned char * collapse_locked,Quadric * vertex_quadrics,const Collapse * collapses,size_t collapse_count,const unsigned int * collapse_order,const unsigned int * remap,const unsigned int * wedge,const unsigned char * vertex_kind,size_t triangle_collapse_goal,float error_goal,float error_limit)790 static size_t performEdgeCollapses(unsigned int* collapse_remap, unsigned char* collapse_locked, Quadric* vertex_quadrics, const Collapse* collapses, size_t collapse_count, const unsigned int* collapse_order, const unsigned int* remap, const unsigned int* wedge, const unsigned char* vertex_kind, size_t triangle_collapse_goal, float error_goal, float error_limit)
791 {
792 size_t edge_collapses = 0;
793 size_t triangle_collapses = 0;
794
795 for (size_t i = 0; i < collapse_count; ++i)
796 {
797 const Collapse& c = collapses[collapse_order[i]];
798
799 if (c.error > error_limit)
800 break;
801
802 if (c.error > error_goal && triangle_collapses > triangle_collapse_goal / 10)
803 break;
804
805 if (triangle_collapses >= triangle_collapse_goal)
806 break;
807
808 unsigned int i0 = c.v0;
809 unsigned int i1 = c.v1;
810
811 unsigned int r0 = remap[i0];
812 unsigned int r1 = remap[i1];
813
814 // we don't collapse vertices that had source or target vertex involved in a collapse
815 // it's important to not move the vertices twice since it complicates the tracking/remapping logic
816 // it's important to not move other vertices towards a moved vertex to preserve error since we don't re-rank collapses mid-pass
817 if (collapse_locked[r0] | collapse_locked[r1])
818 continue;
819
820 assert(collapse_remap[r0] == r0);
821 assert(collapse_remap[r1] == r1);
822
823 quadricAdd(vertex_quadrics[r1], vertex_quadrics[r0]);
824
825 if (vertex_kind[i0] == Kind_Complex)
826 {
827 unsigned int v = i0;
828
829 do
830 {
831 collapse_remap[v] = r1;
832 v = wedge[v];
833 } while (v != i0);
834 }
835 else if (vertex_kind[i0] == Kind_Seam)
836 {
837 // remap v0 to v1 and seam pair of v0 to seam pair of v1
838 unsigned int s0 = wedge[i0];
839 unsigned int s1 = wedge[i1];
840
841 assert(s0 != i0 && s1 != i1);
842 assert(wedge[s0] == i0 && wedge[s1] == i1);
843
844 collapse_remap[i0] = i1;
845 collapse_remap[s0] = s1;
846 }
847 else
848 {
849 assert(wedge[i0] == i0);
850
851 collapse_remap[i0] = i1;
852 }
853
854 collapse_locked[r0] = 1;
855 collapse_locked[r1] = 1;
856
857 // border edges collapse 1 triangle, other edges collapse 2 or more
858 triangle_collapses += (vertex_kind[i0] == Kind_Border) ? 1 : 2;
859 edge_collapses++;
860 }
861
862 return edge_collapses;
863 }
864
remapIndexBuffer(unsigned int * indices,size_t index_count,const unsigned int * collapse_remap)865 static size_t remapIndexBuffer(unsigned int* indices, size_t index_count, const unsigned int* collapse_remap)
866 {
867 size_t write = 0;
868
869 for (size_t i = 0; i < index_count; i += 3)
870 {
871 unsigned int v0 = collapse_remap[indices[i + 0]];
872 unsigned int v1 = collapse_remap[indices[i + 1]];
873 unsigned int v2 = collapse_remap[indices[i + 2]];
874
875 // we never move the vertex twice during a single pass
876 assert(collapse_remap[v0] == v0);
877 assert(collapse_remap[v1] == v1);
878 assert(collapse_remap[v2] == v2);
879
880 if (v0 != v1 && v0 != v2 && v1 != v2)
881 {
882 indices[write + 0] = v0;
883 indices[write + 1] = v1;
884 indices[write + 2] = v2;
885 write += 3;
886 }
887 }
888
889 return write;
890 }
891
remapEdgeLoops(unsigned int * loop,size_t vertex_count,const unsigned int * collapse_remap)892 static void remapEdgeLoops(unsigned int* loop, size_t vertex_count, const unsigned int* collapse_remap)
893 {
894 for (size_t i = 0; i < vertex_count; ++i)
895 {
896 if (loop[i] != ~0u)
897 {
898 unsigned int l = loop[i];
899 unsigned int r = collapse_remap[l];
900
901 // i == r is a special case when the seam edge is collapsed in a direction opposite to where loop goes
902 loop[i] = (i == r) ? loop[l] : r;
903 }
904 }
905 }
906
907 struct CellHasher
908 {
909 const unsigned int* vertex_ids;
910
hashmeshopt::CellHasher911 size_t hash(unsigned int i) const
912 {
913 unsigned int h = vertex_ids[i];
914
915 // MurmurHash2 finalizer
916 h ^= h >> 13;
917 h *= 0x5bd1e995;
918 h ^= h >> 15;
919 return h;
920 }
921
equalmeshopt::CellHasher922 bool equal(unsigned int lhs, unsigned int rhs) const
923 {
924 return vertex_ids[lhs] == vertex_ids[rhs];
925 }
926 };
927
928 struct IdHasher
929 {
hashmeshopt::IdHasher930 size_t hash(unsigned int id) const
931 {
932 unsigned int h = id;
933
934 // MurmurHash2 finalizer
935 h ^= h >> 13;
936 h *= 0x5bd1e995;
937 h ^= h >> 15;
938 return h;
939 }
940
equalmeshopt::IdHasher941 bool equal(unsigned int lhs, unsigned int rhs) const
942 {
943 return lhs == rhs;
944 }
945 };
946
947 struct TriangleHasher
948 {
949 unsigned int* indices;
950
hashmeshopt::TriangleHasher951 size_t hash(unsigned int i) const
952 {
953 const unsigned int* tri = indices + i * 3;
954
955 // Optimized Spatial Hashing for Collision Detection of Deformable Objects
956 return (tri[0] * 73856093) ^ (tri[1] * 19349663) ^ (tri[2] * 83492791);
957 }
958
equalmeshopt::TriangleHasher959 bool equal(unsigned int lhs, unsigned int rhs) const
960 {
961 const unsigned int* lt = indices + lhs * 3;
962 const unsigned int* rt = indices + rhs * 3;
963
964 return lt[0] == rt[0] && lt[1] == rt[1] && lt[2] == rt[2];
965 }
966 };
967
computeVertexIds(unsigned int * vertex_ids,const Vector3 * vertex_positions,size_t vertex_count,int grid_size)968 static void computeVertexIds(unsigned int* vertex_ids, const Vector3* vertex_positions, size_t vertex_count, int grid_size)
969 {
970 assert(grid_size >= 1 && grid_size <= 1024);
971 float cell_scale = float(grid_size - 1);
972
973 for (size_t i = 0; i < vertex_count; ++i)
974 {
975 const Vector3& v = vertex_positions[i];
976
977 int xi = int(v.x * cell_scale + 0.5f);
978 int yi = int(v.y * cell_scale + 0.5f);
979 int zi = int(v.z * cell_scale + 0.5f);
980
981 vertex_ids[i] = (xi << 20) | (yi << 10) | zi;
982 }
983 }
984
countTriangles(const unsigned int * vertex_ids,const unsigned int * indices,size_t index_count)985 static size_t countTriangles(const unsigned int* vertex_ids, const unsigned int* indices, size_t index_count)
986 {
987 size_t result = 0;
988
989 for (size_t i = 0; i < index_count; i += 3)
990 {
991 unsigned int id0 = vertex_ids[indices[i + 0]];
992 unsigned int id1 = vertex_ids[indices[i + 1]];
993 unsigned int id2 = vertex_ids[indices[i + 2]];
994
995 result += (id0 != id1) & (id0 != id2) & (id1 != id2);
996 }
997
998 return result;
999 }
1000
fillVertexCells(unsigned int * table,size_t table_size,unsigned int * vertex_cells,const unsigned int * vertex_ids,size_t vertex_count)1001 static size_t fillVertexCells(unsigned int* table, size_t table_size, unsigned int* vertex_cells, const unsigned int* vertex_ids, size_t vertex_count)
1002 {
1003 CellHasher hasher = {vertex_ids};
1004
1005 memset(table, -1, table_size * sizeof(unsigned int));
1006
1007 size_t result = 0;
1008
1009 for (size_t i = 0; i < vertex_count; ++i)
1010 {
1011 unsigned int* entry = hashLookup2(table, table_size, hasher, unsigned(i), ~0u);
1012
1013 if (*entry == ~0u)
1014 {
1015 *entry = unsigned(i);
1016 vertex_cells[i] = unsigned(result++);
1017 }
1018 else
1019 {
1020 vertex_cells[i] = vertex_cells[*entry];
1021 }
1022 }
1023
1024 return result;
1025 }
1026
countVertexCells(unsigned int * table,size_t table_size,const unsigned int * vertex_ids,size_t vertex_count)1027 static size_t countVertexCells(unsigned int* table, size_t table_size, const unsigned int* vertex_ids, size_t vertex_count)
1028 {
1029 IdHasher hasher;
1030
1031 memset(table, -1, table_size * sizeof(unsigned int));
1032
1033 size_t result = 0;
1034
1035 for (size_t i = 0; i < vertex_count; ++i)
1036 {
1037 unsigned int id = vertex_ids[i];
1038 unsigned int* entry = hashLookup2(table, table_size, hasher, id, ~0u);
1039
1040 result += (*entry == ~0u);
1041 *entry = id;
1042 }
1043
1044 return result;
1045 }
1046
fillCellQuadrics(Quadric * cell_quadrics,const unsigned int * indices,size_t index_count,const Vector3 * vertex_positions,const unsigned int * vertex_cells)1047 static void fillCellQuadrics(Quadric* cell_quadrics, const unsigned int* indices, size_t index_count, const Vector3* vertex_positions, const unsigned int* vertex_cells)
1048 {
1049 for (size_t i = 0; i < index_count; i += 3)
1050 {
1051 unsigned int i0 = indices[i + 0];
1052 unsigned int i1 = indices[i + 1];
1053 unsigned int i2 = indices[i + 2];
1054
1055 unsigned int c0 = vertex_cells[i0];
1056 unsigned int c1 = vertex_cells[i1];
1057 unsigned int c2 = vertex_cells[i2];
1058
1059 bool single_cell = (c0 == c1) & (c0 == c2);
1060
1061 Quadric Q;
1062 quadricFromTriangle(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], single_cell ? 3.f : 1.f);
1063
1064 if (single_cell)
1065 {
1066 quadricAdd(cell_quadrics[c0], Q);
1067 }
1068 else
1069 {
1070 quadricAdd(cell_quadrics[c0], Q);
1071 quadricAdd(cell_quadrics[c1], Q);
1072 quadricAdd(cell_quadrics[c2], Q);
1073 }
1074 }
1075 }
1076
fillCellQuadrics(Quadric * cell_quadrics,const Vector3 * vertex_positions,size_t vertex_count,const unsigned int * vertex_cells)1077 static void fillCellQuadrics(Quadric* cell_quadrics, const Vector3* vertex_positions, size_t vertex_count, const unsigned int* vertex_cells)
1078 {
1079 for (size_t i = 0; i < vertex_count; ++i)
1080 {
1081 unsigned int c = vertex_cells[i];
1082 const Vector3& v = vertex_positions[i];
1083
1084 Quadric Q;
1085 quadricFromPoint(Q, v.x, v.y, v.z, 1.f);
1086
1087 quadricAdd(cell_quadrics[c], Q);
1088 }
1089 }
1090
fillCellRemap(unsigned int * cell_remap,float * cell_errors,size_t cell_count,const unsigned int * vertex_cells,const Quadric * cell_quadrics,const Vector3 * vertex_positions,size_t vertex_count)1091 static void fillCellRemap(unsigned int* cell_remap, float* cell_errors, size_t cell_count, const unsigned int* vertex_cells, const Quadric* cell_quadrics, const Vector3* vertex_positions, size_t vertex_count)
1092 {
1093 memset(cell_remap, -1, cell_count * sizeof(unsigned int));
1094
1095 for (size_t i = 0; i < vertex_count; ++i)
1096 {
1097 unsigned int cell = vertex_cells[i];
1098 float error = quadricError(cell_quadrics[cell], vertex_positions[i]);
1099
1100 if (cell_remap[cell] == ~0u || cell_errors[cell] > error)
1101 {
1102 cell_remap[cell] = unsigned(i);
1103 cell_errors[cell] = error;
1104 }
1105 }
1106 }
1107
filterTriangles(unsigned int * destination,unsigned int * tritable,size_t tritable_size,const unsigned int * indices,size_t index_count,const unsigned int * vertex_cells,const unsigned int * cell_remap)1108 static size_t filterTriangles(unsigned int* destination, unsigned int* tritable, size_t tritable_size, const unsigned int* indices, size_t index_count, const unsigned int* vertex_cells, const unsigned int* cell_remap)
1109 {
1110 TriangleHasher hasher = {destination};
1111
1112 memset(tritable, -1, tritable_size * sizeof(unsigned int));
1113
1114 size_t result = 0;
1115
1116 for (size_t i = 0; i < index_count; i += 3)
1117 {
1118 unsigned int c0 = vertex_cells[indices[i + 0]];
1119 unsigned int c1 = vertex_cells[indices[i + 1]];
1120 unsigned int c2 = vertex_cells[indices[i + 2]];
1121
1122 if (c0 != c1 && c0 != c2 && c1 != c2)
1123 {
1124 unsigned int a = cell_remap[c0];
1125 unsigned int b = cell_remap[c1];
1126 unsigned int c = cell_remap[c2];
1127
1128 if (b < a && b < c)
1129 {
1130 unsigned int t = a;
1131 a = b, b = c, c = t;
1132 }
1133 else if (c < a && c < b)
1134 {
1135 unsigned int t = c;
1136 c = b, b = a, a = t;
1137 }
1138
1139 destination[result * 3 + 0] = a;
1140 destination[result * 3 + 1] = b;
1141 destination[result * 3 + 2] = c;
1142
1143 unsigned int* entry = hashLookup2(tritable, tritable_size, hasher, unsigned(result), ~0u);
1144
1145 if (*entry == ~0u)
1146 *entry = unsigned(result++);
1147 }
1148 }
1149
1150 return result * 3;
1151 }
1152
interpolate(float y,float x0,float y0,float x1,float y1,float x2,float y2)1153 static float interpolate(float y, float x0, float y0, float x1, float y1, float x2, float y2)
1154 {
1155 // three point interpolation from "revenge of interpolation search" paper
1156 float num = (y1 - y) * (x1 - x2) * (x1 - x0) * (y2 - y0);
1157 float den = (y2 - y) * (x1 - x2) * (y0 - y1) + (y0 - y) * (x1 - x0) * (y1 - y2);
1158 return x1 + num / den;
1159 }
1160
1161 } // namespace meshopt
1162
1163 #if TRACE
1164 unsigned char* meshopt_simplifyDebugKind = 0;
1165 unsigned int* meshopt_simplifyDebugLoop = 0;
1166 #endif
1167
meshopt_simplify(unsigned int * destination,const unsigned int * indices,size_t index_count,const float * vertex_positions_data,size_t vertex_count,size_t vertex_positions_stride,size_t target_index_count,float target_error)1168 size_t meshopt_simplify(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride, size_t target_index_count, float target_error)
1169 {
1170 using namespace meshopt;
1171
1172 assert(index_count % 3 == 0);
1173 assert(vertex_positions_stride > 0 && vertex_positions_stride <= 256);
1174 assert(vertex_positions_stride % sizeof(float) == 0);
1175 assert(target_index_count <= index_count);
1176
1177 meshopt_Allocator allocator;
1178
1179 unsigned int* result = destination;
1180
1181 // build adjacency information
1182 EdgeAdjacency adjacency = {};
1183 buildEdgeAdjacency(adjacency, indices, index_count, vertex_count, allocator);
1184
1185 // build position remap that maps each vertex to the one with identical position
1186 unsigned int* remap = allocator.allocate<unsigned int>(vertex_count);
1187 unsigned int* wedge = allocator.allocate<unsigned int>(vertex_count);
1188 buildPositionRemap(remap, wedge, vertex_positions_data, vertex_count, vertex_positions_stride, allocator);
1189
1190 // classify vertices; vertex kind determines collapse rules, see kCanCollapse
1191 unsigned char* vertex_kind = allocator.allocate<unsigned char>(vertex_count);
1192 unsigned int* loop = allocator.allocate<unsigned int>(vertex_count);
1193 classifyVertices(vertex_kind, loop, vertex_count, adjacency, remap, wedge);
1194
1195 #if TRACE
1196 size_t unique_positions = 0;
1197 for (size_t i = 0; i < vertex_count; ++i)
1198 unique_positions += remap[i] == i;
1199
1200 printf("position remap: %d vertices => %d positions\n", int(vertex_count), int(unique_positions));
1201
1202 size_t kinds[Kind_Count] = {};
1203 for (size_t i = 0; i < vertex_count; ++i)
1204 kinds[vertex_kind[i]] += remap[i] == i;
1205
1206 printf("kinds: manifold %d, border %d, seam %d, complex %d, locked %d\n",
1207 int(kinds[Kind_Manifold]), int(kinds[Kind_Border]), int(kinds[Kind_Seam]), int(kinds[Kind_Complex]), int(kinds[Kind_Locked]));
1208 #endif
1209
1210 Vector3* vertex_positions = allocator.allocate<Vector3>(vertex_count);
1211 rescalePositions(vertex_positions, vertex_positions_data, vertex_count, vertex_positions_stride);
1212
1213 Quadric* vertex_quadrics = allocator.allocate<Quadric>(vertex_count);
1214 memset(vertex_quadrics, 0, vertex_count * sizeof(Quadric));
1215
1216 fillFaceQuadrics(vertex_quadrics, indices, index_count, vertex_positions, remap);
1217 fillEdgeQuadrics(vertex_quadrics, indices, index_count, vertex_positions, remap, vertex_kind, loop);
1218
1219 if (result != indices)
1220 memcpy(result, indices, index_count * sizeof(unsigned int));
1221
1222 #if TRACE
1223 size_t pass_count = 0;
1224 float worst_error = 0;
1225 #endif
1226
1227 Collapse* edge_collapses = allocator.allocate<Collapse>(index_count);
1228 unsigned int* collapse_order = allocator.allocate<unsigned int>(index_count);
1229 unsigned int* collapse_remap = allocator.allocate<unsigned int>(vertex_count);
1230 unsigned char* collapse_locked = allocator.allocate<unsigned char>(vertex_count);
1231
1232 size_t result_count = index_count;
1233
1234 // target_error input is linear; we need to adjust it to match quadricError units
1235 float error_limit = target_error * target_error;
1236
1237 while (result_count > target_index_count)
1238 {
1239 size_t edge_collapse_count = pickEdgeCollapses(edge_collapses, result, result_count, remap, vertex_kind, loop);
1240
1241 // no edges can be collapsed any more due to topology restrictions
1242 if (edge_collapse_count == 0)
1243 break;
1244
1245 rankEdgeCollapses(edge_collapses, edge_collapse_count, vertex_positions, vertex_quadrics, remap);
1246
1247 #if TRACE > 1
1248 dumpEdgeCollapses(edge_collapses, edge_collapse_count, vertex_kind);
1249 #endif
1250
1251 sortEdgeCollapses(collapse_order, edge_collapses, edge_collapse_count);
1252
1253 // most collapses remove 2 triangles; use this to establish a bound on the pass in terms of error limit
1254 // note that edge_collapse_goal is an estimate; triangle_collapse_goal will be used to actually limit collapses
1255 size_t triangle_collapse_goal = (result_count - target_index_count) / 3;
1256 size_t edge_collapse_goal = triangle_collapse_goal / 2;
1257
1258 // we limit the error in each pass based on the error of optimal last collapse; since many collapses will be locked
1259 // as they will share vertices with other successfull collapses, we need to increase the acceptable error by this factor
1260 const float kPassErrorBound = 1.5f;
1261
1262 float error_goal = edge_collapse_goal < edge_collapse_count ? edge_collapses[collapse_order[edge_collapse_goal]].error * kPassErrorBound : FLT_MAX;
1263
1264 for (size_t i = 0; i < vertex_count; ++i)
1265 collapse_remap[i] = unsigned(i);
1266
1267 memset(collapse_locked, 0, vertex_count);
1268
1269 size_t collapses = performEdgeCollapses(collapse_remap, collapse_locked, vertex_quadrics, edge_collapses, edge_collapse_count, collapse_order, remap, wedge, vertex_kind, triangle_collapse_goal, error_goal, error_limit);
1270
1271 // no edges can be collapsed any more due to hitting the error limit or triangle collapse limit
1272 if (collapses == 0)
1273 break;
1274
1275 remapEdgeLoops(loop, vertex_count, collapse_remap);
1276
1277 size_t new_count = remapIndexBuffer(result, result_count, collapse_remap);
1278 assert(new_count < result_count);
1279
1280 #if TRACE
1281 float pass_error = 0.f;
1282 for (size_t i = 0; i < edge_collapse_count; ++i)
1283 {
1284 Collapse& c = edge_collapses[collapse_order[i]];
1285
1286 if (collapse_remap[c.v0] == c.v1)
1287 pass_error = c.error;
1288 }
1289
1290 pass_count++;
1291 worst_error = (worst_error < pass_error) ? pass_error : worst_error;
1292
1293 printf("pass %d: triangles: %d -> %d, collapses: %d/%d (goal: %d), error: %e (limit %e goal %e)\n", int(pass_count), int(result_count / 3), int(new_count / 3), int(collapses), int(edge_collapse_count), int(edge_collapse_goal), pass_error, error_limit, error_goal);
1294 #endif
1295
1296 result_count = new_count;
1297 }
1298
1299 #if TRACE
1300 printf("passes: %d, worst error: %e\n", int(pass_count), worst_error);
1301 #endif
1302
1303 #if TRACE > 1
1304 dumpLockedCollapses(result, result_count, vertex_kind);
1305 #endif
1306
1307 #if TRACE
1308 if (meshopt_simplifyDebugKind)
1309 memcpy(meshopt_simplifyDebugKind, vertex_kind, vertex_count);
1310
1311 if (meshopt_simplifyDebugLoop)
1312 memcpy(meshopt_simplifyDebugLoop, loop, vertex_count * sizeof(unsigned int));
1313 #endif
1314
1315 return result_count;
1316 }
1317
meshopt_simplifySloppy(unsigned int * destination,const unsigned int * indices,size_t index_count,const float * vertex_positions_data,size_t vertex_count,size_t vertex_positions_stride,size_t target_index_count)1318 size_t meshopt_simplifySloppy(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride, size_t target_index_count)
1319 {
1320 using namespace meshopt;
1321
1322 assert(index_count % 3 == 0);
1323 assert(vertex_positions_stride > 0 && vertex_positions_stride <= 256);
1324 assert(vertex_positions_stride % sizeof(float) == 0);
1325 assert(target_index_count <= index_count);
1326
1327 // we expect to get ~2 triangles/vertex in the output
1328 size_t target_cell_count = target_index_count / 6;
1329
1330 if (target_cell_count == 0)
1331 return 0;
1332
1333 meshopt_Allocator allocator;
1334
1335 Vector3* vertex_positions = allocator.allocate<Vector3>(vertex_count);
1336 rescalePositions(vertex_positions, vertex_positions_data, vertex_count, vertex_positions_stride);
1337
1338 // find the optimal grid size using guided binary search
1339 #if TRACE
1340 printf("source: %d vertices, %d triangles\n", int(vertex_count), int(index_count / 3));
1341 printf("target: %d cells, %d triangles\n", int(target_cell_count), int(target_index_count / 3));
1342 #endif
1343
1344 unsigned int* vertex_ids = allocator.allocate<unsigned int>(vertex_count);
1345
1346 const int kInterpolationPasses = 5;
1347
1348 // invariant: # of triangles in min_grid <= target_count
1349 int min_grid = 0;
1350 int max_grid = 1025;
1351 size_t min_triangles = 0;
1352 size_t max_triangles = index_count / 3;
1353
1354 // instead of starting in the middle, let's guess as to what the answer might be! triangle count usually grows as a square of grid size...
1355 int next_grid_size = int(sqrtf(float(target_cell_count)) + 0.5f);
1356
1357 for (int pass = 0; pass < 10 + kInterpolationPasses; ++pass)
1358 {
1359 assert(min_triangles < target_index_count / 3);
1360 assert(max_grid - min_grid > 1);
1361
1362 // we clamp the prediction of the grid size to make sure that the search converges
1363 int grid_size = next_grid_size;
1364 grid_size = (grid_size <= min_grid) ? min_grid + 1 : (grid_size >= max_grid) ? max_grid - 1 : grid_size;
1365
1366 computeVertexIds(vertex_ids, vertex_positions, vertex_count, grid_size);
1367 size_t triangles = countTriangles(vertex_ids, indices, index_count);
1368
1369 #if TRACE
1370 printf("pass %d (%s): grid size %d, triangles %d, %s\n",
1371 pass, (pass == 0) ? "guess" : (pass <= kInterpolationPasses) ? "lerp" : "binary",
1372 grid_size, int(triangles),
1373 (triangles <= target_index_count / 3) ? "under" : "over");
1374 #endif
1375
1376 float tip = interpolate(float(target_index_count / 3), float(min_grid), float(min_triangles), float(grid_size), float(triangles), float(max_grid), float(max_triangles));
1377
1378 if (triangles <= target_index_count / 3)
1379 {
1380 min_grid = grid_size;
1381 min_triangles = triangles;
1382 }
1383 else
1384 {
1385 max_grid = grid_size;
1386 max_triangles = triangles;
1387 }
1388
1389 if (triangles == target_index_count / 3 || max_grid - min_grid <= 1)
1390 break;
1391
1392 // we start by using interpolation search - it usually converges faster
1393 // however, interpolation search has a worst case of O(N) so we switch to binary search after a few iterations which converges in O(logN)
1394 next_grid_size = (pass < kInterpolationPasses) ? int(tip + 0.5f) : (min_grid + max_grid) / 2;
1395 }
1396
1397 if (min_triangles == 0)
1398 return 0;
1399
1400 // build vertex->cell association by mapping all vertices with the same quantized position to the same cell
1401 size_t table_size = hashBuckets2(vertex_count);
1402 unsigned int* table = allocator.allocate<unsigned int>(table_size);
1403
1404 unsigned int* vertex_cells = allocator.allocate<unsigned int>(vertex_count);
1405
1406 computeVertexIds(vertex_ids, vertex_positions, vertex_count, min_grid);
1407 size_t cell_count = fillVertexCells(table, table_size, vertex_cells, vertex_ids, vertex_count);
1408
1409 // build a quadric for each target cell
1410 Quadric* cell_quadrics = allocator.allocate<Quadric>(cell_count);
1411 memset(cell_quadrics, 0, cell_count * sizeof(Quadric));
1412
1413 fillCellQuadrics(cell_quadrics, indices, index_count, vertex_positions, vertex_cells);
1414
1415 // for each target cell, find the vertex with the minimal error
1416 unsigned int* cell_remap = allocator.allocate<unsigned int>(cell_count);
1417 float* cell_errors = allocator.allocate<float>(cell_count);
1418
1419 fillCellRemap(cell_remap, cell_errors, cell_count, vertex_cells, cell_quadrics, vertex_positions, vertex_count);
1420
1421 // collapse triangles!
1422 // note that we need to filter out triangles that we've already output because we very frequently generate redundant triangles between cells :(
1423 size_t tritable_size = hashBuckets2(min_triangles);
1424 unsigned int* tritable = allocator.allocate<unsigned int>(tritable_size);
1425
1426 size_t write = filterTriangles(destination, tritable, tritable_size, indices, index_count, vertex_cells, cell_remap);
1427 assert(write <= target_index_count);
1428
1429 #if TRACE
1430 printf("result: %d cells, %d triangles (%d unfiltered)\n", int(cell_count), int(write / 3), int(min_triangles));
1431 #endif
1432
1433 return write;
1434 }
1435
meshopt_simplifyPoints(unsigned int * destination,const float * vertex_positions_data,size_t vertex_count,size_t vertex_positions_stride,size_t target_vertex_count)1436 size_t meshopt_simplifyPoints(unsigned int* destination, const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride, size_t target_vertex_count)
1437 {
1438 using namespace meshopt;
1439
1440 assert(vertex_positions_stride > 0 && vertex_positions_stride <= 256);
1441 assert(vertex_positions_stride % sizeof(float) == 0);
1442 assert(target_vertex_count <= vertex_count);
1443
1444 size_t target_cell_count = target_vertex_count;
1445
1446 if (target_cell_count == 0)
1447 return 0;
1448
1449 meshopt_Allocator allocator;
1450
1451 Vector3* vertex_positions = allocator.allocate<Vector3>(vertex_count);
1452 rescalePositions(vertex_positions, vertex_positions_data, vertex_count, vertex_positions_stride);
1453
1454 // find the optimal grid size using guided binary search
1455 #if TRACE
1456 printf("source: %d vertices\n", int(vertex_count));
1457 printf("target: %d cells\n", int(target_cell_count));
1458 #endif
1459
1460 unsigned int* vertex_ids = allocator.allocate<unsigned int>(vertex_count);
1461
1462 size_t table_size = hashBuckets2(vertex_count);
1463 unsigned int* table = allocator.allocate<unsigned int>(table_size);
1464
1465 const int kInterpolationPasses = 5;
1466
1467 // invariant: # of vertices in min_grid <= target_count
1468 int min_grid = 0;
1469 int max_grid = 1025;
1470 size_t min_vertices = 0;
1471 size_t max_vertices = vertex_count;
1472
1473 // instead of starting in the middle, let's guess as to what the answer might be! triangle count usually grows as a square of grid size...
1474 int next_grid_size = int(sqrtf(float(target_cell_count)) + 0.5f);
1475
1476 for (int pass = 0; pass < 10 + kInterpolationPasses; ++pass)
1477 {
1478 assert(min_vertices < target_vertex_count);
1479 assert(max_grid - min_grid > 1);
1480
1481 // we clamp the prediction of the grid size to make sure that the search converges
1482 int grid_size = next_grid_size;
1483 grid_size = (grid_size <= min_grid) ? min_grid + 1 : (grid_size >= max_grid) ? max_grid - 1 : grid_size;
1484
1485 computeVertexIds(vertex_ids, vertex_positions, vertex_count, grid_size);
1486 size_t vertices = countVertexCells(table, table_size, vertex_ids, vertex_count);
1487
1488 #if TRACE
1489 printf("pass %d (%s): grid size %d, vertices %d, %s\n",
1490 pass, (pass == 0) ? "guess" : (pass <= kInterpolationPasses) ? "lerp" : "binary",
1491 grid_size, int(vertices),
1492 (vertices <= target_vertex_count) ? "under" : "over");
1493 #endif
1494
1495 float tip = interpolate(float(target_vertex_count), float(min_grid), float(min_vertices), float(grid_size), float(vertices), float(max_grid), float(max_vertices));
1496
1497 if (vertices <= target_vertex_count)
1498 {
1499 min_grid = grid_size;
1500 min_vertices = vertices;
1501 }
1502 else
1503 {
1504 max_grid = grid_size;
1505 max_vertices = vertices;
1506 }
1507
1508 if (vertices == target_vertex_count || max_grid - min_grid <= 1)
1509 break;
1510
1511 // we start by using interpolation search - it usually converges faster
1512 // however, interpolation search has a worst case of O(N) so we switch to binary search after a few iterations which converges in O(logN)
1513 next_grid_size = (pass < kInterpolationPasses) ? int(tip + 0.5f) : (min_grid + max_grid) / 2;
1514 }
1515
1516 if (min_vertices == 0)
1517 return 0;
1518
1519 // build vertex->cell association by mapping all vertices with the same quantized position to the same cell
1520 unsigned int* vertex_cells = allocator.allocate<unsigned int>(vertex_count);
1521
1522 computeVertexIds(vertex_ids, vertex_positions, vertex_count, min_grid);
1523 size_t cell_count = fillVertexCells(table, table_size, vertex_cells, vertex_ids, vertex_count);
1524
1525 // build a quadric for each target cell
1526 Quadric* cell_quadrics = allocator.allocate<Quadric>(cell_count);
1527 memset(cell_quadrics, 0, cell_count * sizeof(Quadric));
1528
1529 fillCellQuadrics(cell_quadrics, vertex_positions, vertex_count, vertex_cells);
1530
1531 // for each target cell, find the vertex with the minimal error
1532 unsigned int* cell_remap = allocator.allocate<unsigned int>(cell_count);
1533 float* cell_errors = allocator.allocate<float>(cell_count);
1534
1535 fillCellRemap(cell_remap, cell_errors, cell_count, vertex_cells, cell_quadrics, vertex_positions, vertex_count);
1536
1537 // copy results to the output
1538 assert(cell_count <= target_vertex_count);
1539 memcpy(destination, cell_remap, sizeof(unsigned int) * cell_count);
1540
1541 #if TRACE
1542 printf("result: %d cells\n", int(cell_count));
1543 #endif
1544
1545 return cell_count;
1546 }
1547