zeux / meshoptimizer
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
// This work is based on:
24
// Michael Garland and Paul S. Heckbert. Surface simplification using quadric error metrics. 1997
25
// Michael Garland. Quadric-based polygonal surface simplification. 1999
26
// Peter Lindstrom. Out-of-Core Simplification of Large Polygonal Models. 2000
27
// Matthias Teschner, Bruno Heidelberger, Matthias Mueller, Danat Pomeranets, Markus Gross. Optimized Spatial Hashing for Collision Detection of Deformable Objects. 2003
28
// Peter Van Sandt, Yannis Chronis, Jignesh M. Patel. Efficiently Searching In-Memory Sorted Arrays: Revenge of the Interpolation Search? 2019
29
namespace meshopt
30
{
31

32
struct EdgeAdjacency
33
{
34
	struct Edge
35
	{
36
		unsigned int next;
37
		unsigned int prev;
38
	};
39

40
	unsigned int* counts;
41
	unsigned int* offsets;
42
	Edge* data;
43
};
44

45 3
static void prepareEdgeAdjacency(EdgeAdjacency& adjacency, size_t index_count, size_t vertex_count, meshopt_Allocator& allocator)
46
{
47 3
	adjacency.counts = allocator.allocate<unsigned int>(vertex_count);
48 3
	adjacency.offsets = allocator.allocate<unsigned int>(vertex_count);
49 3
	adjacency.data = allocator.allocate<EdgeAdjacency::Edge>(index_count);
50
}
51

52 3
static void updateEdgeAdjacency(EdgeAdjacency& adjacency, const unsigned int* indices, size_t index_count, size_t vertex_count, const unsigned int* remap)
53
{
54 3
	size_t face_count = index_count / 3;
55

56
	// fill edge counts
57 3
	memset(adjacency.counts, 0, vertex_count * sizeof(unsigned int));
58

59 3
	for (size_t i = 0; i < index_count; ++i)
60
	{
61 3
		unsigned int v = remap ? remap[indices[i]] : indices[i];
62 1
		assert(v < vertex_count);
63

64 3
		adjacency.counts[v]++;
65
	}
66

67
	// fill offset table
68 3
	unsigned int offset = 0;
69

70 3
	for (size_t i = 0; i < vertex_count; ++i)
71
	{
72 3
		adjacency.offsets[i] = offset;
73 3
		offset += adjacency.counts[i];
74
	}
75

76 1
	assert(offset == index_count);
77

78
	// fill edge data
79 3
	for (size_t i = 0; i < face_count; ++i)
80
	{
81 3
		unsigned int a = indices[i * 3 + 0], b = indices[i * 3 + 1], c = indices[i * 3 + 2];
82

83 3
		if (remap)
84
		{
85 3
			a = remap[a];
86 3
			b = remap[b];
87 3
			c = remap[c];
88
		}
89

90 3
		adjacency.data[adjacency.offsets[a]].next = b;
91 3
		adjacency.data[adjacency.offsets[a]].prev = c;
92 3
		adjacency.offsets[a]++;
93

94 3
		adjacency.data[adjacency.offsets[b]].next = c;
95 3
		adjacency.data[adjacency.offsets[b]].prev = a;
96 3
		adjacency.offsets[b]++;
97

98 3
		adjacency.data[adjacency.offsets[c]].next = a;
99 3
		adjacency.data[adjacency.offsets[c]].prev = b;
100 3
		adjacency.offsets[c]++;
101
	}
102

103
	// fix offsets that have been disturbed by the previous pass
104 3
	for (size_t i = 0; i < vertex_count; ++i)
105
	{
106 1
		assert(adjacency.offsets[i] >= adjacency.counts[i]);
107

108 3
		adjacency.offsets[i] -= adjacency.counts[i];
109
	}
110
}
111

112
struct PositionHasher
113
{
114
	const float* vertex_positions;
115
	size_t vertex_stride_float;
116

117 3
	size_t hash(unsigned int index) const
118
	{
119 3
		const unsigned int* key = reinterpret_cast<const unsigned int*>(vertex_positions + index * vertex_stride_float);
120

121
		// scramble bits to make sure that integer coordinates have entropy in lower bits
122 3
		unsigned int x = key[0] ^ (key[0] >> 17);
123 3
		unsigned int y = key[1] ^ (key[1] >> 17);
124 3
		unsigned int z = key[2] ^ (key[2] >> 17);
125

126
		// Optimized Spatial Hashing for Collision Detection of Deformable Objects
127 3
		return (x * 73856093) ^ (y * 19349663) ^ (z * 83492791);
128
	}
129

130 3
	bool equal(unsigned int lhs, unsigned int rhs) const
131
	{
132 3
		return memcmp(vertex_positions + lhs * vertex_stride_float, vertex_positions + rhs * vertex_stride_float, sizeof(float) * 3) == 0;
133
	}
134
};
135

136 3
static size_t hashBuckets2(size_t count)
137
{
138 3
	size_t buckets = 1;
139 3
	while (buckets < count + count / 4)
140 3
		buckets *= 2;
141

142 3
	return buckets;
143
}
144

145
template <typename T, typename Hash>
146 3
static T* hashLookup2(T* table, size_t buckets, const Hash& hash, const T& key, const T& empty)
147
{
148 1
	assert(buckets > 0);
149 1
	assert((buckets & (buckets - 1)) == 0);
150

151 3
	size_t hashmod = buckets - 1;
152 3
	size_t bucket = hash.hash(key) & hashmod;
153

154 3
	for (size_t probe = 0; probe <= hashmod; ++probe)
155
	{
156 3
		T& item = table[bucket];
157

158 3
		if (item == empty)
159 3
			return &item;
160

161 3
		if (hash.equal(item, key))
162 3
			return &item;
163

164
		// hash collision, quadratic probing
165 3
		bucket = (bucket + probe + 1) & hashmod;
166
	}
167

168
	assert(false && "Hash table is full"); // unreachable
169
	return 0;
170
}
171

172 3
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)
173
{
174 3
	PositionHasher hasher = {vertex_positions_data, vertex_positions_stride / sizeof(float)};
175

176 3
	size_t table_size = hashBuckets2(vertex_count);
177 3
	unsigned int* table = allocator.allocate<unsigned int>(table_size);
178 3
	memset(table, -1, table_size * sizeof(unsigned int));
179

180
	// build forward remap: for each vertex, which other (canonical) vertex does it map to?
181
	// we use position equivalence for this, and remap vertices to other existing vertices
182 3
	for (size_t i = 0; i < vertex_count; ++i)
183
	{
184 3
		unsigned int index = unsigned(i);
185 3
		unsigned int* entry = hashLookup2(table, table_size, hasher, index, ~0u);
186

187 3
		if (*entry == ~0u)
188 3
			*entry = index;
189

190 3
		remap[index] = *entry;
191
	}
192

193
	// build wedge table: for each vertex, which other vertex is the next wedge that also maps to the same vertex?
194
	// entries in table form a (cyclic) wedge loop per vertex; for manifold vertices, wedge[i] == remap[i] == i
195 3
	for (size_t i = 0; i < vertex_count; ++i)
196 3
		wedge[i] = unsigned(i);
197

198 3
	for (size_t i = 0; i < vertex_count; ++i)
199 3
		if (remap[i] != i)
200
		{
201 3
			unsigned int r = remap[i];
202

203 3
			wedge[i] = wedge[r];
204 3
			wedge[r] = unsigned(i);
205
		}
206
}
207

208
enum VertexKind
209
{
210
	Kind_Manifold, // not on an attribute seam, not on any boundary
211
	Kind_Border,   // not on an attribute seam, has exactly two open edges
212
	Kind_Seam,     // on an attribute seam with exactly two attribute seam edges
213
	Kind_Complex,  // none of the above; these vertices can move as long as all wedges move to the target vertex
214
	Kind_Locked,   // none of the above; these vertices can't move
215

216
	Kind_Count
217
};
218

219
// manifold vertices can collapse onto anything
220
// border/seam vertices can only be collapsed onto border/seam respectively
221
// complex vertices can collapse onto complex/locked
222
// a rule of thumb is that collapsing kind A into kind B preserves the kind B in the target vertex
223
// for example, while we could collapse Complex into Manifold, this would mean the target vertex isn't Manifold anymore
224
const unsigned char kCanCollapse[Kind_Count][Kind_Count] = {
225
    {1, 1, 1, 1, 1},
226
    {0, 1, 0, 0, 0},
227
    {0, 0, 1, 0, 0},
228
    {0, 0, 0, 1, 1},
229
    {0, 0, 0, 0, 0},
230
};
231

