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