quiet double-promotion warnings, change octree.cpp to use a float (vector accumulated...
[blender.git] / intern / dualcon / intern / octree.cpp
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * Contributor(s): Tao Ju, Nicholas Bishop
19  *
20  * ***** END GPL LICENSE BLOCK *****
21  */
22
23 #ifdef WITH_CXX_GUARDEDALLOC
24 #  include "MEM_guardedalloc.h"
25 #endif
26
27 #include "octree.h"
28 #include <Eigen/Dense>
29 #include <limits>
30 #include <time.h>
31
32 /**
33  * Implementations of Octree member functions.
34  *
35  * @author Tao Ju
36  */
37
38 /* set to non-zero value to enable debugging output */
39 #define DC_DEBUG 0
40
41 #if DC_DEBUG
42 /* enable debug printfs */
43 #define dc_printf printf
44 #else
45 /* disable debug printfs */
46 #define dc_printf(...) do {} while (0)
47 #endif
48
49 Octree::Octree(ModelReader *mr,
50                DualConAllocOutput alloc_output_func,
51                DualConAddVert add_vert_func,
52                DualConAddQuad add_quad_func,
53                DualConFlags flags, DualConMode dualcon_mode, int depth,
54                float threshold, float sharpness)
55         : use_flood_fill(flags & DUALCON_FLOOD_FILL),
56         /* note on `use_manifold':
57
58            After playing around with this option, the only case I could
59            find where this option gives different results is on
60            relatively thin corners. Sometimes along these corners two
61            vertices from separate sides will be placed in the same
62            position, so hole gets filled with a 5-sided face, where two
63            of those vertices are in the same 3D location. If
64            `use_manifold' is disabled, then the modifier doesn't
65            generate two separate vertices so the results end up as all
66            quads.
67
68            Since the results are just as good with all quads, there
69            doesn't seem any reason to allow this to be toggled in
70            Blender. -nicholasbishop
71          */
72         use_manifold(false),
73         hermite_num(sharpness),
74         mode(dualcon_mode),
75         alloc_output(alloc_output_func),
76         add_vert(add_vert_func),
77         add_quad(add_quad_func)
78 {
79         thresh = threshold;
80         reader = mr;
81         dimen = 1 << GRID_DIMENSION;
82         range = reader->getBoundingBox(origin);
83         nodeCount = nodeSpace = 0;
84         maxDepth = depth;
85         mindimen = (dimen >> maxDepth);
86         minshift = (GRID_DIMENSION - maxDepth);
87         buildTable();
88
89         maxTrianglePerCell = 0;
90
91         // Initialize memory
92 #ifdef IN_VERBOSE_MODE
93         dc_printf("Range: %f origin: %f, %f,%f \n", range, origin[0], origin[1], origin[2]);
94         dc_printf("Initialize memory...\n");
95 #endif
96         initMemory();
97         root = (Node *)createInternal(0);
98
99         // Read MC table
100 #ifdef IN_VERBOSE_MODE
101         dc_printf("Reading contour table...\n");
102 #endif
103         cubes = new Cubes();
104
105 }
106
107 Octree::~Octree()
108 {
109         delete cubes;
110         freeMemory();
111 }
112
113 void Octree::scanConvert()
114 {
115         // Scan triangles
116 #if DC_DEBUG
117         clock_t start, finish;
118         start = clock();
119 #endif
120
121         addAllTriangles();
122         resetMinimalEdges();
123         preparePrimalEdgesMask(&root->internal);
124
125 #if DC_DEBUG
126         finish = clock();
127         dc_printf("Time taken: %f seconds                   \n",
128                   (double)(finish - start) / CLOCKS_PER_SEC);
129 #endif
130
131         // Generate signs
132         // Find holes
133 #if DC_DEBUG
134         dc_printf("Patching...\n");
135         start = clock();
136 #endif
137         trace();
138 #if DC_DEBUG
139         finish = clock();
140         dc_printf("Time taken: %f seconds \n",  (double)(finish - start) / CLOCKS_PER_SEC);
141 #ifdef IN_VERBOSE_MODE
142         dc_printf("Holes: %d Average Length: %f Max Length: %d \n", numRings, (float)totRingLengths / (float) numRings, maxRingLength);
143 #endif
144 #endif
145
146         // Check again
147         int tnumRings = numRings;
148         trace();
149 #ifdef IN_VERBOSE_MODE
150         dc_printf("Holes after patching: %d \n", numRings);
151 #endif
152         numRings = tnumRings;
153
154 #if DC_DEBUG
155         dc_printf("Building signs...\n");
156         start = clock();
157 #endif
158         buildSigns();
159 #if DC_DEBUG
160         finish = clock();
161         dc_printf("Time taken: %f seconds \n",  (double)(finish - start) / CLOCKS_PER_SEC);
162 #endif
163
164         if (use_flood_fill) {
165                 /*
166                    start = clock();
167                    floodFill();
168                    // Check again
169                    tnumRings = numRings;
170                    trace();
171                    dc_printf("Holes after filling: %d \n", numRings);
172                    numRings = tnumRings;
173                    buildSigns();
174                    finish = clock();
175                    dc_printf("Time taken: %f seconds \n",       (double)(finish - start) / CLOCKS_PER_SEC);
176                  */
177 #if DC_DEBUG
178                 start = clock();
179                 dc_printf("Removing components...\n");
180 #endif
181                 floodFill();
182                 buildSigns();
183                 //      dc_printf("Checking...\n");
184                 //      floodFill();
185 #if DC_DEBUG
186                 finish = clock();
187                 dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC);
188 #endif
189         }
190
191         // Output
192 #if DC_DEBUG
193         start = clock();
194 #endif
195         writeOut();
196 #if DC_DEBUG
197         finish = clock();
198 #endif
199         // dc_printf("Time taken: %f seconds \n",       (double)(finish - start) / CLOCKS_PER_SEC);
200
201         // Print info
202 #ifdef IN_VERBOSE_MODE
203         printMemUsage();
204 #endif
205 }
206
207 void Octree::initMemory()
208 {
209         leafalloc[0] = new MemoryAllocator<sizeof(LeafNode)>();
210         leafalloc[1] = new MemoryAllocator<sizeof(LeafNode) + sizeof(float) *EDGE_FLOATS>();
211         leafalloc[2] = new MemoryAllocator<sizeof(LeafNode) + sizeof(float) *EDGE_FLOATS * 2>();
212         leafalloc[3] = new MemoryAllocator<sizeof(LeafNode) + sizeof(float) *EDGE_FLOATS * 3>();
213
214         alloc[0] = new MemoryAllocator<sizeof(InternalNode)>();
215         alloc[1] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *)>();
216         alloc[2] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *) * 2>();
217         alloc[3] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *) * 3>();
218         alloc[4] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *) * 4>();
219         alloc[5] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *) * 5>();
220         alloc[6] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *) * 6>();
221         alloc[7] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *) * 7>();
222         alloc[8] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *) * 8>();
223 }
224
225 void Octree::freeMemory()
226 {
227         for (int i = 0; i < 9; i++) {
228                 alloc[i]->destroy();
229                 delete alloc[i];
230         }
231
232         for (int i = 0; i < 4; i++) {
233                 leafalloc[i]->destroy();
234                 delete leafalloc[i];
235         }
236 }
237
238 void Octree::printMemUsage()
239 {
240         int totalbytes = 0;
241         dc_printf("********* Internal nodes: \n");
242         for (int i = 0; i < 9; i++) {
243                 alloc[i]->printInfo();
244
245                 totalbytes += alloc[i]->getAll() * alloc[i]->getBytes();
246         }
247         dc_printf("********* Leaf nodes: \n");
248         int totalLeafs = 0;
249         for (int i = 0; i < 4; i++) {
250                 leafalloc[i]->printInfo();
251
252                 totalbytes += leafalloc[i]->getAll() * leafalloc[i]->getBytes();
253                 totalLeafs += leafalloc[i]->getAllocated();
254         }
255
256         dc_printf("Total allocated bytes on disk: %d \n", totalbytes);
257         dc_printf("Total leaf nodes: %d\n", totalLeafs);
258 }
259
260 void Octree::resetMinimalEdges()
261 {
262         cellProcParity(root, 0, maxDepth);
263 }
264
265 void Octree::addAllTriangles()
266 {
267         Triangle *trian;
268         int count = 0;
269
270 #if DC_DEBUG
271         int total = reader->getNumTriangles();
272         int unitcount = 1000;
273         dc_printf("\nScan converting to depth %d...\n", maxDepth);
274 #endif
275
276         srand(0);
277
278         while ((trian = reader->getNextTriangle()) != NULL) {
279                 // Drop triangles
280                 {
281                         addTriangle(trian, count);
282                 }
283                 delete trian;
284
285                 count++;
286
287 #if DC_DEBUG
288                 if (count % unitcount == 0) {
289                         putchar(13);
290
291                         switch ((count / unitcount) % 4) {
292                                 case 0: dc_printf("-");
293                                         break;
294                                 case 1: dc_printf("/");
295                                         break;
296                                 case 2: dc_printf("|");
297                                         break;
298                                 case 3: dc_printf("\\");
299                                         break;
300                         }
301
302                         float percent = (float) count / total;
303
304                         /*
305                            int totbars = 50;
306                            int bars =(int)(percent * totbars);
307                            for(int i = 0; i < bars; i ++) {
308                             putchar(219);
309                            }
310                            for(i = bars; i < totbars; i ++) {
311                             putchar(176);
312                            }
313                          */
314
315                         dc_printf(" %d triangles: ", count);
316                         dc_printf(" %f%% complete.", 100 * percent);
317                 }
318 #endif
319
320         }
321         putchar(13);
322 }
323
324 /* Prepare a triangle for insertion into the octree; call the other
325    addTriangle() to (recursively) build the octree */
326 void Octree::addTriangle(Triangle *trian, int triind)
327 {
328         int i, j;
329
330         /* Project the triangle's coordinates into the grid */
331         for (i = 0; i < 3; i++) {
332                 for (j = 0; j < 3; j++)
333                         trian->vt[i][j] = dimen * (trian->vt[i][j] - origin[j]) / range;
334         }
335
336         /* Generate projections */
337         int64_t cube[2][3] = {{0, 0, 0}, {dimen, dimen, dimen}};
338         int64_t trig[3][3];
339         for (i = 0; i < 3; i++) {
340                 for (j = 0; j < 3; j++)
341                         trig[i][j] = (int64_t)(trian->vt[i][j]);
342         }
343
344         /* Add triangle to the octree */
345         int64_t errorvec = (int64_t)(0);
346         CubeTriangleIsect *proj = new CubeTriangleIsect(cube, trig, errorvec, triind);
347         root = (Node *)addTriangle(&root->internal, proj, maxDepth);
348
349         delete proj->inherit;
350         delete proj;
351 }
352
353 #if 0
354 static void print_depth(int height, int maxDepth)
355 {
356         for (int i = 0; i < maxDepth - height; i++)
357                 printf("  ");
358 }
359 #endif
360
361 InternalNode *Octree::addTriangle(InternalNode *node, CubeTriangleIsect *p, int height)
362 {
363         int i;
364         const int vertdiff[8][3] = {
365                 {0,  0,  0},
366                 {0,  0,  1},
367                 {0,  1, -1},
368                 {0,  0,  1},
369                 {1, -1, -1},
370                 {0,  0,  1},
371                 {0,  1, -1},
372                 {0,  0,  1}};
373         unsigned char boxmask = p->getBoxMask();
374         CubeTriangleIsect *subp = new CubeTriangleIsect(p);
375         
376         int count = 0;
377         int tempdiff[3] = {0, 0, 0};
378
379         /* Check triangle against each of the input node's children */
380         for (i = 0; i < 8; i++) {
381                 tempdiff[0] += vertdiff[i][0];
382                 tempdiff[1] += vertdiff[i][1];
383                 tempdiff[2] += vertdiff[i][2];
384
385                 /* Quick pruning using bounding box */
386                 if (boxmask & (1 << i)) {
387                         subp->shift(tempdiff);
388                         tempdiff[0] = tempdiff[1] = tempdiff[2] = 0;
389
390                         /* Pruning using intersection test */
391                         if (subp->isIntersecting()) {
392                                 if (!node->has_child(i)) {
393                                         if (height == 1)
394                                                 node = addLeafChild(node, i, count, createLeaf(0));
395                                         else
396                                                 node = addInternalChild(node, i, count, createInternal(0));
397                                 }
398                                 Node *chd = node->get_child(count);
399
400                                 if (node->is_child_leaf(i))
401                                         node->set_child(count, (Node *)updateCell(&chd->leaf, subp));
402                                 else
403                                         node->set_child(count, (Node *)addTriangle(&chd->internal, subp, height - 1));
404                         }
405                 }
406
407                 if (node->has_child(i))
408                         count++;
409         }
410
411         delete subp;
412
413         return node;
414 }
415
416 LeafNode *Octree::updateCell(LeafNode *node, CubeTriangleIsect *p)
417 {
418         int i;
419
420         // Edge connectivity
421         int mask[3] = {0, 4, 8  };
422         int oldc = 0, newc = 0;
423         float offs[3];
424         float a[3], b[3], c[3];
425
426         for (i = 0; i < 3; i++) {
427                 if (!getEdgeParity(node, mask[i])) {
428                         if (p->isIntersectingPrimary(i)) {
429                                 // actualQuads ++;
430                                 setEdge(node, mask[i]);
431                                 offs[newc] = p->getIntersectionPrimary(i);
432                                 a[newc] = (float) p->inherit->norm[0];
433                                 b[newc] = (float) p->inherit->norm[1];
434                                 c[newc] = (float) p->inherit->norm[2];
435                                 newc++;
436                         }
437                 }
438                 else {
439                         offs[newc] = getEdgeOffsetNormal(node, oldc, a[newc], b[newc], c[newc]);
440
441                         oldc++;
442                         newc++;
443                 }
444         }
445
446         if (newc > oldc) {
447                 // New offsets added, update this node
448                 node = updateEdgeOffsetsNormals(node, oldc, newc, offs, a, b, c);
449         }
450
451         return node;
452 }
453
454 void Octree::preparePrimalEdgesMask(InternalNode *node)
455 {
456         int count = 0;
457         for (int i = 0; i < 8; i++) {
458                 if (node->has_child(i)) {
459                         if (node->is_child_leaf(i))
460                                 createPrimalEdgesMask(&node->get_child(count)->leaf);
461                         else
462                                 preparePrimalEdgesMask(&node->get_child(count)->internal);
463
464                         count++;
465                 }
466         }
467 }
468
469 void Octree::trace()
470 {
471         int st[3] = {0, 0, 0, };
472         numRings = 0;
473         totRingLengths = 0;
474         maxRingLength = 0;
475
476         PathList *chdpath = NULL;
477         root = trace(root, st, dimen, maxDepth, chdpath);
478
479         if (chdpath != NULL) {
480                 dc_printf("there are incomplete rings.\n");
481                 printPaths(chdpath);
482         };
483 }
484
485 Node *Octree::trace(Node *newnode, int *st, int len, int depth, PathList *& paths)
486 {
487         len >>= 1;
488         PathList *chdpaths[8];
489         Node *chd[8];
490         int nst[8][3];
491         int i, j;
492
493         // Get children paths
494         int chdleaf[8];
495         newnode->internal.fill_children(chd, chdleaf);
496
497         // int count = 0;
498         for (i = 0; i < 8; i++) {
499                 for (j = 0; j < 3; j++) {
500                         nst[i][j] = st[j] + len * vertmap[i][j];
501                 }
502
503                 if (chd[i] == NULL || newnode->internal.is_child_leaf(i)) {
504                         chdpaths[i] = NULL;
505                 }
506                 else {
507                         trace(chd[i], nst[i], len, depth - 1, chdpaths[i]);
508                 }
509         }
510
511         // Get connectors on the faces
512         PathList *conn[12];
513         Node *nf[2];
514         int lf[2];
515         int df[2] = {depth - 1, depth - 1};
516         int *nstf[2];
517
518         newnode->internal.fill_children(chd, chdleaf);
519
520         for (i = 0; i < 12; i++) {
521                 int c[2] = {cellProcFaceMask[i][0], cellProcFaceMask[i][1]};
522
523                 for (int j = 0; j < 2; j++) {
524                         lf[j] = chdleaf[c[j]];
525                         nf[j] = chd[c[j]];
526                         nstf[j] = nst[c[j]];
527                 }
528
529                 conn[i] = NULL;
530
531                 findPaths((Node **)nf, lf, df, nstf, depth - 1, cellProcFaceMask[i][2], conn[i]);
532
533                 //if(conn[i]) {
534                 //              printPath(conn[i]);
535                 //}
536         }
537
538         // Connect paths
539         PathList *rings = NULL;
540         combinePaths(chdpaths[0], chdpaths[1], conn[8], rings);
541         combinePaths(chdpaths[2], chdpaths[3], conn[9], rings);
542         combinePaths(chdpaths[4], chdpaths[5], conn[10], rings);
543         combinePaths(chdpaths[6], chdpaths[7], conn[11], rings);
544
545         combinePaths(chdpaths[0], chdpaths[2], conn[4], rings);
546         combinePaths(chdpaths[4], chdpaths[6], conn[5], rings);
547         combinePaths(chdpaths[0], NULL, conn[6], rings);
548         combinePaths(chdpaths[4], NULL, conn[7], rings);
549
550         combinePaths(chdpaths[0], chdpaths[4], conn[0], rings);
551         combinePaths(chdpaths[0], NULL, conn[1], rings);
552         combinePaths(chdpaths[0], NULL, conn[2], rings);
553         combinePaths(chdpaths[0], NULL, conn[3], rings);
554
555         // By now, only chdpaths[0] and rings have contents
556
557         // Process rings
558         if (rings) {
559                 // printPath(rings);
560
561                 /* Let's count first */
562                 PathList *trings = rings;
563                 while (trings) {
564                         numRings++;
565                         totRingLengths += trings->length;
566                         if (trings->length > maxRingLength) {
567                                 maxRingLength = trings->length;
568                         }
569                         trings = trings->next;
570                 }
571
572                 // printPath(rings);
573                 newnode = patch(newnode, st, (len << 1), rings);
574         }
575
576         // Return incomplete paths
577         paths = chdpaths[0];
578         return newnode;
579 }
580
581 void Octree::findPaths(Node *node[2], int leaf[2], int depth[2], int *st[2], int maxdep, int dir, PathList *& paths)
582 {
583         if (!(node[0] && node[1])) {
584                 return;
585         }
586
587         if (!(leaf[0] && leaf[1])) {
588                 // Not at the bottom, recur
589
590                 // Fill children nodes
591                 int i, j;
592                 Node *chd[2][8];
593                 int chdleaf[2][8];
594                 int nst[2][8][3];
595
596                 for (j = 0; j < 2; j++) {
597                         if (!leaf[j]) {
598                                 node[j]->internal.fill_children(chd[j], chdleaf[j]);
599
600                                 int len = (dimen >> (maxDepth - depth[j] + 1));
601                                 for (i = 0; i < 8; i++) {
602                                         for (int k = 0; k < 3; k++) {
603                                                 nst[j][i][k] = st[j][k] + len * vertmap[i][k];
604                                         }
605                                 }
606
607                         }
608                 }
609
610                 // 4 face calls
611                 Node *nf[2];
612                 int df[2];
613                 int lf[2];
614                 int *nstf[2];
615                 for (i = 0; i < 4; i++) {
616                         int c[2] = {faceProcFaceMask[dir][i][0], faceProcFaceMask[dir][i][1]};
617                         for (int j = 0; j < 2; j++) {
618                                 if (leaf[j]) {
619                                         lf[j] = leaf[j];
620                                         nf[j] = node[j];
621                                         df[j] = depth[j];
622                                         nstf[j] = st[j];
623                                 }
624                                 else {
625                                         lf[j] = chdleaf[j][c[j]];
626                                         nf[j] = chd[j][c[j]];
627                                         df[j] = depth[j] - 1;
628                                         nstf[j] = nst[j][c[j]];
629                                 }
630                         }
631                         findPaths(nf, lf, df, nstf, maxdep - 1, faceProcFaceMask[dir][i][2], paths);
632                 }
633
634         }
635         else {
636                 // At the bottom, check this face
637                 int ind = (depth[0] == maxdep ? 0 : 1);
638                 int fcind = 2 * dir + (1 - ind);
639                 if (getFaceParity((LeafNode *)node[ind], fcind)) {
640                         // Add into path
641                         PathElement *ele1 = new PathElement;
642                         PathElement *ele2 = new PathElement;
643
644                         ele1->pos[0] = st[0][0];
645                         ele1->pos[1] = st[0][1];
646                         ele1->pos[2] = st[0][2];
647
648                         ele2->pos[0] = st[1][0];
649                         ele2->pos[1] = st[1][1];
650                         ele2->pos[2] = st[1][2];
651
652                         ele1->next = ele2;
653                         ele2->next = NULL;
654
655                         PathList *lst = new PathList;
656                         lst->head = ele1;
657                         lst->tail = ele2;
658                         lst->length = 2;
659                         lst->next = paths;
660                         paths = lst;
661
662                         // int l =(dimen >> maxDepth);
663                 }
664         }
665
666 }
667
668 void Octree::combinePaths(PathList *& list1, PathList *list2, PathList *paths, PathList *& rings)
669 {
670         // Make new list of paths
671         PathList *nlist = NULL;
672
673         // Search for each connectors in paths
674         PathList *tpaths = paths;
675         PathList *tlist, *pre;
676         while (tpaths) {
677                 PathList *singlist = tpaths;
678                 PathList *templist;
679                 tpaths = tpaths->next;
680                 singlist->next = NULL;
681
682                 // Look for hookup in list1
683                 tlist = list1;
684                 pre = NULL;
685                 while (tlist) {
686                         if ((templist = combineSinglePath(list1, pre, tlist, singlist, NULL, singlist)) != NULL) {
687                                 singlist = templist;
688                                 continue;
689                         }
690                         pre = tlist;
691                         tlist = tlist->next;
692                 }
693
694                 // Look for hookup in list2
695                 tlist = list2;
696                 pre = NULL;
697                 while (tlist) {
698                         if ((templist = combineSinglePath(list2, pre, tlist, singlist, NULL, singlist)) != NULL) {
699                                 singlist = templist;
700                                 continue;
701                         }
702                         pre = tlist;
703                         tlist = tlist->next;
704                 }
705
706                 // Look for hookup in nlist
707                 tlist = nlist;
708                 pre = NULL;
709                 while (tlist) {
710                         if ((templist = combineSinglePath(nlist, pre, tlist, singlist, NULL, singlist)) != NULL) {
711                                 singlist = templist;
712                                 continue;
713                         }
714                         pre = tlist;
715                         tlist = tlist->next;
716                 }
717
718                 // Add to nlist or rings
719                 if (isEqual(singlist->head, singlist->tail)) {
720                         PathElement *temp = singlist->head;
721                         singlist->head = temp->next;
722                         delete temp;
723                         singlist->length--;
724                         singlist->tail->next = singlist->head;
725
726                         singlist->next = rings;
727                         rings = singlist;
728                 }
729                 else {
730                         singlist->next = nlist;
731                         nlist = singlist;
732                 }
733
734         }
735
736         // Append list2 and nlist to the end of list1
737         tlist = list1;
738         if (tlist != NULL) {
739                 while (tlist->next != NULL) {
740                         tlist = tlist->next;
741                 }
742                 tlist->next = list2;
743         }
744         else {
745                 tlist = list2;
746                 list1 = list2;
747         }
748
749         if (tlist != NULL) {
750                 while (tlist->next != NULL) {
751                         tlist = tlist->next;
752                 }
753                 tlist->next = nlist;
754         }
755         else {
756                 tlist = nlist;
757                 list1 = nlist;
758         }
759
760 }
761
762 PathList *Octree::combineSinglePath(PathList *& head1, PathList *pre1, PathList *& list1, PathList *& head2, PathList *pre2, PathList *& list2)
763 {
764         if (isEqual(list1->head, list2->head) || isEqual(list1->tail, list2->tail)) {
765                 // Reverse the list
766                 if (list1->length < list2->length) {
767                         // Reverse list1
768                         PathElement *prev = list1->head;
769                         PathElement *next = prev->next;
770                         prev->next = NULL;
771                         while (next != NULL) {
772                                 PathElement *tnext = next->next;
773                                 next->next = prev;
774
775                                 prev = next;
776                                 next = tnext;
777                         }
778
779                         list1->tail = list1->head;
780                         list1->head = prev;
781                 }
782                 else {
783                         // Reverse list2
784                         PathElement *prev = list2->head;
785                         PathElement *next = prev->next;
786                         prev->next = NULL;
787                         while (next != NULL) {
788                                 PathElement *tnext = next->next;
789                                 next->next = prev;
790
791                                 prev = next;
792                                 next = tnext;
793                         }
794
795                         list2->tail = list2->head;
796                         list2->head = prev;
797                 }
798         }
799
800         if (isEqual(list1->head, list2->tail)) {
801
802                 // Easy case
803                 PathElement *temp = list1->head->next;
804                 delete list1->head;
805                 list2->tail->next = temp;
806
807                 PathList *nlist = new PathList;
808                 nlist->length = list1->length + list2->length - 1;
809                 nlist->head = list2->head;
810                 nlist->tail = list1->tail;
811                 nlist->next = NULL;
812
813                 deletePath(head1, pre1, list1);
814                 deletePath(head2, pre2, list2);
815
816                 return nlist;
817         }
818         else if (isEqual(list1->tail, list2->head)) {
819                 // Easy case
820                 PathElement *temp = list2->head->next;
821                 delete list2->head;
822                 list1->tail->next = temp;
823
824                 PathList *nlist = new PathList;
825                 nlist->length = list1->length + list2->length - 1;
826                 nlist->head = list1->head;
827                 nlist->tail = list2->tail;
828                 nlist->next = NULL;
829
830                 deletePath(head1, pre1, list1);
831                 deletePath(head2, pre2, list2);
832
833                 return nlist;
834         }
835
836         return NULL;
837 }
838
839 void Octree::deletePath(PathList *& head, PathList *pre, PathList *& curr)
840 {
841         PathList *temp = curr;
842         curr = temp->next;
843         delete temp;
844
845         if (pre == NULL) {
846                 head = curr;
847         }
848         else {
849                 pre->next = curr;
850         }
851 }
852
853 void Octree::printElement(PathElement *ele)
854 {
855         if (ele != NULL) {
856                 dc_printf("(%d %d %d)", ele->pos[0], ele->pos[1], ele->pos[2]);
857         }
858 }
859
860 void Octree::printPath(PathList *path)
861 {
862         PathElement *n = path->head;
863         int same = 0;
864
865 #if DC_DEBUG
866         int len = (dimen >> maxDepth);
867 #endif
868         while (n && (same == 0 || n != path->head)) {
869                 same++;
870                 dc_printf("(%d %d %d)", n->pos[0] / len, n->pos[1] / len, n->pos[2] / len);
871                 n = n->next;
872         }
873
874         if (n == path->head) {
875                 dc_printf(" Ring!\n");
876         }
877         else {
878                 dc_printf(" %p end!\n", n);
879         }
880 }
881
882 void Octree::printPath(PathElement *path)
883 {
884         PathElement *n = path;
885         int same = 0;
886 #if DC_DEBUG
887         int len = (dimen >> maxDepth);
888 #endif
889         while (n && (same == 0 || n != path)) {
890                 same++;
891                 dc_printf("(%d %d %d)", n->pos[0] / len, n->pos[1] / len, n->pos[2] / len);
892                 n = n->next;
893         }
894
895         if (n == path) {
896                 dc_printf(" Ring!\n");
897         }
898         else {
899                 dc_printf(" %p end!\n", n);
900         }
901
902 }
903
904
905 void Octree::printPaths(PathList *path)
906 {
907         PathList *iter = path;
908         int i = 0;
909         while (iter != NULL) {
910                 dc_printf("Path %d:\n", i);
911                 printPath(iter);
912                 iter = iter->next;
913                 i++;
914         }
915 }
916
917 Node *Octree::patch(Node *newnode, int st[3], int len, PathList *rings)
918 {
919 #ifdef IN_DEBUG_MODE
920         dc_printf("Call to PATCH with rings: \n");
921         printPaths(rings);
922 #endif
923
924         /* Do nothing but couting
925            PathList* tlist = rings;
926            PathList* ttlist;
927            PathElement* telem, * ttelem;
928            while(tlist!= NULL) {
929             // printPath(tlist);
930             numRings ++;
931             totRingLengths += tlist->length;
932             if(tlist->length > maxRingLength) {
933                 maxRingLength = tlist->length;
934             }
935             ttlist = tlist;
936             tlist = tlist->next;
937            }
938            return node;
939          */
940
941
942         /* Pass onto separate calls in each direction */
943         if (len == mindimen) {
944                 dc_printf("Error! should have no list by now.\n");
945                 exit(0);
946         }
947
948         // YZ plane
949         PathList *xlists[2];
950         newnode = patchSplit(newnode, st, len, rings, 0, xlists[0], xlists[1]);
951
952         // XZ plane
953         PathList *ylists[4];
954         newnode = patchSplit(newnode, st, len, xlists[0], 1, ylists[0], ylists[1]);
955         newnode = patchSplit(newnode, st, len, xlists[1], 1, ylists[2], ylists[3]);
956
957         // XY plane
958         PathList *zlists[8];
959         newnode = patchSplit(newnode, st, len, ylists[0], 2, zlists[0], zlists[1]);
960         newnode = patchSplit(newnode, st, len, ylists[1], 2, zlists[2], zlists[3]);
961         newnode = patchSplit(newnode, st, len, ylists[2], 2, zlists[4], zlists[5]);
962         newnode = patchSplit(newnode, st, len, ylists[3], 2, zlists[6], zlists[7]);
963
964         // Recur
965         len >>= 1;
966         int count = 0;
967         for (int i = 0; i < 8; i++) {
968                 if (zlists[i] != NULL) {
969                         int nori[3] = {
970                                 st[0] + len * vertmap[i][0],
971                                 st[1] + len * vertmap[i][1],
972                                 st[2] + len * vertmap[i][2]
973                         };
974                         patch(newnode->internal.get_child(count), nori, len, zlists[i]);
975                 }
976
977                 if (newnode->internal.has_child(i)) {
978                         count++;
979                 }
980         }
981 #ifdef IN_DEBUG_MODE
982         dc_printf("Return from PATCH\n");
983 #endif
984         return newnode;
985
986 }
987
988
989 Node *Octree::patchSplit(Node *newnode, int st[3], int len, PathList *rings,
990                          int dir, PathList *& nrings1, PathList *& nrings2)
991 {
992 #ifdef IN_DEBUG_MODE
993         dc_printf("Call to PATCHSPLIT with direction %d and rings: \n", dir);
994         printPaths(rings);
995 #endif
996
997         nrings1 = NULL;
998         nrings2 = NULL;
999         PathList *tmp;
1000         while (rings != NULL) {
1001                 // Process this ring
1002                 newnode = patchSplitSingle(newnode, st, len, rings->head, dir, nrings1, nrings2);
1003
1004                 // Delete this ring from the group
1005                 tmp = rings;
1006                 rings = rings->next;
1007                 delete tmp;
1008         }
1009
1010 #ifdef IN_DEBUG_MODE
1011         dc_printf("Return from PATCHSPLIT with \n");
1012         dc_printf("Rings gourp 1:\n");
1013         printPaths(nrings1);
1014         dc_printf("Rings group 2:\n");
1015         printPaths(nrings2);
1016 #endif
1017
1018         return newnode;
1019 }
1020
1021 Node *Octree::patchSplitSingle(Node *newnode, int st[3], int len, PathElement *head, int dir, PathList *& nrings1, PathList *& nrings2)
1022 {
1023 #ifdef IN_DEBUG_MODE
1024         dc_printf("Call to PATCHSPLITSINGLE with direction %d and path: \n", dir);
1025         printPath(head);
1026 #endif
1027
1028         if (head == NULL) {
1029 #ifdef IN_DEBUG_MODE
1030                 dc_printf("Return from PATCHSPLITSINGLE with head==NULL.\n");
1031 #endif
1032                 return newnode;
1033         }
1034         else {
1035                 // printPath(head);
1036         }
1037
1038         // Walk along the ring to find pair of intersections
1039         PathElement *pre1 = NULL;
1040         PathElement *pre2 = NULL;
1041         int side = findPair(head, st[dir] + len / 2, dir, pre1, pre2);
1042
1043         /*
1044            if(pre1 == pre2) {
1045             int edgelen =(dimen >> maxDepth);
1046             dc_printf("Location: %d %d %d Direction: %d Reso: %d\n", st[0]/edgelen, st[1]/edgelen, st[2]/edgelen, dir, len/edgelen);
1047             printPath(head);
1048             exit(0);
1049            }
1050          */
1051
1052         if (side) {
1053                 // Entirely on one side
1054                 PathList *nring = new PathList();
1055                 nring->head = head;
1056
1057                 if (side == -1) {
1058                         nring->next = nrings1;
1059                         nrings1 = nring;
1060                 }
1061                 else {
1062                         nring->next = nrings2;
1063                         nrings2 = nring;
1064                 }
1065         }
1066         else {
1067                 // Break into two parts
1068                 PathElement *nxt1 = pre1->next;
1069                 PathElement *nxt2 = pre2->next;
1070                 pre1->next = nxt2;
1071                 pre2->next = nxt1;
1072
1073                 newnode = connectFace(newnode, st, len, dir, pre1, pre2);
1074
1075                 if (isEqual(pre1, pre1->next)) {
1076                         if (pre1 == pre1->next) {
1077                                 delete pre1;
1078                                 pre1 = NULL;
1079                         }
1080                         else {
1081                                 PathElement *temp = pre1->next;
1082                                 pre1->next = temp->next;
1083                                 delete temp;
1084                         }
1085                 }
1086                 if (isEqual(pre2, pre2->next)) {
1087                         if (pre2 == pre2->next) {
1088                                 delete pre2;
1089                                 pre2 = NULL;
1090                         }
1091                         else {
1092                                 PathElement *temp = pre2->next;
1093                                 pre2->next = temp->next;
1094                                 delete temp;
1095                         }
1096                 }
1097
1098                 compressRing(pre1);
1099                 compressRing(pre2);
1100
1101                 // Recur
1102                 newnode = patchSplitSingle(newnode, st, len, pre1, dir, nrings1, nrings2);
1103                 newnode = patchSplitSingle(newnode, st, len, pre2, dir, nrings1, nrings2);
1104
1105         }
1106
1107 #ifdef IN_DEBUG_MODE
1108         dc_printf("Return from PATCHSPLITSINGLE with \n");
1109         dc_printf("Rings gourp 1:\n");
1110         printPaths(nrings1);
1111         dc_printf("Rings group 2:\n");
1112         printPaths(nrings2);
1113 #endif
1114
1115         return newnode;
1116 }
1117
1118 Node *Octree::connectFace(Node *newnode, int st[3], int len, int dir,
1119                           PathElement *f1, PathElement *f2)
1120 {
1121 #ifdef IN_DEBUG_MODE
1122         dc_printf("Call to CONNECTFACE with direction %d and length %d path: \n", dir, len);
1123         dc_printf("Path(low side): \n");
1124         printPath(f1);
1125 //      checkPath(f1);
1126         dc_printf("Path(high side): \n");
1127         printPath(f2);
1128 //      checkPath(f2);
1129 #endif
1130
1131         // Setup 2D
1132         int pos = st[dir] + len / 2;
1133         int xdir = (dir + 1) % 3;
1134         int ydir = (dir + 2) % 3;
1135
1136         // Use existing intersections on f1 and f2
1137         int x1, y1, x2, y2;
1138         float p1, q1, p2, q2;
1139
1140         getFacePoint(f2->next, dir, x1, y1, p1, q1);
1141         getFacePoint(f2, dir, x2, y2, p2, q2);
1142
1143         float dx = x2 + p2 - x1 - p1;
1144         float dy = y2 + q2 - y1 - q1;
1145
1146         // Do adapted Bresenham line drawing
1147         float rx = p1, ry = q1;
1148         int incx = 1, incy = 1;
1149         int lx = x1, ly = y1;
1150         int hx = x2, hy = y2;
1151         int choice;
1152         if (x2 < x1) {
1153                 incx = -1;
1154                 rx = 1 - rx;
1155                 lx = x2;
1156                 hx = x1;
1157         }
1158         if (y2 < y1) {
1159                 incy = -1;
1160                 ry = 1 - ry;
1161                 ly = y2;
1162                 hy = y1;
1163         }
1164
1165         float sx = dx * incx;
1166         float sy = dy * incy;
1167
1168         int ori[3];
1169         ori[dir] = pos / mindimen;
1170         ori[xdir] = x1;
1171         ori[ydir] = y1;
1172         int walkdir;
1173         int inc;
1174         float alpha;
1175
1176         PathElement *curEleN = f1;
1177         PathElement *curEleP = f2->next;
1178         Node *nodeN = NULL, *nodeP = NULL;
1179         LeafNode *curN = locateLeaf(&newnode->internal, len, f1->pos);
1180         LeafNode *curP = locateLeaf(&newnode->internal, len, f2->next->pos);
1181         if (curN == NULL || curP == NULL) {
1182                 exit(0);
1183         }
1184         int stN[3], stP[3];
1185         int lenN, lenP;
1186
1187         /* Unused code, leaving for posterity
1188
1189            float stpt[3], edpt[3];
1190            stpt[dir] = edpt[dir] =(float) pos;
1191            stpt[xdir] =(x1 + p1) * mindimen;
1192            stpt[ydir] =(y1 + q1) * mindimen;
1193            edpt[xdir] =(x2 + p2) * mindimen;
1194            edpt[ydir] =(y2 + q2) * mindimen;
1195          */
1196         while (ori[xdir] != x2 || ori[ydir] != y2) {
1197                 int next;
1198                 if (sy * (1 - rx) > sx * (1 - ry)) {
1199                         choice = 1;
1200                         next = ori[ydir] + incy;
1201                         if (next < ly || next > hy) {
1202                                 choice = 4;
1203                                 next = ori[xdir] + incx;
1204                         }
1205                 }
1206                 else {
1207                         choice = 2;
1208                         next = ori[xdir] + incx;
1209                         if (next < lx || next > hx) {
1210                                 choice = 3;
1211                                 next = ori[ydir] + incy;
1212                         }
1213                 }
1214
1215                 if (choice & 1) {
1216                         ori[ydir] = next;
1217                         if (choice == 1) {
1218                                 rx += (sy == 0 ? 0 : (1 - ry) * sx / sy);
1219                                 ry = 0;
1220                         }
1221
1222                         walkdir = 2;
1223                         inc = incy;
1224                         alpha = x2 < x1 ? 1 - rx : rx;
1225                 }
1226                 else {
1227                         ori[xdir] = next;
1228                         if (choice == 2) {
1229                                 ry += (sx == 0 ? 0 : (1 - rx) * sy / sx);
1230                                 rx = 0;
1231                         }
1232
1233                         walkdir = 1;
1234                         inc = incx;
1235                         alpha = y2 < y1 ? 1 - ry : ry;
1236                 }
1237
1238
1239
1240                 // Get the exact location of the marcher
1241                 int nori[3] = {ori[0] * mindimen, ori[1] * mindimen, ori[2] * mindimen};
1242                 float spt[3] = {(float) nori[0], (float) nori[1], (float) nori[2]};
1243                 spt[(dir + (3 - walkdir)) % 3] += alpha * mindimen;
1244                 if (inc < 0) {
1245                         spt[(dir + walkdir) % 3] += mindimen;
1246                 }
1247
1248                 // dc_printf("new x,y: %d %d\n", ori[xdir] / edgelen, ori[ydir] / edgelen);
1249                 // dc_printf("nori: %d %d %d alpha: %f walkdir: %d\n", nori[0], nori[1], nori[2], alpha, walkdir);
1250                 // dc_printf("%f %f %f\n", spt[0], spt[1], spt[2]);
1251
1252                 // Locate the current cells on both sides
1253                 newnode = locateCell(&newnode->internal, st, len, nori, dir, 1, nodeN, stN, lenN);
1254                 newnode = locateCell(&newnode->internal, st, len, nori, dir, 0, nodeP, stP, lenP);
1255
1256                 updateParent(&newnode->internal, len, st);
1257
1258                 int flag = 0;
1259                 // Add the cells to the rings and fill in the patch
1260                 PathElement *newEleN;
1261                 if (curEleN->pos[0] != stN[0] || curEleN->pos[1] != stN[1] || curEleN->pos[2] != stN[2]) {
1262                         if (curEleN->next->pos[0] != stN[0] || curEleN->next->pos[1] != stN[1] || curEleN->next->pos[2] != stN[2]) {
1263                                 newEleN = new PathElement;
1264                                 newEleN->next = curEleN->next;
1265                                 newEleN->pos[0] = stN[0];
1266                                 newEleN->pos[1] = stN[1];
1267                                 newEleN->pos[2] = stN[2];
1268
1269                                 curEleN->next = newEleN;
1270                         }
1271                         else {
1272                                 newEleN = curEleN->next;
1273                         }
1274                         curN = patchAdjacent(&newnode->internal, len, curEleN->pos, curN,
1275                                              newEleN->pos, (LeafNode *)nodeN, walkdir,
1276                                              inc, dir, 1, alpha);
1277
1278                         curEleN = newEleN;
1279                         flag++;
1280                 }
1281
1282                 PathElement *newEleP;
1283                 if (curEleP->pos[0] != stP[0] || curEleP->pos[1] != stP[1] || curEleP->pos[2] != stP[2]) {
1284                         if (f2->pos[0] != stP[0] || f2->pos[1] != stP[1] || f2->pos[2] != stP[2]) {
1285                                 newEleP = new PathElement;
1286                                 newEleP->next = curEleP;
1287                                 newEleP->pos[0] = stP[0];
1288                                 newEleP->pos[1] = stP[1];
1289                                 newEleP->pos[2] = stP[2];
1290
1291                                 f2->next = newEleP;
1292                         }
1293                         else {
1294                                 newEleP = f2;
1295                         }
1296                         curP = patchAdjacent(&newnode->internal, len, curEleP->pos, curP,
1297                                              newEleP->pos, (LeafNode *)nodeP, walkdir,
1298                                              inc, dir, 0, alpha);
1299
1300
1301
1302                         curEleP = newEleP;
1303                         flag++;
1304                 }
1305
1306
1307                 /*
1308                    if(flag == 0) {
1309                     dc_printf("error: non-synchronized patching! at \n");
1310                    }
1311                  */
1312         }
1313
1314 #ifdef IN_DEBUG_MODE
1315         dc_printf("Return from CONNECTFACE with \n");
1316         dc_printf("Path(low side):\n");
1317         printPath(f1);
1318         checkPath(f1);
1319         dc_printf("Path(high side):\n");
1320         printPath(f2);
1321         checkPath(f2);
1322 #endif
1323
1324
1325         return newnode;
1326 }
1327
1328 LeafNode *Octree::patchAdjacent(InternalNode *node, int len, int st1[3],
1329                                 LeafNode *leaf1, int st2[3], LeafNode *leaf2,
1330                                 int walkdir, int inc, int dir, int side,
1331                                 float alpha)
1332 {
1333 #ifdef IN_DEBUG_MODE
1334         dc_printf("Before patching.\n");
1335         printInfo(st1);
1336         printInfo(st2);
1337         dc_printf("-----------------%d %d %d; %d %d %d\n", st1[0], st2[1], st1[2], st2[0], st2[1], st2[2]);
1338 #endif
1339
1340         // Get edge index on each leaf
1341         int edgedir = (dir + (3 - walkdir)) % 3;
1342         int incdir = (dir + walkdir) % 3;
1343         int ind1 = (edgedir == 1 ? (dir + 3 - edgedir) % 3 - 1 : 2 - (dir + 3 - edgedir) % 3);
1344         int ind2 = (edgedir == 1 ? (incdir + 3 - edgedir) % 3 - 1 : 2 - (incdir + 3 - edgedir) % 3);
1345
1346         int eind1 = ((edgedir << 2) | (side << ind1) | ((inc > 0 ? 1 : 0) << ind2));
1347         int eind2 = ((edgedir << 2) | (side << ind1) | ((inc > 0 ? 0 : 1) << ind2));
1348
1349 #ifdef IN_DEBUG_MODE
1350         dc_printf("Index 1: %d Alpha 1: %f Index 2: %d Alpha 2: %f\n", eind1, alpha, eind2, alpha);
1351         /*
1352            if(alpha < 0 || alpha > 1) {
1353             dc_printf("Index 1: %d Alpha 1: %f Index 2: %d Alpha 2: %f\n", eind1, alpha, eind2, alpha);
1354             printInfo(st1);
1355             printInfo(st2);
1356            }
1357          */
1358 #endif
1359
1360         // Flip edge parity
1361         LeafNode *nleaf1 = flipEdge(leaf1, eind1, alpha);
1362         LeafNode *nleaf2 = flipEdge(leaf2, eind2, alpha);
1363
1364         // Update parent link
1365         updateParent(node, len, st1, nleaf1);
1366         updateParent(node, len, st2, nleaf2);
1367         // updateParent(nleaf1, mindimen, st1);
1368         // updateParent(nleaf2, mindimen, st2);
1369
1370         /*
1371            float m[3];
1372            dc_printf("Adding new point: %f %f %f\n", spt[0], spt[1], spt[2]);
1373            getMinimizer(leaf1, m);
1374            dc_printf("Cell %d now has minimizer %f %f %f\n", leaf1, m[0], m[1], m[2]);
1375            getMinimizer(leaf2, m);
1376            dc_printf("Cell %d now has minimizer %f %f %f\n", leaf2, m[0], m[1], m[2]);
1377          */
1378
1379 #ifdef IN_DEBUG_MODE
1380         dc_printf("After patching.\n");
1381         printInfo(st1);
1382         printInfo(st2);
1383 #endif
1384         return nleaf2;
1385 }
1386
1387 Node *Octree::locateCell(InternalNode *node, int st[3], int len, int ori[3], int dir, int side, Node *& rleaf, int rst[3], int& rlen)
1388 {
1389 #ifdef IN_DEBUG_MODE
1390         // dc_printf("Call to LOCATECELL with node ");
1391         // printNode(node);
1392 #endif
1393
1394         int i;
1395         len >>= 1;
1396         int ind = 0;
1397         for (i = 0; i < 3; i++) {
1398                 ind <<= 1;
1399                 if (i == dir && side == 1) {
1400                         ind |= (ori[i] <= (st[i] + len) ? 0 : 1);
1401                 }
1402                 else {
1403                         ind |= (ori[i] < (st[i] + len) ? 0 : 1);
1404                 }
1405         }
1406
1407 #ifdef IN_DEBUG_MODE
1408         // dc_printf("In LOCATECELL index of ori(%d %d %d) with dir %d side %d in st(%d %d %d, %d) is: %d\n",
1409         //      ori[0], ori[1], ori[2], dir, side, st[0], st[1], st[2], len, ind);
1410 #endif
1411
1412         rst[0] = st[0] + vertmap[ind][0] * len;
1413         rst[1] = st[1] + vertmap[ind][1] * len;
1414         rst[2] = st[2] + vertmap[ind][2] * len;
1415
1416         if (node->has_child(ind)) {
1417                 int count = node->get_child_count(ind);
1418                 Node *chd = node->get_child(count);
1419                 if (node->is_child_leaf(ind)) {
1420                         rleaf = chd;
1421                         rlen = len;
1422                 }
1423                 else {
1424                         // Recur
1425                         node->set_child(count, locateCell(&chd->internal, rst, len, ori, dir, side, rleaf, rst, rlen));
1426                 }
1427         }
1428         else {
1429                 // Create a new child here
1430                 if (len == mindimen) {
1431                         LeafNode *chd = createLeaf(0);
1432                         node = addChild(node, ind, (Node *)chd, 1);
1433                         rleaf = (Node *)chd;
1434                         rlen = len;
1435                 }
1436                 else {
1437                         // Subdivide the empty cube
1438                         InternalNode *chd = createInternal(0);
1439                         node = addChild(node, ind,
1440                                         locateCell(chd, rst, len, ori, dir, side, rleaf, rst, rlen), 0);
1441                 }
1442         }
1443
1444 #ifdef IN_DEBUG_MODE
1445         // dc_printf("Return from LOCATECELL with node ");
1446         // printNode(newnode);
1447 #endif
1448         return (Node *)node;
1449 }
1450
1451 void Octree::checkElement(PathElement *ele)
1452 {
1453         /*
1454            if(ele != NULL && locateLeafCheck(ele->pos) != ele->node) {
1455             dc_printf("Screwed! at pos: %d %d %d\n", ele->pos[0]>>minshift, ele->pos[1]>>minshift, ele->pos[2]>>minshift);
1456             exit(0);
1457            }
1458          */
1459 }
1460
1461 void Octree::checkPath(PathElement *path)
1462 {
1463         PathElement *n = path;
1464         int same = 0;
1465         while (n && (same == 0 || n != path)) {
1466                 same++;
1467                 checkElement(n);
1468                 n = n->next;
1469         }
1470
1471 }
1472
1473 void Octree::testFacePoint(PathElement *e1, PathElement *e2)
1474 {
1475         int i;
1476         PathElement *e = NULL;
1477         for (i = 0; i < 3; i++) {
1478                 if (e1->pos[i] != e2->pos[i]) {
1479                         if (e1->pos[i] < e2->pos[i]) {
1480                                 e = e2;
1481                         }
1482                         else {
1483                                 e = e1;
1484                         }
1485                         break;
1486                 }
1487         }
1488
1489         int x, y;
1490         float p, q;
1491         dc_printf("Test.");
1492         getFacePoint(e, i, x, y, p, q);
1493 }
1494
1495 void Octree::getFacePoint(PathElement *leaf, int dir, int& x, int& y, float& p, float& q)
1496 {
1497         // Find average intersections
1498         float avg[3] = {0, 0, 0};
1499         float off[3];
1500         int num = 0, num2 = 0;
1501
1502         LeafNode *leafnode = locateLeaf(leaf->pos);
1503         for (int i = 0; i < 4; i++) {
1504                 int edgeind = faceMap[dir * 2][i];
1505                 int nst[3];
1506                 for (int j = 0; j < 3; j++) {
1507                         nst[j] = leaf->pos[j] + mindimen * vertmap[edgemap[edgeind][0]][j];
1508                 }
1509
1510                 if (getEdgeIntersectionByIndex(nst, edgeind / 4, off, 1)) {
1511                         avg[0] += off[0];
1512                         avg[1] += off[1];
1513                         avg[2] += off[2];
1514                         num++;
1515                 }
1516                 if (getEdgeParity(leafnode, edgeind)) {
1517                         num2++;
1518                 }
1519         }
1520         if (num == 0) {
1521                 dc_printf("Wrong! dir: %d pos: %d %d %d num: %d\n", dir, leaf->pos[0] >> minshift, leaf->pos[1] >> minshift, leaf->pos[2] >> minshift, num2);
1522                 avg[0] = (float) leaf->pos[0];
1523                 avg[1] = (float) leaf->pos[1];
1524                 avg[2] = (float) leaf->pos[2];
1525         }
1526         else {
1527
1528                 avg[0] /= num;
1529                 avg[1] /= num;
1530                 avg[2] /= num;
1531
1532                 //avg[0] =(float) leaf->pos[0];
1533                 //avg[1] =(float) leaf->pos[1];
1534                 //avg[2] =(float) leaf->pos[2];
1535         }
1536
1537         int xdir = (dir + 1) % 3;
1538         int ydir = (dir + 2) % 3;
1539
1540         float xf = avg[xdir];
1541         float yf = avg[ydir];
1542
1543 #ifdef IN_DEBUG_MODE
1544         // Is it outside?
1545         // PathElement* leaf = leaf1->len < leaf2->len ? leaf1 : leaf2;
1546         /*
1547            float* m =(leaf == leaf1 ? m1 : m2);
1548            if(xf < leaf->pos[xdir] ||
1549              yf < leaf->pos[ydir] ||
1550              xf > leaf->pos[xdir] + leaf->len ||
1551              yf > leaf->pos[ydir] + leaf->len) {
1552             dc_printf("Outside cube(%d %d %d), %d : %d %d %f %f\n", leaf->pos[0], leaf->pos[1], leaf->pos[2], leaf->len,
1553                             pos, dir, xf, yf);
1554
1555             // For now, snap to cell
1556             xf = m[xdir];
1557             yf = m[ydir];
1558            }
1559          */
1560
1561         /*
1562            if(alpha < 0 || alpha > 1 ||
1563              xf < leaf->pos[xdir] || xf > leaf->pos[xdir] + leaf->len ||
1564              yf < leaf->pos[ydir] || yf > leaf->pos[ydir] + leaf->len) {
1565             dc_printf("Alpha: %f Address: %d and %d\n", alpha, leaf1->node, leaf2->node);
1566             dc_printf("GETFACEPOINT result:(%d %d %d) %d min:(%f %f %f);(%d %d %d) %d min:(%f %f %f).\n",
1567                 leaf1->pos[0], leaf1->pos[1], leaf1->pos[2], leaf1->len, m1[0], m1[1], m1[2],
1568                 leaf2->pos[0], leaf2->pos[1], leaf2->pos[2], leaf2->len, m2[0], m2[1], m2[2]);
1569             dc_printf("Face point at dir %d pos %d: %f %f\n", dir, pos, xf, yf);
1570            }
1571          */
1572 #endif
1573
1574
1575         // Get the integer and float part
1576         x = ((leaf->pos[xdir]) >> minshift);
1577         y = ((leaf->pos[ydir]) >> minshift);
1578
1579         p = (xf - leaf->pos[xdir]) / mindimen;
1580         q = (yf - leaf->pos[ydir]) / mindimen;
1581
1582
1583 #ifdef IN_DEBUG_MODE
1584         dc_printf("Face point at dir %d : %f %f\n", dir, xf, yf);
1585 #endif
1586 }
1587
1588 int Octree::findPair(PathElement *head, int pos, int dir, PathElement *& pre1, PathElement *& pre2)
1589 {
1590         int side = getSide(head, pos, dir);
1591         PathElement *cur = head;
1592         PathElement *anchor;
1593         PathElement *ppre1, *ppre2;
1594
1595         // Start from this face, find a pair
1596         anchor = cur;
1597         ppre1 = cur;
1598         cur = cur->next;
1599         while (cur != anchor && (getSide(cur, pos, dir) == side)) {
1600                 ppre1 = cur;
1601                 cur = cur->next;
1602         }
1603         if (cur == anchor) {
1604                 // No pair found
1605                 return side;
1606         }
1607
1608         side = getSide(cur, pos, dir);
1609         ppre2 = cur;
1610         cur = cur->next;
1611         while (getSide(cur, pos, dir) == side) {
1612                 ppre2 = cur;
1613                 cur = cur->next;
1614         }
1615
1616
1617         // Switch pre1 and pre2 if we start from the higher side
1618         if (side == -1) {
1619                 cur = ppre1;
1620                 ppre1 = ppre2;
1621                 ppre2 = cur;
1622         }
1623
1624         pre1 = ppre1;
1625         pre2 = ppre2;
1626
1627         return 0;
1628 }
1629
1630 int Octree::getSide(PathElement *e, int pos, int dir)
1631 {
1632         return (e->pos[dir] < pos ? -1 : 1);
1633 }
1634
1635 int Octree::isEqual(PathElement *e1, PathElement *e2)
1636 {
1637         return (e1->pos[0] == e2->pos[0] && e1->pos[1] == e2->pos[1] && e1->pos[2] == e2->pos[2]);
1638 }
1639
1640 void Octree::compressRing(PathElement *& ring)
1641 {
1642         if (ring == NULL) {
1643                 return;
1644         }
1645 #ifdef IN_DEBUG_MODE
1646         dc_printf("Call to COMPRESSRING with path: \n");
1647         printPath(ring);
1648 #endif
1649
1650         PathElement *cur = ring->next->next;
1651         PathElement *pre = ring->next;
1652         PathElement *prepre = ring;
1653         PathElement *anchor = prepre;
1654
1655         do {
1656                 while (isEqual(cur, prepre)) {
1657                         // Delete
1658                         if (cur == prepre) {
1659                                 // The ring has shrinked to a point
1660                                 delete pre;
1661                                 delete cur;
1662                                 anchor = NULL;
1663                                 break;
1664                         }
1665                         else {
1666                                 prepre->next = cur->next;
1667                                 delete pre;
1668                                 delete cur;
1669                                 pre = prepre->next;
1670                                 cur = pre->next;
1671                                 anchor = prepre;
1672                         }
1673                 }
1674
1675                 if (anchor == NULL) {
1676                         break;
1677                 }
1678
1679                 prepre = pre;
1680                 pre = cur;
1681                 cur = cur->next;
1682         } while (prepre != anchor);
1683
1684         ring = anchor;
1685
1686 #ifdef IN_DEBUG_MODE
1687         dc_printf("Return from COMPRESSRING with path: \n");
1688         printPath(ring);
1689 #endif
1690 }
1691
1692 void Octree::buildSigns()
1693 {
1694         // First build a lookup table
1695         // dc_printf("Building up look up table...\n");
1696         int size = 1 << 12;
1697         unsigned char table[1 << 12];
1698         for (int i = 0; i < size; i++) {
1699                 table[i] = 0;
1700         }
1701         for (int i = 0; i < 256; i++) {
1702                 int ind = 0;
1703                 for (int j = 11; j >= 0; j--) {
1704                         ind <<= 1;
1705                         if (((i >> edgemap[j][0]) & 1) ^ ((i >> edgemap[j][1]) & 1)) {
1706                                 ind |= 1;
1707                         }
1708                 }
1709
1710                 table[ind] = i;
1711         }
1712
1713         // Next, traverse the grid
1714         int sg = 1;
1715         int cube[8];
1716         buildSigns(table, root, 0, sg, cube);
1717 }
1718
1719 void Octree::buildSigns(unsigned char table[], Node *node, int isLeaf, int sg, int rvalue[8])
1720 {
1721         if (node == NULL) {
1722                 for (int i = 0; i < 8; i++) {
1723                         rvalue[i] = sg;
1724                 }
1725                 return;
1726         }
1727
1728         if (isLeaf == 0) {
1729                 // Internal node
1730                 Node *chd[8];
1731                 int leaf[8];
1732                 node->internal.fill_children(chd, leaf);
1733
1734                 // Get the signs at the corners of the first cube
1735                 rvalue[0] = sg;
1736                 int oris[8];
1737                 buildSigns(table, chd[0], leaf[0], sg, oris);
1738
1739                 // Get the rest
1740                 int cube[8];
1741                 for (int i = 1; i < 8; i++) {
1742                         buildSigns(table, chd[i], leaf[i], oris[i], cube);
1743                         rvalue[i] = cube[i];
1744                 }
1745
1746         }
1747         else {
1748                 // Leaf node
1749                 generateSigns(&node->leaf, table, sg);
1750
1751                 for (int i = 0; i < 8; i++) {
1752                         rvalue[i] = getSign(&node->leaf, i);
1753                 }
1754         }
1755 }
1756
1757 void Octree::floodFill()
1758 {
1759         // int threshold =(int)((dimen/mindimen) *(dimen/mindimen) * 0.5f);
1760         int st[3] = {0, 0, 0};
1761
1762         // First, check for largest component
1763         // size stored in -threshold
1764         clearProcessBits(root, maxDepth);
1765         int threshold = floodFill(root, st, dimen, maxDepth, 0);
1766
1767         // Next remove
1768         dc_printf("Largest component: %d\n", threshold);
1769         threshold *= thresh;
1770         dc_printf("Removing all components smaller than %d\n", threshold);
1771
1772         int st2[3] = {0, 0, 0};
1773         clearProcessBits(root, maxDepth);
1774         floodFill(root, st2, dimen, maxDepth, threshold);
1775
1776 }
1777
1778 void Octree::clearProcessBits(Node *node, int height)
1779 {
1780         int i;
1781
1782         if (height == 0) {
1783                 // Leaf cell,
1784                 for (i = 0; i < 12; i++) {
1785                         setOutProcess(&node->leaf, i);
1786                 }
1787         }
1788         else {
1789                 // Internal cell, recur
1790                 int count = 0;
1791                 for (i = 0; i < 8; i++) {
1792                         if (node->internal.has_child(i)) {
1793                                 clearProcessBits(node->internal.get_child(count), height - 1);
1794                                 count++;
1795                         }
1796                 }
1797         }
1798 }
1799
1800 int Octree::floodFill(LeafNode *leaf, int st[3], int len, int height, int threshold)
1801 {
1802         int i, j;
1803         int maxtotal = 0;
1804
1805         // Leaf cell,
1806         int par, inp;
1807
1808         // Test if the leaf has intersection edges
1809         for (i = 0; i < 12; i++) {
1810                 par = getEdgeParity(leaf, i);
1811                 inp = isInProcess(leaf, i);
1812
1813                 if (par == 1 && inp == 0) {
1814                         // Intersection edge, hasn't been processed
1815                         // Let's start filling
1816                         GridQueue *queue = new GridQueue();
1817                         int total = 1;
1818
1819                         // Set to in process
1820                         int mst[3];
1821                         mst[0] = st[0] + vertmap[edgemap[i][0]][0] * len;
1822                         mst[1] = st[1] + vertmap[edgemap[i][0]][1] * len;
1823                         mst[2] = st[2] + vertmap[edgemap[i][0]][2] * len;
1824                         int mdir = i / 4;
1825                         setInProcessAll(mst, mdir);
1826
1827                         // Put this edge into queue
1828                         queue->pushQueue(mst, mdir);
1829
1830                         // Queue processing
1831                         int nst[3], dir;
1832                         while (queue->popQueue(nst, dir) == 1) {
1833                                 // dc_printf("nst: %d %d %d, dir: %d\n", nst[0]/mindimen, nst[1]/mindimen, nst[2]/mindimen, dir);
1834                                 // locations
1835                                 int stMask[3][3] = {
1836                                         {0, 0 - len, 0 - len},
1837                                         {0 - len, 0, 0 - len},
1838                                         {0 - len, 0 - len, 0}
1839                                 };
1840                                 int cst[2][3];
1841                                 for (j = 0; j < 3; j++) {
1842                                         cst[0][j] = nst[j];
1843                                         cst[1][j] = nst[j] + stMask[dir][j];
1844                                 }
1845
1846                                 // cells
1847                                 LeafNode *cs[2];
1848                                 for (j = 0; j < 2; j++) {
1849                                         cs[j] = locateLeaf(cst[j]);
1850                                 }
1851
1852                                 // Middle sign
1853                                 int s = getSign(cs[0], 0);
1854
1855                                 // Masks
1856                                 int fcCells[4] = {1, 0, 1, 0};
1857                                 int fcEdges[3][4][3] = {
1858                                         {{9, 2, 11}, {8, 1, 10}, {5, 1, 7}, {4, 2, 6}},
1859                                         {{10, 6, 11}, {8, 5, 9}, {1, 5, 3}, {0, 6, 2}},
1860                                         {{6, 10, 7}, {4, 9, 5}, {2, 9, 3}, {0, 10, 1}}
1861                                 };
1862
1863                                 // Search for neighboring connected intersection edges
1864                                 for (int find = 0; find < 4; find++) {
1865                                         int cind = fcCells[find];
1866                                         int eind, edge;
1867                                         if (s == 0) {
1868                                                 // Original order
1869                                                 for (eind = 0; eind < 3; eind++) {
1870                                                         edge = fcEdges[dir][find][eind];
1871                                                         if (getEdgeParity(cs[cind], edge) == 1) {
1872                                                                 break;
1873                                                         }
1874                                                 }
1875                                         }
1876                                         else {
1877                                                 // Inverse order
1878                                                 for (eind = 2; eind >= 0; eind--) {
1879                                                         edge = fcEdges[dir][find][eind];
1880                                                         if (getEdgeParity(cs[cind], edge) == 1) {
1881                                                                 break;
1882                                                         }
1883                                                 }
1884                                         }
1885
1886                                         if (eind == 3 || eind == -1) {
1887                                                 dc_printf("Wrong! this is not a consistent sign. %d\n", eind);
1888                                         }
1889                                         else {
1890                                                 int est[3];
1891                                                 est[0] = cst[cind][0] + vertmap[edgemap[edge][0]][0] * len;
1892                                                 est[1] = cst[cind][1] + vertmap[edgemap[edge][0]][1] * len;
1893                                                 est[2] = cst[cind][2] + vertmap[edgemap[edge][0]][2] * len;
1894                                                 int edir = edge / 4;
1895
1896                                                 if (isInProcess(cs[cind], edge) == 0) {
1897                                                         setInProcessAll(est, edir);
1898                                                         queue->pushQueue(est, edir);
1899                                                         // dc_printf("Pushed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir);
1900                                                         total++;
1901                                                 }
1902                                         }
1903                                 }
1904                         }
1905
1906                         dc_printf("Size of component: %d ", total);
1907
1908                         if (threshold == 0) {
1909                                 // Measuring stage
1910                                 if (total > maxtotal) {
1911                                         maxtotal = total;
1912                                 }
1913                                 dc_printf(".\n");
1914                                 delete queue;
1915                                 continue;
1916                         }
1917
1918                         if (total >= threshold) {
1919                                 dc_printf("Maintained.\n");
1920                                 delete queue;
1921                                 continue;
1922                         }
1923                         dc_printf("Less then %d, removing...\n", threshold);
1924
1925                         // We have to remove this noise
1926
1927                         // Flip parity
1928                         // setOutProcessAll(mst, mdir);
1929                         flipParityAll(mst, mdir);
1930
1931                         // Put this edge into queue
1932                         queue->pushQueue(mst, mdir);
1933
1934                         // Queue processing
1935                         while (queue->popQueue(nst, dir) == 1) {
1936                                 // dc_printf("nst: %d %d %d, dir: %d\n", nst[0]/mindimen, nst[1]/mindimen, nst[2]/mindimen, dir);
1937                                 // locations
1938                                 int stMask[3][3] = {
1939                                         {0, 0 - len, 0 - len},
1940                                         {0 - len, 0, 0 - len},
1941                                         {0 - len, 0 - len, 0}
1942                                 };
1943                                 int cst[2][3];
1944                                 for (j = 0; j < 3; j++) {
1945                                         cst[0][j] = nst[j];
1946                                         cst[1][j] = nst[j] + stMask[dir][j];
1947                                 }
1948
1949                                 // cells
1950                                 LeafNode *cs[2];
1951                                 for (j = 0; j < 2; j++)
1952                                         cs[j] = locateLeaf(cst[j]);
1953
1954                                 // Middle sign
1955                                 int s = getSign(cs[0], 0);
1956
1957                                 // Masks
1958                                 int fcCells[4] = {1, 0, 1, 0};
1959                                 int fcEdges[3][4][3] = {
1960                                         {{9, 2, 11}, {8, 1, 10}, {5, 1, 7}, {4, 2, 6}},
1961                                         {{10, 6, 11}, {8, 5, 9}, {1, 5, 3}, {0, 6, 2}},
1962                                         {{6, 10, 7}, {4, 9, 5}, {2, 9, 3}, {0, 10, 1}}
1963                                 };
1964
1965                                 // Search for neighboring connected intersection edges
1966                                 for (int find = 0; find < 4; find++) {
1967                                         int cind = fcCells[find];
1968                                         int eind, edge;
1969                                         if (s == 0) {
1970                                                 // Original order
1971                                                 for (eind = 0; eind < 3; eind++) {
1972                                                         edge = fcEdges[dir][find][eind];
1973                                                         if (isInProcess(cs[cind], edge) == 1) {
1974                                                                 break;
1975                                                         }
1976                                                 }
1977                                         }
1978                                         else {
1979                                                 // Inverse order
1980                                                 for (eind = 2; eind >= 0; eind--) {
1981                                                         edge = fcEdges[dir][find][eind];
1982                                                         if (isInProcess(cs[cind], edge) == 1) {
1983                                                                 break;
1984                                                         }
1985                                                 }
1986                                         }
1987
1988                                         if (eind == 3 || eind == -1) {
1989                                                 dc_printf("Wrong! this is not a consistent sign. %d\n", eind);
1990                                         }
1991                                         else {
1992                                                 int est[3];
1993                                                 est[0] = cst[cind][0] + vertmap[edgemap[edge][0]][0] * len;
1994                                                 est[1] = cst[cind][1] + vertmap[edgemap[edge][0]][1] * len;
1995                                                 est[2] = cst[cind][2] + vertmap[edgemap[edge][0]][2] * len;
1996                                                 int edir = edge / 4;
1997
1998                                                 if (getEdgeParity(cs[cind], edge) == 1) {
1999                                                         flipParityAll(est, edir);
2000                                                         queue->pushQueue(est, edir);
2001                                                         // dc_printf("Pushed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir);
2002                                                         total++;
2003                                                 }
2004                                         }
2005                                 }
2006                         }
2007
2008                         delete queue;
2009                 }
2010         }
2011
2012         return maxtotal;
2013 }
2014
2015 int Octree::floodFill(Node *node, int st[3], int len, int height, int threshold)
2016 {
2017         int i;
2018         int maxtotal = 0;
2019
2020         if (height == 0) {
2021                 maxtotal = floodFill(&node->leaf, st, len, height, threshold);
2022         }
2023         else {
2024                 // Internal cell, recur
2025                 int count = 0;
2026                 len >>= 1;
2027                 for (i = 0; i < 8; i++) {
2028                         if (node->internal.has_child(i)) {
2029                                 int nst[3];
2030                                 nst[0] = st[0] + vertmap[i][0] * len;
2031                                 nst[1] = st[1] + vertmap[i][1] * len;
2032                                 nst[2] = st[2] + vertmap[i][2] * len;
2033
2034                                 int d = floodFill(node->internal.get_child(count), nst, len, height - 1, threshold);
2035                                 if (d > maxtotal) {
2036                                         maxtotal = d;
2037                                 }
2038                                 count++;
2039                         }
2040                 }
2041         }
2042
2043
2044         return maxtotal;
2045
2046 }
2047
2048 void Octree::writeOut()
2049 {
2050         int numQuads = 0;
2051         int numVertices = 0;
2052         int numEdges = 0;
2053
2054         countIntersection(root, maxDepth, numQuads, numVertices, numEdges);
2055
2056         dc_printf("Vertices counted: %d Polys counted: %d \n", numVertices, numQuads);
2057         output_mesh = alloc_output(numVertices, numQuads);
2058         int offset = 0;
2059         int st[3] = {0, 0, 0};
2060
2061         // First, output vertices
2062         offset = 0;
2063         actualVerts = 0;
2064         actualQuads = 0;
2065
2066         generateMinimizer(root, st, dimen, maxDepth, offset);
2067         cellProcContour(root, 0, maxDepth);
2068         dc_printf("Vertices written: %d Quads written: %d \n", offset, actualQuads);
2069 }
2070
2071 void Octree::countIntersection(Node *node, int height, int& nedge, int& ncell, int& nface)
2072 {
2073         if (height > 0) {
2074                 int total = node->internal.get_num_children();
2075                 for (int i = 0; i < total; i++) {
2076                         countIntersection(node->internal.get_child(i), height - 1, nedge, ncell, nface);
2077                 }
2078         }
2079         else {
2080                 nedge += getNumEdges2(&node->leaf);
2081
2082                 int smask = getSignMask(&node->leaf);
2083
2084                 if (use_manifold) {
2085                         int comps = manifold_table[smask].comps;
2086                         ncell += comps;
2087                 }
2088                 else {
2089                         if (smask > 0 && smask < 255) {
2090                                 ncell++;
2091                         }
2092                 }
2093
2094                 for (int i = 0; i < 3; i++) {
2095                         if (getFaceEdgeNum(&node->leaf, i * 2)) {
2096                                 nface++;
2097                         }
2098                 }
2099         }
2100 }
2101
2102 /* from http://eigen.tuxfamily.org/bz/show_bug.cgi?id=257 */
2103 template<typename _Matrix_Type_>
2104 void pseudoInverse(const _Matrix_Type_ &a,
2105                    _Matrix_Type_ &result,
2106                    double epsilon = std::numeric_limits<typename _Matrix_Type_::Scalar>::epsilon())
2107 {
2108         Eigen::JacobiSVD< _Matrix_Type_ > svd = a.jacobiSvd(Eigen::ComputeFullU |
2109                                                             Eigen::ComputeFullV);
2110
2111         typename _Matrix_Type_::Scalar tolerance = epsilon * std::max(a.cols(),
2112                                                                       a.rows()) *
2113                                                    svd.singularValues().array().abs().maxCoeff();
2114
2115         result = svd.matrixV() *
2116                  _Matrix_Type_((svd.singularValues().array().abs() >
2117                                 tolerance).select(svd.singularValues().
2118                                                   array().inverse(), 0)).asDiagonal() *
2119                  svd.matrixU().adjoint();
2120 }
2121
2122 static void solve_least_squares(const float halfA[], const float b[],
2123                                 const float midpoint[], float rvalue[])
2124 {
2125         /* calculate pseudo-inverse */
2126         Eigen::MatrixXf A(3, 3), pinv(3, 3);
2127         A << halfA[0], halfA[1], halfA[2],
2128             halfA[1], halfA[3], halfA[4],
2129             halfA[2], halfA[4], halfA[5];
2130         pseudoInverse(A, pinv);
2131
2132         Eigen::Vector3f b2(b), mp(midpoint), result;
2133         b2 = b2 + A * -mp;
2134         result = pinv * b2 + mp;
2135
2136         for (int i = 0; i < 3; i++)
2137                 rvalue[i] = result(i);
2138 }
2139
2140 static void mass_point(float mp[3], const float pts[12][3], const int parity[12])
2141 {
2142         int ec = 0;
2143
2144         for (int i = 0; i < 12; i++) {
2145                 if (parity[i]) {
2146                         const float *p = pts[i];
2147
2148                         mp[0] += p[0];
2149                         mp[1] += p[1];
2150                         mp[2] += p[2];
2151
2152                         ec++;
2153                 }
2154         }
2155
2156         if (ec == 0) {
2157                 return;
2158         }
2159         mp[0] /= ec;
2160         mp[1] /= ec;
2161         mp[2] /= ec;
2162 }
2163
2164 static void minimize(float rvalue[3], float mp[3], const float pts[12][3],
2165                                          const float norms[12][3], const int parity[12])
2166 {
2167         float ata[6] = {0, 0, 0, 0, 0, 0};
2168         float atb[3] = {0, 0, 0};
2169         int ec = 0;
2170
2171         for (int i = 0; i < 12; i++) {
2172                 // if(getEdgeParity(leaf, i))
2173                 if (parity[i]) {
2174                         const float *norm = norms[i];
2175                         const float *p = pts[i];
2176
2177                         // QEF
2178                         ata[0] += (float)(norm[0] * norm[0]);
2179                         ata[1] += (float)(norm[0] * norm[1]);
2180                         ata[2] += (float)(norm[0] * norm[2]);
2181                         ata[3] += (float)(norm[1] * norm[1]);
2182                         ata[4] += (float)(norm[1] * norm[2]);
2183                         ata[5] += (float)(norm[2] * norm[2]);
2184
2185                         const float pn = p[0] * norm[0] + p[1] * norm[1] + p[2] * norm[2];
2186
2187                         atb[0] += (float)(norm[0] * pn);
2188                         atb[1] += (float)(norm[1] * pn);
2189                         atb[2] += (float)(norm[2] * pn);
2190
2191                         // Minimizer
2192                         mp[0] += p[0];
2193                         mp[1] += p[1];
2194                         mp[2] += p[2];
2195
2196                         ec++;
2197                 }
2198         }
2199
2200         if (ec == 0) {
2201                 return;
2202         }
2203         mp[0] /= ec;
2204         mp[1] /= ec;
2205         mp[2] /= ec;
2206
2207         // Solve least squares
2208         solve_least_squares(ata, atb, mp, rvalue);
2209 }
2210
2211 void Octree::computeMinimizer(const LeafNode *leaf, int st[3], int len,
2212                                                           float rvalue[3]) const
2213 {
2214         // First, gather all edge intersections
2215         float pts[12][3], norms[12][3];
2216         int parity[12];
2217         fillEdgeIntersections(leaf, st, len, pts, norms, parity);
2218
2219         switch (mode) {
2220                 case DUALCON_CENTROID:
2221                         rvalue[0] = (float) st[0] + len / 2;
2222                         rvalue[1] = (float) st[1] + len / 2;
2223                         rvalue[2] = (float) st[2] + len / 2;
2224                         break;
2225
2226                 case DUALCON_MASS_POINT:
2227                         rvalue[0] = rvalue[1] = rvalue[2] = 0;
2228                         mass_point(rvalue, pts, parity);
2229                         break;
2230
2231                 default: {
2232                         // Sharp features */
2233                         
2234                         // construct QEF and minimizer
2235                         float mp[3] = {0, 0, 0};
2236                         minimize(rvalue, mp, pts, norms, parity);
2237
2238                         /* Restraining the location of the minimizer */
2239                         float nh1 = hermite_num * len;
2240                         float nh2 = (1 + hermite_num) * len;
2241
2242                         if (rvalue[0] < st[0] - nh1 ||
2243                                 rvalue[1] < st[1] - nh1 ||
2244                                 rvalue[2] < st[2] - nh1 ||
2245
2246                                 rvalue[0] > st[0] + nh2 ||
2247                                 rvalue[1] > st[1] + nh2 ||
2248                                 rvalue[2] > st[2] + nh2)
2249                         {
2250                                 // Use mass point instead
2251                                 rvalue[0] = mp[0];
2252                                 rvalue[1] = mp[1];
2253                                 rvalue[2] = mp[2];
2254                         }
2255                         break;
2256                 }
2257         }
2258 }
2259
2260 void Octree::generateMinimizer(Node *node, int st[3], int len, int height, int& offset)
2261 {
2262         int i, j;
2263
2264         if (height == 0) {
2265                 // Leaf cell, generate
2266
2267                 // First, find minimizer
2268                 float rvalue[3];
2269                 rvalue[0] = (float) st[0] + len / 2;
2270                 rvalue[1] = (float) st[1] + len / 2;
2271                 rvalue[2] = (float) st[2] + len / 2;
2272                 computeMinimizer(&node->leaf, st, len, rvalue);
2273
2274                 // Update
2275                 //float fnst[3];
2276                 for (j = 0; j < 3; j++) {
2277                         rvalue[j] = rvalue[j] * range / dimen + origin[j];
2278                         //fnst[j] = st[j] * range / dimen + origin[j];
2279                 }
2280
2281                 int mult = 0, smask = getSignMask(&node->leaf);
2282
2283                 if (use_manifold) {
2284                         mult = manifold_table[smask].comps;
2285                 }
2286                 else {
2287                         if (smask > 0 && smask < 255) {
2288                                 mult = 1;
2289                         }
2290                 }
2291
2292                 for (j = 0; j < mult; j++) {
2293                         add_vert(output_mesh, rvalue);
2294                 }
2295
2296                 // Store the index
2297                 setMinimizerIndex(&node->leaf, offset);
2298
2299                 offset += mult;
2300         }
2301         else {
2302                 // Internal cell, recur
2303                 int count = 0;
2304                 len >>= 1;
2305                 for (i = 0; i < 8; i++) {
2306                         if (node->internal.has_child(i)) {
2307                                 int nst[3];
2308                                 nst[0] = st[0] + vertmap[i][0] * len;
2309                                 nst[1] = st[1] + vertmap[i][1] * len;
2310                                 nst[2] = st[2] + vertmap[i][2] * len;
2311
2312                                 generateMinimizer(node->internal.get_child(count),
2313                                                   nst, len, height - 1, offset);
2314                                 count++;
2315                         }
2316                 }
2317         }
2318 }
2319
2320 void Octree::processEdgeWrite(Node *node[4], int depth[4], int maxdep, int dir)
2321 {
2322         //int color = 0;
2323
2324         int i = 3;
2325         {
2326                 if (getEdgeParity((LeafNode *)(node[i]), processEdgeMask[dir][i])) {
2327                         int flip = 0;
2328                         int edgeind = processEdgeMask[dir][i];
2329                         if (getSign((LeafNode *)node[i], edgemap[edgeind][1]) > 0) {
2330                                 flip = 1;
2331                         }
2332
2333                         int num = 0;
2334                         {
2335                                 int ind[8];
2336                                 if (use_manifold) {
2337                                         int vind[2];
2338                                         int seq[4] = {0, 1, 3, 2};
2339                                         for (int k = 0; k < 4; k++) {
2340                                                 getMinimizerIndices((LeafNode *)(node[seq[k]]), processEdgeMask[dir][seq[k]], vind);
2341                                                 ind[num] = vind[0];
2342                                                 num++;
2343
2344                                                 if (vind[1] != -1) {
2345                                                         ind[num] = vind[1];
2346                                                         num++;
2347                                                         if (flip == 0) {
2348                                                                 ind[num - 1] = vind[0];
2349                                                                 ind[num - 2] = vind[1];
2350                                                         }
2351                                                 }
2352                                         }
2353
2354                                         /* we don't use the manifold option, but if it is
2355                                            ever enabled again note that it can output
2356                                            non-quads */
2357                                 }
2358                                 else {
2359                                         if (flip) {
2360                                                 ind[0] = getMinimizerIndex((LeafNode *)(node[2]));
2361                                                 ind[1] = getMinimizerIndex((LeafNode *)(node[3]));
2362                                                 ind[2] = getMinimizerIndex((LeafNode *)(node[1]));
2363                                                 ind[3] = getMinimizerIndex((LeafNode *)(node[0]));
2364                                         }
2365                                         else {
2366                                                 ind[0] = getMinimizerIndex((LeafNode *)(node[0]));
2367                                                 ind[1] = getMinimizerIndex((LeafNode *)(node[1]));
2368                                                 ind[2] = getMinimizerIndex((LeafNode *)(node[3]));
2369                                                 ind[3] = getMinimizerIndex((LeafNode *)(node[2]));
2370                                         }
2371
2372                                         add_quad(output_mesh, ind);
2373                                 }
2374                         }
2375                         return;
2376                 }
2377                 else {
2378                         return;
2379                 }
2380         }
2381 }
2382
2383
2384 void Octree::edgeProcContour(Node *node[4], int leaf[4], int depth[4], int maxdep, int dir)
2385 {
2386         if (!(node[0] && node[1] && node[2] && node[3])) {
2387                 return;
2388         }
2389         if (leaf[0] && leaf[1] && leaf[2] && leaf[3]) {
2390                 processEdgeWrite(node, depth, maxdep, dir);
2391         }
2392         else {
2393                 int i, j;
2394                 Node *chd[4][8];
2395                 for (j = 0; j < 4; j++) {
2396                         for (i = 0; i < 8; i++) {
2397                                 chd[j][i] = ((!leaf[j]) && node[j]->internal.has_child(i)) ?
2398                                             node[j]->internal.get_child(
2399                                                      node[j]->internal.get_child_count(i)) : NULL;
2400                         }
2401                 }
2402
2403                 // 2 edge calls
2404                 Node *ne[4];
2405                 int le[4];
2406                 int de[4];
2407                 for (i = 0; i < 2; i++) {
2408                         int c[4] = {edgeProcEdgeMask[dir][i][0],
2409                                         edgeProcEdgeMask[dir][i][1],
2410                                         edgeProcEdgeMask[dir][i][2],
2411                                         edgeProcEdgeMask[dir][i][3]};
2412
2413                         for (int j = 0; j < 4; j++) {
2414                                 if (leaf[j]) {
2415                                         le[j] = leaf[j];
2416                                         ne[j] = node[j];
2417                                         de[j] = depth[j];
2418                                 }
2419                                 else {
2420                                         le[j] = node[j]->internal.is_child_leaf(c[j]);
2421                                         ne[j] = chd[j][c[j]];
2422                                         de[j] = depth[j] - 1;
2423                                 }
2424                         }
2425
2426                         edgeProcContour(ne, le, de, maxdep - 1, edgeProcEdgeMask[dir][i][4]);
2427                 }
2428
2429         }
2430 }
2431
2432 void Octree::faceProcContour(Node *node[2], int leaf[2], int depth[2], int maxdep, int dir)
2433 {
2434         if (!(node[0] && node[1])) {
2435                 return;
2436         }
2437
2438         if (!(leaf[0] && leaf[1])) {
2439                 int i, j;
2440                 // Fill children nodes
2441                 Node *chd[2][8];
2442                 for (j = 0; j < 2; j++) {
2443                         for (i = 0; i < 8; i++) {
2444                                 chd[j][i] = ((!leaf[j]) && node[j]->internal.has_child(i)) ?
2445                                             node[j]->internal.get_child(
2446                                                      node[j]->internal.get_child_count(i)) : NULL;
2447                         }
2448                 }
2449
2450                 // 4 face calls
2451                 Node *nf[2];
2452                 int df[2];
2453                 int lf[2];
2454                 for (i = 0; i < 4; i++) {
2455                         int c[2] = {faceProcFaceMask[dir][i][0], faceProcFaceMask[dir][i][1]};
2456                         for (int j = 0; j < 2; j++) {
2457                                 if (leaf[j]) {
2458                                         lf[j] = leaf[j];
2459                                         nf[j] = node[j];
2460                                         df[j] = depth[j];
2461                                 }
2462                                 else {
2463                                         lf[j] = node[j]->internal.is_child_leaf(c[j]);
2464                                         nf[j] = chd[j][c[j]];
2465                                         df[j] = depth[j] - 1;
2466                                 }
2467                         }
2468                         faceProcContour(nf, lf, df, maxdep - 1, faceProcFaceMask[dir][i][2]);
2469                 }
2470
2471                 // 4 edge calls
2472                 int orders[2][4] = {{0, 0, 1, 1}, {0, 1, 0, 1}};
2473                 Node *ne[4];
2474                 int le[4];
2475                 int de[4];
2476
2477                 for (i = 0; i < 4; i++) {
2478                         int c[4] = {faceProcEdgeMask[dir][i][1], faceProcEdgeMask[dir][i][2],
2479                                         faceProcEdgeMask[dir][i][3], faceProcEdgeMask[dir][i][4]};
2480                         int *order = orders[faceProcEdgeMask[dir][i][0]];
2481
2482                         for (int j = 0; j < 4; j++) {
2483                                 if (leaf[order[j]]) {
2484                                         le[j] = leaf[order[j]];
2485                                         ne[j] = node[order[j]];
2486                                         de[j] = depth[order[j]];
2487                                 }
2488                                 else {
2489                                         le[j] = node[order[j]]->internal.is_child_leaf(c[j]);
2490                                         ne[j] = chd[order[j]][c[j]];
2491                                         de[j] = depth[order[j]] - 1;
2492                                 }
2493                         }
2494
2495                         edgeProcContour(ne, le, de, maxdep - 1, faceProcEdgeMask[dir][i][5]);
2496                 }
2497         }
2498 }
2499
2500
2501 void Octree::cellProcContour(Node *node, int leaf, int depth)
2502 {
2503         if (node == NULL) {
2504                 return;
2505         }
2506
2507         if (!leaf) {
2508                 int i;
2509
2510                 // Fill children nodes
2511                 Node *chd[8];
2512                 for (i = 0; i < 8; i++) {
2513                         chd[i] = ((!leaf) && node->internal.has_child(i)) ?
2514                                  node->internal.get_child(node->internal.get_child_count(i)) : NULL;
2515                 }
2516
2517                 // 8 Cell calls
2518                 for (i = 0; i < 8; i++) {
2519                         cellProcContour(chd[i], node->internal.is_child_leaf(i), depth - 1);
2520                 }
2521
2522                 // 12 face calls
2523                 Node *nf[2];
2524                 int lf[2];
2525                 int df[2] = {depth - 1, depth - 1};
2526                 for (i = 0; i < 12; i++) {
2527                         int c[2] = {cellProcFaceMask[i][0], cellProcFaceMask[i][1]};
2528
2529                         lf[0] = node->internal.is_child_leaf(c[0]);
2530                         lf[1] = node->internal.is_child_leaf(c[1]);
2531
2532                         nf[0] = chd[c[0]];
2533                         nf[1] = chd[c[1]];
2534
2535                         faceProcContour(nf, lf, df, depth - 1, cellProcFaceMask[i][2]);
2536                 }
2537
2538                 // 6 edge calls
2539                 Node *ne[4];
2540                 int le[4];
2541                 int de[4] = {depth - 1, depth - 1, depth - 1, depth - 1};
2542                 for (i = 0; i < 6; i++) {
2543                         int c[4] = {cellProcEdgeMask[i][0], cellProcEdgeMask[i][1], cellProcEdgeMask[i][2], cellProcEdgeMask[i][3]};
2544
2545                         for (int j = 0; j < 4; j++) {
2546                                 le[j] = node->internal.is_child_leaf(c[j]);
2547                                 ne[j] = chd[c[j]];
2548                         }
2549
2550                         edgeProcContour(ne, le, de, depth - 1, cellProcEdgeMask[i][4]);
2551                 }
2552         }
2553
2554 }
2555
2556 void Octree::processEdgeParity(LeafNode *node[4], int depth[4], int maxdep, int dir)
2557 {
2558         int con = 0;
2559         for (int i = 0; i < 4; i++) {
2560                 // Minimal cell
2561                 // if(op == 0)
2562                 {
2563                         if (getEdgeParity(node[i], processEdgeMask[dir][i])) {
2564                                 con = 1;
2565                                 break;
2566                         }
2567                 }
2568         }
2569
2570         if (con == 1) {
2571                 for (int i = 0; i < 4; i++) {
2572                         setEdge(node[i], processEdgeMask[dir][i]);
2573                 }
2574         }
2575
2576 }
2577
2578 void Octree::edgeProcParity(Node *node[4], int leaf[4], int depth[4], int maxdep, int dir)
2579 {
2580         if (!(node[0] && node[1] && node[2] && node[3])) {
2581                 return;
2582         }
2583         if (leaf[0] && leaf[1] && leaf[2] && leaf[3]) {
2584                 processEdgeParity((LeafNode **)node, depth, maxdep, dir);
2585         }
2586         else {
2587                 int i, j;
2588                 Node *chd[4][8];
2589                 for (j = 0; j < 4; j++) {
2590                         for (i = 0; i < 8; i++) {
2591                                 chd[j][i] = ((!leaf[j]) && node[j]->internal.has_child(i)) ?
2592                                             node[j]->internal.get_child( node[j]->internal.get_child_count(i)) : NULL;
2593                         }
2594                 }
2595
2596                 // 2 edge calls
2597                 Node *ne[4];
2598                 int le[4];
2599                 int de[4];
2600                 for (i = 0; i < 2; i++) {
2601                         int c[4] = {edgeProcEdgeMask[dir][i][0],
2602                                         edgeProcEdgeMask[dir][i][1],
2603                                         edgeProcEdgeMask[dir][i][2],
2604                                         edgeProcEdgeMask[dir][i][3]};
2605
2606                         // int allleaf = 1;
2607                         for (int j = 0; j < 4; j++) {
2608
2609                                 if (leaf[j]) {
2610                                         le[j] = leaf[j];
2611                                         ne[j] = node[j];
2612                                         de[j] = depth[j];
2613                                 }
2614                                 else {
2615                                         le[j] = node[j]->internal.is_child_leaf(c[j]);
2616                                         ne[j] = chd[j][c[j]];
2617                                         de[j] = depth[j] - 1;
2618
2619                                 }
2620
2621                         }
2622
2623                         edgeProcParity(ne, le, de, maxdep - 1, edgeProcEdgeMask[dir][i][4]);
2624                 }
2625
2626         }
2627 }
2628
2629 void Octree::faceProcParity(Node *node[2], int leaf[2], int depth[2], int maxdep, int dir)
2630 {
2631         if (!(node[0] && node[1])) {
2632                 return;
2633         }
2634
2635         if (!(leaf[0] && leaf[1])) {
2636                 int i, j;
2637                 // Fill children nodes
2638                 Node *chd[2][8];
2639                 for (j = 0; j < 2; j++) {
2640                         for (i = 0; i < 8; i++) {
2641                                 chd[j][i] = ((!leaf[j]) && node[j]->internal.has_child(i)) ?
2642                                             node[j]->internal.get_child(
2643                                                      node[j]->internal.get_child_count(i)) : NULL;
2644                         }
2645                 }
2646
2647                 // 4 face calls
2648                 Node *nf[2];
2649                 int df[2];
2650                 int lf[2];
2651                 for (i = 0; i < 4; i++) {
2652                         int c[2] = {faceProcFaceMask[dir][i][0], faceProcFaceMask[dir][i][1]};
2653                         for (int j = 0; j < 2; j++) {
2654                                 if (leaf[j]) {
2655                                         lf[j] = leaf[j];
2656                                         nf[j] = node[j];
2657                                         df[j] = depth[j];
2658                                 }
2659                                 else {
2660                                         lf[j] = node[j]->internal.is_child_leaf(c[j]);
2661                                         nf[j] = chd[j][c[j]];
2662                                         df[j] = depth[j] - 1;
2663                                 }
2664                         }
2665                         faceProcParity(nf, lf, df, maxdep - 1, faceProcFaceMask[dir][i][2]);
2666                 }
2667
2668                 // 4 edge calls
2669                 int orders[2][4] = {{0, 0, 1, 1}, {0, 1, 0, 1}};
2670                 Node *ne[4];
2671                 int le[4];
2672                 int de[4];
2673
2674                 for (i = 0; i < 4; i++) {
2675                         int c[4] = {faceProcEdgeMask[dir][i][1], faceProcEdgeMask[dir][i][2],
2676                                         faceProcEdgeMask[dir][i][3], faceProcEdgeMask[dir][i][4]};
2677                         int *order = orders[faceProcEdgeMask[dir][i][0]];
2678
2679                         for (int j = 0; j < 4; j++) {
2680                                 if (leaf[order[j]]) {
2681                                         le[j] = leaf[order[j]];
2682                                         ne[j] = node[order[j]];
2683                                         de[j] = depth[order[j]];
2684                                 }
2685                                 else {
2686                                         le[j] = node[order[j]]->internal.is_child_leaf(c[j]);
2687                                         ne[j] = chd[order[j]][c[j]];
2688                                         de[j] = depth[order[j]] - 1;
2689                                 }
2690                         }
2691
2692                         edgeProcParity(ne, le, de, maxdep - 1, faceProcEdgeMask[dir][i][5]);
2693                 }
2694         }
2695 }
2696
2697
2698 void Octree::cellProcParity(Node *node, int leaf, int depth)
2699 {
2700         if (node == NULL) {
2701                 return;
2702         }
2703
2704         if (!leaf) {
2705                 int i;
2706
2707                 // Fill children nodes
2708                 Node *chd[8];
2709                 for (i = 0; i < 8; i++) {
2710                         chd[i] = ((!leaf) && node->internal.has_child(i)) ?
2711                                  node->internal.get_child(node->internal.get_child_count(i)) : NULL;
2712                 }
2713
2714                 // 8 Cell calls
2715                 for (i = 0; i < 8; i++) {
2716                         cellProcParity(chd[i], node->internal.is_child_leaf(i), depth - 1);
2717                 }
2718
2719                 // 12 face calls
2720                 Node *nf[2];
2721                 int lf[2];
2722                 int df[2] = {depth - 1, depth - 1};
2723                 for (i = 0; i < 12; i++) {
2724                         int c[2] = {cellProcFaceMask[i][0], cellProcFaceMask[i][1]};
2725
2726                         lf[0] = node->internal.is_child_leaf(c[0]);
2727                         lf[1] = node->internal.is_child_leaf(c[1]);
2728
2729                         nf[0] = chd[c[0]];
2730                         nf[1] = chd[c[1]];
2731
2732                         faceProcParity(nf, lf, df, depth - 1, cellProcFaceMask[i][2]);
2733                 }
2734
2735                 // 6 edge calls
2736                 Node *ne[4];
2737                 int le[4];
2738                 int de[4] = {depth - 1, depth - 1, depth - 1, depth - 1};
2739                 for (i = 0; i < 6; i++) {
2740                         int c[4] = {cellProcEdgeMask[i][0], cellProcEdgeMask[i][1], cellProcEdgeMask[i][2], cellProcEdgeMask[i][3]};
2741
2742                         for (int j = 0; j < 4; j++) {
2743                                 le[j] = node->internal.is_child_leaf(c[j]);
2744                                 ne[j] = chd[c[j]];
2745                         }
2746
2747                         edgeProcParity(ne, le, de, depth - 1, cellProcEdgeMask[i][4]);
2748                 }
2749         }
2750
2751 }
2752
2753 /* definitions for global arrays */
2754 const int edgemask[3] = {5, 3, 6};
2755
2756 const int faceMap[6][4] = {
2757         {4, 8, 5, 9},
2758         {6, 10, 7, 11},
2759         {0, 8, 1, 10},
2760         {2, 9, 3, 11},
2761         {0, 4, 2, 6},
2762         {1, 5, 3, 7}
2763 };
2764
2765 const int cellProcFaceMask[12][3] = {
2766         {0, 4, 0},
2767         {1, 5, 0},
2768         {2, 6, 0},
2769         {3, 7, 0},
2770         {0, 2, 1},
2771         {4, 6, 1},
2772         {1, 3, 1},
2773         {5, 7, 1},
2774         {0, 1, 2},
2775         {2, 3, 2},
2776         {4, 5, 2},
2777         {6, 7, 2}
2778 };
2779
2780 const int cellProcEdgeMask[6][5] = {
2781         {0, 1, 2, 3, 0},
2782         {4, 5, 6, 7, 0},
2783         {0, 4, 1, 5, 1},
2784         {2, 6, 3, 7, 1},
2785         {0, 2, 4, 6, 2},
2786         {1, 3, 5, 7, 2}
2787 };
2788
2789 const int faceProcFaceMask[3][4][3] = {
2790         {{4, 0, 0},
2791          {5, 1, 0},
2792          {6, 2, 0},
2793          {7, 3, 0}},
2794         {{2, 0, 1},
2795          {6, 4, 1},
2796          {3, 1, 1},
2797          {7, 5, 1}},
2798         {{1, 0, 2},
2799          {3, 2, 2},
2800          {5, 4, 2},
2801          {7, 6, 2}}
2802 };
2803
2804 const int faceProcEdgeMask[3][4][6] = {
2805         {{1, 4, 0, 5, 1, 1},
2806          {1, 6, 2, 7, 3, 1},
2807          {0, 4, 6, 0, 2, 2},
2808          {0, 5, 7, 1, 3, 2}},
2809         {{0, 2, 3, 0, 1, 0},
2810          {0, 6, 7, 4, 5, 0},
2811          {1, 2, 0, 6, 4, 2},
2812          {1, 3, 1, 7, 5, 2}},
2813         {{1, 1, 0, 3, 2, 0},
2814          {1, 5, 4, 7, 6, 0},
2815          {0, 1, 5, 0, 4, 1},
2816          {0, 3, 7, 2, 6, 1}}
2817 };
2818
2819 const int edgeProcEdgeMask[3][2][5] = {
2820         {{3, 2, 1, 0, 0},
2821          {7, 6, 5, 4, 0}},
2822         {{5, 1, 4, 0, 1},
2823          {7, 3, 6, 2, 1}},
2824         {{6, 4, 2, 0, 2},
2825          {7, 5, 3, 1, 2}},
2826 };
2827
2828 const int processEdgeMask[3][4] = {
2829         {3, 2, 1, 0},
2830         {7, 5, 6, 4},
2831         {11, 10, 9, 8}
2832 };
2833
2834 const int dirCell[3][4][3] = {
2835         {{0, -1, -1},
2836          {0, -1, 0},
2837          {0, 0, -1},
2838          {0, 0, 0}},
2839         {{-1, 0, -1},
2840          {-1, 0, 0},
2841          {0, 0, -1},
2842          {0, 0, 0}},
2843         {{-1, -1, 0},
2844          {-1, 0, 0},
2845          {0, -1, 0},
2846          {0, 0, 0}}
2847 };
2848
2849 const int dirEdge[3][4] = {
2850         {3, 2, 1, 0},
2851         {7, 6, 5, 4},
2852         {11, 10, 9, 8}
2853 };
2854
2855 int InternalNode::numChildrenTable[256];
2856 int InternalNode::childrenCountTable[256][8];
2857 int InternalNode::childrenIndexTable[256][8];