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