232
// if a vertex is manifold or seam, adjoining edges are guaranteed to have an opposite edge
233
// note that for seam edges, the opposite edge isn't present in the attribute-based topology
234
// but is present if you consider a position-only mesh variant
235
const unsigned char kHasOpposite[Kind_Count][Kind_Count] = {
236
    {1, 1, 1, 0, 1},
237
    {1, 0, 1, 0, 0},
238
    {1, 1, 1, 0, 1},
239
    {0, 0, 0, 0, 0},
240
    {1, 0, 1, 0, 0},
241
};
242

243 3
static bool hasEdge(const EdgeAdjacency& adjacency, unsigned int a, unsigned int b)
244
{
245 3
	unsigned int count = adjacency.counts[a];
246 3
	const EdgeAdjacency::Edge* edges = adjacency.data + adjacency.offsets[a];
247

248 3
	for (size_t i = 0; i < count; ++i)
249 3
		if (edges[i].next == b)
250 3
			return true;
251

252 3
	return false;
253
}
254

255 3
static 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)
256
{
257 3
	memset(loop, -1, vertex_count * sizeof(unsigned int));
258 3
	memset(loopback, -1, vertex_count * sizeof(unsigned int));
259

260
	// incoming & outgoing open edges: ~0u if no open edges, i if there are more than 1
261
	// note that this is the same data as required in loop[] arrays; loop[] data is only valid for border/seam
262
	// but here it's okay to fill the data out for other types of vertices as well
263 3
	unsigned int* openinc = loopback;
264 3
	unsigned int* openout = loop;
265

266 3
	for (size_t i = 0; i < vertex_count; ++i)
267
	{
268 3
		unsigned int vertex = unsigned(i);
269

270 3
		unsigned int count = adjacency.counts[vertex];
271 3
		const EdgeAdjacency::Edge* edges = adjacency.data + adjacency.offsets[vertex];
272

273 3
		for (size_t j = 0; j < count; ++j)
274
		{
275 3
			unsigned int target = edges[j].next;
276

277 3
			if (!hasEdge(adjacency, target, vertex))
278
			{
279 3
				openinc[target] = (openinc[target] == ~0u) ? vertex : target;
280 3
				openout[vertex] = (openout[vertex] == ~0u) ? target : vertex;
281
			}
282
		}
283
	}
284

285
#if TRACE
286
	size_t stats[4] = {};
287
#endif
288

289 3
	for (size_t i = 0; i < vertex_count; ++i)
290
	{
291 3
		if (remap[i] == i)
292
		{
293 3
			if (wedge[i] == i)
294
			{
295
				// no attribute seam, need to check if it's manifold
296 3
				unsigned int openi = openinc[i], openo = openout[i];
297

298
				// note: we classify any vertices with no open edges as manifold
299
				// this is technically incorrect - if 4 triangles share an edge, we'll classify vertices as manifold
300
				// it's unclear if this is a problem in practice
301 3
				if (openi == ~0u && openo == ~0u)
302
				{
303 3
					result[i] = Kind_Manifold;
304
				}
305 3
				else if (openi != i && openo != i)
306
				{
307 3
					result[i] = Kind_Border;
308
				}
309
				else
310
				{
311 3
					result[i] = Kind_Locked;
312
					TRACESTATS(0);
313
				}
314
			}
315 3
			else if (wedge[wedge[i]] == i)
316
			{
317
				// attribute seam; need to distinguish between Seam and Locked
318 3
				unsigned int w = wedge[i];
319 3
				unsigned int openiv = openinc[i], openov = openout[i];
320 3
				unsigned int openiw = openinc[w], openow = openout[w];
321

322
				// seam should have one open half-edge for each vertex, and the edges need to "connect" - point to the same vertex post-remap
323 3
				if (openiv != ~0u && openiv != i && openov != ~0u && openov != i &&
324 3
				    openiw != ~0u && openiw != w && openow != ~0u && openow != w)
325
				{
326 3
					if (remap[openiv] == remap[openow] && remap[openov] == remap[openiw])
327
					{
328 3
						result[i] = Kind_Seam;
329
					}
330
					else
331
					{
332 3
						result[i] = Kind_Locked;
333
						TRACESTATS(1);
334
					}
335
				}
336
				else
337
				{
338 3
					result[i] = Kind_Locked;
339
					TRACESTATS(2);
340
				}
341
			}
342
			else
343
			{
344
				// more than one vertex maps to this one; we don't have classification available
345 3
				result[i] = Kind_Locked;
346
				TRACESTATS(3);
347
			}
348
		}
349
		else
350
		{
351 1
			assert(remap[i] < i);
352

353 3
			result[i] = result[remap[i]];
354
		}
355
	}
356

357
#if TRACE
358
	printf("locked: many open edges %d, disconnected seam %d, many seam edges %d, many wedges %d\n",
359
	    int(stats[0]), int(stats[1]), int(stats[2]), int(stats[3]));
360
#endif
361
}
362

363
struct Vector3
364
{
365
	float x, y, z;
366
};
367

368 3
static float rescalePositions(Vector3* result, const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride)
369
{
370 3
	size_t vertex_stride_float = vertex_positions_stride / sizeof(float);
371

372 3
	float minv[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
373 3
	float maxv[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX};
374

375 3
	for (size_t i = 0; i < vertex_count; ++i)
376
	{
377 3
		const float* v = vertex_positions_data + i * vertex_stride_float;
378

379 3
		if (result)
380
		{
381 3
			result[i].x = v[0];
382 3
			result[i].y = v[1];
383 3
			result[i].z = v[2];
384
		}
385

386 3
		for (int j = 0; j < 3; ++j)
387
		{
388 3
			float vj = v[j];
389

390 3
			minv[j] = minv[j] > vj ? vj : minv[j];
391 3
			maxv[j] = maxv[j] < vj ? vj : maxv[j];
392
		}
393
	}
394

395 3
	float extent = 0.f;
396

397 1
	extent = (maxv[0] - minv[0]) < extent ? extent : (maxv[0] - minv[0]);
398 3
	extent = (maxv[1] - minv[1]) < extent ? extent : (maxv[1] - minv[1]);
399 3
	extent = (maxv[2] - minv[2]) < extent ? extent : (maxv[2] - minv[2]);
400

401 3
	if (result)
402
	{
403 3
		float scale = extent == 0 ? 0.f : 1.f / extent;
404

405 3
		for (size_t i = 0; i < vertex_count; ++i)
406
		{
407 3
			result[i].x = (result[i].x - minv[0]) * scale;
408 3
			result[i].y = (result[i].y - minv[1]) * scale;
409 3
			result[i].z = (result[i].z - minv[2]) * scale;
410
		}
411
	}
412

413 3
	return extent;
414
}
415

416
struct Quadric
417
{
418
	float a00, a11, a22;
419
	float a10, a20, a21;
420
	float b0, b1, b2, c;
421
	float w;
422
};
423

424
struct Collapse
425
{
426
	unsigned int v0;
427
	unsigned int v1;
428

429
	union
430
	{
431
		unsigned int bidi;
432
		float error;
433
		unsigned int errorui;
434
	};
435
};
436

437 3
static float normalize(Vector3& v)
438
{
439 3
	float length = sqrtf(v.x * v.x + v.y * v.y + v.z * v.z);
440

441 3
	if (length > 0)
442
	{
443 3
		v.x /= length;
444 3
		v.y /= length;
445 3
		v.z /= length;
446
	}
447

448 3
	return length;
449
}
450

451 3
static void quadricAdd(Quadric& Q, const Quadric& R)
452
{
453 3
	Q.a00 += R.a00;
454 3
	Q.a11 += R.a11;
455 3
	Q.a22 += R.a22;
456 3
	Q.a10 += R.a10;
457 3
	Q.a20 += R.a20;
458 3
	Q.a21 += R.a21;
459 3
	Q.b0 += R.b0;
460 3
	Q.b1 += R.b1;
461 3
	Q.b2 += R.b2;
462 3
	Q.c += R.c;
463 3
	Q.w += R.w;
464
}
465

466 3
static float quadricError(const Quadric& Q, const Vector3& v)
467
{
468 3
	float rx = Q.b0;
469 3
	float ry = Q.b1;
470 3
	float rz = Q.b2;
471

472 3
	rx += Q.a10 * v.y;
473 3
	ry += Q.a21 * v.z;
474 3
	rz += Q.a20 * v.x;
475

476 3
	rx *= 2;
477 3
	ry *= 2;
478 3
	rz *= 2;
479

480 3
	rx += Q.a00 * v.x;
481 3
	ry += Q.a11 * v.y;
482 3
	rz += Q.a22 * v.z;
483

484 3
	float r = Q.c;
485 3
	r += rx * v.x;
486 3
	r += ry * v.y;
487 3
	r += rz * v.z;
488

489 1
	float s = Q.w == 0.f ? 0.f : 1.f / Q.w;
490

491 3
	return fabsf(r) * s;
492
}
493

494 3
static void quadricFromPlane(Quadric& Q, float a, float b, float c, float d, float w)
495
{
496 3
	float aw = a * w;
497 3
	float bw = b * w;
498 3
	float cw = c * w;
499 3
	float dw = d * w;
500

501 3
	Q.a00 = a * aw;
502 3
	Q.a11 = b * bw;
503 3
	Q.a22 = c * cw;
504 3
	Q.a10 = a * bw;
505 3
	Q.a20 = a * cw;
506 3
	Q.a21 = b * cw;
507 3
	Q.b0 = a * dw;
508 3
	Q.b1 = b * dw;
509 3
	Q.b2 = c * dw;
510 3
	Q.c = d * dw;
511 3
	Q.w = w;
512
}
513

514 3
static void quadricFromPoint(Quadric& Q, float x, float y, float z, float w)
515
{
516
	// we need to encode (x - X) ^ 2 + (y - Y)^2 + (z - Z)^2 into the quadric
517 3
	Q.a00 = w;
518 3
	Q.a11 = w;
519 3
	Q.a22 = w;
520 3
	Q.a10 = 0.f;
521 3
	Q.a20 = 0.f;
522 3
	Q.a21 = 0.f;
523 3
	Q.b0 = -2.f * x * w;
524 3
	Q.b1 = -2.f * y * w;
525 3
	Q.b2 = -2.f * z * w;
526 3
	Q.c = (x * x + y * y + z * z) * w;
527 3
	Q.w = w;
528
}
529

530 3
static void quadricFromTriangle(Quadric& Q, const Vector3& p0, const Vector3& p1, const Vector3& p2, float weight)
531
{
532 3
	Vector3 p10 = {p1.x - p0.x, p1.y - p0.y, p1.z - p0.z};
533 3
	Vector3 p20 = {p2.x - p0.x, p2.y - p0.y, p2.z - p0.z};
534

535
	// normal = cross(p1 - p0, p2 - p0)
536 3
	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};
537 3
	float area = normalize(normal);
538

539 3
	float distance = normal.x * p0.x + normal.y * p0.y + normal.z * p0.z;
540

541
	// we use sqrtf(area) so that the error is scaled linearly; this tends to improve silhouettes
542 3
	quadricFromPlane(Q, normal.x, normal.y, normal.z, -distance, sqrtf(area) * weight);
543
}
544

