octree.h

Go to the documentation of this file.

/* SPDX-FileCopyrightText: 2011-2023 Blender Authors
 *
 * SPDX-License-Identifier: GPL-2.0-or-later */
 
#ifndef __OCTREE_H__
#define __OCTREE_H__
 
#include "GeoCommon.h"
#include "MemoryAllocator.h"
#include "ModelReader.h"
#include "Projections.h"
#include "Queue.h"
#include "cubes.h"
#include "dualcon.h"
#include "manifold_table.h"
#include <cassert>
#include <cstring>
#include <math.h>
#include <stdio.h>
 
/* Global defines */
/* Uncomment to see debug information */
// #define IN_DEBUG_MODE
 
/* Uncomment to see more output messages during repair */
// #define IN_VERBOSE_MODE
 
/* Set scan convert params */
 
#define EDGE_FLOATS 4
 
union Node;
struct LeafNode;
 
struct InternalNode {
  /* Initialized in Octree::BuildTable */
  static int numChildrenTable[256];
  static int childrenCountTable[256][8];
  static int childrenIndexTable[256][8];
 
  /* Bit N indicates whether child N exists or not */
  unsigned char has_child_bitfield;
  /* Bit N indicates whether child N is a leaf or not */
  unsigned char child_is_leaf_bitfield;
 
  /* Can have up to eight children */
  Node *children[0];
 
  int is_child_leaf(int index) const
  {
    return (child_is_leaf_bitfield >> index) & 1;
  }
 
  int has_child(int index) const
  {
    return (has_child_bitfield >> index) & 1;
  }
 
  Node *get_child(int count)
  {
    return children[count];
  }
 
  const Node *get_child(int count) const
  {
    return children[count];
  }
 
  int get_num_children() const
  {
    return numChildrenTable[has_child_bitfield];
  }
 
  int get_child_count(int index) const
  {
    return childrenCountTable[has_child_bitfield][index];
  }
  int get_child_index(int count)
  {
    return childrenIndexTable[has_child_bitfield][count];
  }
  const int *get_child_counts() const
  {
    return childrenCountTable[has_child_bitfield];
  }
 
  void fill_children(Node *children[8], int leaf[8])
  {
    int count = 0;
    for (int i = 0; i < 8; i++) {
      leaf[i] = is_child_leaf(i);
      if (has_child(i)) {
        children[i] = get_child(count);
        count++;
      }
      else {
        children[i] = NULL;
        leaf[i] = 0;
      }
    }
  }
 
  void set_child(int count, Node *chd)
  {
    children[count] = chd;
  }
  void set_internal_child(int index, int count, InternalNode *chd)
  {
    set_child(count, (Node *)chd);
    has_child_bitfield |= (1 << index);
  }
  void set_leaf_child(int index, int count, LeafNode *chd)
  {
    set_child(count, (Node *)chd);
    has_child_bitfield |= (1 << index);
    child_is_leaf_bitfield |= (1 << index);
  }
};
 
struct LeafNode /* TODO: remove this attribute once everything is fixed */ {
  unsigned short edge_parity : 12;
  unsigned short primary_edge_intersections : 3;
 
  /* XXX: maybe actually unused? */
  unsigned short in_process : 1;
 
  /* bitfield */
  char signs;
 
  int minimizer_index;
 
  unsigned short flood_fill;
 
  float edge_intersections[0];
};
 
/* Doesn't get directly allocated anywhere, just used for passing
   pointers to nodes that could be internal or leaf. */
union Node {
  InternalNode internal;
  LeafNode leaf;
};
 
/* Global variables */
extern const int edgemask[3];
extern const int faceMap[6][4];
extern const int cellProcFaceMask[12][3];
extern const int cellProcEdgeMask[6][5];
extern const int faceProcFaceMask[3][4][3];
extern const int edgeProcEdgeMask[3][2][5];
extern const int faceProcEdgeMask[3][4][6];
extern const int processEdgeMask[3][4];
extern const int dirCell[3][4][3];
extern const int dirEdge[3][4];
 
struct PathElement {
  /* Origin */
  int pos[3];
 
  /* link */
  PathElement *next;
};
 
struct PathList {
  /* Head */
  PathElement *head;
  PathElement *tail;
 
  /* Length of the list */
  int length;
 
  /* Next list */
  PathList *next;
};
 
class Octree {
 public:
  /* Public members */
 
  VirtualMemoryAllocator *alloc[9];
  VirtualMemoryAllocator *leafalloc[4];
 
  Node *root;
 
  ModelReader *reader;
 
  Cubes *cubes;
 
  int dimen;
  int mindimen, minshift;
 
  int maxDepth;
 
  float origin[3];
  float range;
 
  int nodeCount;
  int nodeSpace;
  int nodeCounts[9];
 
  int actualQuads, actualVerts;
 
  PathList *ringList;
 
  int maxTrianglePerCell;
  int outType; /* 0 for OFF, 1 for PLY, 2 for VOL */
 
  /* For flood filling */
  int use_flood_fill;
  float thresh;
 
