code cleanup: spelling
[blender.git] / intern / dualcon / intern / octree.h
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 #ifndef OCTREE_H
24 #define OCTREE_H
25
26 #include <cassert>
27 #include <cstring>
28 #include <stdio.h>
29 #include <math.h>
30 #include "GeoCommon.h"
31 #include "Projections.h"
32 #include "ModelReader.h"
33 #include "MemoryAllocator.h"
34 #include "cubes.h"
35 #include "Queue.h"
36 #include "manifold_table.h"
37 #include "dualcon.h"
38
39 /**
40  * Main class and structures for scan-convertion, sign-generation,
41  * and surface reconstruction.
42  *
43  * @author Tao Ju
44  */
45
46
47 /* Global defines */
48 // Uncomment to see debug information
49 // #define IN_DEBUG_MODE
50
51 // Uncomment to see more output messages during repair
52 // #define IN_VERBOSE_MODE
53
54 /* Set scan convert params */
55
56 #define EDGE_FLOATS 4
57
58 union Node;
59
60 struct InternalNode {
61         /* Treat as bitfield, bit N indicates whether child N exists or not */
62         unsigned char has_child;
63         /* Treat as bitfield, bit N indicates whether child N is a leaf or not */
64         unsigned char child_is_leaf;
65
66         /* Can have up to eight children */
67         Node *children[0];
68
69         /// Test if child is leaf
70         int is_child_leaf(int index) const
71         {
72                 return (child_is_leaf >> index) & 1;
73         }
74 };
75
76
77 /**
78  * Bits order
79  *
80  * Leaf node:
81  * Byte 0,1(0-11): edge parity
82  * Byte 1(4,5,6): mask of primary edges intersections stored
83  * Byte 1(7): in flood fill mode, whether the cell is in process
84  * Byte 2(0-8): signs
85  * Byte 3,4: in coloring mode, the mask for edges
86  * Byte 5: edge intersections(4 bytes per inter, or 12 bytes if USE_HERMIT)
87  */
88 struct LeafNode /* TODO: remove this attribute once everything is fixed */ {
89         unsigned short edge_parity : 12;
90         unsigned short primary_edge_intersections : 3;
91
92         /* XXX: maybe actually unused? */
93         unsigned short in_process : 1;
94
95         /* bitfield */
96         char signs;
97
98         int minimizer_index;
99
100         unsigned short flood_fill;
101
102         float edge_intersections[0];
103 };
104
105 /* Doesn't get directly allocated anywhere, just used for passing
106    pointers to nodes that could be internal or leaf. */
107 union Node {
108         InternalNode internal;
109         LeafNode leaf;
110 };
111
112 /* Global variables */
113 extern const int edgemask[3];
114 extern const int faceMap[6][4];
115 extern const int cellProcFaceMask[12][3];
116 extern const int cellProcEdgeMask[6][5];
117 extern const int faceProcFaceMask[3][4][3];
118 extern const int edgeProcEdgeMask[3][2][5];
119 extern const int faceProcEdgeMask[3][4][6];
120 extern const int processEdgeMask[3][4];
121 extern const int dirCell[3][4][3];
122 extern const int dirEdge[3][4];
123
124 /**
125  * Structures for detecting/patching open cycles on the dual surface
126  */
127
128 struct PathElement {
129         // Origin
130         int pos[3];
131
132         // link
133         PathElement *next;
134 };
135
136 struct PathList {
137         // Head
138         PathElement *head;
139         PathElement *tail;
140
141         // Length of the list
142         int length;
143
144         // Next list
145         PathList *next;
146 };
147
148
149 /**
150  * Class for building and processing an octree
151  */
152 class Octree
153 {
154  public:
155         /* Public members */
156
157         /// Memory allocators
158         VirtualMemoryAllocator *alloc[9];
159         VirtualMemoryAllocator *leafalloc[4];
160
161         /// Root node
162         Node *root;
163
164         /// Model reader
165         ModelReader *reader;
166
167         /// Marching cubes table
168         Cubes *cubes;
169
170         /// Length of grid
171         int dimen;
172         int mindimen, minshift;
173
174         /// Maximum depth
175         int maxDepth;
176
177         /// The lower corner of the bounding box and the size
178         float origin[3];
179         float range;
180
181         /// Counting information
182         int nodeCount;
183         int nodeSpace;
184         int nodeCounts[9];
185
186         int actualQuads, actualVerts;
187
188         PathList *ringList;
189
190         int maxTrianglePerCell;
191         int outType;     // 0 for OFF, 1 for PLY, 2 for VOL
192
193         // For flood filling
194         int use_flood_fill;
195         float thresh;
196
197         int use_manifold;
198
199         float hermite_num;
200
201         DualConMode mode;
202
203  public:
204         /**
205          * Construtor
206          */
207         Octree(ModelReader *mr,
208                    DualConAllocOutput alloc_output_func,
209                    DualConAddVert add_vert_func,
210                    DualConAddQuad add_quad_func,
211                    DualConFlags flags, DualConMode mode, int depth,
212                    float threshold, float hermite_num);
213
214         /**
215          * Destructor
216          */
217         ~Octree();
218
219         /**
220          * Scan convert
221          */
222         void scanConvert();
223
224         void *getOutputMesh() {
225                 return output_mesh;
226         }
227
228  private:
229         /* Helper functions */
230
231         /**
232          * Initialize memory allocators
233          */
234         void initMemory();
235
236         /**
237          * Release memory
238          */
239         void freeMemory();
240
241         /**
242          * Print memory usage
243          */
244         void printMemUsage();
245
246
247         /**
248          * Methods to set / restore minimum edges
249          */
250         void resetMinimalEdges();
251
252         void cellProcParity(Node *node, int leaf, int depth);
253         void faceProcParity(Node * node[2], int leaf[2], int depth[2], int maxdep, int dir);
254         void edgeProcParity(Node * node[4], int leaf[4], int depth[4], int maxdep, int dir);
255
256         void processEdgeParity(LeafNode * node[4], int depths[4], int maxdep, int dir);
257
258         /**
259          * Add triangles to the tree
260          */
261         void addAllTriangles();
262         void addTriangle(Triangle *trian, int triind);
263         InternalNode *addTriangle(InternalNode *node, CubeTriangleIsect *p, int height);
264
265         /**
266          * Method to update minimizer in a cell: update edge intersections instead
267          */
268         LeafNode *updateCell(LeafNode *node, CubeTriangleIsect *p);
269
270         /* Routines to detect and patch holes */
271         int numRings;
272         int totRingLengths;
273         int maxRingLength;
274
275         /**
276          * Entry routine.
277          */
278         void trace();
279         /**
280          * Trace the given node, find patches and fill them in
281          */
282         Node *trace(Node *node, int *st, int len, int depth, PathList *& paths);
283         /**
284          * Look for path on the face and add to paths
285          */
286         void findPaths(Node * node[2], int leaf[2], int depth[2], int *st[2], int maxdep, int dir, PathList * &paths);
287         /**
288          * Combine two list1 and list2 into list1 using connecting paths list3,
289          * while closed paths are appended to rings
290          */
291         void combinePaths(PathList *& list1, PathList *list2, PathList *paths, PathList *& rings);
292         /**
293          * Helper function: combine current paths in list1 and list2 to a single path and append to list3
294          */
295         PathList *combineSinglePath(PathList *& head1, PathList *pre1, PathList *& list1, PathList *& head2, PathList *pre2, PathList *& list2);
296
297         /**
298          * Functions to patch rings in a node
299          */
300         Node *patch(Node * node, int st[3], int len, PathList * rings);
301         Node *patchSplit(Node * node, int st[3], int len, PathList * rings, int dir, PathList * &nrings1, PathList * &nrings2);
302         Node *patchSplitSingle(Node * node, int st[3], int len, PathElement * head, int dir, PathList * &nrings1, PathList * &nrings2);
303         Node *connectFace(Node * node, int st[3], int len, int dir, PathElement * f1, PathElement * f2);
304         Node *locateCell(InternalNode * node, int st[3], int len, int ori[3], int dir, int side, Node * &rleaf, int rst[3], int& rlen);
305         void compressRing(PathElement *& ring);
306         void getFacePoint(PathElement *leaf, int dir, int& x, int& y, float& p, float& q);
307         LeafNode *patchAdjacent(InternalNode * node, int len, int st1[3], LeafNode * leaf1, int st2[3], LeafNode * leaf2, int walkdir, int inc, int dir, int side, float alpha);
308         int findPair(PathElement *head, int pos, int dir, PathElement *& pre1, PathElement *& pre2);
309         int getSide(PathElement *e, int pos, int dir);
310         int isEqual(PathElement *e1, PathElement *e2);
311         void preparePrimalEdgesMask(InternalNode *node);
312         void testFacePoint(PathElement *e1, PathElement *e2);
313
314         /**
315          * Path-related functions
316          */
317         void deletePath(PathList *& head, PathList *pre, PathList *& curr);
318         void printPath(PathList *path);
319         void printPath(PathElement *path);
320         void printElement(PathElement *ele);
321         void printPaths(PathList *path);
322         void checkElement(PathElement *ele);
323         void checkPath(PathElement *path);
324
325
326         /**
327          * Routines to build signs to create a partitioned volume
328          *(after patching rings)
329          */
330         void buildSigns();
331         void buildSigns(unsigned char table[], Node * node, int isLeaf, int sg, int rvalue[8]);
332
333         /************************************************************************/
334         /* To remove disconnected components */
335         /************************************************************************/
336         void floodFill();
337         void clearProcessBits(Node *node, int height);
338         int floodFill(LeafNode * leaf, int st[3], int len, int height, int threshold);
339         int floodFill(Node * node, int st[3], int len, int height, int threshold);
340
341         /**
342          * Write out polygon file
343          */
344         void writeOut();
345
346         void countIntersection(Node *node, int height, int& nedge, int& ncell, int& nface);
347         void generateMinimizer(Node * node, int st[3], int len, int height, int& offset);
348         void computeMinimizer(LeafNode * leaf, int st[3], int len, float rvalue[3]);
349         /**
350          * Traversal functions to generate polygon model
351          * op: 0 for counting, 1 for writing OBJ, 2 for writing OFF, 3 for writing PLY
352          */
353         void cellProcContour(Node *node, int leaf, int depth);
354         void faceProcContour(Node * node[2], int leaf[2], int depth[2], int maxdep, int dir);
355         void edgeProcContour(Node * node[4], int leaf[4], int depth[4], int maxdep, int dir);
356         void processEdgeWrite(Node * node[4], int depths[4], int maxdep, int dir);
357
358         /* output callbacks/data */
359         DualConAllocOutput alloc_output;
360         DualConAddVert add_vert;
361         DualConAddQuad add_quad;
362         void *output_mesh;
363
364  private:
365         /************ Operators for all nodes ************/
366
367         /// Lookup table
368         int numChildrenTable[256];
369         int childrenCountTable[256][8];
370         int childrenIndexTable[256][8];
371         int numEdgeTable[8];
372         int edgeCountTable[8][3];
373
374         /// Build up lookup table
375         void buildTable()
376         {
377                 for (int i = 0; i < 256; i++) {
378                         numChildrenTable[i] = 0;
379                         int count = 0;
380                         for (int j = 0; j < 8; j++) {
381                                 numChildrenTable[i] += ((i >> j) & 1);
382                                 childrenCountTable[i][j] = count;
383                                 childrenIndexTable[i][count] = j;
384                                 count += ((i >> j) & 1);
385                         }
386                 }
387
388                 for (int i = 0; i < 8; i++) {
389                         numEdgeTable[i] = 0;
390                         int count = 0;
391                         for (int j = 0; j < 3; j++) {
392                                 numEdgeTable[i] += ((i >> j) & 1);
393                                 edgeCountTable[i][j] = count;
394                                 count += ((i >> j) & 1);
395                         }
396                 }
397         }
398
399         int getSign(Node *node, int height, int index)
400         {
401                 if (height == 0) {
402                         return getSign(&node->leaf, index);
403                 }
404                 else {
405                         if (hasChild(&node->internal, index)) {
406                                 return getSign(getChild(&node->internal, getChildCount(&node->internal, index)),
407                                                            height - 1,
408                                                            index);
409                         }
410                         else {
411                                 return getSign(getChild(&node->internal, 0),
412                                                            height - 1,
413                                                            7 - getChildIndex(&node->internal, 0));
414                         }
415                 }
416         }
417
418         /************ Operators for leaf nodes ************/
419
420         void printInfo(int st[3])
421         {
422                 printf("INFO AT: %d %d %d\n", st[0] >> minshift, st[1] >> minshift, st[2] >> minshift);
423                 LeafNode *leaf = (LeafNode *)locateLeafCheck(st);
424                 if (leaf)
425                         printInfo(leaf);
426                 else
427                         printf("Leaf not exists!\n");
428         }
429
430         void printInfo(const LeafNode *leaf)
431         {
432                 /*
433                   printf("Edge mask: ");
434                   for(int i = 0; i < 12; i ++)
435                   {
436                   printf("%d ", getEdgeParity(leaf, i));
437                   }
438                   printf("\n");
439                   printf("Stored edge mask: ");
440                   for(i = 0; i < 3; i ++)
441                   {
442                   printf("%d ", getStoredEdgesParity(leaf, i));
443                   }
444                   printf("\n");
445                 */
446                 printf("Sign mask: ");
447                 for (int i = 0; i < 8; i++) {
448                         printf("%d ", getSign(leaf, i));
449                 }
450                 printf("\n");
451
452         }
453
454         /// Retrieve signs
455         int getSign(const LeafNode *leaf, int index)
456         {
457                 return ((leaf->signs >> index) & 1);
458         }
459
460         /// Set sign
461         void setSign(LeafNode *leaf, int index)
462         {
463                 leaf->signs |= (1 << index);
464         }
465
466         void setSign(LeafNode *leaf, int index, int sign)
467         {
468                 leaf->signs &= (~(1 << index));
469                 leaf->signs |= ((sign & 1) << index);
470         }
471
472         int getSignMask(const LeafNode *leaf)
473         {
474                 return leaf->signs;
475         }
476
477         void setInProcessAll(int st[3], int dir)
478         {
479                 int nst[3], eind;
480                 for (int i = 0; i < 4; i++) {
481                         nst[0] = st[0] + dirCell[dir][i][0] * mindimen;
482                         nst[1] = st[1] + dirCell[dir][i][1] * mindimen;
483                         nst[2] = st[2] + dirCell[dir][i][2] * mindimen;
484                         eind = dirEdge[dir][i];
485
486                         LeafNode *cell = locateLeafCheck(nst);
487                         assert(cell);
488
489                         setInProcess(cell, eind);
490                 }
491         }
492
493         void flipParityAll(int st[3], int dir)
494         {
495                 int nst[3], eind;
496                 for (int i = 0; i < 4; i++) {
497                         nst[0] = st[0] + dirCell[dir][i][0] * mindimen;
498                         nst[1] = st[1] + dirCell[dir][i][1] * mindimen;
499                         nst[2] = st[2] + dirCell[dir][i][2] * mindimen;
500                         eind = dirEdge[dir][i];
501
502                         LeafNode *cell = locateLeaf(nst);
503                         flipEdge(cell, eind);
504                 }
505         }
506
507         void setInProcess(LeafNode *leaf, int eind)
508         {
509                 assert(eind >= 0 && eind <= 11);
510
511                 leaf->flood_fill |= (1 << eind);
512         }
513
514         void setOutProcess(LeafNode *leaf, int eind)
515         {
516                 assert(eind >= 0 && eind <= 11);
517
518                 leaf->flood_fill &= ~(1 << eind);
519         }
520
521         int isInProcess(LeafNode *leaf, int eind)
522         {
523                 assert(eind >= 0 && eind <= 11);
524
525                 return (leaf->flood_fill >> eind) & 1;
526         }
527
528         /// Generate signs at the corners from the edge parity
529         void generateSigns(LeafNode *leaf, unsigned char table[], int start)
530         {
531                 leaf->signs = table[leaf->edge_parity];
532
533                 if ((start ^ leaf->signs) & 1) {
534                         leaf->signs = ~(leaf->signs);
535                 }
536         }
537
538         /// Get edge parity
539         int getEdgeParity(LeafNode *leaf, int index)
540         {
541                 assert(index >= 0 && index <= 11);
542
543                 return (leaf->edge_parity >> index) & 1;
544         }
545
546         /// Get edge parity on a face
547         int getFaceParity(LeafNode *leaf, int index)
548         {
549                 int a = getEdgeParity(leaf, faceMap[index][0]) +
550                 getEdgeParity(leaf, faceMap[index][1]) +
551                 getEdgeParity(leaf, faceMap[index][2]) +
552                 getEdgeParity(leaf, faceMap[index][3]);
553                 return (a & 1);
554         }
555         int getFaceEdgeNum(LeafNode *leaf, int index)
556         {
557                 int a = getEdgeParity(leaf, faceMap[index][0]) +
558                 getEdgeParity(leaf, faceMap[index][1]) +
559                 getEdgeParity(leaf, faceMap[index][2]) +
560                 getEdgeParity(leaf, faceMap[index][3]);
561                 return a;
562         }
563
564         /// Set edge parity
565         void flipEdge(LeafNode *leaf, int index)
566         {
567                 assert(index >= 0 && index <= 11);
568
569                 leaf->edge_parity ^= (1 << index);
570         }
571
572         /// Set 1
573         void setEdge(LeafNode *leaf, int index)
574         {
575                 assert(index >= 0 && index <= 11);
576
577                 leaf->edge_parity |= (1 << index);
578         }
579
580         /// Set 0
581         void resetEdge(LeafNode *leaf, int index)
582         {
583                 assert(index >= 0 && index <= 11);
584
585                 leaf->edge_parity &= ~(1 << index);
586         }
587
588         /// Flipping with a new intersection offset
589         void createPrimalEdgesMask(LeafNode *leaf)
590         {
591                 leaf->primary_edge_intersections = getPrimalEdgesMask2(leaf);
592         }
593
594         void setStoredEdgesParity(LeafNode *leaf, int pindex)
595         {
596                 assert(pindex <= 2 && pindex >= 0);
597
598                 leaf->primary_edge_intersections |= (1 << pindex);
599         }
600         int getStoredEdgesParity(LeafNode *leaf, int pindex)
601         {
602                 assert(pindex <= 2 && pindex >= 0);
603
604                 return (leaf->primary_edge_intersections >> pindex) & 1;
605         }
606
607         LeafNode *flipEdge(LeafNode *leaf, int index, float alpha)
608         {
609                 flipEdge(leaf, index);
610
611                 if ((index & 3) == 0) {
612                         int ind = index / 4;
613                         if (getEdgeParity(leaf, index) && !getStoredEdgesParity(leaf, ind)) {
614                                 // Create a new node
615                                 int num = getNumEdges(leaf) + 1;
616                                 setStoredEdgesParity(leaf, ind);
617                                 int count = getEdgeCount(leaf, ind);
618                                 LeafNode *nleaf = createLeaf(num);
619                                 *nleaf = *leaf;
620
621                                 setEdgeOffset(nleaf, alpha, count);
622
623                                 if (num > 1) {
624                                         float *pts = leaf->edge_intersections;
625                                         float *npts = nleaf->edge_intersections;
626                                         for (int i = 0; i < count; i++) {
627                                                 for (int j = 0; j < EDGE_FLOATS; j++) {
628                                                         npts[i * EDGE_FLOATS + j] = pts[i * EDGE_FLOATS + j];
629                                                 }
630                                         }
631                                         for (int i = count + 1; i < num; i++) {
632                                                 for (int j = 0; j < EDGE_FLOATS; j++) {
633                                                         npts[i * EDGE_FLOATS + j] = pts[(i - 1) * EDGE_FLOATS + j];
634                                                 }
635                                         }
636                                 }
637
638
639                                 removeLeaf(num - 1, (LeafNode *)leaf);
640                                 leaf = nleaf;
641                         }
642                 }
643
644                 return leaf;
645         }
646
647         /// Update parent link
648         void updateParent(InternalNode *node, int len, int st[3], LeafNode *leaf)
649         {
650                 // First, locate the parent
651                 int count;
652                 InternalNode *parent = locateParent(node, len, st, count);
653
654                 // Update
655                 setChild(parent, count, (Node *)leaf);
656         }
657
658         void updateParent(InternalNode *node, int len, int st[3])
659         {
660                 if (len == dimen) {
661                         root = (Node *)node;
662                         return;
663                 }
664
665                 // First, locate the parent
666                 int count;
667                 InternalNode *parent = locateParent(len, st, count);
668
669                 // UPdate
670                 setChild(parent, count, (Node *)node);
671         }
672
673         /// Find edge intersection on a given edge
674         int getEdgeIntersectionByIndex(int st[3], int index, float pt[3], int check)
675         {
676                 // First, locat the leaf
677                 LeafNode *leaf;
678                 if (check) {
679                         leaf = locateLeafCheck(st);
680                 }
681                 else {
682                         leaf = locateLeaf(st);
683                 }
684
685                 if (leaf && getStoredEdgesParity(leaf, index)) {
686                         float off = getEdgeOffset(leaf, getEdgeCount(leaf, index));
687                         pt[0] = (float) st[0];
688                         pt[1] = (float) st[1];
689                         pt[2] = (float) st[2];
690                         pt[index] += off * mindimen;
691
692                         return 1;
693                 }
694                 else {
695                         return 0;
696                 }
697         }
698
699         /// Retrieve number of edges intersected
700         int getPrimalEdgesMask(LeafNode *leaf)
701         {
702                 return leaf->primary_edge_intersections;
703         }
704
705         int getPrimalEdgesMask2(LeafNode *leaf)
706         {
707                 return (((leaf->edge_parity &   0x1) >> 0) |
708                                 ((leaf->edge_parity &  0x10) >> 3) |
709                                 ((leaf->edge_parity & 0x100) >> 6));
710         }
711
712         /// Get the count for a primary edge
713         int getEdgeCount(LeafNode *leaf, int index)
714         {
715                 return edgeCountTable[getPrimalEdgesMask(leaf)][index];
716         }
717         int getNumEdges(LeafNode *leaf)
718         {
719                 return numEdgeTable[getPrimalEdgesMask(leaf)];
720         }
721
722         int getNumEdges2(LeafNode *leaf)
723         {
724                 return numEdgeTable[getPrimalEdgesMask2(leaf)];
725         }
726
727         /// Set edge intersection
728         void setEdgeOffset(LeafNode *leaf, float pt, int count)
729         {
730                 float *pts = leaf->edge_intersections;
731                 pts[EDGE_FLOATS * count] = pt;
732                 pts[EDGE_FLOATS * count + 1] = 0;
733                 pts[EDGE_FLOATS * count + 2] = 0;
734                 pts[EDGE_FLOATS * count + 3] = 0;
735         }
736
737         /// Set multiple edge intersections
738         void setEdgeOffsets(LeafNode *leaf, float pt[3], int len)
739         {
740                 float *pts = leaf->edge_intersections;
741                 for (int i = 0; i < len; i++) {
742                         pts[i] = pt[i];
743                 }
744         }
745
746         /// Retrieve edge intersection
747         float getEdgeOffset(LeafNode *leaf, int count)
748         {
749                 return leaf->edge_intersections[4 * count];
750         }
751
752         /// Update method
753         LeafNode *updateEdgeOffsets(LeafNode *leaf, int oldlen, int newlen, float offs[3])
754         {
755                 // First, create a new leaf node
756                 LeafNode *nleaf = createLeaf(newlen);
757                 *nleaf = *leaf;
758
759                 // Next, fill in the offsets
760                 setEdgeOffsets(nleaf, offs, newlen);
761
762                 // Finally, delete the old leaf
763                 removeLeaf(oldlen, leaf);
764
765                 return nleaf;
766         }
767
768         /// Set minimizer index
769         void setMinimizerIndex(LeafNode *leaf, int index)
770         {
771                 leaf->minimizer_index = index;
772         }
773
774         /// Get minimizer index
775         int getMinimizerIndex(LeafNode *leaf)
776         {
777                 return leaf->minimizer_index;
778         }
779
780         int getMinimizerIndex(LeafNode *leaf, int eind)
781         {
782                 int add = manifold_table[getSignMask(leaf)].pairs[eind][0] - 1;
783                 assert(add >= 0);
784                 return leaf->minimizer_index + add;
785         }
786
787         void getMinimizerIndices(LeafNode *leaf, int eind, int inds[2])
788         {
789                 const int *add = manifold_table[getSignMask(leaf)].pairs[eind];
790                 inds[0] = leaf->minimizer_index + add[0] - 1;
791                 if (add[0] == add[1]) {
792                         inds[1] = -1;
793                 }
794                 else {
795                         inds[1] = leaf->minimizer_index + add[1] - 1;
796                 }
797         }
798
799
800         /// Set edge intersection
801         void setEdgeOffsetNormal(LeafNode *leaf, float pt, float a, float b, float c, int count)
802         {
803                 float *pts = leaf->edge_intersections;
804                 pts[4 * count] = pt;
805                 pts[4 * count + 1] = a;
806                 pts[4 * count + 2] = b;
807                 pts[4 * count + 3] = c;
808         }
809
810         float getEdgeOffsetNormal(LeafNode *leaf, int count, float& a, float& b, float& c)
811         {
812                 float *pts = leaf->edge_intersections;
813                 a = pts[4 * count + 1];
814                 b = pts[4 * count + 2];
815                 c = pts[4 * count + 3];
816                 return pts[4 * count];
817         }
818
819         /// Set multiple edge intersections
820         void setEdgeOffsetsNormals(LeafNode *leaf, const float pt[],
821                                                            const float a[], const float b[],
822                                                            const float c[], int len)
823         {
824                 float *pts = leaf->edge_intersections;
825                 for (int i = 0; i < len; i++) {
826                         if (pt[i] > 1 || pt[i] < 0) {
827                                 printf("\noffset: %f\n", pt[i]);
828                         }
829                         pts[i * 4] = pt[i];
830                         pts[i * 4 + 1] = a[i];
831                         pts[i * 4 + 2] = b[i];
832                         pts[i * 4 + 3] = c[i];
833                 }
834         }
835
836         /// Retrieve complete edge intersection
837         void getEdgeIntersectionByIndex(LeafNode *leaf, int index, int st[3], int len, float pt[3], float nm[3])
838         {
839                 int count = getEdgeCount(leaf, index);
840                 float *pts = leaf->edge_intersections;
841
842                 float off = pts[4 * count];
843
844                 pt[0] = (float) st[0];
845                 pt[1] = (float) st[1];
846                 pt[2] = (float) st[2];
847                 pt[index] += (off * len);
848
849                 nm[0] = pts[4 * count + 1];
850                 nm[1] = pts[4 * count + 2];
851                 nm[2] = pts[4 * count + 3];
852         }
853
854         float getEdgeOffsetNormalByIndex(LeafNode *leaf, int index, float nm[3])
855         {
856                 int count = getEdgeCount(leaf, index);
857                 float *pts = leaf->edge_intersections;
858
859                 float off = pts[4 * count];
860
861                 nm[0] = pts[4 * count + 1];
862                 nm[1] = pts[4 * count + 2];
863                 nm[2] = pts[4 * count + 3];
864
865                 return off;
866         }
867
868         void fillEdgeIntersections(LeafNode *leaf, int st[3], int len, float pts[12][3], float norms[12][3])
869         {
870                 int i;
871                 // int stt[3] = {0, 0, 0};
872
873                 // The three primal edges are easy
874                 int pmask[3] = {0, 4, 8};
875                 for (i = 0; i < 3; i++) {
876                         if (getEdgeParity(leaf, pmask[i])) {
877                                 // getEdgeIntersectionByIndex(leaf, i, stt, 1, pts[pmask[i]], norms[pmask[i]]);
878                                 getEdgeIntersectionByIndex(leaf, i, st, len, pts[pmask[i]], norms[pmask[i]]);
879                         }
880                 }
881
882                 // 3 face adjacent cubes
883                 int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}};
884                 int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}};
885                 for (i = 0; i < 3; i++) {
886                         int e1 = getEdgeParity(leaf, fmask[i][0]);
887                         int e2 = getEdgeParity(leaf, fmask[i][1]);
888                         if (e1 || e2) {
889                                 int nst[3] = {st[0], st[1], st[2]};
890                                 nst[i] += len;
891                                 // int nstt[3] = {0, 0, 0};
892                                 // nstt[i] += 1;
893                                 LeafNode *node = locateLeaf(nst);
894
895                                 if (e1) {
896                                         // getEdgeIntersectionByIndex(node, femask[i][0], nstt, 1, pts[fmask[i][0]], norms[fmask[i][0]]);
897                                         getEdgeIntersectionByIndex(node, femask[i][0], nst, len, pts[fmask[i][0]], norms[fmask[i][0]]);
898                                 }
899                                 if (e2) {
900                                         // getEdgeIntersectionByIndex(node, femask[i][1], nstt, 1, pts[fmask[i][1]], norms[fmask[i][1]]);
901                                         getEdgeIntersectionByIndex(node, femask[i][1], nst, len, pts[fmask[i][1]], norms[fmask[i][1]]);
902                                 }
903                         }
904                 }
905
906                 // 3 edge adjacent cubes
907                 int emask[3] = {3, 7, 11};
908                 int eemask[3] = {0, 1, 2};
909                 for (i = 0; i < 3; i++) {
910                         if (getEdgeParity(leaf, emask[i])) {
911                                 int nst[3] = {st[0] + len, st[1] + len, st[2] + len};
912                                 nst[i] -= len;
913                                 // int nstt[3] = {1, 1, 1};
914                                 // nstt[i] -= 1;
915                                 LeafNode *node = locateLeaf(nst);
916
917                                 // getEdgeIntersectionByIndex(node, eemask[i], nstt, 1, pts[emask[i]], norms[emask[i]]);
918                                 getEdgeIntersectionByIndex(node, eemask[i], nst, len, pts[emask[i]], norms[emask[i]]);
919                         }
920                 }
921         }
922
923
924         void fillEdgeIntersections(LeafNode *leaf, int st[3], int len, float pts[12][3], float norms[12][3], int parity[12])
925         {
926                 int i;
927                 for (i = 0; i < 12; i++) {
928                         parity[i] = 0;
929                 }
930                 // int stt[3] = {0, 0, 0};
931
932                 // The three primal edges are easy
933                 int pmask[3] = {0, 4, 8};
934                 for (i = 0; i < 3; i++) {
935                         if (getStoredEdgesParity(leaf, i)) {
936                                 // getEdgeIntersectionByIndex(leaf, i, stt, 1, pts[pmask[i]], norms[pmask[i]]);
937                                 getEdgeIntersectionByIndex(leaf, i, st, len, pts[pmask[i]], norms[pmask[i]]);
938                                 parity[pmask[i]] = 1;
939                         }
940                 }
941
942                 // 3 face adjacent cubes
943                 int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}};
944                 int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}};
945                 for (i = 0; i < 3; i++) {
946                         {
947                                 int nst[3] = {st[0], st[1], st[2]};
948                                 nst[i] += len;
949                                 // int nstt[3] = {0, 0, 0};
950                                 // nstt[i] += 1;
951                                 LeafNode *node = locateLeafCheck(nst);
952                                 if (node == NULL) {
953                                         continue;
954                                 }
955
956                                 int e1 = getStoredEdgesParity(node, femask[i][0]);
957                                 int e2 = getStoredEdgesParity(node, femask[i][1]);
958
959                                 if (e1) {
960                                         // getEdgeIntersectionByIndex(node, femask[i][0], nstt, 1, pts[fmask[i][0]], norms[fmask[i][0]]);
961                                         getEdgeIntersectionByIndex(node, femask[i][0], nst, len, pts[fmask[i][0]], norms[fmask[i][0]]);
962                                         parity[fmask[i][0]] = 1;
963                                 }
964                                 if (e2) {
965                                         // getEdgeIntersectionByIndex(node, femask[i][1], nstt, 1, pts[fmask[i][1]], norms[fmask[i][1]]);
966                                         getEdgeIntersectionByIndex(node, femask[i][1], nst, len, pts[fmask[i][1]], norms[fmask[i][1]]);
967                                         parity[fmask[i][1]] = 1;
968                                 }
969                         }
970                 }
971
972                 // 3 edge adjacent cubes
973                 int emask[3] = {3, 7, 11};
974                 int eemask[3] = {0, 1, 2};
975                 for (i = 0; i < 3; i++) {
976                         //                      if(getEdgeParity(leaf, emask[i]))
977                         {
978                                 int nst[3] = {st[0] + len, st[1] + len, st[2] + len};
979                                 nst[i] -= len;
980                                 // int nstt[3] = {1, 1, 1};
981                                 // nstt[i] -= 1;
982                                 LeafNode *node = locateLeafCheck(nst);
983                                 if (node == NULL) {
984                                         continue;
985                                 }
986
987                                 if (getStoredEdgesParity(node, eemask[i])) {
988                                         // getEdgeIntersectionByIndex(node, eemask[i], nstt, 1, pts[emask[i]], norms[emask[i]]);
989                                         getEdgeIntersectionByIndex(node, eemask[i], nst, len, pts[emask[i]], norms[emask[i]]);
990                                         parity[emask[i]] = 1;
991                                 }
992                         }
993                 }
994         }
995
996         void fillEdgeOffsetsNormals(LeafNode *leaf, int st[3], int len, float pts[12], float norms[12][3], int parity[12])
997         {
998                 int i;
999                 for (i = 0; i < 12; i++) {
1000                         parity[i] = 0;
1001                 }
1002                 // int stt[3] = {0, 0, 0};
1003
1004                 // The three primal edges are easy
1005                 int pmask[3] = {0, 4, 8};
1006                 for (i = 0; i < 3; i++) {
1007                         if (getStoredEdgesParity(leaf, i)) {
1008                                 pts[pmask[i]] = getEdgeOffsetNormalByIndex(leaf, i, norms[pmask[i]]);
1009                                 parity[pmask[i]] = 1;
1010                         }
1011                 }
1012
1013                 // 3 face adjacent cubes
1014                 int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}};
1015                 int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}};
1016                 for (i = 0; i < 3; i++) {
1017                         {
1018                                 int nst[3] = {st[0], st[1], st[2]};
1019                                 nst[i] += len;
1020                                 // int nstt[3] = {0, 0, 0};
1021                                 // nstt[i] += 1;
1022                                 LeafNode *node = locateLeafCheck(nst);
1023                                 if (node == NULL) {
1024                                         continue;
1025                                 }
1026
1027                                 int e1 = getStoredEdgesParity(node, femask[i][0]);
1028                                 int e2 = getStoredEdgesParity(node, femask[i][1]);
1029
1030                                 if (e1) {
1031                                         pts[fmask[i][0]] = getEdgeOffsetNormalByIndex(node, femask[i][0], norms[fmask[i][0]]);
1032                                         parity[fmask[i][0]] = 1;
1033                                 }
1034                                 if (e2) {
1035                                         pts[fmask[i][1]] = getEdgeOffsetNormalByIndex(node, femask[i][1], norms[fmask[i][1]]);
1036                                         parity[fmask[i][1]] = 1;
1037                                 }
1038                         }
1039                 }
1040
1041                 // 3 edge adjacent cubes
1042                 int emask[3] = {3, 7, 11};
1043                 int eemask[3] = {0, 1, 2};
1044                 for (i = 0; i < 3; i++) {
1045                         //                      if(getEdgeParity(leaf, emask[i]))
1046                         {
1047                                 int nst[3] = {st[0] + len, st[1] + len, st[2] + len};
1048                                 nst[i] -= len;
1049                                 // int nstt[3] = {1, 1, 1};
1050                                 // nstt[i] -= 1;
1051                                 LeafNode *node = locateLeafCheck(nst);
1052                                 if (node == NULL) {
1053                                         continue;
1054                                 }
1055
1056                                 if (getStoredEdgesParity(node, eemask[i])) {
1057                                         pts[emask[i]] = getEdgeOffsetNormalByIndex(node, eemask[i], norms[emask[i]]);
1058                                         parity[emask[i]] = 1;
1059                                 }
1060                         }
1061                 }
1062         }
1063
1064
1065         /// Update method
1066         LeafNode *updateEdgeOffsetsNormals(LeafNode *leaf, int oldlen, int newlen, float offs[3], float a[3], float b[3], float c[3])
1067         {
1068                 // First, create a new leaf node
1069                 LeafNode *nleaf = createLeaf(newlen);
1070                 *nleaf = *leaf;
1071
1072                 // Next, fill in the offsets
1073                 setEdgeOffsetsNormals(nleaf, offs, a, b, c, newlen);
1074
1075                 // Finally, delete the old leaf
1076                 removeLeaf(oldlen, leaf);
1077
1078                 return nleaf;
1079         }
1080
1081         /// Locate a leaf
1082         /// WARNING: assuming this leaf already exists!
1083
1084         LeafNode *locateLeaf(int st[3])
1085         {
1086                 Node *node = (Node *)root;
1087                 for (int i = GRID_DIMENSION - 1; i > GRID_DIMENSION - maxDepth - 1; i--) {
1088                         int index = (((st[0] >> i) & 1) << 2) |
1089                                 (((st[1] >> i) & 1) << 1) |
1090                                 (((st[2] >> i) & 1));
1091                         node = getChild(&node->internal, getChildCount(&node->internal, index));
1092                 }
1093
1094                 return &node->leaf;
1095         }
1096
1097         LeafNode *locateLeaf(InternalNode *parent, int len, int st[3])
1098         {
1099                 Node *node = (Node *)parent;
1100                 int index;
1101                 for (int i = len / 2; i >= mindimen; i >>= 1) {
1102                         index = (((st[0] & i) ? 4 : 0) |
1103                                          ((st[1] & i) ? 2 : 0) |
1104                                          ((st[2] & i) ? 1 : 0));
1105                         node = getChild(&node->internal,
1106                                                         getChildCount(&node->internal, index));
1107                 }
1108
1109                 return &node->leaf;
1110         }
1111
1112         LeafNode *locateLeafCheck(int st[3])
1113         {
1114                 Node *node = (Node *)root;
1115                 for (int i = GRID_DIMENSION - 1; i > GRID_DIMENSION - maxDepth - 1; i--) {
1116                         int index = (((st[0] >> i) & 1) << 2) |
1117                                 (((st[1] >> i) & 1) << 1) |
1118                                 (((st[2] >> i) & 1));
1119                         if (!hasChild(&node->internal, index)) {
1120                                 return NULL;
1121                         }
1122                         node = getChild(&node->internal, getChildCount(&node->internal, index));
1123                 }
1124
1125                 return &node->leaf;
1126         }
1127
1128         InternalNode *locateParent(int len, int st[3], int& count)
1129         {
1130                 InternalNode *node = (InternalNode *)root;
1131                 InternalNode *pre = NULL;
1132                 int index = 0;
1133                 for (int i = dimen / 2; i >= len; i >>= 1) {
1134                         index = (((st[0] & i) ? 4 : 0) |
1135                                          ((st[1] & i) ? 2 : 0) |
1136                                          ((st[2] & i) ? 1 : 0));
1137                         pre = node;
1138                         node = &getChild(node, getChildCount(node, index))->internal;
1139                 }
1140
1141                 count = getChildCount(pre, index);
1142                 return pre;
1143         }
1144
1145         InternalNode *locateParent(InternalNode *parent, int len, int st[3], int& count)
1146         {
1147                 InternalNode *node = parent;
1148                 InternalNode *pre = NULL;
1149                 int index = 0;
1150                 for (int i = len / 2; i >= mindimen; i >>= 1) {
1151                         index = (((st[0] & i) ? 4 : 0) |
1152                                          ((st[1] & i) ? 2 : 0) |
1153                                          ((st[2] & i) ? 1 : 0));
1154                         pre = node;
1155                         node = (InternalNode *)getChild(node, getChildCount(node, index));
1156                 }
1157
1158                 count = getChildCount(pre, index);
1159                 return pre;
1160         }
1161
1162         /************ Operators for internal nodes ************/
1163
1164         /// If child index exists
1165         int hasChild(InternalNode *node, int index)
1166         {
1167                 return (node->has_child >> index) & 1;
1168         }
1169
1170         /// Get the pointer to child index
1171         Node *getChild(InternalNode *node, int count)
1172         {
1173                 return node->children[count];
1174         };
1175
1176         /// Get total number of children
1177         int getNumChildren(InternalNode *node)
1178         {
1179                 return numChildrenTable[node->has_child];
1180         }
1181
1182         /// Get the count of children
1183         int getChildCount(InternalNode *node, int index)
1184         {
1185                 return childrenCountTable[node->has_child][index];
1186         }
1187         int getChildIndex(InternalNode *node, int count)
1188         {
1189                 return childrenIndexTable[node->has_child][count];
1190         }
1191         int *getChildCounts(InternalNode *node)
1192         {
1193                 return childrenCountTable[node->has_child];
1194         }
1195
1196         /// Get all children
1197         void fillChildren(InternalNode *node, Node *children[8], int leaf[8])
1198         {
1199                 int count = 0;
1200                 for (int i = 0; i < 8; i++) {
1201                         leaf[i] = node->is_child_leaf(i);
1202                         if (hasChild(node, i)) {
1203                                 children[i] = getChild(node, count);
1204                                 count++;
1205                         }
1206                         else {
1207                                 children[i] = NULL;
1208                                 leaf[i] = 0;
1209                         }
1210                 }
1211         }
1212
1213         /// Sets the child pointer
1214         void setChild(InternalNode *node, int count, Node *chd)
1215         {
1216                 node->children[count] = chd;
1217         }
1218         void setInternalChild(InternalNode *node, int index, int count, InternalNode *chd)
1219         {
1220                 setChild(node, count, (Node *)chd);
1221                 node->has_child |= (1 << index);
1222         }
1223         void setLeafChild(InternalNode *node, int index, int count, LeafNode *chd)
1224         {
1225                 setChild(node, count, (Node *)chd);
1226                 node->has_child |= (1 << index);
1227                 node->child_is_leaf |= (1 << index);
1228         }
1229
1230         /// Add a kid to an existing internal node
1231         /// Fix me: can we do this without wasting memory ?
1232         /// Fixed: using variable memory
1233         InternalNode *addChild(InternalNode *node, int index, Node *child, int aLeaf)
1234         {
1235                 // Create new internal node
1236                 int num = getNumChildren(node);
1237                 InternalNode *rnode = createInternal(num + 1);
1238
1239                 // Establish children
1240                 int i;
1241                 int count1 = 0, count2 = 0;
1242                 for (i = 0; i < 8; i++) {
1243                         if (i == index) {
1244                                 if (aLeaf) {
1245                                         setLeafChild(rnode, i, count2, &child->leaf);
1246                                 }
1247                                 else {
1248                                         setInternalChild(rnode, i, count2, &child->internal);
1249                                 }
1250                                 count2++;
1251                         }
1252                         else if (hasChild(node, i)) {
1253                                 if (node->is_child_leaf(i)) {
1254                                         setLeafChild(rnode, i, count2, &getChild(node, count1)->leaf);
1255                                 }
1256                                 else {
1257                                         setInternalChild(rnode, i, count2, &getChild(node, count1)->internal);
1258                                 }
1259                                 count1++;
1260                                 count2++;
1261                         }
1262                 }
1263
1264                 removeInternal(num, node);
1265                 return rnode;
1266         }
1267
1268         /// Allocate a node
1269         InternalNode *createInternal(int length)
1270         {
1271                 InternalNode *inode = (InternalNode *)alloc[length]->allocate();
1272                 inode->has_child = 0;
1273                 inode->child_is_leaf = 0;
1274                 return inode;
1275         }
1276
1277         LeafNode *createLeaf(int length)
1278         {
1279                 assert(length <= 3);
1280
1281                 LeafNode *lnode = (LeafNode *)leafalloc[length]->allocate();
1282                 lnode->edge_parity = 0;
1283                 lnode->primary_edge_intersections = 0;
1284                 lnode->signs = 0;
1285
1286                 return lnode;
1287         }
1288
1289         void removeInternal(int num, InternalNode *node)
1290         {
1291                 alloc[num]->deallocate(node);
1292         }
1293
1294         void removeLeaf(int num, LeafNode *leaf)
1295         {
1296                 assert(num >= 0 && num <= 3);
1297                 leafalloc[num]->deallocate(leaf);
1298         }
1299
1300         /// Add a leaf (by creating a new par node with the leaf added)
1301         InternalNode *addLeafChild(InternalNode *par, int index, int count,
1302                                                            LeafNode *leaf)
1303         {
1304                 int num = getNumChildren(par) + 1;
1305                 InternalNode *npar = createInternal(num);
1306                 *npar = *par;
1307
1308                 if (num == 1) {
1309                         setLeafChild(npar, index, 0, leaf);
1310                 }
1311                 else {
1312                         int i;
1313                         for (i = 0; i < count; i++) {
1314                                 setChild(npar, i, getChild(par, i));
1315                         }
1316                         setLeafChild(npar, index, count, leaf);
1317                         for (i = count + 1; i < num; i++) {
1318                                 setChild(npar, i, getChild(par, i - 1));
1319                         }
1320                 }
1321
1322                 removeInternal(num - 1, par);
1323                 return npar;
1324         }
1325
1326         InternalNode *addInternalChild(InternalNode *par, int index, int count,
1327                                                                    InternalNode *node)
1328         {
1329                 int num = getNumChildren(par) + 1;
1330                 InternalNode *npar = createInternal(num);
1331                 *npar = *par;
1332
1333                 if (num == 1) {
1334                         setInternalChild(npar, index, 0, node);
1335                 }
1336                 else {
1337                         int i;
1338                         for (i = 0; i < count; i++) {
1339                                 setChild(npar, i, getChild(par, i));
1340                         }
1341                         setInternalChild(npar, index, count, node);
1342                         for (i = count + 1; i < num; i++) {
1343                                 setChild(npar, i, getChild(par, i - 1));
1344                         }
1345                 }
1346
1347                 removeInternal(num - 1, par);
1348                 return npar;
1349         }
1350 };
1351
1352 #endif