545 3
static void quadricFromTriangleEdge(Quadric& Q, const Vector3& p0, const Vector3& p1, const Vector3& p2, float weight)
546
{
547 3
	Vector3 p10 = {p1.x - p0.x, p1.y - p0.y, p1.z - p0.z};
548 3
	float length = normalize(p10);
549

550
	// p20p = length of projection of p2-p0 onto normalize(p1 - p0)
551 3
	Vector3 p20 = {p2.x - p0.x, p2.y - p0.y, p2.z - p0.z};
552 3
	float p20p = p20.x * p10.x + p20.y * p10.y + p20.z * p10.z;
553

554
	// normal = altitude of triangle from point p2 onto edge p1-p0
555 3
	Vector3 normal = {p20.x - p10.x * p20p, p20.y - p10.y * p20p, p20.z - p10.z * p20p};
556 3
	normalize(normal);
557

558 3
	float distance = normal.x * p0.x + normal.y * p0.y + normal.z * p0.z;
559

560
	// note: the weight is scaled linearly with edge length; this has to match the triangle weight
561 3
	quadricFromPlane(Q, normal.x, normal.y, normal.z, -distance, length * weight);
562
}
563

564 3
static void fillFaceQuadrics(Quadric* vertex_quadrics, const unsigned int* indices, size_t index_count, const Vector3* vertex_positions, const unsigned int* remap)
565
{
566 3
	for (size_t i = 0; i < index_count; i += 3)
567
	{
568 3
		unsigned int i0 = indices[i + 0];
569 3
		unsigned int i1 = indices[i + 1];
570 3
		unsigned int i2 = indices[i + 2];
571

572
		Quadric Q;
573 3
		quadricFromTriangle(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], 1.f);
574

575 3
		quadricAdd(vertex_quadrics[remap[i0]], Q);
576 3
		quadricAdd(vertex_quadrics[remap[i1]], Q);
577 3
		quadricAdd(vertex_quadrics[remap[i2]], Q);
578
	}
579
}
580

581 3
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, const unsigned int* loopback)
582
{
583 3
	for (size_t i = 0; i < index_count; i += 3)
584
	{
585
		static const int next[3] = {1, 2, 0};
586

587 3
		for (int e = 0; e < 3; ++e)
588
		{
589 3
			unsigned int i0 = indices[i + e];
590 3
			unsigned int i1 = indices[i + next[e]];
591

592 3
			unsigned char k0 = vertex_kind[i0];
593 3
			unsigned char k1 = vertex_kind[i1];
594

595
			// check that either i0 or i1 are border/seam and are on the same edge loop
596
			// note that we need to add the error even for edged that connect e.g. border & locked
597
			// if we don't do that, the adjacent border->border edge won't have correct errors for corners
598 3
			if (k0 != Kind_Border && k0 != Kind_Seam && k1 != Kind_Border && k1 != Kind_Seam)
599 3
				continue;
600

601 3
			if ((k0 == Kind_Border || k0 == Kind_Seam) && loop[i0] != i1)
602 3
				continue;
603

604 3
			if ((k1 == Kind_Border || k1 == Kind_Seam) && loopback[i1] != i0)
605 3
				continue;
606

607
			// seam edges should occur twice (i0->i1 and i1->i0) - skip redundant edges
608 3
			if (kHasOpposite[k0][k1] && remap[i1] > remap[i0])
609 3
				continue;
610

611 3
			unsigned int i2 = indices[i + next[next[e]]];
612

613
			// we try hard to maintain border edge geometry; seam edges can move more freely
614
			// due to topological restrictions on collapses, seam quadrics slightly improves collapse structure but aren't critical
615 3
			const float kEdgeWeightSeam = 1.f;
616 3
			const float kEdgeWeightBorder = 10.f;
617

618 3
			float edgeWeight = (k0 == Kind_Border || k1 == Kind_Border) ? kEdgeWeightBorder : kEdgeWeightSeam;
619

620
			Quadric Q;
621 3
			quadricFromTriangleEdge(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], edgeWeight);
622

623 3
			quadricAdd(vertex_quadrics[remap[i0]], Q);
624 3
			quadricAdd(vertex_quadrics[remap[i1]], Q);
625
		}
626
	}
627
}
628

629
// does triangle ABC flip when C is replaced with D?
630 3
static bool hasTriangleFlip(const Vector3& a, const Vector3& b, const Vector3& c, const Vector3& d)
631
{
632 3
	Vector3 eb = {b.x - a.x, b.y - a.y, b.z - a.z};
633 3
	Vector3 ec = {c.x - a.x, c.y - a.y, c.z - a.z};
634 3
	Vector3 ed = {d.x - a.x, d.y - a.y, d.z - a.z};
635

636 3
	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};
637 3
	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};
638

639 3
	return nbc.x * nbd.x + nbc.y * nbd.y + nbc.z * nbd.z < 0;
640
}
641

642 3
static bool hasTriangleFlips(const EdgeAdjacency& adjacency, const Vector3* vertex_positions, const unsigned int* collapse_remap, unsigned int i0, unsigned int i1)
643
{
644 1
	assert(collapse_remap[i0] == i0);
645 1
	assert(collapse_remap[i1] == i1);
646

647 3
	const Vector3& v0 = vertex_positions[i0];
648 3
	const Vector3& v1 = vertex_positions[i1];
649

650 3
	const EdgeAdjacency::Edge* edges = &adjacency.data[adjacency.offsets[i0]];
651 3
	size_t count = adjacency.counts[i0];
652

653 3
	for (size_t i = 0; i < count; ++i)
654
	{
655 3
		unsigned int a = collapse_remap[edges[i].next];
656 3
		unsigned int b = collapse_remap[edges[i].prev];
657

658
		// skip triangles that get collapsed
659
		// note: this is mathematically redundant as if either of these is true, the dot product in hasTriangleFlip should be 0
660 3
		if (a == i1 || b == i1)
661 3
			continue;
662

663
		// early-out when at least one triangle flips due to a collapse
664 3
		if (hasTriangleFlip(vertex_positions[a], vertex_positions[b], v0, v1))
665 3
			return true;
666
	}
667

668 3
	return false;
669
}
670