  int use_manifold;
 
  float hermite_num;
 
  DualConMode mode;
 
 public:
  Octree(ModelReader *mr,
         DualConAllocOutput alloc_output_func,
         DualConAddVert add_vert_func,
         DualConAddQuad add_quad_func,
         DualConFlags flags,
         DualConMode mode,
         int depth,
         float threshold,
         float hermite_num);
 
  ~Octree();
 
  void scanConvert();
 
  void *getOutputMesh()
  {
    return output_mesh;
  }
 
 private:
  /* Helper functions */
 
  void initMemory();
 
  void freeMemory();
 
  void printMemUsage();
 
  void resetMinimalEdges();
 
  void cellProcParity(Node *node, int leaf, int depth);
  void faceProcParity(Node *node[2], int leaf[2], int depth[2], int maxdep, int dir);
  void edgeProcParity(Node *node[4], int leaf[4], int depth[4], int maxdep, int dir);
 
  void processEdgeParity(LeafNode *node[4], int depths[4], int maxdep, int dir);
 
  void addAllTriangles();
  void addTriangle(Triangle *trian, int triind);
  InternalNode *addTriangle(InternalNode *node, CubeTriangleIsect *p, int height);
 
  LeafNode *updateCell(LeafNode *node, CubeTriangleIsect *p);
 
  /* Routines to detect and patch holes */
  int numRings;
  int totRingLengths;
  int maxRingLength;
 
  void trace();
  Node *trace(Node *node, int *st, int len, int depth, PathList *&paths);
  void findPaths(
      Node *node[2], int leaf[2], int depth[2], int *st[2], int maxdep, int dir, PathList *&paths);
  void combinePaths(PathList *&list1, PathList *list2, PathList *paths, PathList *&rings);
  PathList *combineSinglePath(PathList *&head1,
                              PathList *pre1,
                              PathList *&list1,
                              PathList *&head2,
                              PathList *pre2,
                              PathList *&list2);
 
