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#if TRACE
18#define TRACESTATS(i) stats[i]++;
19#else
20#define TRACESTATS(i) (void)0
21#endif
22
23#define ATTRIBUTES 3
24
25// This work is based on:
26// Michael Garland and Paul S. Heckbert. Surface simplification using quadric error metrics. 1997
27// Michael Garland. Quadric-based polygonal surface simplification. 1999
28// Peter Lindstrom. Out-of-Core Simplification of Large Polygonal Models. 2000
29// Matthias Teschner, Bruno Heidelberger, Matthias Mueller, Danat Pomeranets, Markus Gross. Optimized Spatial Hashing for Collision Detection of Deformable Objects. 2003
30// Peter Van Sandt, Yannis Chronis, Jignesh M. Patel. Efficiently Searching In-Memory Sorted Arrays: Revenge of the Interpolation Search? 2019
31namespace meshopt
32{
33
34struct EdgeAdjacency
35{
36 struct Edge
37 {
38 unsigned int next;
39 unsigned int prev;
40 };
41
42 unsigned int* counts;
43 unsigned int* offsets;
44 Edge* data;
45};
46
47static void prepareEdgeAdjacency(EdgeAdjacency& adjacency, size_t index_count, size_t vertex_count, meshopt_Allocator& allocator)
48{
49 adjacency.counts = allocator.allocate<unsigned int>(vertex_count);
50 adjacency.offsets = allocator.allocate<unsigned int>(vertex_count);
51 adjacency.data = allocator.allocate<EdgeAdjacency::Edge>(index_count);
52}
53
54static void updateEdgeAdjacency(EdgeAdjacency& adjacency, const unsigned int* indices, size_t index_count, size_t vertex_count, const unsigned int* remap)
55{
56 size_t face_count = index_count / 3;
57
58 // fill edge counts
59 memset(adjacency.counts, 0, vertex_count * sizeof(unsigned int));
60
61 for (size_t i = 0; i < index_count; ++i)
62 {
63 unsigned int v = remap ? remap[indices[i]] : indices[i];
64 assert(v < vertex_count);
65
66 adjacency.counts[v]++;
67 }
68
69 // fill offset table
70 unsigned int offset = 0;
71
72 for (size_t i = 0; i < vertex_count; ++i)
73 {
74 adjacency.offsets[i] = offset;
75 offset += adjacency.counts[i];
76 }
77
78 assert(offset == index_count);
79
80 // fill edge data
81 for (size_t i = 0; i < face_count; ++i)
82 {
83 unsigned int a = indices[i * 3 + 0], b = indices[i * 3 + 1], c = indices[i * 3 + 2];
84
85 if (remap)
86 {
87 a = remap[a];
88 b = remap[b];
89 c = remap[c];
90 }
91
92 adjacency.data[adjacency.offsets[a]].next = b;
93 adjacency.data[adjacency.offsets[a]].prev = c;
94 adjacency.offsets[a]++;
95
96 adjacency.data[adjacency.offsets[b]].next = c;
97 adjacency.data[adjacency.offsets[b]].prev = a;
98 adjacency.offsets[b]++;
99
100 adjacency.data[adjacency.offsets[c]].next = a;
101 adjacency.data[adjacency.offsets[c]].prev = b;
102 adjacency.offsets[c]++;
103 }
104
105 // fix offsets that have been disturbed by the previous pass
106 for (size_t i = 0; i < vertex_count; ++i)
107 {
108 assert(adjacency.offsets[i] >= adjacency.counts[i]);
109
110 adjacency.offsets[i] -= adjacency.counts[i];
111 }
112}
113
114struct PositionHasher
115{
116 const float* vertex_positions;
117 size_t vertex_stride_float;
118
119 size_t hash(unsigned int index) const
120 {
121 const unsigned int* key = reinterpret_cast<const unsigned int*>(vertex_positions + index * vertex_stride_float);
122
123 // scramble bits to make sure that integer coordinates have entropy in lower bits
124 unsigned int x = key[0] ^ (key[0] >> 17);
125 unsigned int y = key[1] ^ (key[1] >> 17);
126 unsigned int z = key[2] ^ (key[2] >> 17);
127
128 // Optimized Spatial Hashing for Collision Detection of Deformable Objects
129 return (x * 73856093) ^ (y * 19349663) ^ (z * 83492791);
130 }
131
132 bool equal(unsigned int lhs, unsigned int rhs) const
133 {
134 return memcmp(vertex_positions + lhs * vertex_stride_float, vertex_positions + rhs * vertex_stride_float, sizeof(float) * 3) == 0;
135 }
136};
137
138static size_t hashBuckets2(size_t count)
139{
140 size_t buckets = 1;
141 while (buckets < count + count / 4)
142 buckets *= 2;
143
144 return buckets;
145}
146
147template <typename T, typename Hash>
148static T* hashLookup2(T* table, size_t buckets, const Hash& hash, const T& key, const T& empty)
149{
150 assert(buckets > 0);
151 assert((buckets & (buckets - 1)) == 0);
152
153 size_t hashmod = buckets - 1;
154 size_t bucket = hash.hash(key) & hashmod;
155
156 for (size_t probe = 0; probe <= hashmod; ++probe)
157 {
158 T& item = table[bucket];
159
160 if (item == empty)
161 return &item;
162
163 if (hash.equal(item, key))
164 return &item;
165
166 // hash collision, quadratic probing
167 bucket = (bucket + probe + 1) & hashmod;
168 }
169
170 assert(false && "Hash table is full"); // unreachable
171 return 0;
172}
173
174static 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)
175{
176 PositionHasher hasher = {vertex_positions_data, vertex_positions_stride / sizeof(float)};
177
178 size_t table_size = hashBuckets2(vertex_count);
179 unsigned int* table = allocator.allocate<unsigned int>(table_size);
180 memset(table, -1, table_size * sizeof(unsigned int));
181
182 // build forward remap: for each vertex, which other (canonical) vertex does it map to?
183 // we use position equivalence for this, and remap vertices to other existing vertices
184 for (size_t i = 0; i < vertex_count; ++i)
185 {
186 unsigned int index = unsigned(i);
187 unsigned int* entry = hashLookup2(table, table_size, hasher, index, ~0u);
188
189 if (*entry == ~0u)
190 *entry = index;
191
192 remap[index] = *entry;
193 }
194
195 // build wedge table: for each vertex, which other vertex is the next wedge that also maps to the same vertex?
196 // entries in table form a (cyclic) wedge loop per vertex; for manifold vertices, wedge[i] == remap[i] == i
197 for (size_t i = 0; i < vertex_count; ++i)
198 wedge[i] = unsigned(i);
199
200 for (size_t i = 0; i < vertex_count; ++i)
201 if (remap[i] != i)
202 {
203 unsigned int r = remap[i];
204
205 wedge[i] = wedge[r];
206 wedge[r] = unsigned(i);
207 }
208}
209
210enum VertexKind
211{
212 Kind_Manifold, // not on an attribute seam, not on any boundary
213 Kind_Border, // not on an attribute seam, has exactly two open edges
214 Kind_Seam, // on an attribute seam with exactly two attribute seam edges
215 Kind_Complex, // none of the above; these vertices can move as long as all wedges move to the target vertex
216 Kind_Locked, // none of the above; these vertices can't move
217
218 Kind_Count
219};
220
221// manifold vertices can collapse onto anything
222// border/seam vertices can only be collapsed onto border/seam respectively
223// complex vertices can collapse onto complex/locked
224// a rule of thumb is that collapsing kind A into kind B preserves the kind B in the target vertex
225// for example, while we could collapse Complex into Manifold, this would mean the target vertex isn't Manifold anymore
226const unsigned char kCanCollapse[Kind_Count][Kind_Count] = {
227 {1, 1, 1, 1, 1},
228 {0, 1, 0, 0, 0},
229 {0, 0, 1, 0, 0},
230 {0, 0, 0, 1, 1},
231 {0, 0, 0, 0, 0},
232};
233
234// if a vertex is manifold or seam, adjoining edges are guaranteed to have an opposite edge
235// note that for seam edges, the opposite edge isn't present in the attribute-based topology
236// but is present if you consider a position-only mesh variant
237const unsigned char kHasOpposite[Kind_Count][Kind_Count] = {
238 {1, 1, 1, 0, 1},
239 {1, 0, 1, 0, 0},
240 {1, 1, 1, 0, 1},
241 {0, 0, 0, 0, 0},
242 {1, 0, 1, 0, 0},
243};
244
245static bool hasEdge(const EdgeAdjacency& adjacency, unsigned int a, unsigned int b)
246{
247 unsigned int count = adjacency.counts[a];
248 const EdgeAdjacency::Edge* edges = adjacency.data + adjacency.offsets[a];
249
250 for (size_t i = 0; i < count; ++i)
251 if (edges[i].next == b)
252 return true;
253
254 return false;
255}
256
257static void classifyVertices(unsigned char* result, unsigned int* loop, unsigned int* loopback, size_t vertex_count, const EdgeAdjacency& adjacency, const unsigned int* remap, const unsigned int* wedge, unsigned int options)
258{
259 memset(loop, -1, vertex_count * sizeof(unsigned int));
260 memset(loopback, -1, vertex_count * sizeof(unsigned int));
261
262 // incoming & outgoing open edges: ~0u if no open edges, i if there are more than 1
263 // note that this is the same data as required in loop[] arrays; loop[] data is only valid for border/seam
264 // but here it's okay to fill the data out for other types of vertices as well
265 unsigned int* openinc = loopback;
266 unsigned int* openout = loop;
267
268 for (size_t i = 0; i < vertex_count; ++i)
269 {
270 unsigned int vertex = unsigned(i);
271
272 unsigned int count = adjacency.counts[vertex];
273 const EdgeAdjacency::Edge* edges = adjacency.data + adjacency.offsets[vertex];
274
275 for (size_t j = 0; j < count; ++j)
276 {
277 unsigned int target = edges[j].next;
278
279 if (target == vertex)
280 {
281 // degenerate triangles have two distinct edges instead of three, and the self edge
282 // is bi-directional by definition; this can break border/seam classification by "closing"
283 // the open edge from another triangle and falsely marking the vertex as manifold
284 // instead we mark the vertex as having >1 open edges which turns it into locked/complex
285 openinc[vertex] = openout[vertex] = vertex;
286 }
287 else if (!hasEdge(adjacency, target, vertex))
288 {
289 openinc[target] = (openinc[target] == ~0u) ? vertex : target;
290 openout[vertex] = (openout[vertex] == ~0u) ? target : vertex;
291 }
292 }
293 }
294
295#if TRACE
296 size_t stats[4] = {};
297#endif
298
299 for (size_t i = 0; i < vertex_count; ++i)
300 {
301 if (remap[i] == i)
302 {
303 if (wedge[i] == i)
304 {
305 // no attribute seam, need to check if it's manifold
306 unsigned int openi = openinc[i], openo = openout[i];
307
308 // note: we classify any vertices with no open edges as manifold
309 // this is technically incorrect - if 4 triangles share an edge, we'll classify vertices as manifold
310 // it's unclear if this is a problem in practice
311 if (openi == ~0u && openo == ~0u)
312 {
313 result[i] = Kind_Manifold;
314 }
315 else if (openi != i && openo != i)
316 {
317 result[i] = Kind_Border;
318 }
319 else
320 {
321 result[i] = Kind_Locked;
322 TRACESTATS(0);
323 }
324 }
325 else if (wedge[wedge[i]] == i)
326 {
327 // attribute seam; need to distinguish between Seam and Locked
328 unsigned int w = wedge[i];
329 unsigned int openiv = openinc[i], openov = openout[i];
330 unsigned int openiw = openinc[w], openow = openout[w];
331
332 // seam should have one open half-edge for each vertex, and the edges need to "connect" - point to the same vertex post-remap
333 if (openiv != ~0u && openiv != i && openov != ~0u && openov != i &&
334 openiw != ~0u && openiw != w && openow != ~0u && openow != w)
335 {
336 if (remap[openiv] == remap[openow] && remap[openov] == remap[openiw])
337 {
338 result[i] = Kind_Seam;
339 }
340 else
341 {
342 result[i] = Kind_Locked;
343 TRACESTATS(1);
344 }
345 }
346 else
347 {
348 result[i] = Kind_Locked;
349 TRACESTATS(2);
350 }
351 }
352 else
353 {
354 // more than one vertex maps to this one; we don't have classification available
355 result[i] = Kind_Locked;
356 TRACESTATS(3);
357 }
358 }
359 else
360 {
361 assert(remap[i] < i);
362
363 result[i] = result[remap[i]];
364 }
365 }
366
367 if (options & meshopt_SimplifyLockBorder)
368 for (size_t i = 0; i < vertex_count; ++i)
369 if (result[i] == Kind_Border)
370 result[i] = Kind_Locked;
371
372#if TRACE
373 printf("locked: many open edges %d, disconnected seam %d, many seam edges %d, many wedges %d\n",
374 int(stats[0]), int(stats[1]), int(stats[2]), int(stats[3]));
375#endif
376}
377
378struct Vector3
379{
380 float x, y, z;
381
382#if ATTRIBUTES
383 float a[ATTRIBUTES];
384#endif
385};
386
387static float rescalePositions(Vector3* result, const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride)
388{
389 size_t vertex_stride_float = vertex_positions_stride / sizeof(float);
390
391 float minv[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
392 float maxv[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX};
393
394 for (size_t i = 0; i < vertex_count; ++i)
395 {
396 const float* v = vertex_positions_data + i * vertex_stride_float;
397
398 if (result)
399 {
400 result[i].x = v[0];
401 result[i].y = v[1];
402 result[i].z = v[2];
403 }
404
405 for (int j = 0; j < 3; ++j)
406 {
407 float vj = v[j];
408
409 minv[j] = minv[j] > vj ? vj : minv[j];
410 maxv[j] = maxv[j] < vj ? vj : maxv[j];
411 }
412 }
413
414 float extent = 0.f;
415
416 extent = (maxv[0] - minv[0]) < extent ? extent : (maxv[0] - minv[0]);
417 extent = (maxv[1] - minv[1]) < extent ? extent : (maxv[1] - minv[1]);
418 extent = (maxv[2] - minv[2]) < extent ? extent : (maxv[2] - minv[2]);
419
420 if (result)
421 {
422 float scale = extent == 0 ? 0.f : 1.f / extent;
423
424 for (size_t i = 0; i < vertex_count; ++i)
425 {
426 result[i].x = (result[i].x - minv[0]) * scale;
427 result[i].y = (result[i].y - minv[1]) * scale;
428 result[i].z = (result[i].z - minv[2]) * scale;
429 }
430 }
431
432 return extent;
433}
434
435struct Quadric
436{
437 float a00, a11, a22;
438 float a10, a20, a21;
439 float b0, b1, b2, c;
440 float w;
441
442#if ATTRIBUTES
443 float gx[ATTRIBUTES];
444 float gy[ATTRIBUTES];
445 float gz[ATTRIBUTES];
446 float gw[ATTRIBUTES];
447#endif
448};
449
450struct Collapse
451{
452 unsigned int v0;
453 unsigned int v1;
454
455 union
456 {
457 unsigned int bidi;
458 float error;
459 unsigned int errorui;
460 };
461 float distance_error;
462};
463
464static float normalize(Vector3& v)
465{
466 float length = sqrtf(v.x * v.x + v.y * v.y + v.z * v.z);
467
468 if (length > 0)
469 {
470 v.x /= length;
471 v.y /= length;
472 v.z /= length;
473 }
474
475 return length;
476}
477
478static void quadricAdd(Quadric& Q, const Quadric& R)
479{
480 Q.a00 += R.a00;
481 Q.a11 += R.a11;
482 Q.a22 += R.a22;
483 Q.a10 += R.a10;
484 Q.a20 += R.a20;
485 Q.a21 += R.a21;
486 Q.b0 += R.b0;
487 Q.b1 += R.b1;
488 Q.b2 += R.b2;
489 Q.c += R.c;
490 Q.w += R.w;
491
492#if ATTRIBUTES
493 for (int k = 0; k < ATTRIBUTES; ++k)
494 {
495 Q.gx[k] += R.gx[k];
496 Q.gy[k] += R.gy[k];
497 Q.gz[k] += R.gz[k];
498 Q.gw[k] += R.gw[k];
499 }
500#endif
501}
502
503static float quadricError(const Quadric& Q, const Vector3& v)
504{
505 float rx = Q.b0;
506 float ry = Q.b1;
507 float rz = Q.b2;
508
509 rx += Q.a10 * v.y;
510 ry += Q.a21 * v.z;
511 rz += Q.a20 * v.x;
512
513 rx *= 2;
514 ry *= 2;
515 rz *= 2;
516
517 rx += Q.a00 * v.x;
518 ry += Q.a11 * v.y;
519 rz += Q.a22 * v.z;
520
521 float r = Q.c;
522 r += rx * v.x;
523 r += ry * v.y;
524 r += rz * v.z;
525
526#if ATTRIBUTES
527 // see quadricUpdateAttributes for general derivation; here we need to add the parts of (eval(pos) - attr)^2 that depend on attr
528 for (int k = 0; k < ATTRIBUTES; ++k)
529 {
530 float a = v.a[k];
531
532 r += a * a * Q.w;
533 r -= 2 * a * (v.x * Q.gx[k] + v.y * Q.gy[k] + v.z * Q.gz[k] + Q.gw[k]);
534 }
535#endif
536
537 float s = Q.w == 0.f ? 0.f : 1.f / Q.w;
538
539 return fabsf(r) * s;
540}
541
542static float quadricErrorNoAttributes(const Quadric& Q, const Vector3& v)
543{
544 float rx = Q.b0;
545 float ry = Q.b1;
546 float rz = Q.b2;
547
548 rx += Q.a10 * v.y;
549 ry += Q.a21 * v.z;
550 rz += Q.a20 * v.x;
551
552 rx *= 2;
553 ry *= 2;
554 rz *= 2;
555
556 rx += Q.a00 * v.x;
557 ry += Q.a11 * v.y;
558 rz += Q.a22 * v.z;
559
560 float r = Q.c;
561 r += rx * v.x;
562 r += ry * v.y;
563 r += rz * v.z;
564
565 float s = Q.w == 0.f ? 0.f : 1.f / Q.w;
566
567 return fabsf(r) * s;
568}
569
570static void quadricFromPlane(Quadric& Q, float a, float b, float c, float d, float w)
571{
572 float aw = a * w;
573 float bw = b * w;
574 float cw = c * w;
575 float dw = d * w;
576
577 Q.a00 = a * aw;
578 Q.a11 = b * bw;
579 Q.a22 = c * cw;
580 Q.a10 = a * bw;
581 Q.a20 = a * cw;
582 Q.a21 = b * cw;
583 Q.b0 = a * dw;
584 Q.b1 = b * dw;
585 Q.b2 = c * dw;
586 Q.c = d * dw;
587 Q.w = w;
588
589#if ATTRIBUTES
590 memset(Q.gx, 0, sizeof(Q.gx));
591 memset(Q.gy, 0, sizeof(Q.gy));
592 memset(Q.gz, 0, sizeof(Q.gz));
593 memset(Q.gw, 0, sizeof(Q.gw));
594#endif
595}
596
597static void quadricFromPoint(Quadric& Q, float x, float y, float z, float w)
598{
599 // we need to encode (x - X) ^ 2 + (y - Y)^2 + (z - Z)^2 into the quadric
600 Q.a00 = w;
601 Q.a11 = w;
602 Q.a22 = w;
603 Q.a10 = 0.f;
604 Q.a20 = 0.f;
605 Q.a21 = 0.f;
606 Q.b0 = -2.f * x * w;
607 Q.b1 = -2.f * y * w;
608 Q.b2 = -2.f * z * w;
609 Q.c = (x * x + y * y + z * z) * w;
610 Q.w = w;
611}
612
613static void quadricFromTriangle(Quadric& Q, const Vector3& p0, const Vector3& p1, const Vector3& p2, float weight)
614{
615 Vector3 p10 = {p1.x - p0.x, p1.y - p0.y, p1.z - p0.z};
616 Vector3 p20 = {p2.x - p0.x, p2.y - p0.y, p2.z - p0.z};
617
618 // normal = cross(p1 - p0, p2 - p0)
619 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};
620 float area = normalize(normal);
621
622 float distance = normal.x * p0.x + normal.y * p0.y + normal.z * p0.z;
623
624 // we use sqrtf(area) so that the error is scaled linearly; this tends to improve silhouettes
625 quadricFromPlane(Q, normal.x, normal.y, normal.z, -distance, sqrtf(area) * weight);
626}
627
628static void quadricFromTriangleEdge(Quadric& Q, const Vector3& p0, const Vector3& p1, const Vector3& p2, float weight)
629{
630 Vector3 p10 = {p1.x - p0.x, p1.y - p0.y, p1.z - p0.z};
631 float length = normalize(p10);
632
633 // p20p = length of projection of p2-p0 onto normalize(p1 - p0)
634 Vector3 p20 = {p2.x - p0.x, p2.y - p0.y, p2.z - p0.z};
635 float p20p = p20.x * p10.x + p20.y * p10.y + p20.z * p10.z;
636
637 // normal = altitude of triangle from point p2 onto edge p1-p0
638 Vector3 normal = {p20.x - p10.x * p20p, p20.y - p10.y * p20p, p20.z - p10.z * p20p};
639 normalize(normal);
640
641 float distance = normal.x * p0.x + normal.y * p0.y + normal.z * p0.z;
642
643 // note: the weight is scaled linearly with edge length; this has to match the triangle weight
644 quadricFromPlane(Q, normal.x, normal.y, normal.z, -distance, length * weight);
645}
646
647#if ATTRIBUTES
648static void quadricUpdateAttributes(Quadric& Q, const Vector3& p0, const Vector3& p1, const Vector3& p2, float w)
649{
650 // for each attribute we want to encode the following function into the quadric:
651 // (eval(pos) - attr)^2
652 // where eval(pos) interpolates attribute across the triangle like so:
653 // eval(pos) = pos.x * gx + pos.y * gy + pos.z * gz + gw
654 // where gx/gy/gz/gw are gradients
655 Vector3 p10 = {p1.x - p0.x, p1.y - p0.y, p1.z - p0.z};
656 Vector3 p20 = {p2.x - p0.x, p2.y - p0.y, p2.z - p0.z};
657
658 // we compute gradients using barycentric coordinates; barycentric coordinates can be computed as follows:
659 // v = (d11 * d20 - d01 * d21) / denom
660 // w = (d00 * d21 - d01 * d20) / denom
661 // u = 1 - v - w
662 // here v0, v1 are triangle edge vectors, v2 is a vector from point to triangle corner, and dij = dot(vi, vj)
663 const Vector3& v0 = p10;
664 const Vector3& v1 = p20;
665 float d00 = v0.x * v0.x + v0.y * v0.y + v0.z * v0.z;
666 float d01 = v0.x * v1.x + v0.y * v1.y + v0.z * v1.z;
667 float d11 = v1.x * v1.x + v1.y * v1.y + v1.z * v1.z;
668 float denom = d00 * d11 - d01 * d01;
669 float denomr = denom == 0 ? 0.f : 1.f / denom;
670
671 // precompute gradient factors
672 // these are derived by directly computing derivative of eval(pos) = a0 * u + a1 * v + a2 * w and factoring out common factors that are shared between attributes
673 float gx1 = (d11 * v0.x - d01 * v1.x) * denomr;
674 float gx2 = (d00 * v1.x - d01 * v0.x) * denomr;
675 float gy1 = (d11 * v0.y - d01 * v1.y) * denomr;
676 float gy2 = (d00 * v1.y - d01 * v0.y) * denomr;
677 float gz1 = (d11 * v0.z - d01 * v1.z) * denomr;
678 float gz2 = (d00 * v1.z - d01 * v0.z) * denomr;
679
680 for (int k = 0; k < ATTRIBUTES; ++k)
681 {
682 float a0 = p0.a[k], a1 = p1.a[k], a2 = p2.a[k];
683
684 // compute gradient of eval(pos) for x/y/z/w
685 // the formulas below are obtained by directly computing derivative of eval(pos) = a0 * u + a1 * v + a2 * w
686 float gx = gx1 * (a1 - a0) + gx2 * (a2 - a0);
687 float gy = gy1 * (a1 - a0) + gy2 * (a2 - a0);
688 float gz = gz1 * (a1 - a0) + gz2 * (a2 - a0);
689 float gw = a0 - p0.x * gx - p0.y * gy - p0.z * gz;
690
691 // quadric encodes (eval(pos)-attr)^2; this means that the resulting expansion needs to compute, for example, pos.x * pos.y * K
692 // since quadrics already encode factors for pos.x * pos.y, we can accumulate almost everything in basic quadric fields
693 Q.a00 += w * (gx * gx);
694 Q.a11 += w * (gy * gy);
695 Q.a22 += w * (gz * gz);
696
697 Q.a10 += w * (gy * gx);
698 Q.a20 += w * (gz * gx);
699 Q.a21 += w * (gz * gy);
700
701 Q.b0 += w * (gx * gw);
702 Q.b1 += w * (gy * gw);
703 Q.b2 += w * (gz * gw);
704
705 Q.c += w * (gw * gw);
706
707 // the only remaining sum components are ones that depend on attr; these will be addded during error evaluation, see quadricError
708 Q.gx[k] = w * gx;
709 Q.gy[k] = w * gy;
710 Q.gz[k] = w * gz;
711 Q.gw[k] = w * gw;
712
713#if TRACE > 2
714 printf("attr%d: %e %e %e\n",
715 k,
716 (gx * p0.x + gy * p0.y + gz * p0.z + gw - a0),
717 (gx * p1.x + gy * p1.y + gz * p1.z + gw - a1),
718 (gx * p2.x + gy * p2.y + gz * p2.z + gw - a2)
719 );
720#endif
721 }
722}
723#endif
724
725static void fillFaceQuadrics(Quadric* vertex_quadrics, Quadric* vertex_no_attrib_quadrics, const unsigned int* indices, size_t index_count, const Vector3* vertex_positions, const unsigned int* remap)
726{
727 for (size_t i = 0; i < index_count; i += 3)
728 {
729 unsigned int i0 = indices[i + 0];
730 unsigned int i1 = indices[i + 1];
731 unsigned int i2 = indices[i + 2];
732
733 Quadric Q;
734 quadricFromTriangle(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], 1.f);
735 quadricAdd(vertex_no_attrib_quadrics[remap[i0]], Q);
736 quadricAdd(vertex_no_attrib_quadrics[remap[i1]], Q);
737 quadricAdd(vertex_no_attrib_quadrics[remap[i2]], Q);
738
739#if ATTRIBUTES
740 quadricUpdateAttributes(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], Q.w);
741#endif
742 quadricAdd(vertex_quadrics[remap[i0]], Q);
743 quadricAdd(vertex_quadrics[remap[i1]], Q);
744 quadricAdd(vertex_quadrics[remap[i2]], Q);
745 }
746}
747
748static void fillEdgeQuadrics(Quadric* vertex_quadrics, Quadric* vertex_no_attrib_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, const unsigned int* loopback)
749{
750 for (size_t i = 0; i < index_count; i += 3)
751 {
752 static const int next[3] = {1, 2, 0};
753
754 for (int e = 0; e < 3; ++e)
755 {
756 unsigned int i0 = indices[i + e];
757 unsigned int i1 = indices[i + next[e]];
758
759 unsigned char k0 = vertex_kind[i0];
760 unsigned char k1 = vertex_kind[i1];
761
762 // check that either i0 or i1 are border/seam and are on the same edge loop
763 // note that we need to add the error even for edged that connect e.g. border & locked
764 // if we don't do that, the adjacent border->border edge won't have correct errors for corners
765 if (k0 != Kind_Border && k0 != Kind_Seam && k1 != Kind_Border && k1 != Kind_Seam)
766 continue;
767
768 if ((k0 == Kind_Border || k0 == Kind_Seam) && loop[i0] != i1)
769 continue;
770
771 if ((k1 == Kind_Border || k1 == Kind_Seam) && loopback[i1] != i0)
772 continue;
773
774 // seam edges should occur twice (i0->i1 and i1->i0) - skip redundant edges
775 if (kHasOpposite[k0][k1] && remap[i1] > remap[i0])
776 continue;
777
778 unsigned int i2 = indices[i + next[next[e]]];
779
780 // we try hard to maintain border edge geometry; seam edges can move more freely
781 // due to topological restrictions on collapses, seam quadrics slightly improves collapse structure but aren't critical
782 const float kEdgeWeightSeam = 1.f;
783 const float kEdgeWeightBorder = 10.f;
784
785 float edgeWeight = (k0 == Kind_Border || k1 == Kind_Border) ? kEdgeWeightBorder : kEdgeWeightSeam;
786
787 Quadric Q;
788 quadricFromTriangleEdge(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], edgeWeight);
789
790 quadricAdd(vertex_quadrics[remap[i0]], Q);
791 quadricAdd(vertex_quadrics[remap[i1]], Q);
792
793 quadricAdd(vertex_no_attrib_quadrics[remap[i0]], Q);
794 quadricAdd(vertex_no_attrib_quadrics[remap[i1]], Q);
795 }
796 }
797}
798
799// does triangle ABC flip when C is replaced with D?
800static bool hasTriangleFlip(const Vector3& a, const Vector3& b, const Vector3& c, const Vector3& d)
801{
802 Vector3 eb = {b.x - a.x, b.y - a.y, b.z - a.z};
803 Vector3 ec = {c.x - a.x, c.y - a.y, c.z - a.z};
804 Vector3 ed = {d.x - a.x, d.y - a.y, d.z - a.z};
805
806 Vector3 nbc = {eb.y * ec.z - eb.z * ec.y, eb.z * ec.x - eb.x * ec.z, eb.x * ec.y - eb.y * ec.x};
807 Vector3 nbd = {eb.y * ed.z - eb.z * ed.y, eb.z * ed.x - eb.x * ed.z, eb.x * ed.y - eb.y * ed.x};
808
809 return nbc.x * nbd.x + nbc.y * nbd.y + nbc.z * nbd.z < 0;
810}
811
812static bool hasTriangleFlips(const EdgeAdjacency& adjacency, const Vector3* vertex_positions, const unsigned int* collapse_remap, unsigned int i0, unsigned int i1)
813{
814 assert(collapse_remap[i0] == i0);
815 assert(collapse_remap[i1] == i1);
816
817 const Vector3& v0 = vertex_positions[i0];
818 const Vector3& v1 = vertex_positions[i1];
819
820 const EdgeAdjacency::Edge* edges = &adjacency.data[adjacency.offsets[i0]];
821 size_t count = adjacency.counts[i0];
822
823 for (size_t i = 0; i < count; ++i)
824 {
825 unsigned int a = collapse_remap[edges[i].next];
826 unsigned int b = collapse_remap[edges[i].prev];
827
828 // skip triangles that get collapsed
829 // note: this is mathematically redundant as if either of these is true, the dot product in hasTriangleFlip should be 0
830 if (a == i1 || b == i1)
831 continue;
832
833 // early-out when at least one triangle flips due to a collapse
834 if (hasTriangleFlip(vertex_positions[a], vertex_positions[b], v0, v1))
835 return true;
836 }
837
838 return false;
839}
840
841static 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)
842{
843 size_t collapse_count = 0;
844
845 for (size_t i = 0; i < index_count; i += 3)
846 {
847 static const int next[3] = {1, 2, 0};
848
849 for (int e = 0; e < 3; ++e)
850 {
851 unsigned int i0 = indices[i + e];
852 unsigned int i1 = indices[i + next[e]];
853
854 // this can happen either when input has a zero-length edge, or when we perform collapses for complex
855 // topology w/seams and collapse a manifold vertex that connects to both wedges onto one of them
856 // we leave edges like this alone since they may be important for preserving mesh integrity
857 if (remap[i0] == remap[i1])
858 continue;
859
860 unsigned char k0 = vertex_kind[i0];
861 unsigned char k1 = vertex_kind[i1];
862
863 // the edge has to be collapsible in at least one direction
864 if (!(kCanCollapse[k0][k1] | kCanCollapse[k1][k0]))
865 continue;
866
867 // manifold and seam edges should occur twice (i0->i1 and i1->i0) - skip redundant edges
868 if (kHasOpposite[k0][k1] && remap[i1] > remap[i0])
869 continue;
870
871 // two vertices are on a border or a seam, but there's no direct edge between them
872 // this indicates that they belong to two different edge loops and we should not collapse this edge
873 // loop[] tracks half edges so we only need to check i0->i1
874 if (k0 == k1 && (k0 == Kind_Border || k0 == Kind_Seam) && loop[i0] != i1)
875 continue;
876
877 // edge can be collapsed in either direction - we will pick the one with minimum error
878 // note: we evaluate error later during collapse ranking, here we just tag the edge as bidirectional
879 if (kCanCollapse[k0][k1] & kCanCollapse[k1][k0])
880 {
881 Collapse c = {i0, i1, {/* bidi= */ 1}};
882 collapses[collapse_count++] = c;
883 }
884 else
885 {
886 // edge can only be collapsed in one direction
887 unsigned int e0 = kCanCollapse[k0][k1] ? i0 : i1;
888 unsigned int e1 = kCanCollapse[k0][k1] ? i1 : i0;
889
890 Collapse c = {e0, e1, {/* bidi= */ 0}};
891 collapses[collapse_count++] = c;
892 }
893 }
894 }
895
896 return collapse_count;
897}
898
899static void rankEdgeCollapses(Collapse* collapses, size_t collapse_count, const Vector3* vertex_positions, const Quadric* vertex_quadrics, const Quadric* vertex_no_attrib_quadrics, const unsigned int* remap)
900{
901 for (size_t i = 0; i < collapse_count; ++i)
902 {
903 Collapse& c = collapses[i];
904
905 unsigned int i0 = c.v0;
906 unsigned int i1 = c.v1;
907
908 // most edges are bidirectional which means we need to evaluate errors for two collapses
909 // to keep this code branchless we just use the same edge for unidirectional edges
910 unsigned int j0 = c.bidi ? i1 : i0;
911 unsigned int j1 = c.bidi ? i0 : i1;
912
913 const Quadric& qi = vertex_quadrics[remap[i0]];
914 const Quadric& qj = vertex_quadrics[remap[j0]];
915
916 float ei = quadricError(qi, vertex_positions[i1]);
917 float ej = quadricError(qj, vertex_positions[j1]);
918
919 const Quadric& naqi = vertex_no_attrib_quadrics[remap[i0]];
920 const Quadric& naqj = vertex_no_attrib_quadrics[remap[j0]];
921
922 // pick edge direction with minimal error
923 c.v0 = ei <= ej ? i0 : j0;
924 c.v1 = ei <= ej ? i1 : j1;
925 c.error = ei <= ej ? ei : ej;
926 c.distance_error = ei <= ej ? quadricErrorNoAttributes(naqi, vertex_positions[i1]) : quadricErrorNoAttributes(naqj, vertex_positions[j1]);
927 }
928}
929
930#if TRACE > 1
931static void dumpEdgeCollapses(const Collapse* collapses, size_t collapse_count, const unsigned char* vertex_kind)
932{
933 size_t ckinds[Kind_Count][Kind_Count] = {};
934 float cerrors[Kind_Count][Kind_Count] = {};
935
936 for (int k0 = 0; k0 < Kind_Count; ++k0)
937 for (int k1 = 0; k1 < Kind_Count; ++k1)
938 cerrors[k0][k1] = FLT_MAX;
939
940 for (size_t i = 0; i < collapse_count; ++i)
941 {
942 unsigned int i0 = collapses[i].v0;
943 unsigned int i1 = collapses[i].v1;
944
945 unsigned char k0 = vertex_kind[i0];
946 unsigned char k1 = vertex_kind[i1];
947
948 ckinds[k0][k1]++;
949 cerrors[k0][k1] = (collapses[i].error < cerrors[k0][k1]) ? collapses[i].error : cerrors[k0][k1];
950 }
951
952 for (int k0 = 0; k0 < Kind_Count; ++k0)
953 for (int k1 = 0; k1 < Kind_Count; ++k1)
954 if (ckinds[k0][k1])
955 printf("collapses %d -> %d: %d, min error %e\n", k0, k1, int(ckinds[k0][k1]), ckinds[k0][k1] ? sqrtf(cerrors[k0][k1]) : 0.f);
956}
957
958static void dumpLockedCollapses(const unsigned int* indices, size_t index_count, const unsigned char* vertex_kind)
959{
960 size_t locked_collapses[Kind_Count][Kind_Count] = {};
961
962 for (size_t i = 0; i < index_count; i += 3)
963 {
964 static const int next[3] = {1, 2, 0};
965
966 for (int e = 0; e < 3; ++e)
967 {
968 unsigned int i0 = indices[i + e];
969 unsigned int i1 = indices[i + next[e]];
970
971 unsigned char k0 = vertex_kind[i0];
972 unsigned char k1 = vertex_kind[i1];
973
974 locked_collapses[k0][k1] += !kCanCollapse[k0][k1] && !kCanCollapse[k1][k0];
975 }
976 }
977
978 for (int k0 = 0; k0 < Kind_Count; ++k0)
979 for (int k1 = 0; k1 < Kind_Count; ++k1)
980 if (locked_collapses[k0][k1])
981 printf("locked collapses %d -> %d: %d\n", k0, k1, int(locked_collapses[k0][k1]));
982}
983#endif
984
985static void sortEdgeCollapses(unsigned int* sort_order, const Collapse* collapses, size_t collapse_count)
986{
987 const int sort_bits = 11;
988
989 // fill histogram for counting sort
990 unsigned int histogram[1 << sort_bits];
991 memset(histogram, 0, sizeof(histogram));
992
993 for (size_t i = 0; i < collapse_count; ++i)
994 {
995 // skip sign bit since error is non-negative
996 unsigned int key = (collapses[i].errorui << 1) >> (32 - sort_bits);
997
998 histogram[key]++;
999 }
1000
1001 // compute offsets based on histogram data
1002 size_t histogram_sum = 0;
1003
1004 for (size_t i = 0; i < 1 << sort_bits; ++i)
1005 {
1006 size_t count = histogram[i];
1007 histogram[i] = unsigned(histogram_sum);
1008 histogram_sum += count;
1009 }
1010
1011 assert(histogram_sum == collapse_count);
1012
1013 // compute sort order based on offsets
1014 for (size_t i = 0; i < collapse_count; ++i)
1015 {
1016 // skip sign bit since error is non-negative
1017 unsigned int key = (collapses[i].errorui << 1) >> (32 - sort_bits);
1018
1019 sort_order[histogram[key]++] = unsigned(i);
1020 }
1021}
1022
1023static size_t performEdgeCollapses(unsigned int* collapse_remap, unsigned char* collapse_locked, Quadric* vertex_quadrics, Quadric* vertex_no_attrib_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, const Vector3* vertex_positions, const EdgeAdjacency& adjacency, size_t triangle_collapse_goal, float error_limit, float& result_error)
1024{
1025 size_t edge_collapses = 0;
1026 size_t triangle_collapses = 0;
1027
1028 // most collapses remove 2 triangles; use this to establish a bound on the pass in terms of error limit
1029 // note that edge_collapse_goal is an estimate; triangle_collapse_goal will be used to actually limit collapses
1030 size_t edge_collapse_goal = triangle_collapse_goal / 2;
1031
1032#if TRACE
1033 size_t stats[4] = {};
1034#endif
1035
1036 for (size_t i = 0; i < collapse_count; ++i)
1037 {
1038 const Collapse& c = collapses[collapse_order[i]];
1039
1040 TRACESTATS(0);
1041
1042 if (c.error > error_limit)
1043 break;
1044
1045 if (triangle_collapses >= triangle_collapse_goal)
1046 break;
1047
1048 // we limit the error in each pass based on the error of optimal last collapse; since many collapses will be locked
1049 // as they will share vertices with other successfull collapses, we need to increase the acceptable error by some factor
1050 float error_goal = edge_collapse_goal < collapse_count ? 1.5f * collapses[collapse_order[edge_collapse_goal]].error : FLT_MAX;
1051
1052 // on average, each collapse is expected to lock 6 other collapses; to avoid degenerate passes on meshes with odd
1053 // topology, we only abort if we got over 1/6 collapses accordingly.
1054 if (c.error > error_goal && triangle_collapses > triangle_collapse_goal / 6)
1055 break;
1056
1057 unsigned int i0 = c.v0;
1058 unsigned int i1 = c.v1;
1059
1060 unsigned int r0 = remap[i0];
1061 unsigned int r1 = remap[i1];
1062
1063 // we don't collapse vertices that had source or target vertex involved in a collapse
1064 // it's important to not move the vertices twice since it complicates the tracking/remapping logic
1065 // it's important to not move other vertices towards a moved vertex to preserve error since we don't re-rank collapses mid-pass
1066 if (collapse_locked[r0] | collapse_locked[r1])
1067 {
1068 TRACESTATS(1);
1069 continue;
1070 }
1071
1072 if (hasTriangleFlips(adjacency, vertex_positions, collapse_remap, r0, r1))
1073 {
1074 // adjust collapse goal since this collapse is invalid and shouldn't factor into error goal
1075 edge_collapse_goal++;
1076
1077 TRACESTATS(2);
1078 continue;
1079 }
1080
1081 assert(collapse_remap[r0] == r0);
1082 assert(collapse_remap[r1] == r1);
1083
1084 quadricAdd(vertex_quadrics[r1], vertex_quadrics[r0]);
1085 quadricAdd(vertex_no_attrib_quadrics[r1], vertex_no_attrib_quadrics[r0]);
1086
1087 if (vertex_kind[i0] == Kind_Complex)
1088 {
1089 unsigned int v = i0;
1090
1091 do
1092 {
1093 collapse_remap[v] = r1;
1094 v = wedge[v];
1095 } while (v != i0);
1096 }
1097 else if (vertex_kind[i0] == Kind_Seam)
1098 {
1099 // remap v0 to v1 and seam pair of v0 to seam pair of v1
1100 unsigned int s0 = wedge[i0];
1101 unsigned int s1 = wedge[i1];
1102
1103 assert(s0 != i0 && s1 != i1);
1104 assert(wedge[s0] == i0 && wedge[s1] == i1);
1105
1106 collapse_remap[i0] = i1;
1107 collapse_remap[s0] = s1;
1108 }
1109 else
1110 {
1111 assert(wedge[i0] == i0);
1112
1113 collapse_remap[i0] = i1;
1114 }
1115
1116 collapse_locked[r0] = 1;
1117 collapse_locked[r1] = 1;
1118
1119 // border edges collapse 1 triangle, other edges collapse 2 or more
1120 triangle_collapses += (vertex_kind[i0] == Kind_Border) ? 1 : 2;
1121 edge_collapses++;
1122
1123 result_error = result_error < c.distance_error ? c.distance_error : result_error;
1124 }
1125
1126#if TRACE
1127 float error_goal_perfect = edge_collapse_goal < collapse_count ? collapses[collapse_order[edge_collapse_goal]].error : 0.f;
1128
1129 printf("removed %d triangles, error %e (goal %e); evaluated %d/%d collapses (done %d, skipped %d, invalid %d)\n",
1130 int(triangle_collapses), sqrtf(result_error), sqrtf(error_goal_perfect),
1131 int(stats[0]), int(collapse_count), int(edge_collapses), int(stats[1]), int(stats[2]));
1132#endif
1133
1134 return edge_collapses;
1135}
1136
1137static size_t remapIndexBuffer(unsigned int* indices, size_t index_count, const unsigned int* collapse_remap)
1138{
1139 size_t write = 0;
1140
1141 for (size_t i = 0; i < index_count; i += 3)
1142 {
1143 unsigned int v0 = collapse_remap[indices[i + 0]];
1144 unsigned int v1 = collapse_remap[indices[i + 1]];
1145 unsigned int v2 = collapse_remap[indices[i + 2]];
1146
1147 // we never move the vertex twice during a single pass
1148 assert(collapse_remap[v0] == v0);
1149 assert(collapse_remap[v1] == v1);
1150 assert(collapse_remap[v2] == v2);
1151
1152 if (v0 != v1 && v0 != v2 && v1 != v2)
1153 {
1154 indices[write + 0] = v0;
1155 indices[write + 1] = v1;
1156 indices[write + 2] = v2;
1157 write += 3;
1158 }
1159 }
1160
1161 return write;
1162}
1163
1164static void remapEdgeLoops(unsigned int* loop, size_t vertex_count, const unsigned int* collapse_remap)
1165{
1166 for (size_t i = 0; i < vertex_count; ++i)
1167 {
1168 if (loop[i] != ~0u)
1169 {
1170 unsigned int l = loop[i];
1171 unsigned int r = collapse_remap[l];
1172
1173 // i == r is a special case when the seam edge is collapsed in a direction opposite to where loop goes
1174 loop[i] = (i == r) ? loop[l] : r;
1175 }
1176 }
1177}
1178
1179struct CellHasher
1180{
1181 const unsigned int* vertex_ids;
1182
1183 size_t hash(unsigned int i) const
1184 {
1185 unsigned int h = vertex_ids[i];
1186
1187 // MurmurHash2 finalizer
1188 h ^= h >> 13;
1189 h *= 0x5bd1e995;
1190 h ^= h >> 15;
1191 return h;
1192 }
1193
1194 bool equal(unsigned int lhs, unsigned int rhs) const
1195 {
1196 return vertex_ids[lhs] == vertex_ids[rhs];
1197 }
1198};
1199
1200struct IdHasher
1201{
1202 size_t hash(unsigned int id) const
1203 {
1204 unsigned int h = id;
1205
1206 // MurmurHash2 finalizer
1207 h ^= h >> 13;
1208 h *= 0x5bd1e995;
1209 h ^= h >> 15;
1210 return h;
1211 }
1212
1213 bool equal(unsigned int lhs, unsigned int rhs) const
1214 {
1215 return lhs == rhs;
1216 }
1217};
1218
1219struct TriangleHasher
1220{
1221 const unsigned int* indices;
1222
1223 size_t hash(unsigned int i) const
1224 {
1225 const unsigned int* tri = indices + i * 3;
1226
1227 // Optimized Spatial Hashing for Collision Detection of Deformable Objects
1228 return (tri[0] * 73856093) ^ (tri[1] * 19349663) ^ (tri[2] * 83492791);
1229 }
1230
1231 bool equal(unsigned int lhs, unsigned int rhs) const
1232 {
1233 const unsigned int* lt = indices + lhs * 3;
1234 const unsigned int* rt = indices + rhs * 3;
1235
1236 return lt[0] == rt[0] && lt[1] == rt[1] && lt[2] == rt[2];
1237 }
1238};
1239
1240static void computeVertexIds(unsigned int* vertex_ids, const Vector3* vertex_positions, size_t vertex_count, int grid_size)
1241{
1242 assert(grid_size >= 1 && grid_size <= 1024);
1243 float cell_scale = float(grid_size - 1);
1244
1245 for (size_t i = 0; i < vertex_count; ++i)
1246 {
1247 const Vector3& v = vertex_positions[i];
1248
1249 int xi = int(v.x * cell_scale + 0.5f);
1250 int yi = int(v.y * cell_scale + 0.5f);
1251 int zi = int(v.z * cell_scale + 0.5f);
1252
1253 vertex_ids[i] = (xi << 20) | (yi << 10) | zi;
1254 }
1255}
1256
1257static size_t countTriangles(const unsigned int* vertex_ids, const unsigned int* indices, size_t index_count)
1258{
1259 size_t result = 0;
1260
1261 for (size_t i = 0; i < index_count; i += 3)
1262 {
1263 unsigned int id0 = vertex_ids[indices[i + 0]];
1264 unsigned int id1 = vertex_ids[indices[i + 1]];
1265 unsigned int id2 = vertex_ids[indices[i + 2]];
1266
1267 result += (id0 != id1) & (id0 != id2) & (id1 != id2);
1268 }
1269
1270 return result;
1271}
1272
1273static size_t fillVertexCells(unsigned int* table, size_t table_size, unsigned int* vertex_cells, const unsigned int* vertex_ids, size_t vertex_count)
1274{
1275 CellHasher hasher = {vertex_ids};
1276
1277 memset(table, -1, table_size * sizeof(unsigned int));
1278
1279 size_t result = 0;
1280
1281 for (size_t i = 0; i < vertex_count; ++i)
1282 {
1283 unsigned int* entry = hashLookup2(table, table_size, hasher, unsigned(i), ~0u);
1284
1285 if (*entry == ~0u)
1286 {
1287 *entry = unsigned(i);
1288 vertex_cells[i] = unsigned(result++);
1289 }
1290 else
1291 {
1292 vertex_cells[i] = vertex_cells[*entry];
1293 }
1294 }
1295
1296 return result;
1297}
1298
1299static size_t countVertexCells(unsigned int* table, size_t table_size, const unsigned int* vertex_ids, size_t vertex_count)
1300{
1301 IdHasher hasher;
1302
1303 memset(table, -1, table_size * sizeof(unsigned int));
1304
1305 size_t result = 0;
1306
1307 for (size_t i = 0; i < vertex_count; ++i)
1308 {
1309 unsigned int id = vertex_ids[i];
1310 unsigned int* entry = hashLookup2(table, table_size, hasher, id, ~0u);
1311
1312 result += (*entry == ~0u);
1313 *entry = id;
1314 }
1315
1316 return result;
1317}
1318
1319static void fillCellQuadrics(Quadric* cell_quadrics, const unsigned int* indices, size_t index_count, const Vector3* vertex_positions, const unsigned int* vertex_cells)
1320{
1321 for (size_t i = 0; i < index_count; i += 3)
1322 {
1323 unsigned int i0 = indices[i + 0];
1324 unsigned int i1 = indices[i + 1];
1325 unsigned int i2 = indices[i + 2];
1326
1327 unsigned int c0 = vertex_cells[i0];
1328 unsigned int c1 = vertex_cells[i1];
1329 unsigned int c2 = vertex_cells[i2];
1330
1331 bool single_cell = (c0 == c1) & (c0 == c2);
1332
1333 Quadric Q;
1334 quadricFromTriangle(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], single_cell ? 3.f : 1.f);
1335
1336 if (single_cell)
1337 {
1338 quadricAdd(cell_quadrics[c0], Q);
1339 }
1340 else
1341 {
1342 quadricAdd(cell_quadrics[c0], Q);
1343 quadricAdd(cell_quadrics[c1], Q);
1344 quadricAdd(cell_quadrics[c2], Q);
1345 }
1346 }
1347}
1348
1349static void fillCellQuadrics(Quadric* cell_quadrics, const Vector3* vertex_positions, size_t vertex_count, const unsigned int* vertex_cells)
1350{
1351 for (size_t i = 0; i < vertex_count; ++i)
1352 {
1353 unsigned int c = vertex_cells[i];
1354 const Vector3& v = vertex_positions[i];
1355
1356 Quadric Q;
1357 quadricFromPoint(Q, v.x, v.y, v.z, 1.f);
1358
1359 quadricAdd(cell_quadrics[c], Q);
1360 }
1361}
1362
1363static 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)
1364{
1365 memset(cell_remap, -1, cell_count * sizeof(unsigned int));
1366
1367 for (size_t i = 0; i < vertex_count; ++i)
1368 {
1369 unsigned int cell = vertex_cells[i];
1370 float error = quadricError(cell_quadrics[cell], vertex_positions[i]);
1371
1372 if (cell_remap[cell] == ~0u || cell_errors[cell] > error)
1373 {
1374 cell_remap[cell] = unsigned(i);
1375 cell_errors[cell] = error;
1376 }
1377 }
1378}
1379
1380static 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)
1381{
1382 TriangleHasher hasher = {destination};
1383
1384 memset(tritable, -1, tritable_size * sizeof(unsigned int));
1385
1386 size_t result = 0;
1387
1388 for (size_t i = 0; i < index_count; i += 3)
1389 {
1390 unsigned int c0 = vertex_cells[indices[i + 0]];
1391 unsigned int c1 = vertex_cells[indices[i + 1]];
1392 unsigned int c2 = vertex_cells[indices[i + 2]];
1393
1394 if (c0 != c1 && c0 != c2 && c1 != c2)
1395 {
1396 unsigned int a = cell_remap[c0];
1397 unsigned int b = cell_remap[c1];
1398 unsigned int c = cell_remap[c2];
1399
1400 if (b < a && b < c)
1401 {
1402 unsigned int t = a;
1403 a = b, b = c, c = t;
1404 }
1405 else if (c < a && c < b)
1406 {
1407 unsigned int t = c;
1408 c = b, b = a, a = t;
1409 }
1410
1411 destination[result * 3 + 0] = a;
1412 destination[result * 3 + 1] = b;
1413 destination[result * 3 + 2] = c;
1414
1415 unsigned int* entry = hashLookup2(tritable, tritable_size, hasher, unsigned(result), ~0u);
1416
1417 if (*entry == ~0u)
1418 *entry = unsigned(result++);
1419 }
1420 }
1421
1422 return result * 3;
1423}
1424
1425static float interpolate(float y, float x0, float y0, float x1, float y1, float x2, float y2)
1426{
1427 // three point interpolation from "revenge of interpolation search" paper
1428 float num = (y1 - y) * (x1 - x2) * (x1 - x0) * (y2 - y0);
1429 float den = (y2 - y) * (x1 - x2) * (y0 - y1) + (y0 - y) * (x1 - x0) * (y1 - y2);
1430 return x1 + num / den;
1431}
1432
1433} // namespace meshopt
1434
1435#ifndef NDEBUG
1436// Note: this is only exposed for debug visualization purposes; do *not* use these in debug builds
1437MESHOPTIMIZER_API unsigned char* meshopt_simplifyDebugKind = 0;
1438MESHOPTIMIZER_API unsigned int* meshopt_simplifyDebugLoop = 0;
1439MESHOPTIMIZER_API unsigned int* meshopt_simplifyDebugLoopBack = 0;
1440#endif
1441
1442size_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, unsigned int options, float* out_result_error)
1443{
1444 return meshopt_simplifyWithAttributes(destination, indices, index_count, vertex_positions_data, vertex_count, vertex_positions_stride, target_index_count, target_error, options, out_result_error, 0, 0, 0);
1445}
1446
1447size_t meshopt_simplifyWithAttributes(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_data, size_t vertex_count, size_t vertex_stride, size_t target_index_count, float target_error, unsigned int options, float* out_result_error, const float* attributes, const float* attribute_weights, size_t attribute_count)
1448{
1449 using namespace meshopt;
1450
1451 assert(index_count % 3 == 0);
1452 assert(vertex_stride >= 12 && vertex_stride <= 256);
1453 assert(vertex_stride % sizeof(float) == 0);
1454 assert(target_index_count <= index_count);
1455 assert((options & ~(meshopt_SimplifyLockBorder)) == 0);
1456 assert(attribute_count <= ATTRIBUTES);
1457
1458 meshopt_Allocator allocator;
1459
1460 unsigned int* result = destination;
1461
1462 // build adjacency information
1463 EdgeAdjacency adjacency = {};
1464 prepareEdgeAdjacency(adjacency, index_count, vertex_count, allocator);
1465 updateEdgeAdjacency(adjacency, indices, index_count, vertex_count, NULL);
1466
1467 // build position remap that maps each vertex to the one with identical position
1468 unsigned int* remap = allocator.allocate<unsigned int>(vertex_count);
1469 unsigned int* wedge = allocator.allocate<unsigned int>(vertex_count);
1470 buildPositionRemap(remap, wedge, vertex_data, vertex_count, vertex_stride, allocator);
1471
1472 // classify vertices; vertex kind determines collapse rules, see kCanCollapse
1473 unsigned char* vertex_kind = allocator.allocate<unsigned char>(vertex_count);
1474 unsigned int* loop = allocator.allocate<unsigned int>(vertex_count);
1475 unsigned int* loopback = allocator.allocate<unsigned int>(vertex_count);
1476 classifyVertices(vertex_kind, loop, loopback, vertex_count, adjacency, remap, wedge, options);
1477
1478#if TRACE
1479 size_t unique_positions = 0;
1480 for (size_t i = 0; i < vertex_count; ++i)
1481 unique_positions += remap[i] == i;
1482
1483 printf("position remap: %d vertices => %d positions\n", int(vertex_count), int(unique_positions));
1484
1485 size_t kinds[Kind_Count] = {};
1486 for (size_t i = 0; i < vertex_count; ++i)
1487 kinds[vertex_kind[i]] += remap[i] == i;
1488
1489 printf("kinds: manifold %d, border %d, seam %d, complex %d, locked %d\n",
1490 int(kinds[Kind_Manifold]), int(kinds[Kind_Border]), int(kinds[Kind_Seam]), int(kinds[Kind_Complex]), int(kinds[Kind_Locked]));
1491#endif
1492
1493 Vector3* vertex_positions = allocator.allocate<Vector3>(vertex_count);
1494 rescalePositions(vertex_positions, vertex_data, vertex_count, vertex_stride);
1495
1496#if ATTRIBUTES
1497 for (size_t i = 0; i < vertex_count; ++i)
1498 {
1499 memset(vertex_positions[i].a, 0, sizeof(vertex_positions[i].a));
1500
1501 for (size_t k = 0; k < attribute_count; ++k)
1502 {
1503 float a = attributes[i * attribute_count + k];
1504
1505 vertex_positions[i].a[k] = a * attribute_weights[k];
1506 }
1507 }
1508#endif
1509
1510 Quadric* vertex_quadrics = allocator.allocate<Quadric>(vertex_count);
1511 memset(vertex_quadrics, 0, vertex_count * sizeof(Quadric));
1512 Quadric* vertex_no_attrib_quadrics = allocator.allocate<Quadric>(vertex_count);
1513 memset(vertex_no_attrib_quadrics, 0, vertex_count * sizeof(Quadric));
1514
1515 fillFaceQuadrics(vertex_quadrics, vertex_no_attrib_quadrics, indices, index_count, vertex_positions, remap);
1516 fillEdgeQuadrics(vertex_quadrics, vertex_no_attrib_quadrics, indices, index_count, vertex_positions, remap, vertex_kind, loop, loopback);
1517
1518 if (result != indices)
1519 memcpy(result, indices, index_count * sizeof(unsigned int));
1520
1521#if TRACE
1522 size_t pass_count = 0;
1523#endif
1524
1525 Collapse* edge_collapses = allocator.allocate<Collapse>(index_count);
1526 unsigned int* collapse_order = allocator.allocate<unsigned int>(index_count);
1527 unsigned int* collapse_remap = allocator.allocate<unsigned int>(vertex_count);
1528 unsigned char* collapse_locked = allocator.allocate<unsigned char>(vertex_count);
1529
1530 size_t result_count = index_count;
1531 float result_error = 0;
1532
1533 // target_error input is linear; we need to adjust it to match quadricError units
1534 float error_limit = target_error * target_error;
1535
1536 while (result_count > target_index_count)
1537 {
1538 // note: throughout the simplification process adjacency structure reflects welded topology for result-in-progress
1539 updateEdgeAdjacency(adjacency, result, result_count, vertex_count, remap);
1540
1541 size_t edge_collapse_count = pickEdgeCollapses(edge_collapses, result, result_count, remap, vertex_kind, loop);
1542
1543 // no edges can be collapsed any more due to topology restrictions
1544 if (edge_collapse_count == 0)
1545 break;
1546
1547 rankEdgeCollapses(edge_collapses, edge_collapse_count, vertex_positions, vertex_quadrics, vertex_no_attrib_quadrics, remap);
1548
1549#if TRACE > 1
1550 dumpEdgeCollapses(edge_collapses, edge_collapse_count, vertex_kind);
1551#endif
1552
1553 sortEdgeCollapses(collapse_order, edge_collapses, edge_collapse_count);
1554
1555 size_t triangle_collapse_goal = (result_count - target_index_count) / 3;
1556
1557 for (size_t i = 0; i < vertex_count; ++i)
1558 collapse_remap[i] = unsigned(i);
1559
1560 memset(collapse_locked, 0, vertex_count);
1561
1562#if TRACE
1563 printf("pass %d: ", int(pass_count++));
1564#endif
1565
1566 size_t collapses = performEdgeCollapses(collapse_remap, collapse_locked, vertex_quadrics, vertex_no_attrib_quadrics, edge_collapses, edge_collapse_count, collapse_order, remap, wedge, vertex_kind, vertex_positions, adjacency, triangle_collapse_goal, error_limit, result_error);
1567
1568 // no edges can be collapsed any more due to hitting the error limit or triangle collapse limit
1569 if (collapses == 0)
1570 break;
1571
1572 remapEdgeLoops(loop, vertex_count, collapse_remap);
1573 remapEdgeLoops(loopback, vertex_count, collapse_remap);
1574
1575 size_t new_count = remapIndexBuffer(result, result_count, collapse_remap);
1576 assert(new_count < result_count);
1577
1578 result_count = new_count;
1579 }
1580
1581#if TRACE
1582 printf("result: %d triangles, error: %e; total %d passes\n", int(result_count), sqrtf(result_error), int(pass_count));
1583#endif
1584
1585#if TRACE > 1
1586 dumpLockedCollapses(result, result_count, vertex_kind);
1587#endif
1588
1589#ifndef NDEBUG
1590 if (meshopt_simplifyDebugKind)
1591 memcpy(meshopt_simplifyDebugKind, vertex_kind, vertex_count);
1592
1593 if (meshopt_simplifyDebugLoop)
1594 memcpy(meshopt_simplifyDebugLoop, loop, vertex_count * sizeof(unsigned int));
1595
1596 if (meshopt_simplifyDebugLoopBack)
1597 memcpy(meshopt_simplifyDebugLoopBack, loopback, vertex_count * sizeof(unsigned int));
1598#endif
1599
1600 // result_error is quadratic; we need to remap it back to linear
1601 if (out_result_error)
1602 {
1603 *out_result_error = sqrtf(result_error);
1604 }
1605
1606 return result_count;
1607}
1608
1609size_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, float target_error, float* out_result_error)
1610{
1611 using namespace meshopt;
1612
1613 assert(index_count % 3 == 0);
1614 assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
1615 assert(vertex_positions_stride % sizeof(float) == 0);
1616 assert(target_index_count <= index_count);
1617
1618 // we expect to get ~2 triangles/vertex in the output
1619 size_t target_cell_count = target_index_count / 6;
1620
1621 meshopt_Allocator allocator;
1622
1623 Vector3* vertex_positions = allocator.allocate<Vector3>(vertex_count);
1624 rescalePositions(vertex_positions, vertex_positions_data, vertex_count, vertex_positions_stride);
1625
1626 // find the optimal grid size using guided binary search
1627#if TRACE
1628 printf("source: %d vertices, %d triangles\n", int(vertex_count), int(index_count / 3));
1629 printf("target: %d cells, %d triangles\n", int(target_cell_count), int(target_index_count / 3));
1630#endif
1631
1632 unsigned int* vertex_ids = allocator.allocate<unsigned int>(vertex_count);
1633
1634 const int kInterpolationPasses = 5;
1635
1636 // invariant: # of triangles in min_grid <= target_count
1637 int min_grid = int(1.f / (target_error < 1e-3f ? 1e-3f : target_error));
1638 int max_grid = 1025;
1639 size_t min_triangles = 0;
1640 size_t max_triangles = index_count / 3;
1641
1642 // when we're error-limited, we compute the triangle count for the min. size; this accelerates convergence and provides the correct answer when we can't use a larger grid
1643 if (min_grid > 1)
1644 {
1645 computeVertexIds(vertex_ids, vertex_positions, vertex_count, min_grid);
1646 min_triangles = countTriangles(vertex_ids, indices, index_count);
1647 }
1648
1649 // 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...
1650 int next_grid_size = int(sqrtf(float(target_cell_count)) + 0.5f);
1651
1652 for (int pass = 0; pass < 10 + kInterpolationPasses; ++pass)
1653 {
1654 if (min_triangles >= target_index_count / 3 || max_grid - min_grid <= 1)
1655 break;
1656
1657 // we clamp the prediction of the grid size to make sure that the search converges
1658 int grid_size = next_grid_size;
1659 grid_size = (grid_size <= min_grid) ? min_grid + 1 : (grid_size >= max_grid) ? max_grid - 1 : grid_size;
1660
1661 computeVertexIds(vertex_ids, vertex_positions, vertex_count, grid_size);
1662 size_t triangles = countTriangles(vertex_ids, indices, index_count);
1663
1664#if TRACE
1665 printf("pass %d (%s): grid size %d, triangles %d, %s\n",
1666 pass, (pass == 0) ? "guess" : (pass <= kInterpolationPasses) ? "lerp" : "binary",
1667 grid_size, int(triangles),
1668 (triangles <= target_index_count / 3) ? "under" : "over");
1669#endif
1670
1671 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));
1672
1673 if (triangles <= target_index_count / 3)
1674 {
1675 min_grid = grid_size;
1676 min_triangles = triangles;
1677 }
1678 else
1679 {
1680 max_grid = grid_size;
1681 max_triangles = triangles;
1682 }
1683
1684 // we start by using interpolation search - it usually converges faster
1685 // 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)
1686 next_grid_size = (pass < kInterpolationPasses) ? int(tip + 0.5f) : (min_grid + max_grid) / 2;
1687 }
1688
1689 if (min_triangles == 0)
1690 {
1691 if (out_result_error)
1692 *out_result_error = 1.f;
1693
1694 return 0;
1695 }
1696
1697 // build vertex->cell association by mapping all vertices with the same quantized position to the same cell
1698 size_t table_size = hashBuckets2(vertex_count);
1699 unsigned int* table = allocator.allocate<unsigned int>(table_size);
1700
1701 unsigned int* vertex_cells = allocator.allocate<unsigned int>(vertex_count);
1702
1703 computeVertexIds(vertex_ids, vertex_positions, vertex_count, min_grid);
1704 size_t cell_count = fillVertexCells(table, table_size, vertex_cells, vertex_ids, vertex_count);
1705
1706 // build a quadric for each target cell
1707 Quadric* cell_quadrics = allocator.allocate<Quadric>(cell_count);
1708 memset(cell_quadrics, 0, cell_count * sizeof(Quadric));
1709
1710 fillCellQuadrics(cell_quadrics, indices, index_count, vertex_positions, vertex_cells);
1711
1712 // for each target cell, find the vertex with the minimal error
1713 unsigned int* cell_remap = allocator.allocate<unsigned int>(cell_count);
1714 float* cell_errors = allocator.allocate<float>(cell_count);
1715
1716 fillCellRemap(cell_remap, cell_errors, cell_count, vertex_cells, cell_quadrics, vertex_positions, vertex_count);
1717
1718 // compute error
1719 float result_error = 0.f;
1720
1721 for (size_t i = 0; i < cell_count; ++i)
1722 result_error = result_error < cell_errors[i] ? cell_errors[i] : result_error;
1723
1724 // collapse triangles!
1725 // note that we need to filter out triangles that we've already output because we very frequently generate redundant triangles between cells :(
1726 size_t tritable_size = hashBuckets2(min_triangles);
1727 unsigned int* tritable = allocator.allocate<unsigned int>(tritable_size);
1728
1729 size_t write = filterTriangles(destination, tritable, tritable_size, indices, index_count, vertex_cells, cell_remap);
1730
1731#if TRACE
1732 printf("result: %d cells, %d triangles (%d unfiltered), error %e\n", int(cell_count), int(write / 3), int(min_triangles), sqrtf(result_error));
1733#endif
1734
1735 if (out_result_error)
1736 *out_result_error = sqrtf(result_error);
1737
1738 return write;
1739}
1740
1741size_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)
1742{
1743 using namespace meshopt;
1744
1745 assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
1746 assert(vertex_positions_stride % sizeof(float) == 0);
1747 assert(target_vertex_count <= vertex_count);
1748
1749 size_t target_cell_count = target_vertex_count;
1750
1751 if (target_cell_count == 0)
1752 return 0;
1753
1754 meshopt_Allocator allocator;
1755
1756 Vector3* vertex_positions = allocator.allocate<Vector3>(vertex_count);
1757 rescalePositions(vertex_positions, vertex_positions_data, vertex_count, vertex_positions_stride);
1758
1759 // find the optimal grid size using guided binary search
1760#if TRACE
1761 printf("source: %d vertices\n", int(vertex_count));
1762 printf("target: %d cells\n", int(target_cell_count));
1763#endif
1764
1765 unsigned int* vertex_ids = allocator.allocate<unsigned int>(vertex_count);
1766
1767 size_t table_size = hashBuckets2(vertex_count);
1768 unsigned int* table = allocator.allocate<unsigned int>(table_size);
1769
1770 const int kInterpolationPasses = 5;
1771
1772 // invariant: # of vertices in min_grid <= target_count
1773 int min_grid = 0;
1774 int max_grid = 1025;
1775 size_t min_vertices = 0;
1776 size_t max_vertices = vertex_count;
1777
1778 // 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...
1779 int next_grid_size = int(sqrtf(float(target_cell_count)) + 0.5f);
1780
1781 for (int pass = 0; pass < 10 + kInterpolationPasses; ++pass)
1782 {
1783 assert(min_vertices < target_vertex_count);
1784 assert(max_grid - min_grid > 1);
1785
1786 // we clamp the prediction of the grid size to make sure that the search converges
1787 int grid_size = next_grid_size;
1788 grid_size = (grid_size <= min_grid) ? min_grid + 1 : (grid_size >= max_grid) ? max_grid - 1 : grid_size;
1789
1790 computeVertexIds(vertex_ids, vertex_positions, vertex_count, grid_size);
1791 size_t vertices = countVertexCells(table, table_size, vertex_ids, vertex_count);
1792
1793#if TRACE
1794 printf("pass %d (%s): grid size %d, vertices %d, %s\n",
1795 pass, (pass == 0) ? "guess" : (pass <= kInterpolationPasses) ? "lerp" : "binary",
1796 grid_size, int(vertices),
1797 (vertices <= target_vertex_count) ? "under" : "over");
1798#endif
1799
1800 float tip = interpolate(float(target_vertex_count), float(min_grid), float(min_vertices), float(grid_size), float(vertices), float(max_grid), float(max_vertices));
1801
1802 if (vertices <= target_vertex_count)
1803 {
1804 min_grid = grid_size;
1805 min_vertices = vertices;
1806 }
1807 else
1808 {
1809 max_grid = grid_size;
1810 max_vertices = vertices;
1811 }
1812
1813 if (vertices == target_vertex_count || max_grid - min_grid <= 1)
1814 break;
1815
1816 // we start by using interpolation search - it usually converges faster
1817 // 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)
1818 next_grid_size = (pass < kInterpolationPasses) ? int(tip + 0.5f) : (min_grid + max_grid) / 2;
1819 }
1820
1821 if (min_vertices == 0)
1822 return 0;
1823
1824 // build vertex->cell association by mapping all vertices with the same quantized position to the same cell
1825 unsigned int* vertex_cells = allocator.allocate<unsigned int>(vertex_count);
1826
1827 computeVertexIds(vertex_ids, vertex_positions, vertex_count, min_grid);
1828 size_t cell_count = fillVertexCells(table, table_size, vertex_cells, vertex_ids, vertex_count);
1829
1830 // build a quadric for each target cell
1831 Quadric* cell_quadrics = allocator.allocate<Quadric>(cell_count);
1832 memset(cell_quadrics, 0, cell_count * sizeof(Quadric));
1833
1834 fillCellQuadrics(cell_quadrics, vertex_positions, vertex_count, vertex_cells);
1835
1836 // for each target cell, find the vertex with the minimal error
1837 unsigned int* cell_remap = allocator.allocate<unsigned int>(cell_count);
1838 float* cell_errors = allocator.allocate<float>(cell_count);
1839
1840 fillCellRemap(cell_remap, cell_errors, cell_count, vertex_cells, cell_quadrics, vertex_positions, vertex_count);
1841
1842 // copy results to the output
1843 assert(cell_count <= target_vertex_count);
1844 memcpy(destination, cell_remap, sizeof(unsigned int) * cell_count);
1845
1846#if TRACE
1847 printf("result: %d cells\n", int(cell_count));
1848#endif
1849
1850 return cell_count;
1851}
1852
1853float meshopt_simplifyScale(const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride)
1854{
1855 using namespace meshopt;
1856
1857 assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
1858 assert(vertex_positions_stride % sizeof(float) == 0);
1859
1860 float extent = rescalePositions(NULL, vertex_positions, vertex_count, vertex_positions_stride);
1861
1862 return extent;
1863}
1864