671 3
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)
672
{
673 3
	size_t collapse_count = 0;
674

675 3
	for (size_t i = 0; i < index_count; i += 3)
676
	{
677
		static const int next[3] = {1, 2, 0};
678

679 3
		for (int e = 0; e < 3; ++e)
680
		{
681 3
			unsigned int i0 = indices[i + e];
682 3
			unsigned int i1 = indices[i + next[e]];
683

684
			// this can happen either when input has a zero-length edge, or when we perform collapses for complex
685
			// topology w/seams and collapse a manifold vertex that connects to both wedges onto one of them
686
			// we leave edges like this alone since they may be important for preserving mesh integrity
687 1
			if (remap[i0] == remap[i1])
688 0
				continue;
689

690 3
			unsigned char k0 = vertex_kind[i0];
691 3
			unsigned char k1 = vertex_kind[i1];
692

693
			// the edge has to be collapsible in at least one direction
694 3
			if (!(kCanCollapse[k0][k1] | kCanCollapse[k1][k0]))
695 3
				continue;
696

697
			// manifold and seam edges should occur twice (i0->i1 and i1->i0) - skip redundant edges
698 3
			if (kHasOpposite[k0][k1] && remap[i1] > remap[i0])
699 3
				continue;
700

701
			// two vertices are on a border or a seam, but there's no direct edge between them
702
			// this indicates that they belong to two different edge loops and we should not collapse this edge
703
			// loop[] tracks half edges so we only need to check i0->i1
704 3
			if (k0 == k1 && (k0 == Kind_Border || k0 == Kind_Seam) && loop[i0] != i1)
705 3
				continue;
706

707
			// edge can be collapsed in either direction - we will pick the one with minimum error
708
			// note: we evaluate error later during collapse ranking, here we just tag the edge as bidirectional
709 3
			if (kCanCollapse[k0][k1] & kCanCollapse[k1][k0])
710
			{
711 3
				Collapse c = {i0, i1, {/* bidi= */ 1}};
712 3
				collapses[collapse_count++] = c;
713
			}
714
			else
715
			{
716
				// edge can only be collapsed in one direction
717 3
				unsigned int e0 = kCanCollapse[k0][k1] ? i0 : i1;
718 3
				unsigned int e1 = kCanCollapse[k0][k1] ? i1 : i0;
719

720 3
				Collapse c = {e0, e1, {/* bidi= */ 0}};
721 3
				collapses[collapse_count++] = c;
722
			}
723
		}
724
	}
725

726 3
	return collapse_count;
727
}
728

729 3
static void rankEdgeCollapses(Collapse* collapses, size_t collapse_count, const Vector3* vertex_positions, const Quadric* vertex_quadrics, const unsigned int* remap)
730
{
731 3
	for (size_t i = 0; i < collapse_count; ++i)
732
	{
733 3
		Collapse& c = collapses[i];
734

735 3
		unsigned int i0 = c.v0;
736 3
		unsigned int i1 = c.v1;
737

738
		// most edges are bidirectional which means we need to evaluate errors for two collapses
739
		// to keep this code branchless we just use the same edge for unidirectional edges
740 3
		unsigned int j0 = c.bidi ? i1 : i0;
741 3
		unsigned int j1 = c.bidi ? i0 : i1;
742

743 3
		const Quadric& qi = vertex_quadrics[remap[i0]];
744 3
		const Quadric& qj = vertex_quadrics[remap[j0]];
745

746 3
		float ei = quadricError(qi, vertex_positions[i1]);
747 3
		float ej = quadricError(qj, vertex_positions[j1]);
748

749
		// pick edge direction with minimal error
750 3
		c.v0 = ei <= ej ? i0 : j0;
751 3
		c.v1 = ei <= ej ? i1 : j1;
752 3
		c.error = ei <= ej ? ei : ej;
753
	}
754
}
755

756
#if TRACE > 1
757
static void dumpEdgeCollapses(const Collapse* collapses, size_t collapse_count, const unsigned char* vertex_kind)
758
{
759
	size_t ckinds[Kind_Count][Kind_Count] = {};
760
	float cerrors[Kind_Count][Kind_Count] = {};
761

762
	for (int k0 = 0; k0 < Kind_Count; ++k0)
763
		for (int k1 = 0; k1 < Kind_Count; ++k1)
764
			cerrors[k0][k1] = FLT_MAX;
765

766
	for (size_t i = 0; i < collapse_count; ++i)
767
	{
768
		unsigned int i0 = collapses[i].v0;
769
		unsigned int i1 = collapses[i].v1;
770

771
		unsigned char k0 = vertex_kind[i0];
772
		unsigned char k1 = vertex_kind[i1];
773

774
		ckinds[k0][k1]++;
775
		cerrors[k0][k1] = (collapses[i].error < cerrors[k0][k1]) ? collapses[i].error : cerrors[k0][k1];
776
	}
777

778
	for (int k0 = 0; k0 < Kind_Count; ++k0)
779
		for (int k1 = 0; k1 < Kind_Count; ++k1)
780
			if (ckinds[k0][k1])
781
				printf("collapses %d -> %d: %d, min error %e\n", k0, k1, int(ckinds[k0][k1]), ckinds[k0][k1] ? sqrtf(cerrors[k0][k1]) : 0.f);
782
}
783

784
static void dumpLockedCollapses(const unsigned int* indices, size_t index_count, const unsigned char* vertex_kind)
785
{
786
	size_t locked_collapses[Kind_Count][Kind_Count] = {};
787

788
	for (size_t i = 0; i < index_count; i += 3)
789
	{
790
		static const int next[3] = {1, 2, 0};
791

792
		for (int e = 0; e < 3; ++e)
793
		{
794
			unsigned int i0 = indices[i + e];
795
			unsigned int i1 = indices[i + next[e]];
796

797
			unsigned char k0 = vertex_kind[i0];
798
			unsigned char k1 = vertex_kind[i1];
799

800
			locked_collapses[k0][k1] += !kCanCollapse[k0][k1] && !kCanCollapse[k1][k0];
801
		}
802
	}
803

804
	for (int k0 = 0; k0 < Kind_Count; ++k0)
805
		for (int k1 = 0; k1 < Kind_Count; ++k1)
806
			if (locked_collapses[k0][k1])
807
				printf("locked collapses %d -> %d: %d\n", k0, k1, int(locked_collapses[k0][k1]));
808
}
809
#endif
810

811 3
static void sortEdgeCollapses(unsigned int* sort_order, const Collapse* collapses, size_t collapse_count)
812
{
813 3
	const int sort_bits = 11;
814

815
	// fill histogram for counting sort
816
	unsigned int histogram[1 << sort_bits];
817 3
	memset(histogram, 0, sizeof(histogram));
818

819 3
	for (size_t i = 0; i < collapse_count; ++i)
820
	{
821
		// skip sign bit since error is non-negative
822 3
		unsigned int key = (collapses[i].errorui << 1) >> (32 - sort_bits);
823

824 3
		histogram[key]++;
825
	}
826

827
	// compute offsets based on histogram data
828 3
	size_t histogram_sum = 0;
829

830 3
	for (size_t i = 0; i < 1 << sort_bits; ++i)
831
	{
832 3
		size_t count = histogram[i];
833 3
		histogram[i] = unsigned(histogram_sum);
834 3
		histogram_sum += count;
835
	}
836

837 1
	assert(histogram_sum == collapse_count);
838

839
	// compute sort order based on offsets
840 3
	for (size_t i = 0; i < collapse_count; ++i)
841
	{
842
		// skip sign bit since error is non-negative
843 3
		unsigned int key = (collapses[i].errorui << 1) >> (32 - sort_bits);
844

845 3
		sort_order[histogram[key]++] = unsigned(i);
846
	}
847
}
848