  Node *patch(Node *node, int st[3], int len, PathList *rings);
  Node *patchSplit(Node *node,
                   int st[3],
                   int len,
                   PathList *rings,
                   int dir,
                   PathList *&nrings1,
                   PathList *&nrings2);
  Node *patchSplitSingle(Node *node,
                         int st[3],
                         int len,
                         PathElement *head,
                         int dir,
                         PathList *&nrings1,
                         PathList *&nrings2);
  Node *connectFace(Node *node, int st[3], int len, int dir, PathElement *f1, PathElement *f2);
  Node *locateCell(InternalNode *node,
                   int st[3],
                   int len,
                   int ori[3],
                   int dir,
                   int side,
                   Node *&rleaf,
                   int rst[3],
                   int &rlen);
  void compressRing(PathElement *&ring);
  void getFacePoint(PathElement *leaf, int dir, int &x, int &y, float &p, float &q);
  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);
  int findPair(PathElement *head, int pos, int dir, PathElement *&pre1, PathElement *&pre2);
  int getSide(PathElement *e, int pos, int dir);
  int isEqual(PathElement *e1, PathElement *e2);
  void preparePrimalEdgesMask(InternalNode *node);
  void testFacePoint(PathElement *e1, PathElement *e2);
 
  void deletePath(PathList *&head, PathList *pre, PathList *&curr);
  void printPath(PathList *path);
  void printPath(PathElement *path);
  void printElement(PathElement *ele);
  void printPaths(PathList *path);
  void checkElement(PathElement *ele);
  void checkPath(PathElement *path);
 
  void buildSigns();
  void buildSigns(unsigned char table[], Node *node, int isLeaf, int sg, int rvalue[8]);
 
  /************************************************************************/
  /* To remove disconnected components */
  /************************************************************************/
  void floodFill();
  void clearProcessBits(Node *node, int height);
  int floodFill(LeafNode *leaf, int st[3], int len, int height, int threshold);
  int floodFill(Node *node, int st[3], int len, int height, int threshold);
 
  void writeOut();
 
  void countIntersection(Node *node, int height, int &nedge, int &ncell, int &nface);
  void generateMinimizer(Node *node, int st[3], int len, int height, int &offset);
  void computeMinimizer(const LeafNode *leaf, int st[3], int len, float rvalue[3]) const;
  void cellProcContour(Node *node, int leaf, int depth);
  void faceProcContour(Node *node[2], int leaf[2], int depth[2], int maxdep, int dir);
  void edgeProcContour(Node *node[4], int leaf[4], int depth[4], int maxdep, int dir);
  void processEdgeWrite(Node *node[4], int depths[4], int maxdep, int dir);
 
  /* output callbacks/data */
  DualConAllocOutput alloc_output;
  DualConAddVert add_vert;
  DualConAddQuad add_quad;
  void *output_mesh;
 
 private:
  /************ Operators for all nodes ************/
 
  int numEdgeTable[8];
  int edgeCountTable[8][3];
 
  void buildTable()
  {
    for (int i = 0; i < 256; i++) {
      InternalNode::numChildrenTable[i] = 0;
      int count = 0;
      for (int j = 0; j < 8; j++) {
        InternalNode::numChildrenTable[i] += ((i >> j) & 1);
        InternalNode::childrenCountTable[i][j] = count;
        InternalNode::childrenIndexTable[i][count] = j;
        count += ((i >> j) & 1);
      }
    }
 
    for (int i = 0; i < 8; i++) {
      numEdgeTable[i] = 0;
      int count = 0;
      for (int j = 0; j < 3; j++) {
        numEdgeTable[i] += ((i >> j) & 1);
        edgeCountTable[i][j] = count;
        count += ((i >> j) & 1);
      }
    }
  }
 
  int getSign(Node *node, int height, int index)
  {
    if (height == 0) {
      return getSign(&node->leaf, index);
    }
    else {
      if (node->internal.has_child(index)) {
        return getSign(
            node->internal.get_child(node->internal.get_child_count(index)), height - 1, index);
      }
      else {
        return getSign(
            node->internal.get_child(0), height - 1, 7 - node->internal.get_child_index(0));
      }
    }
  }
 
  /************ Operators for leaf nodes ************/
 
  void printInfo(int st[3])
  {
    printf("INFO AT: %d %d %d\n", st[0] >> minshift, st[1] >> minshift, st[2] >> minshift);
    LeafNode *leaf = (LeafNode *)locateLeafCheck(st);
    if (leaf) {
      printInfo(leaf);
    }
    else {
      printf("Leaf not exists!\n");
    }
  }
 
  void printInfo(const LeafNode *leaf)
  {
    /*
      printf("Edge mask: ");
      for(int i = 0; i < 12; i ++)
      {
      printf("%d ", getEdgeParity(leaf, i));
      }
      printf("\n");
      printf("Stored edge mask: ");
      for(i = 0; i < 3; i ++)
      {
      printf("%d ", getStoredEdgesParity(leaf, i));
      }
      printf("\n");
    */
    printf("Sign mask: ");
    for (int i = 0; i < 8; i++) {
      printf("%d ", getSign(leaf, i));
    }
    printf("\n");
  }
 
  int getSign(const LeafNode *leaf, int index)
  {
    return ((leaf->signs >> index) & 1);
  }
 
  void setSign(LeafNode *leaf, int index)
  {
    leaf->signs |= (1 << index);
  }
 
  void setSign(LeafNode *leaf, int index, int sign)
  {
    leaf->signs &= (~(1 << index));
    leaf->signs |= ((sign & 1) << index);
  }
 
  int getSignMask(const LeafNode *leaf) const
  {
    return leaf->signs;
  }
 
  void setInProcessAll(int st[3], int dir)
  {
    int nst[3], eind;
    for (int i = 0; i < 4; i++) {
      nst[0] = st[0] + dirCell[dir][i][0] * mindimen;
      nst[1] = st[1] + dirCell[dir][i][1] * mindimen;
      nst[2] = st[2] + dirCell[dir][i][2] * mindimen;
      eind = dirEdge[dir][i];
 
      LeafNode *cell = locateLeafCheck(nst);
      assert(cell);
 
      setInProcess(cell, eind);
    }
  }
 
  void flipParityAll(int st[3], int dir)
  {
    int nst[3], eind;
    for (int i = 0; i < 4; i++) {
      nst[0] = st[0] + dirCell[dir][i][0] * mindimen;
      nst[1] = st[1] + dirCell[dir][i][1] * mindimen;
      nst[2] = st[2] + dirCell[dir][i][2] * mindimen;
      eind = dirEdge[dir][i];
 
      LeafNode *cell = locateLeaf(nst);
      flipEdge(cell, eind);
    }
  }
 
  void setInProcess(LeafNode *leaf, int eind)
  {
    assert(eind >= 0 && eind <= 11);
 
    leaf->flood_fill |= (1 << eind);
  }
 
  void setOutProcess(LeafNode *leaf, int eind)
  {
    assert(eind >= 0 && eind <= 11);
 
    leaf->flood_fill &= ~(1 << eind);
  }
 
  int isInProcess(LeafNode *leaf, int eind)
  {
    assert(eind >= 0 && eind <= 11);
 
    return (leaf->flood_fill >> eind) & 1;
  }
 
  void generateSigns(LeafNode *leaf, unsigned char table[], int start)
  {
    leaf->signs = table[leaf->edge_parity];
 
    if ((start ^ leaf->signs) & 1) {
      leaf->signs = ~(leaf->signs);
    }
  }
 
  int getEdgeParity(const LeafNode *leaf, int index) const
  {
    assert(index >= 0 && index <= 11);
 
    return (leaf->edge_parity >> index) & 1;
  }
 
  int getFaceParity(LeafNode *leaf, int index)
  {
    int a = getEdgeParity(leaf, faceMap[index][0]) + getEdgeParity(leaf, faceMap[index][1]) +
            getEdgeParity(leaf, faceMap[index][2]) + getEdgeParity(leaf, faceMap[index][3]);
    return (a & 1);
  }
  int getFaceEdgeNum(LeafNode *leaf, int index)
  {
    int a = getEdgeParity(leaf, faceMap[index][0]) + getEdgeParity(leaf, faceMap[index][1]) +
            getEdgeParity(leaf, faceMap[index][2]) + getEdgeParity(leaf, faceMap[index][3]);
    return a;
  }
 
  void flipEdge(LeafNode *leaf, int index)
  {
    assert(index >= 0 && index <= 11);
 
    leaf->edge_parity ^= (1 << index);
  }
 
  void setEdge(LeafNode *leaf, int index)
  {
    assert(index >= 0 && index <= 11);
 
    leaf->edge_parity |= (1 << index);
  }
 
  void resetEdge(LeafNode *leaf, int index)
  {
    assert(index >= 0 && index <= 11);
 
    leaf->edge_parity &= ~(1 << index);
  }
 
  void createPrimalEdgesMask(LeafNode *leaf)
  {
    leaf->primary_edge_intersections = getPrimalEdgesMask2(leaf);
  }
 
  void setStoredEdgesParity(LeafNode *leaf, int pindex)
  {
    assert(pindex <= 2 && pindex >= 0);
 
    leaf->primary_edge_intersections |= (1 << pindex);
  }
 
  int getStoredEdgesParity(const LeafNode *leaf, int pindex) const
  {
    assert(pindex <= 2 && pindex >= 0);
 
    return (leaf->primary_edge_intersections >> pindex) & 1;
  }
 
  LeafNode *flipEdge(LeafNode *leaf, int index, float alpha)
  {
    flipEdge(leaf, index);
 
    if ((index & 3) == 0) {
      int ind = index / 4;
      if (getEdgeParity(leaf, index) && !getStoredEdgesParity(leaf, ind)) {
        /* Create a new node */
        int num = getNumEdges(leaf) + 1;
        setStoredEdgesParity(leaf, ind);
        int count = getEdgeCount(leaf, ind);
        LeafNode *nleaf = createLeaf(num);
        *nleaf = *leaf;
 
        setEdgeOffset(nleaf, alpha, count);
 
        if (num > 1) {
          float *pts = leaf->edge_intersections;
          float *npts = nleaf->edge_intersections;
          for (int i = 0; i < count; i++) {
            for (int j = 0; j < EDGE_FLOATS; j++) {
              npts[i * EDGE_FLOATS + j] = pts[i * EDGE_FLOATS + j];
            }
          }
          for (int i = count + 1; i < num; i++) {
            for (int j = 0; j < EDGE_FLOATS; j++) {
              npts[i * EDGE_FLOATS + j] = pts[(i - 1) * EDGE_FLOATS + j];
            }
          }
        }
 
        removeLeaf(num - 1, (LeafNode *)leaf);
        leaf = nleaf;
      }
    }
 
    return leaf;
  }
 
  void updateParent(InternalNode *node, int len, int st[3], LeafNode *leaf)
  {
    /* First, locate the parent */
    int count;
    InternalNode *parent = locateParent(node, len, st, count);
 
    /* Update */
    parent->set_child(count, (Node *)leaf);
  }
 
  void updateParent(InternalNode *node, int len, int st[3])
  {
    if (len == dimen) {
      root = (Node *)node;
      return;
    }
 
    /* First, locate the parent */
    int count;
    InternalNode *parent = locateParent(len, st, count);
 
    /* Update */
    parent->set_child(count, (Node *)node);
  }
 
  int getEdgeIntersectionByIndex(int st[3], int index, float pt[3], int check) const
  {
    /* First, locat the leaf */
    const LeafNode *leaf;
    if (check) {
      leaf = locateLeafCheck(st);
    }
    else {
      leaf = locateLeaf(st);
    }
 
    if (leaf && getStoredEdgesParity(leaf, index)) {
      float off = getEdgeOffset(leaf, getEdgeCount(leaf, index));
      pt[0] = (float)st[0];
      pt[1] = (float)st[1];
      pt[2] = (float)st[2];
      pt[index] += off * mindimen;
 
      return 1;
    }
    else {
      return 0;
    }
  }
 
  int getPrimalEdgesMask(const LeafNode *leaf) const
  {
    return leaf->primary_edge_intersections;
  }
 
  int getPrimalEdgesMask2(LeafNode *leaf)
  {
    return (((leaf->edge_parity & 0x1) >> 0) | ((leaf->edge_parity & 0x10) >> 3) |
            ((leaf->edge_parity & 0x100) >> 6));
  }
 
  int getEdgeCount(const LeafNode *leaf, int index) const
  {
    return edgeCountTable[getPrimalEdgesMask(leaf)][index];
  }
  int getNumEdges(LeafNode *leaf)
  {
    return numEdgeTable[getPrimalEdgesMask(leaf)];
  }
 
  int getNumEdges2(LeafNode *leaf)
  {
    return numEdgeTable[getPrimalEdgesMask2(leaf)];
  }
 
  void setEdgeOffset(LeafNode *leaf, float pt, int count)
  {
    float *pts = leaf->edge_intersections;
    pts[EDGE_FLOATS * count] = pt;
    pts[EDGE_FLOATS * count + 1] = 0;
    pts[EDGE_FLOATS * count + 2] = 0;
    pts[EDGE_FLOATS * count + 3] = 0;
  }
 
  void setEdgeOffsets(LeafNode *leaf, float pt[3], int len)
  {
    float *pts = leaf->edge_intersections;
    for (int i = 0; i < len; i++) {
      pts[i] = pt[i];
    }
  }
 
  float getEdgeOffset(const LeafNode *leaf, int count) const
  {
    return leaf->edge_intersections[4 * count];
  }
 
  LeafNode *updateEdgeOffsets(LeafNode *leaf, int oldlen, int newlen, float offs[3])
  {
    /* First, create a new leaf node */
    LeafNode *nleaf = createLeaf(newlen);
    *nleaf = *leaf;
 
    /* Next, fill in the offsets */
    setEdgeOffsets(nleaf, offs, newlen);
 
    /* Finally, delete the old leaf */
    removeLeaf(oldlen, leaf);
 
    return nleaf;
  }
 
  void setMinimizerIndex(LeafNode *leaf, int index)
  {
    leaf->minimizer_index = index;
  }
 
  int getMinimizerIndex(LeafNode *leaf)
  {
    return leaf->minimizer_index;
  }
 
  int getMinimizerIndex(LeafNode *leaf, int eind)
  {
    int add = manifold_table[getSignMask(leaf)].pairs[eind][0] - 1;
    assert(add >= 0);
    return leaf->minimizer_index + add;
  }
 
  void getMinimizerIndices(LeafNode *leaf, int eind, int inds[2])
  {
    const int *add = manifold_table[getSignMask(leaf)].pairs[eind];
    inds[0] = leaf->minimizer_index + add[0] - 1;
    if (add[0] == add[1]) {
      inds[1] = -1;
    }
    else {
      inds[1] = leaf->minimizer_index + add[1] - 1;
    }
  }
 
  void setEdgeOffsetNormal(LeafNode *leaf, float pt, float a, float b, float c, int count)
  {
    float *pts = leaf->edge_intersections;
    pts[4 * count] = pt;
    pts[4 * count + 1] = a;
    pts[4 * count + 2] = b;
    pts[4 * count + 3] = c;
  }
 
  float getEdgeOffsetNormal(LeafNode *leaf, int count, float &a, float &b, float &c)
  {
    float *pts = leaf->edge_intersections;
    a = pts[4 * count + 1];
    b = pts[4 * count + 2];
    c = pts[4 * count + 3];
    return pts[4 * count];
  }
 
  void setEdgeOffsetsNormals(
      LeafNode *leaf, const float pt[], const float a[], const float b[], const float c[], int len)
  {
    float *pts = leaf->edge_intersections;
    for (int i = 0; i < len; i++) {
      if (pt[i] > 1 || pt[i] < 0) {
        printf("\noffset: %f\n", pt[i]);
      }
      pts[i * 4] = pt[i];
      pts[i * 4 + 1] = a[i];
      pts[i * 4 + 2] = b[i];
      pts[i * 4 + 3] = c[i];
    }
  }
 
  void getEdgeIntersectionByIndex(
      const LeafNode *leaf, int index, int st[3], int len, float pt[3], float nm[3]) const
  {
    int count = getEdgeCount(leaf, index);
    const float *pts = leaf->edge_intersections;
 
    float off = pts[4 * count];
 
    pt[0] = (float)st[0];
    pt[1] = (float)st[1];
    pt[2] = (float)st[2];
    pt[index] += (off * len);
 
    nm[0] = pts[4 * count + 1];
    nm[1] = pts[4 * count + 2];
    nm[2] = pts[4 * count + 3];
  }
 
  float getEdgeOffsetNormalByIndex(LeafNode *leaf, int index, float nm[3])
  {
    int count = getEdgeCount(leaf, index);
    float *pts = leaf->edge_intersections;
 
    float off = pts[4 * count];
 
    nm[0] = pts[4 * count + 1];
    nm[1] = pts[4 * count + 2];
    nm[2] = pts[4 * count + 3];
 
    return off;
  }
 
  void fillEdgeIntersections(
      const LeafNode *leaf, int st[3], int len, float pts[12][3], float norms[12][3]) const
  {
    int i;
    // int stt[3] = {0, 0, 0};
 
    // The three primal edges are easy
    int pmask[3] = {0, 4, 8};
    for (i = 0; i < 3; i++) {
      if (getEdgeParity(leaf, pmask[i])) {
        // getEdgeIntersectionByIndex(leaf, i, stt, 1, pts[pmask[i]], norms[pmask[i]]);
        getEdgeIntersectionByIndex(leaf, i, st, len, pts[pmask[i]], norms[pmask[i]]);
      }
    }
 
    // 3 face adjacent cubes
    int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}};
    int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}};
    for (i = 0; i < 3; i++) {
      int e1 = getEdgeParity(leaf, fmask[i][0]);
      int e2 = getEdgeParity(leaf, fmask[i][1]);
      if (e1 || e2) {
        int nst[3] = {st[0], st[1], st[2]};
        nst[i] += len;
        // int nstt[3] = {0, 0, 0};
        // nstt[i] += 1;
        const LeafNode *node = locateLeaf(nst);
 
        if (e1) {
          // getEdgeIntersectionByIndex(node, femask[i][0], nstt, 1, pts[fmask[i][0]],
          // norms[fmask[i][0]]);
          getEdgeIntersectionByIndex(
              node, femask[i][0], nst, len, pts[fmask[i][0]], norms[fmask[i][0]]);
        }
        if (e2) {
          // getEdgeIntersectionByIndex(node, femask[i][1], nstt, 1, pts[fmask[i][1]],
          // norms[fmask[i][1]]);
          getEdgeIntersectionByIndex(
              node, femask[i][1], nst, len, pts[fmask[i][1]], norms[fmask[i][1]]);
        }
      }
    }
 
    // 3 edge adjacent cubes
    int emask[3] = {3, 7, 11};
    int eemask[3] = {0, 1, 2};
    for (i = 0; i < 3; i++) {
      if (getEdgeParity(leaf, emask[i])) {
        int nst[3] = {st[0] + len, st[1] + len, st[2] + len};
        nst[i] -= len;
        // int nstt[3] = {1, 1, 1};
        // nstt[i] -= 1;
        const LeafNode *node = locateLeaf(nst);
 
        // getEdgeIntersectionByIndex(node, eemask[i], nstt, 1, pts[emask[i]], norms[emask[i]]);
        getEdgeIntersectionByIndex(node, eemask[i], nst, len, pts[emask[i]], norms[emask[i]]);
      }
    }
  }
 
  void fillEdgeIntersections(const LeafNode *leaf,
                             int st[3],
                             int len,
                             float pts[12][3],
                             float norms[12][3],
                             int parity[12]) const
  {
    int i;
    for (i = 0; i < 12; i++) {
      parity[i] = 0;
    }
    // int stt[3] = {0, 0, 0};
 
    // The three primal edges are easy
    int pmask[3] = {0, 4, 8};
    for (i = 0; i < 3; i++) {
      if (getStoredEdgesParity(leaf, i)) {
        // getEdgeIntersectionByIndex(leaf, i, stt, 1, pts[pmask[i]], norms[pmask[i]]);
        getEdgeIntersectionByIndex(leaf, i, st, len, pts[pmask[i]], norms[pmask[i]]);
        parity[pmask[i]] = 1;
      }
    }
 
    // 3 face adjacent cubes
    int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}};
    int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}};
    for (i = 0; i < 3; i++) {
      {
        int nst[3] = {st[0], st[1], st[2]};
        nst[i] += len;
        // int nstt[3] = {0, 0, 0};
        // nstt[i] += 1;
        const LeafNode *node = locateLeafCheck(nst);
        if (node == NULL) {
          continue;
        }
 
        int e1 = getStoredEdgesParity(node, femask[i][0]);
        int e2 = getStoredEdgesParity(node, femask[i][1]);
 
        if (e1) {
          // getEdgeIntersectionByIndex(node, femask[i][0], nstt, 1, pts[fmask[i][0]],
          // norms[fmask[i][0]]);
          getEdgeIntersectionByIndex(
              node, femask[i][0], nst, len, pts[fmask[i][0]], norms[fmask[i][0]]);
          parity[fmask[i][0]] = 1;
        }
        if (e2) {
          // getEdgeIntersectionByIndex(node, femask[i][1], nstt, 1, pts[fmask[i][1]],
          // norms[fmask[i][1]]);
          getEdgeIntersectionByIndex(
              node, femask[i][1], nst, len, pts[fmask[i][1]], norms[fmask[i][1]]);
          parity[fmask[i][1]] = 1;
        }
      }
    }
 
    // 3 edge adjacent cubes
    int emask[3] = {3, 7, 11};
    int eemask[3] = {0, 1, 2};
    for (i = 0; i < 3; i++) {
      //          if(getEdgeParity(leaf, emask[i]))
      {
        int nst[3] = {st[0] + len, st[1] + len, st[2] + len};
        nst[i] -= len;
        // int nstt[3] = {1, 1, 1};
        // nstt[i] -= 1;
        const LeafNode *node = locateLeafCheck(nst);
        if (node == NULL) {
          continue;
        }
 
        if (getStoredEdgesParity(node, eemask[i])) {
          // getEdgeIntersectionByIndex(node, eemask[i], nstt, 1, pts[emask[i]], norms[emask[i]]);
          getEdgeIntersectionByIndex(node, eemask[i], nst, len, pts[emask[i]], norms[emask[i]]);
          parity[emask[i]] = 1;
        }
      }
    }
  }
 
  void fillEdgeOffsetsNormals(
      LeafNode *leaf, int st[3], int len, float pts[12], float norms[12][3], int parity[12])
  {
    int i;
    for (i = 0; i < 12; i++) {
      parity[i] = 0;
    }
    // int stt[3] = {0, 0, 0};
 
    // The three primal edges are easy
    int pmask[3] = {0, 4, 8};
    for (i = 0; i < 3; i++) {
      if (getStoredEdgesParity(leaf, i)) {
        pts[pmask[i]] = getEdgeOffsetNormalByIndex(leaf, i, norms[pmask[i]]);
        parity[pmask[i]] = 1;
      }
    }
 
    // 3 face adjacent cubes
    int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}};
    int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}};
    for (i = 0; i < 3; i++) {
      {
        int nst[3] = {st[0], st[1], st[2]};
        nst[i] += len;
        // int nstt[3] = {0, 0, 0};
        // nstt[i] += 1;
        LeafNode *node = locateLeafCheck(nst);
        if (node == NULL) {
          continue;
        }
 
        int e1 = getStoredEdgesParity(node, femask[i][0]);
        int e2 = getStoredEdgesParity(node, femask[i][1]);
 
        if (e1) {
          pts[fmask[i][0]] = getEdgeOffsetNormalByIndex(node, femask[i][0], norms[fmask[i][0]]);
          parity[fmask[i][0]] = 1;
        }
        if (e2) {
          pts[fmask[i][1]] = getEdgeOffsetNormalByIndex(node, femask[i][1], norms[fmask[i][1]]);
          parity[fmask[i][1]] = 1;
        }
      }
    }
 
    // 3 edge adjacent cubes
    int emask[3] = {3, 7, 11};
    int eemask[3] = {0, 1, 2};
    for (i = 0; i < 3; i++) {
      //          if(getEdgeParity(leaf, emask[i]))
      {
        int nst[3] = {st[0] + len, st[1] + len, st[2] + len};
        nst[i] -= len;
        // int nstt[3] = {1, 1, 1};
        // nstt[i] -= 1;
        LeafNode *node = locateLeafCheck(nst);
        if (node == NULL) {
          continue;
        }
 
        if (getStoredEdgesParity(node, eemask[i])) {
          pts[emask[i]] = getEdgeOffsetNormalByIndex(node, eemask[i], norms[emask[i]]);
          parity[emask[i]] = 1;
        }
      }
    }
  }
 
  LeafNode *updateEdgeOffsetsNormals(
      LeafNode *leaf, int oldlen, int newlen, float offs[3], float a[3], float b[3], float c[3])
  {
    /* First, create a new leaf node */
    LeafNode *nleaf = createLeaf(newlen);
    *nleaf = *leaf;
 
    /* Next, fill in the offsets */
    setEdgeOffsetsNormals(nleaf, offs, a, b, c, newlen);
 
    /* Finally, delete the old leaf */
    removeLeaf(oldlen, leaf);
 
    return nleaf;
  }
 
  LeafNode *locateLeaf(int st[3])
  {
    Node *node = (Node *)root;
    for (int i = GRID_DIMENSION - 1; i > GRID_DIMENSION - maxDepth - 1; i--) {
      int index = (((st[0] >> i) & 1) << 2) | (((st[1] >> i) & 1) << 1) | (((st[2] >> i) & 1));
      node = node->internal.get_child(node->internal.get_child_count(index));
    }
 
    return &node->leaf;
  }
 
  const LeafNode *locateLeaf(int st[3]) const
  {
    const Node *node = root;
    for (int i = GRID_DIMENSION - 1; i > GRID_DIMENSION - maxDepth - 1; i--) {
      int index = (((st[0] >> i) & 1) << 2) | (((st[1] >> i) & 1) << 1) | (((st[2] >> i) & 1));
      node = node->internal.get_child(node->internal.get_child_count(index));
    }
 
    return &node->leaf;
  }
 
  LeafNode *locateLeaf(InternalNode *parent, int len, int st[3])
  {
    Node *node = (Node *)parent;
    int index;
    for (int i = len / 2; i >= mindimen; i >>= 1) {
      index = (((st[0] & i) ? 4 : 0) | ((st[1] & i) ? 2 : 0) | ((st[2] & i) ? 1 : 0));
      node = node->internal.get_child(node->internal.get_child_count(index));
    }
 
    return &node->leaf;
  }
 
  LeafNode *locateLeafCheck(int st[3])
  {
    Node *node = (Node *)root;
    for (int i = GRID_DIMENSION - 1; i > GRID_DIMENSION - maxDepth - 1; i--) {
      int index = (((st[0] >> i) & 1) << 2) | (((st[1] >> i) & 1) << 1) | (((st[2] >> i) & 1));
      if (!node->internal.has_child(index)) {
        return NULL;
      }
      node = node->internal.get_child(node->internal.get_child_count(index));
    }
 
    return &node->leaf;
  }
 
  const LeafNode *locateLeafCheck(int st[3]) const
  {
    const Node *node = root;
    for (int i = GRID_DIMENSION - 1; i > GRID_DIMENSION - maxDepth - 1; i--) {
      int index = (((st[0] >> i) & 1) << 2) | (((st[1] >> i) & 1) << 1) | (((st[2] >> i) & 1));
      if (!node->internal.has_child(index)) {
        return NULL;
      }
      node = node->internal.get_child(node->internal.get_child_count(index));
    }
 
    return &node->leaf;
  }
 
  InternalNode *locateParent(int len, int st[3], int &count)
  {
    InternalNode *node = (InternalNode *)root;
    InternalNode *pre = NULL;
    int index = 0;
    for (int i = dimen / 2; i >= len; i >>= 1) {
      index = (((st[0] & i) ? 4 : 0) | ((st[1] & i) ? 2 : 0) | ((st[2] & i) ? 1 : 0));
      pre = node;
      node = &node->get_child(node->get_child_count(index))->internal;
    }
 
    count = pre->get_child_count(index);
    return pre;
  }
 
  InternalNode *locateParent(InternalNode *parent, int len, int st[3], int &count)
  {
    InternalNode *node = parent;
    InternalNode *pre = NULL;
    int index = 0;
    for (int i = len / 2; i >= mindimen; i >>= 1) {
      index = (((st[0] & i) ? 4 : 0) | ((st[1] & i) ? 2 : 0) | ((st[2] & i) ? 1 : 0));
      pre = node;
      node = &node->get_child(node->get_child_count(index))->internal;
    }
 
    count = pre->get_child_count(index);
    return pre;
  }
 
  /************ Operators for internal nodes ************/
 
  InternalNode *addChild(InternalNode *node, int index, Node *child, int aLeaf)
  {
    // Create new internal node
    int num = node->get_num_children();
    InternalNode *rnode = createInternal(num + 1);
 
    // Establish children
    int i;
    int count1 = 0, count2 = 0;
    for (i = 0; i < 8; i++) {
      if (i == index) {
        if (aLeaf) {
          rnode->set_leaf_child(i, count2, &child->leaf);
        }
        else {
          rnode->set_internal_child(i, count2, &child->internal);
        }
        count2++;
      }
      else if (node->has_child(i)) {
        if (node->is_child_leaf(i)) {
          rnode->set_leaf_child(i, count2, &node->get_child(count1)->leaf);
        }
        else {
          rnode->set_internal_child(i, count2, &node->get_child(count1)->internal);
        }
        count1++;
        count2++;
      }
    }
 
    removeInternal(num, node);
    return rnode;
  }
 
  InternalNode *createInternal(int length)
  {
    InternalNode *inode = (InternalNode *)alloc[length]->allocate();
    inode->has_child_bitfield = 0;
    inode->child_is_leaf_bitfield = 0;
    return inode;
  }
 
  LeafNode *createLeaf(int length)
  {
    assert(length <= 3);
 
    LeafNode *lnode = (LeafNode *)leafalloc[length]->allocate();
    lnode->edge_parity = 0;
    lnode->primary_edge_intersections = 0;
    lnode->signs = 0;
 
    return lnode;
  }
 
  void removeInternal(int num, InternalNode *node)
  {
    alloc[num]->deallocate(node);
  }
 
  void removeLeaf(int num, LeafNode *leaf)
  {
    assert(num >= 0 && num <= 3);
    leafalloc[num]->deallocate(leaf);
  }
 
  InternalNode *addLeafChild(InternalNode *par, int index, int count, LeafNode *leaf)
  {
    int num = par->get_num_children() + 1;
    InternalNode *npar = createInternal(num);
    *npar = *par;
 
    if (num == 1) {
      npar->set_leaf_child(index, 0, leaf);
    }
    else {
      int i;
      for (i = 0; i < count; i++) {
        npar->set_child(i, par->get_child(i));
      }
      npar->set_leaf_child(index, count, leaf);
      for (i = count + 1; i < num; i++) {
        npar->set_child(i, par->get_child(i - 1));
      }
    }
 
    removeInternal(num - 1, par);
    return npar;
  }
 
  InternalNode *addInternalChild(InternalNode *par, int index, int count, InternalNode *node)
  {
    int num = par->get_num_children() + 1;
    InternalNode *npar = createInternal(num);
    *npar = *par;
 
    if (num == 1) {
      npar->set_internal_child(index, 0, node);
    }
    else {
      int i;
      for (i = 0; i < count; i++) {
        npar->set_child(i, par->get_child(i));
      }
      npar->set_internal_child(index, count, node);
      for (i = count + 1; i < num; i++) {
        npar->set_child(i, par->get_child(i - 1));
      }
    }
 
    removeInternal(num - 1, par);
    return npar;
  }
 
#ifdef WITH_CXX_GUARDEDALLOC
  MEM_CXX_CLASS_ALLOC_FUNCS("DUALCON:Octree")
#endif
};
 
#endif /* __OCTREE_H__ */