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