849 3
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, const Vector3* vertex_positions, const EdgeAdjacency& adjacency, size_t triangle_collapse_goal, float error_limit, float& result_error)
850
{
851 3
	size_t edge_collapses = 0;
852 3
	size_t triangle_collapses = 0;
853

854
	// most collapses remove 2 triangles; use this to establish a bound on the pass in terms of error limit
855
	// note that edge_collapse_goal is an estimate; triangle_collapse_goal will be used to actually limit collapses
856 3
	size_t edge_collapse_goal = triangle_collapse_goal / 2;
857

858
#if TRACE
859
	size_t stats[4] = {};
860
#endif
861

862 3
	for (size_t i = 0; i < collapse_count; ++i)
863
	{
864 3
		const Collapse& c = collapses[collapse_order[i]];
865

866
		TRACESTATS(0);
867

868 3
		if (c.error > error_limit)
869 3
			break;
870

871 3
		if (triangle_collapses >= triangle_collapse_goal)
872 3
			break;
873

874
		// we limit the error in each pass based on the error of optimal last collapse; since many collapses will be locked
875
		// as they will share vertices with other successfull collapses, we need to increase the acceptable error by some factor
876 1
		float error_goal = edge_collapse_goal < collapse_count ? 1.5f * collapses[collapse_order[edge_collapse_goal]].error : FLT_MAX;
877

878
		// on average, each collapse is expected to lock 6 other collapses; to avoid degenerate passes on meshes with odd
879
		// topology, we only abort if we got over 1/6 collapses accordingly.
880 3
		if (c.error > error_goal && triangle_collapses > triangle_collapse_goal / 6)
881 3
			break;
882

883 3
		unsigned int i0 = c.v0;
884 3
		unsigned int i1 = c.v1;
885

886 3
		unsigned int r0 = remap[i0];
887 3
		unsigned int r1 = remap[i1];
888

889
		// we don't collapse vertices that had source or target vertex involved in a collapse
890
		// it's important to not move the vertices twice since it complicates the tracking/remapping logic
891
		// it's important to not move other vertices towards a moved vertex to preserve error since we don't re-rank collapses mid-pass
892 3
		if (collapse_locked[r0] | collapse_locked[r1])
893
		{
894
			TRACESTATS(1);
895 3
			continue;
896
		}
897

898 3
		if (hasTriangleFlips(adjacency, vertex_positions, collapse_remap, r0, r1))
899
		{
900
			// adjust collapse goal since this collapse is invalid and shouldn't factor into error goal
901 3
			edge_collapse_goal++;
902

903
			TRACESTATS(2);
904 3
			continue;
905
		}
906

907 1
		assert(collapse_remap[r0] == r0);
908 1
		assert(collapse_remap[r1] == r1);
909

910 3
		quadricAdd(vertex_quadrics[r1], vertex_quadrics[r0]);
911

912 3
		if (vertex_kind[i0] == Kind_Complex)
913
		{
914 0
			unsigned int v = i0;
915

916 0
			do
917
			{
918 0
				collapse_remap[v] = r1;
919 0
				v = wedge[v];
920 0
			} while (v != i0);
921
		}
922 3
		else if (vertex_kind[i0] == Kind_Seam)
923
		{
924
			// remap v0 to v1 and seam pair of v0 to seam pair of v1
925 3
			unsigned int s0 = wedge[i0];
926 3
			unsigned int s1 = wedge[i1];
927

928 1
			assert(s0 != i0 && s1 != i1);
929 1
			assert(wedge[s0] == i0 && wedge[s1] == i1);
930

931 3
			collapse_remap[i0] = i1;
932 3
			collapse_remap[s0] = s1;
933
		}
934
		else
935
		{
936 1
			assert(wedge[i0] == i0);
937

938 3
			collapse_remap[i0] = i1;
939
		}
940

941 3
		collapse_locked[r0] = 1;
942 3
		collapse_locked[r1] = 1;
943

944
		// border edges collapse 1 triangle, other edges collapse 2 or more
945 3
		triangle_collapses += (vertex_kind[i0] == Kind_Border) ? 1 : 2;
946 3
		edge_collapses++;
947

948 3
		result_error = result_error < c.error ? c.error : result_error;
949
	}
950

951
#if TRACE
952
	float error_goal_perfect = edge_collapse_goal < collapse_count ? collapses[collapse_order[edge_collapse_goal]].error : 0.f;
953

954
	printf("removed %d triangles, error %e (goal %e); evaluated %d/%d collapses (done %d, skipped %d, invalid %d)\n",
955
	    int(triangle_collapses), sqrtf(result_error), sqrtf(error_goal_perfect),
956
	    int(stats[0]), int(collapse_count), int(edge_collapses), int(stats[1]), int(stats[2]));
957
#endif
958

959 3
	return edge_collapses;
960
}
961

962 3
static size_t remapIndexBuffer(unsigned int* indices, size_t index_count, const unsigned int* collapse_remap)
963
{
964 3
	size_t write = 0;
965

966 3
	for (size_t i = 0; i < index_count; i += 3)
967
	{
968 3
		unsigned int v0 = collapse_remap[indices[i + 0]];
969 3
		unsigned int v1 = collapse_remap[indices[i + 1]];
970 3
		unsigned int v2 = collapse_remap[indices[i + 2]];
971

972
		// we never move the vertex twice during a single pass
973 1
		assert(collapse_remap[v0] == v0);
974 1
		assert(collapse_remap[v1] == v1);
975 1
		assert(collapse_remap[v2] == v2);
976

977 3
		if (v0 != v1 && v0 != v2 && v1 != v2)
978
		{
979 3
			indices[write + 0] = v0;
980 3
			indices[write + 1] = v1;
981 3
			indices[write + 2] = v2;
982 3
			write += 3;
983
		}
984
	}
985

986 3
	return write;
987
}
988

989 3
static void remapEdgeLoops(unsigned int* loop, size_t vertex_count, const unsigned int* collapse_remap)
990
{
991 3
	for (size_t i = 0; i < vertex_count; ++i)
992
	{
993 3
		if (loop[i] != ~0u)
994
		{
995 3
			unsigned int l = loop[i];
996 3
			unsigned int r = collapse_remap[l];
997

998
			// i == r is a special case when the seam edge is collapsed in a direction opposite to where loop goes
999 3
			loop[i] = (i == r) ? loop[l] : r;
1000
		}
1001
	}
1002
}
1003

1004
struct CellHasher
1005
{
1006
	const unsigned int* vertex_ids;
1007

1008 3
	size_t hash(unsigned int i) const
1009
	{
1010 3
		unsigned int h = vertex_ids[i];
1011

1012
		// MurmurHash2 finalizer
1013 3
		h ^= h >> 13;
1014 3
		h *= 0x5bd1e995;
1015 3
		h ^= h >> 15;
1016 3
		return h;
1017
	}
1018

1019 3
	bool equal(unsigned int lhs, unsigned int rhs) const
1020
	{
1021 3
		return vertex_ids[lhs] == vertex_ids[rhs];
1022
	}
1023
};
1024

1025
struct IdHasher
1026
{
1027 3
	size_t hash(unsigned int id) const
1028
	{
1029 3
		unsigned int h = id;
1030

1031
		// MurmurHash2 finalizer
1032 3
		h ^= h >> 13;
1033 3
		h *= 0x5bd1e995;
1034 3
		h ^= h >> 15;
1035 3
		return h;
1036
	}
1037

1038 3
	bool equal(unsigned int lhs, unsigned int rhs) const
1039
	{
1040 3
		return lhs == rhs;
1041
	}
1042
};
1043

1044
struct TriangleHasher
1045
{
1046
	const unsigned int* indices;
1047

1048 3
	size_t hash(unsigned int i) const
1049
	{
1050 3
		const unsigned int* tri = indices + i * 3;
1051

1052
		// Optimized Spatial Hashing for Collision Detection of Deformable Objects
1053 3
		return (tri[0] * 73856093) ^ (tri[1] * 19349663) ^ (tri[2] * 83492791);
1054
	}
1055

1056 3
	bool equal(unsigned int lhs, unsigned int rhs) const
1057
	{
1058 3
		const unsigned int* lt = indices + lhs * 3;
1059 3
		const unsigned int* rt = indices + rhs * 3;
1060

1061 3
		return lt[0] == rt[0] && lt[1] == rt[1] && lt[2] == rt[2];
1062
	}
1063
};
1064

1065 3
static void computeVertexIds(unsigned int* vertex_ids, const Vector3* vertex_positions, size_t vertex_count, int grid_size)
1066
{
1067 1
	assert(grid_size >= 1 && grid_size <= 1024);
1068 3
	float cell_scale = float(grid_size - 1);
1069

1070 3
	for (size_t i = 0; i < vertex_count; ++i)
1071
	{
1072 3
		const Vector3& v = vertex_positions[i];
1073

1074 3
		int xi = int(v.x * cell_scale + 0.5f);
1075 3
		int yi = int(v.y * cell_scale + 0.5f);
1076 3
		int zi = int(v.z * cell_scale + 0.5f);
1077

1078 3
		vertex_ids[i] = (xi << 20) | (yi << 10) | zi;
1079
	}
1080
}
1081

1082 3
static size_t countTriangles(const unsigned int* vertex_ids, const unsigned int* indices, size_t index_count)
1083
{
1084 3
	size_t result = 0;
1085

1086 3
	for (size_t i = 0; i < index_count; i += 3)
1087
	{
1088 3
		unsigned int id0 = vertex_ids[indices[i + 0]];
1089 3
		unsigned int id1 = vertex_ids[indices[i + 1]];
1090 3
		unsigned int id2 = vertex_ids[indices[i + 2]];
1091

1092 3
		result += (id0 != id1) & (id0 != id2) & (id1 != id2);
1093
	}
1094

1095 3
	return result;
1096
}
1097

1098 3
static size_t fillVertexCells(unsigned int* table, size_t table_size, unsigned int* vertex_cells, const unsigned int* vertex_ids, size_t vertex_count)
1099
{
1100 3
	CellHasher hasher = {vertex_ids};
1101

1102 3
	memset(table, -1, table_size * sizeof(unsigned int));
1103

1104 3
	size_t result = 0;
1105

1106 3
	for (size_t i = 0; i < vertex_count; ++i)
1107
	{
1108 3
		unsigned int* entry = hashLookup2(table, table_size, hasher, unsigned(i), ~0u);
1109

1110 3
		if (*entry == ~0u)
1111
		{
1112 3
			*entry = unsigned(i);
1113 3
			vertex_cells[i] = unsigned(result++);
1114
		}
1115
		else
1116
		{
1117 3
			vertex_cells[i] = vertex_cells[*entry];
1118
		}
1119
	}
1120

1121 3
	return result;
1122
}
1123

1124 3
static size_t countVertexCells(unsigned int* table, size_t table_size, const unsigned int* vertex_ids, size_t vertex_count)
1125
{
1126
	IdHasher hasher;
1127

1128 3
	memset(table, -1, table_size * sizeof(unsigned int));
1129

1130 3
	size_t result = 0;
1131

1132 3
	for (size_t i = 0; i < vertex_count; ++i)
1133
	{
1134 3
		unsigned int id = vertex_ids[i];
1135 3
		unsigned int* entry = hashLookup2(table, table_size, hasher, id, ~0u);
1136

1137 3
		result += (*entry == ~0u);
1138 3
		*entry = id;
1139
	}
1140

1141 3
	return result;
1142
}
1143

1144 3
static void fillCellQuadrics(Quadric* cell_quadrics, const unsigned int* indices, size_t index_count, const Vector3* vertex_positions, const unsigned int* vertex_cells)
1145
{
1146 3
	for (size_t i = 0; i < index_count; i += 3)
1147
	{
1148 3
		unsigned int i0 = indices[i + 0];
1149 3
		unsigned int i1 = indices[i + 1];
1150 3
		unsigned int i2 = indices[i + 2];
1151

1152 3
		unsigned int c0 = vertex_cells[i0];
1153 3
		unsigned int c1 = vertex_cells[i1];
1154 3
		unsigned int c2 = vertex_cells[i2];
1155

1156 3
		bool single_cell = (c0 == c1) & (c0 == c2);
1157

1158
		Quadric Q;
1159 3
		quadricFromTriangle(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], single_cell ? 3.f : 1.f);
1160

1161 3
		if (single_cell)
1162
		{
1163 3
			quadricAdd(cell_quadrics[c0], Q);
1164
		}
1165
		else
1166
		{
1167 3
			quadricAdd(cell_quadrics[c0], Q);
1168 3
			quadricAdd(cell_quadrics[c1], Q);
1169 3
			quadricAdd(cell_quadrics[c2], Q);
1170
		}
1171
	}
1172
}
1173

1174 3
static void fillCellQuadrics(Quadric* cell_quadrics, const Vector3* vertex_positions, size_t vertex_count, const unsigned int* vertex_cells)
1175
{
1176 3
	for (size_t i = 0; i < vertex_count; ++i)
1177
	{
1178 3
		unsigned int c = vertex_cells[i];
1179 3
		const Vector3& v = vertex_positions[i];
1180

1181
		Quadric Q;
1182 3
		quadricFromPoint(Q, v.x, v.y, v.z, 1.f);
1183

1184 3
		quadricAdd(cell_quadrics[c], Q);
1185
	}
1186
}
1187

1188 3
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)
1189
{
1190 3
	memset(cell_remap, -1, cell_count * sizeof(unsigned int));
1191

1192 3
	for (size_t i = 0; i < vertex_count; ++i)
1193
	{
1194 3
		unsigned int cell = vertex_cells[i];
1195 3
		float error = quadricError(cell_quadrics[cell], vertex_positions[i]);
1196

1197 3
		if (cell_remap[cell] == ~0u || cell_errors[cell] > error)
1198
		{
1199 3
			cell_remap[cell] = unsigned(i);
1200 3
			cell_errors[cell] = error;
1201
		}
1202
	}
1203
}
1204

1205 3
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)
1206
{
1207 3
	TriangleHasher hasher = {destination};
1208

1209 3
	memset(tritable, -1, tritable_size * sizeof(unsigned int));
1210

1211 3
	size_t result = 0;
1212

1213 3
	for (size_t i = 0; i < index_count; i += 3)
1214
	{
1215 3
		unsigned int c0 = vertex_cells[indices[i + 0]];
1216 3
		unsigned int c1 = vertex_cells[indices[i + 1]];
1217 3
		unsigned int c2 = vertex_cells[indices[i + 2]];
1218

1219 3
		if (c0 != c1 && c0 != c2 && c1 != c2)
1220
		{
1221 3
			unsigned int a = cell_remap[c0];
1222 3
			unsigned int b = cell_remap[c1];
1223 3
			unsigned int c = cell_remap[c2];
1224

1225 3
			if (b < a && b < c)
1226
			{
1227 3
				unsigned int t = a;
1228 3
				a = b, b = c, c = t;
1229
			}
1230 3
			else if (c < a && c < b)
1231
			{
1232 3
				unsigned int t = c;
1233 3
				c = b, b = a, a = t;
1234
			}
1235

1236 3
			destination[result * 3 + 0] = a;
1237 3
			destination[result * 3 + 1] = b;
1238 3
			destination[result * 3 + 2] = c;
1239

1240 3
			unsigned int* entry = hashLookup2(tritable, tritable_size, hasher, unsigned(result), ~0u);
1241

1242 3
			if (*entry == ~0u)
1243 3
				*entry = unsigned(result++);
1244
		}
1245
	}
1246

1247 3
	return result * 3;
1248
}
1249

1250 3
static float interpolate(float y, float x0, float y0, float x1, float y1, float x2, float y2)
1251
{
1252
	// three point interpolation from "revenge of interpolation search" paper
1253 3
	float num = (y1 - y) * (x1 - x2) * (x1 - x0) * (y2 - y0);
1254 3
	float den = (y2 - y) * (x1 - x2) * (y0 - y1) + (y0 - y) * (x1 - x0) * (y1 - y2);
1255 3
	return x1 + num / den;
1256
}
1257

1258
} // namespace meshopt
1259

1260
#ifndef NDEBUG
1261
// Note: this is only exposed for debug visualization purposes; do *not* use these in debug builds
1262
MESHOPTIMIZER_API unsigned char* meshopt_simplifyDebugKind = 0;
1263
MESHOPTIMIZER_API unsigned int* meshopt_simplifyDebugLoop = 0;
1264
MESHOPTIMIZER_API unsigned int* meshopt_simplifyDebugLoopBack = 0;
1265
#endif
1266

1267 3
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, float* out_result_error)
1268
{
1269
	using namespace meshopt;
1270

1271 1
	assert(index_count % 3 == 0);
1272 1
	assert(vertex_positions_stride > 0 && vertex_positions_stride <= 256);
1273 1
	assert(vertex_positions_stride % sizeof(float) == 0);
1274 1
	assert(target_index_count <= index_count);
1275

1276 3
	meshopt_Allocator allocator;
1277

1278 3
	unsigned int* result = destination;
1279

1280
	// build adjacency information
1281 3
	EdgeAdjacency adjacency = {};
1282 3
	prepareEdgeAdjacency(adjacency, index_count, vertex_count, allocator);
1283 3
	updateEdgeAdjacency(adjacency, indices, index_count, vertex_count, NULL);
1284

1285
	// build position remap that maps each vertex to the one with identical position
1286 3
	unsigned int* remap = allocator.allocate<unsigned int>(vertex_count);
1287 3
	unsigned int* wedge = allocator.allocate<unsigned int>(vertex_count);
1288 3
	buildPositionRemap(remap, wedge, vertex_positions_data, vertex_count, vertex_positions_stride, allocator);
1289

1290
	// classify vertices; vertex kind determines collapse rules, see kCanCollapse
1291 3
	unsigned char* vertex_kind = allocator.allocate<unsigned char>(vertex_count);
1292 3
	unsigned int* loop = allocator.allocate<unsigned int>(vertex_count);
1293 3
	unsigned int* loopback = allocator.allocate<unsigned int>(vertex_count);
1294 3
	classifyVertices(vertex_kind, loop, loopback, vertex_count, adjacency, remap, wedge);
1295

1296
#if TRACE
1297
	size_t unique_positions = 0;
1298
	for (size_t i = 0; i < vertex_count; ++i)
1299
		unique_positions += remap[i] == i;
1300

1301
	printf("position remap: %d vertices => %d positions\n", int(vertex_count), int(unique_positions));
1302

1303
	size_t kinds[Kind_Count] = {};
1304
	for (size_t i = 0; i < vertex_count; ++i)
1305
		kinds[vertex_kind[i]] += remap[i] == i;
1306

1307
	printf("kinds: manifold %d, border %d, seam %d, complex %d, locked %d\n",
1308
	    int(kinds[Kind_Manifold]), int(kinds[Kind_Border]), int(kinds[Kind_Seam]), int(kinds[Kind_Complex]), int(kinds[Kind_Locked]));
1309
#endif
1310

1311 3
	Vector3* vertex_positions = allocator.allocate<Vector3>(vertex_count);
1312 3
	rescalePositions(vertex_positions, vertex_positions_data, vertex_count, vertex_positions_stride);
1313

1314 3
	Quadric* vertex_quadrics = allocator.allocate<Quadric>(vertex_count);
1315 3
	memset(vertex_quadrics, 0, vertex_count * sizeof(Quadric));
1316

1317 3
	fillFaceQuadrics(vertex_quadrics, indices, index_count, vertex_positions, remap);
1318 3
	fillEdgeQuadrics(vertex_quadrics, indices, index_count, vertex_positions, remap, vertex_kind, loop, loopback);
1319

1320 3
	if (result != indices)
1321 3
		memcpy(result, indices, index_count * sizeof(unsigned int));
1322

1323
#if TRACE
1324
	size_t pass_count = 0;
1325
#endif
1326

1327 3
	Collapse* edge_collapses = allocator.allocate<Collapse>(index_count);
1328 3
	unsigned int* collapse_order = allocator.allocate<unsigned int>(index_count);
1329 3
	unsigned int* collapse_remap = allocator.allocate<unsigned int>(vertex_count);
1330 3
	unsigned char* collapse_locked = allocator.allocate<unsigned char>(vertex_count);
1331

1332 3
	size_t result_count = index_count;
1333 3
	float result_error = 0;
1334

1335
	// target_error input is linear; we need to adjust it to match quadricError units
1336 3
	float error_limit = target_error * target_error;
1337

1338 3
	while (result_count > target_index_count)
1339
	{
1340
		// note: throughout the simplification process adjacency structure reflects welded topology for result-in-progress
1341 3
		updateEdgeAdjacency(adjacency, result, result_count, vertex_count, remap);
1342

1343 3
		size_t edge_collapse_count = pickEdgeCollapses(edge_collapses, result, result_count, remap, vertex_kind, loop);
1344

1345
		// no edges can be collapsed any more due to topology restrictions
1346 3
		if (edge_collapse_count == 0)
1347 3
			break;
1348

1349 3
		rankEdgeCollapses(edge_collapses, edge_collapse_count, vertex_positions, vertex_quadrics, remap);
1350

1351
#if TRACE > 1
1352
		dumpEdgeCollapses(edge_collapses, edge_collapse_count, vertex_kind);
1353
#endif
1354

1355 3
		sortEdgeCollapses(collapse_order, edge_collapses, edge_collapse_count);
1356

1357 3
		size_t triangle_collapse_goal = (result_count - target_index_count) / 3;
1358

1359 3
		for (size_t i = 0; i < vertex_count; ++i)
1360 3
			collapse_remap[i] = unsigned(i);
1361

1362 3
		memset(collapse_locked, 0, vertex_count);
1363

1364
#if TRACE
1365
		printf("pass %d: ", int(pass_count++));
1366
#endif
1367

1368 3
		size_t collapses = performEdgeCollapses(collapse_remap, collapse_locked, vertex_quadrics, edge_collapses, edge_collapse_count, collapse_order, remap, wedge, vertex_kind, vertex_positions, adjacency, triangle_collapse_goal, error_limit, result_error);
1369

1370
		// no edges can be collapsed any more due to hitting the error limit or triangle collapse limit
1371 3
		if (collapses == 0)
1372 3
			break;
1373

1374 3
		remapEdgeLoops(loop, vertex_count, collapse_remap);
1375 3
		remapEdgeLoops(loopback, vertex_count, collapse_remap);
1376

1377 3
		size_t new_count = remapIndexBuffer(result, result_count, collapse_remap);
1378 1
		assert(new_count < result_count);
1379

1380 3
		result_count = new_count;
1381
	}
1382

1383
#if TRACE
1384
	printf("result: %d triangles, error: %e; total %d passes\n", int(result_count), sqrtf(result_error), int(pass_count));
1385
#endif
1386

1387
#if TRACE > 1
1388
	dumpLockedCollapses(result, result_count, vertex_kind);
1389
#endif
1390

1391
#ifndef NDEBUG
1392 3
	if (meshopt_simplifyDebugKind)
1393 0
		memcpy(meshopt_simplifyDebugKind, vertex_kind, vertex_count);
1394

1395 3
	if (meshopt_simplifyDebugLoop)
1396 0
		memcpy(meshopt_simplifyDebugLoop, loop, vertex_count * sizeof(unsigned int));
1397

1398 3
	if (meshopt_simplifyDebugLoopBack)
1399 0
		memcpy(meshopt_simplifyDebugLoopBack, loopback, vertex_count * sizeof(unsigned int));
1400
#endif
1401

1402
	// result_error is quadratic; we need to remap it back to linear
1403 3
	if (out_result_error)
1404 3
		*out_result_error = sqrtf(result_error);
1405

1406 3
	return result_count;
1407
}
1408

1409 3
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, float target_error, float* out_result_error)
1410
{
1411
	using namespace meshopt;
1412

1413 1
	assert(index_count % 3 == 0);
1414 1
	assert(vertex_positions_stride > 0 && vertex_positions_stride <= 256);
1415 1
	assert(vertex_positions_stride % sizeof(float) == 0);
1416 1
	assert(target_index_count <= index_count);
1417

1418
	// we expect to get ~2 triangles/vertex in the output
1419 3
	size_t target_cell_count = target_index_count / 6;
1420

1421 3
	meshopt_Allocator allocator;
1422

1423 3
	Vector3* vertex_positions = allocator.allocate<Vector3>(vertex_count);
1424 3
	rescalePositions(vertex_positions, vertex_positions_data, vertex_count, vertex_positions_stride);
1425

1426
	// find the optimal grid size using guided binary search
1427
#if TRACE
1428
	printf("source: %d vertices, %d triangles\n", int(vertex_count), int(index_count / 3));
1429
	printf("target: %d cells, %d triangles\n", int(target_cell_count), int(target_index_count / 3));
1430
#endif
1431

1432 3
	unsigned int* vertex_ids = allocator.allocate<unsigned int>(vertex_count);
1433

1434 3
	const int kInterpolationPasses = 5;
1435

1436
	// invariant: # of triangles in min_grid <= target_count
1437 3
	int min_grid = int(1.f / (target_error < 1e-3f ? 1e-3f : target_error));
1438 3
	int max_grid = 1025;
1439 3
	size_t min_triangles = 0;
1440 3
	size_t max_triangles = index_count / 3;
1441

1442
	// 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
1443 3
	if (min_grid > 1)
1444
	{
1445 3
		computeVertexIds(vertex_ids, vertex_positions, vertex_count, min_grid);
1446 3
		min_triangles = countTriangles(vertex_ids, indices, index_count);
1447
	}
1448

1449
	// 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...
1450 3
	int next_grid_size = int(sqrtf(float(target_cell_count)) + 0.5f);
1451

1452 3
	for (int pass = 0; pass < 10 + kInterpolationPasses; ++pass)
1453
	{
1454 3
		if (min_triangles >= target_index_count / 3 || max_grid - min_grid <= 1)
1455 1
			break;
1456

1457
		// we clamp the prediction of the grid size to make sure that the search converges
1458 3
		int grid_size = next_grid_size;
1459 3
		grid_size = (grid_size <= min_grid) ? min_grid + 1 : (grid_size >= max_grid) ? max_grid - 1 : grid_size;
1460

1461 3
		computeVertexIds(vertex_ids, vertex_positions, vertex_count, grid_size);
1462 3
		size_t triangles = countTriangles(vertex_ids, indices, index_count);
1463

1464
#if TRACE
1465
		printf("pass %d (%s): grid size %d, triangles %d, %s\n",
1466
		    pass, (pass == 0) ? "guess" : (pass <= kInterpolationPasses) ? "lerp" : "binary",
1467
		    grid_size, int(triangles),
1468
		    (triangles <= target_index_count / 3) ? "under" : "over");
1469
#endif
1470

1471 3
		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));
1472

1473 3
		if (triangles <= target_index_count / 3)
1474
		{
1475 3
			min_grid = grid_size;
1476 3
			min_triangles = triangles;
1477
		}
1478
		else
1479
		{
1480 3
			max_grid = grid_size;
1481 3
			max_triangles = triangles;
1482
		}
1483

1484
		// we start by using interpolation search - it usually converges faster
1485
		// 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)
1486 1
		next_grid_size = (pass < kInterpolationPasses) ? int(tip + 0.5f) : (min_grid + max_grid) / 2;
1487
	}
1488

1489 3
	if (min_triangles == 0)
1490
	{
1491 3
		if (out_result_error)
1492 0
			*out_result_error = 1.f;
1493

1494 3
		return 0;
1495
	}
1496

1497
	// build vertex->cell association by mapping all vertices with the same quantized position to the same cell
1498 3
	size_t table_size = hashBuckets2(vertex_count);
1499 3
	unsigned int* table = allocator.allocate<unsigned int>(table_size);
1500

1501 3
	unsigned int* vertex_cells = allocator.allocate<unsigned int>(vertex_count);
1502

1503 3
	computeVertexIds(vertex_ids, vertex_positions, vertex_count, min_grid);
1504 3
	size_t cell_count = fillVertexCells(table, table_size, vertex_cells, vertex_ids, vertex_count);
1505

1506
	// build a quadric for each target cell
1507 3
	Quadric* cell_quadrics = allocator.allocate<Quadric>(cell_count);
1508 3
	memset(cell_quadrics, 0, cell_count * sizeof(Quadric));
1509

1510 3
	fillCellQuadrics(cell_quadrics, indices, index_count, vertex_positions, vertex_cells);
1511

1512
	// for each target cell, find the vertex with the minimal error
1513 3
	unsigned int* cell_remap = allocator.allocate<unsigned int>(cell_count);
1514 3
	float* cell_errors = allocator.allocate<float>(cell_count);
1515

1516 3
	fillCellRemap(cell_remap, cell_errors, cell_count, vertex_cells, cell_quadrics, vertex_positions, vertex_count);
1517

1518
	// compute error
1519 3
	float result_error = 0.f;
1520

1521 3
	for (size_t i = 0; i < cell_count; ++i)
1522 3
		result_error = result_error < cell_errors[i] ? cell_errors[i] : result_error;
1523

1524
	// collapse triangles!
1525
	// note that we need to filter out triangles that we've already output because we very frequently generate redundant triangles between cells :(
1526 3
	size_t tritable_size = hashBuckets2(min_triangles);
1527 3
	unsigned int* tritable = allocator.allocate<unsigned int>(tritable_size);
1528

1529 3
	size_t write = filterTriangles(destination, tritable, tritable_size, indices, index_count, vertex_cells, cell_remap);
1530

1531
#if TRACE
1532
	printf("result: %d cells, %d triangles (%d unfiltered), error %e\n", int(cell_count), int(write / 3), int(min_triangles), sqrtf(result_error));
1533
#endif
1534

1535 3
	if (out_result_error)
1536 3
		*out_result_error = sqrtf(result_error);
1537

1538 3
	return write;
1539
}
1540

1541 3
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)
1542
{
1543
	using namespace meshopt;
1544

1545 1
	assert(vertex_positions_stride > 0 && vertex_positions_stride <= 256);
1546 1
	assert(vertex_positions_stride % sizeof(float) == 0);
1547 1
	assert(target_vertex_count <= vertex_count);
1548

1549 3
	size_t target_cell_count = target_vertex_count;
1550

1551 3
	if (target_cell_count == 0)
1552 3
		return 0;
1553

1554 3
	meshopt_Allocator allocator;
1555

1556 3
	Vector3* vertex_positions = allocator.allocate<Vector3>(vertex_count);
1557 3
	rescalePositions(vertex_positions, vertex_positions_data, vertex_count, vertex_positions_stride);
1558

1559
	// find the optimal grid size using guided binary search
1560
#if TRACE
1561
	printf("source: %d vertices\n", int(vertex_count));
1562
	printf("target: %d cells\n", int(target_cell_count));
1563
#endif
1564

1565 3
	unsigned int* vertex_ids = allocator.allocate<unsigned int>(vertex_count);
1566

1567 3
	size_t table_size = hashBuckets2(vertex_count);
1568 3
	unsigned int* table = allocator.allocate<unsigned int>(table_size);
1569

1570 3
	const int kInterpolationPasses = 5;
1571

1572
	// invariant: # of vertices in min_grid <= target_count
1573 3
	int min_grid = 0;
1574 3
	int max_grid = 1025;
1575 3
	size_t min_vertices = 0;
1576 3
	size_t max_vertices = vertex_count;
1577

1578
	// 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...
1579 3
	int next_grid_size = int(sqrtf(float(target_cell_count)) + 0.5f);
1580

1581 3
	for (int pass = 0; pass < 10 + kInterpolationPasses; ++pass)
1582
	{
1583 1
		assert(min_vertices < target_vertex_count);
1584 1
		assert(max_grid - min_grid > 1);
1585

1586
		// we clamp the prediction of the grid size to make sure that the search converges
1587 3
		int grid_size = next_grid_size;
1588 1
		grid_size = (grid_size <= min_grid) ? min_grid + 1 : (grid_size >= max_grid) ? max_grid - 1 : grid_size;
1589

1590 3
		computeVertexIds(vertex_ids, vertex_positions, vertex_count, grid_size);
1591 3
		size_t vertices = countVertexCells(table, table_size, vertex_ids, vertex_count);
1592

1593
#if TRACE
1594
		printf("pass %d (%s): grid size %d, vertices %d, %s\n",
1595
		    pass, (pass == 0) ? "guess" : (pass <= kInterpolationPasses) ? "lerp" : "binary",
1596
		    grid_size, int(vertices),
1597
		    (vertices <= target_vertex_count) ? "under" : "over");
1598
#endif
1599

1600 3
		float tip = interpolate(float(target_vertex_count), float(min_grid), float(min_vertices), float(grid_size), float(vertices), float(max_grid), float(max_vertices));
1601

1602 3
		if (vertices <= target_vertex_count)
1603
		{
1604 3
			min_grid = grid_size;
1605 3
			min_vertices = vertices;
1606
		}
1607
		else
1608
		{
1609 3
			max_grid = grid_size;
1610 3
			max_vertices = vertices;
1611
		}
1612

1613 3
		if (vertices == target_vertex_count || max_grid - min_grid <= 1)
1614 1
			break;
1615

1616
		// we start by using interpolation search - it usually converges faster
1617
		// 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)
1618 1
		next_grid_size = (pass < kInterpolationPasses) ? int(tip + 0.5f) : (min_grid + max_grid) / 2;
1619
	}
1620

1621 3
	if (min_vertices == 0)
1622 0
		return 0;
1623

1624
	// build vertex->cell association by mapping all vertices with the same quantized position to the same cell
1625 3
	unsigned int* vertex_cells = allocator.allocate<unsigned int>(vertex_count);
1626

1627 3
	computeVertexIds(vertex_ids, vertex_positions, vertex_count, min_grid);
1628 3
	size_t cell_count = fillVertexCells(table, table_size, vertex_cells, vertex_ids, vertex_count);
1629

1630
	// build a quadric for each target cell
1631 3
	Quadric* cell_quadrics = allocator.allocate<Quadric>(cell_count);
1632 3
	memset(cell_quadrics, 0, cell_count * sizeof(Quadric));
1633

1634 3
	fillCellQuadrics(cell_quadrics, vertex_positions, vertex_count, vertex_cells);
1635

1636
	// for each target cell, find the vertex with the minimal error
1637 3
	unsigned int* cell_remap = allocator.allocate<unsigned int>(cell_count);
1638 3
	float* cell_errors = allocator.allocate<float>(cell_count);
1639

1640 3
	fillCellRemap(cell_remap, cell_errors, cell_count, vertex_cells, cell_quadrics, vertex_positions, vertex_count);
1641

1642
	// copy results to the output
1643 1
	assert(cell_count <= target_vertex_count);
1644 3
	memcpy(destination, cell_remap, sizeof(unsigned int) * cell_count);
1645

1646
#if TRACE
1647
	printf("result: %d cells\n", int(cell_count));
1648
#endif
1649

1650 3
	return cell_count;
1651
}
1652

1653 3
float meshopt_simplifyScale(const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride)
1654
{
1655
	using namespace meshopt;
1656

1657 1
	assert(vertex_positions_stride > 0 && vertex_positions_stride <= 256);
1658 1
	assert(vertex_positions_stride % sizeof(float) == 0);
1659

1660 3
	float extent = rescalePositions(NULL, vertex_positions, vertex_count, vertex_positions_stride);
1661

1662 3
	return extent;
1663
}

Read our documentation on viewing source code .

Loading