octree.h

Go to the documentation of this file.

/*
 * ***** BEGIN GPL LICENSE BLOCK *****
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 *
 * Contributor(s): Tao Ju
 *
 * ***** END GPL LICENSE BLOCK *****
 */

#ifndef OCTREE_H
#define OCTREE_H

#include <stdio.h>
#include <math.h>
#include "GeoCommon.h"
#include "Projections.h"
#include "ModelReader.h"
#include "MemoryAllocator.h"
#include "cubes.h"
#include "Queue.h"
#include "manifold_table.h"
#include "dualcon.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 */
// Uncomment to use Dual Contouring on Hermit representation
// for better sharp feature reproduction, but more mem is required
// The number indicates how far do we allow the minimizer to shoot
// outside the cell
#define USE_HERMIT 1.0f

#ifdef USE_HERMIT
//#define CINDY
#endif


//#define TESTMANIFOLD


/* Set output options */
// Comment out to prevent writing output files
#define OUTPUT_REPAIRED


/* Set node bytes */
#ifdef USE_HERMIT
#define EDGE_BYTES 16
#define EDGE_FLOATS 4
#else
#define EDGE_BYTES 4
#define EDGE_FLOATS 1
#endif

#define CINDY_BYTES 0

/*#define LEAF_EXTRA_BYTES FLOOD_BYTES + CINDY_BYTES

#ifdef USE_HERMIT
#define LEAF_NODE_BYTES 7 + LEAF_EXTRA_BYTES
#else
#define LEAF_NODE_BYTES 3 + LEAF_EXTRA_BYTES
#endif*/

#define INTERNAL_NODE_BYTES 2
#define POINTER_BYTES 8
#define FLOOD_FILL_BYTES 2

#define signtype short
#define nodetype int
#define numtype int

/* 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 ] ;

    UCHAR* 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;

    // testing
    FILE* testfout ;

    float hermite_num;

    DualConMode mode;

    int leaf_node_bytes;
    int leaf_extra_bytes;
    int flood_bytes;

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 ( UCHAR* node, int leaf, int depth ) ;
    void faceProcParity ( UCHAR* node[2], int leaf[2], int depth[2], int maxdep, int dir ) ;
    void edgeProcParity ( UCHAR* node[4], int leaf[4], int depth[4], int maxdep, int dir ) ;

    void processEdgeParity ( UCHAR* node[4], int depths[4], int maxdep, int dir ) ;

    void addTrian ( );
    void addTrian ( Triangle* trian, int triind );
    UCHAR* addTrian ( UCHAR* node, Projections* p, int height );

    UCHAR* updateCell( UCHAR* node, Projections* p ) ;

    /* Routines to detect and patch holes */
    int numRings ;
    int totRingLengths ;
    int maxRingLength ;

    void trace ( ) ;
    UCHAR* trace ( UCHAR* node, int* st, int len, int depth, PathList*& paths ) ;
    void findPaths ( UCHAR* 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 ) ;
    
    UCHAR* patch ( UCHAR* node, int st[3], int len, PathList* rings ) ;
    UCHAR* patchSplit ( UCHAR* node, int st[3], int len, PathList* rings, int dir, PathList*& nrings1, PathList*& nrings2 ) ;
    UCHAR* patchSplitSingle ( UCHAR* node, int st[3], int len, PathElement* head, int dir, PathList*& nrings1, PathList*& nrings2 ) ;
    UCHAR* connectFace ( UCHAR* node, int st[3], int len, int dir, PathElement* f1, PathElement* f2 ) ;
    UCHAR* locateCell( UCHAR* node, int st[3], int len, int ori[3], int dir, int side, UCHAR*& rleaf, int rst[3], int& rlen ) ;
    void compressRing ( PathElement*& ring ) ;
    void getFacePoint( PathElement* leaf, int dir, int& x, int& y, float& p, float& q ) ;
    UCHAR* patchAdjacent( UCHAR* node, int len, int st1[3], UCHAR* leaf1, int st2[3], UCHAR* 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( UCHAR* 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[], UCHAR* node, int isLeaf, int sg, int rvalue[8] ) ;

    /************************************************************************/
    /* To remove disconnected components */
    /************************************************************************/
    void floodFill( ) ;
    void clearProcessBits( UCHAR* node, int height ) ;
    int floodFill( UCHAR* node, int st[3], int len, int height, int threshold ) ;

    void writeOut();
    void writeOFF ( char* fname ) ;
    void writePLY ( char* fname ) ;
    void writeOpenEdges( FILE* fout ) ;
    void writeAllEdges( FILE* fout, int mode ) ;
    void writeAllEdges( FILE* fout, UCHAR* node, int height, int st[3], int len, int mode ) ;
    
    void writeOctree( char* fname ) ;
    void writeOctree( FILE* fout, UCHAR* node, int depth ) ;
#ifdef USE_HERMIT
    void writeOctreeGeom( char* fname ) ;
    void writeOctreeGeom( FILE* fout, UCHAR* node, int st[3], int len, int depth ) ;
#endif
#ifdef USE_HERMIT
    void writeDCF ( char* fname ) ;
    void writeDCF ( FILE* fout, UCHAR* node, int height, int st[3], int len ) ;
    void countEdges ( UCHAR* node, int height, int& nedge, int mode ) ;
    void countIntersection( UCHAR* node, int height, int& nedge, int& ncell, int& nface ) ;
    void generateMinimizer( UCHAR* node, int st[3], int len, int height, int& offset ) ;
    void computeMinimizer( UCHAR* leaf, int st[3], int len, float rvalue[3] ) ;
    void cellProcContour ( UCHAR* node, int leaf, int depth ) ;
    void faceProcContour ( UCHAR* node[2], int leaf[2], int depth[2], int maxdep, int dir ) ;
    void edgeProcContour ( UCHAR* node[4], int leaf[4], int depth[4], int maxdep, int dir ) ;
    void processEdgeWrite ( UCHAR* node[4], int depths[4], int maxdep, int dir ) ;
#else
    void countIntersection( UCHAR* node, int height, int& nquad, int& nvert ) ;
    void writeVertex( UCHAR* node, int st[3], int len, int height, int& offset, FILE* fout ) ;
    void writeQuad( UCHAR* node, int st[3], int len, int height, FILE* fout ) ;
#endif

    void writeModel( char* fname ) ;

    /************************************************************************/
    /* Write out original vertex tags */
    /************************************************************************/
#ifdef CINDY
    void writeTags( char* fname ) ;
    void readVertices( ) ;
    void readVertex(  UCHAR* node, int st[3], int len, int height, float v[3], int index ) ;
    void outputTags( UCHAR* node, int height, FILE* fout ) ;
    void clearCindyBits( UCHAR* node, int height ) ;
#endif

    /* output callbacks/data */
    DualConAllocOutput alloc_output;
    DualConAddVert add_vert;
    DualConAddQuad add_quad;
    void *output_mesh;
    
private:
    /************ Operators for all nodes ************/


    int numChildrenTable[ 256 ] ;
    int childrenCountTable[ 256 ][ 8 ] ;
    int childrenIndexTable[ 256 ][ 8 ] ;
    int numEdgeTable[ 8 ] ;
    int edgeCountTable[ 8 ][ 3 ] ;

    void buildTable ( )
    {
        for ( int i = 0 ; i < 256 ; i ++ )
        {
            numChildrenTable[ i ] = 0 ;
            int count = 0 ;
            for ( int j = 0 ; j < 8 ; j ++ )
            {
                numChildrenTable[ i ] += ( ( i >> j ) & 1 ) ;
                childrenCountTable[ i ][ j ] = count ;
                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( UCHAR* node, int height, int index )
    {
        if ( height == 0 )
        {
            return getSign( node, index ) ;
        }
        else
        {
            if ( hasChild( node, index ) )
            {
                return getSign( getChild( node, getChildCount( node, index ) ), height - 1, index ) ;
            }
            else
            {
                return getSign( getChild( node, 0 ), height - 1, 7 - getChildIndex( node, 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 ) ;
        UCHAR* leaf = locateLeafCheck( st ) ;
        if ( leaf == NULL )
        {
            printf("Leaf not exists!\n") ;
        }
        else
        {
            printInfo( leaf ) ;
        }
    }

    void printInfo( UCHAR* 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 ( UCHAR* leaf, int index )
    {
        return (( leaf[2] >> index ) & 1 );     
    };

    void setSign ( UCHAR* leaf, int index ) 
    {
        leaf[2] |= ( 1 << index ) ;
    };

    void setSign ( UCHAR* leaf, int index, int sign ) 
    {
        leaf[2] &= ( ~ ( 1 << index ) ) ;
        leaf[2] |= ( ( sign & 1 ) << index ) ;
    };

    int getSignMask( UCHAR* leaf )
    {
        return leaf[2] ;
    }

    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] ;

            UCHAR* cell = locateLeafCheck( nst ) ;
            if ( cell == NULL )
            {
                printf("Wrong!\n") ;
            }
            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] ;

            UCHAR* cell = locateLeaf( nst ) ;
            flipEdge( cell, eind ) ;
        }
    }

    void setInProcess( UCHAR* leaf, int eind )
    {
        // leaf[1] |= ( 1 << 7 ) ;
        ( (USHORT*) (leaf + leaf_node_bytes - (flood_bytes + CINDY_BYTES)))[0] |= ( 1 << eind ) ;
    }
    void setOutProcess( UCHAR* leaf, int eind )
    {
        // leaf[1] &= ( ~ ( 1 << 7 ) ) ;
        ( (USHORT*) (leaf + leaf_node_bytes - (flood_bytes + CINDY_BYTES)))[0] &= ( ~ ( 1 << eind ) ) ;
    }

    int isInProcess( UCHAR* leaf, int eind )
    {
        //int a = ( ( leaf[1] >> 7 ) & 1 ) ;
        int a = ( ( ( (USHORT*) (leaf + leaf_node_bytes - (flood_bytes + CINDY_BYTES)))[0] >> eind ) & 1 ) ;
        return a ;
    }

#ifndef USE_HERMIT

    void setEdgeIntersectionIndex( UCHAR* leaf, int count, int index )
    {
        ((int *)( leaf + leaf_node_bytes ))[ count ] = index ;
    }

    int getEdgeIntersectionIndex( UCHAR* leaf, int count )
    {
        return  ((int *)( leaf + leaf_node_bytes ))[ count ] ;
    }

    void fillEdgeIntersectionIndices( UCHAR* leaf, int st[3], int len, int inds[12] )
    {
        int i ;

        // The three primal edges are easy
        int pmask[3] = { 0, 4, 8 } ;
        for ( i = 0 ; i < 3 ; i ++ )
        {
            if ( getEdgeParity( leaf, pmask[i] ) )
            {
                inds[pmask[i]] = getEdgeIntersectionIndex( leaf, getEdgeCount( leaf, 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 ;
                UCHAR* node = locateLeaf( nst ) ;
                
                if ( e1 )
                {
                    inds[ fmask[i][0] ] = getEdgeIntersectionIndex( node, getEdgeCount( node, femask[i][0] ) ) ;
                }
                if ( e2 )
                {
                    inds[ fmask[i][1] ] = getEdgeIntersectionIndex( node, getEdgeCount( node, femask[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 ;
                UCHAR* node = locateLeaf( nst ) ;
                
                inds[ emask[i] ] = getEdgeIntersectionIndex( node, getEdgeCount( node, eemask[i] ) ) ;
            }
        }
    }


#endif
    
    void generateSigns ( UCHAR* leaf, UCHAR table[], int start )
    {
        leaf[2] = table[ ( ((USHORT *) leaf)[ 0 ] ) & ( ( 1 << 12 ) - 1 ) ] ; 

        if ( ( start ^ leaf[2] ) & 1 ) 
        {
            leaf[2] = ~ ( leaf[2] ) ;
        }
    }

    int getEdgeParity( UCHAR* leaf, int index ) 
    {
        int a = ( ( ((USHORT *) leaf)[ 0 ] >> index ) & 1 ) ;
        return  a ;
    };

    int getFaceParity ( UCHAR* 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 ( UCHAR* 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( UCHAR* leaf, int index ) 
    {
        ((USHORT *) leaf)[ 0 ] ^= ( 1 << index ) ;  
    };
    void setEdge( UCHAR* leaf, int index ) 
    {
        ((USHORT *) leaf)[ 0 ] |= ( 1 << index ) ;  
    };
    void resetEdge( UCHAR* leaf, int index ) 
    {
        ((USHORT *) leaf)[ 0 ] &= ( ~ ( 1 << index ) ) ;    
    };

    void createPrimalEdgesMask( UCHAR* leaf )
    {
        int mask = (( leaf[0] & 1 ) | ( (leaf[0] >> 3) & 2 ) | ( (leaf[1] & 1) << 2 ) ) ;
        leaf[1] |= ( mask << 4 ) ;

    }

    void setStoredEdgesParity( UCHAR* leaf, int pindex )
    {
        leaf[1] |= ( 1 << ( 4 + pindex ) ) ;
    }
    int getStoredEdgesParity( UCHAR* leaf, int pindex )
    {
        return ( ( leaf[1] >> ( 4 + pindex ) ) & 1 ) ;
    }

    UCHAR* flipEdge( UCHAR* 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 ) ;
                UCHAR* nleaf = createLeaf( num ) ;
                for ( int i = 0 ; i < leaf_node_bytes ; i ++ )
                {
                    nleaf[i] = leaf[i] ;
                }

                setEdgeOffset( nleaf, alpha, count ) ;

                if ( num > 1 )
                {
                    float * pts = ( float * ) ( leaf + leaf_node_bytes ) ;
                    float * npts = ( float * ) ( nleaf + leaf_node_bytes ) ;
                    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, leaf ) ;
                leaf = nleaf ;
            }
        }

        return leaf ;
    };

    void updateParent( UCHAR* node, int len, int st[3], UCHAR* leaf ) 
    {
        // First, locate the parent
        int count ;
        UCHAR* parent = locateParent( node, len, st, count ) ;

        // UPdate
        setChild( parent, count, leaf ) ;
    }

    void updateParent( UCHAR* node, int len, int st[3] ) 
    {
        if ( len == dimen )
        {
            root = node ;
            return ;
        }

        // First, locate the parent
        int count ;
        UCHAR* parent = locateParent( len, st, count ) ;

        // UPdate
        setChild( parent, count, node ) ;
    }

    int getEdgeIntersectionByIndex( int st[3], int index, float pt[3], int check )
    {
        // First, locat the leaf
        UCHAR* 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( UCHAR* leaf )
    {
        // return (( leaf[0] & 1 ) | ( (leaf[0] >> 3) & 2 ) | ( (leaf[1] & 1) << 2 ) ) ;
        return ( ( leaf[1] >> 4 ) & 7 ) ;
    }

    int getPrimalEdgesMask2( UCHAR* leaf )
    {
        return (( leaf[0] & 1 ) | ( (leaf[0] >> 3) & 2 ) | ( (leaf[1] & 1) << 2 ) ) ;
    }

    int getEdgeCount( UCHAR* leaf, int index )
    {
        return edgeCountTable[ getPrimalEdgesMask( leaf ) ][ index ] ;
    }
    int getNumEdges( UCHAR* leaf )
    {
        return numEdgeTable[ getPrimalEdgesMask( leaf ) ] ;
    }

    int getNumEdges2( UCHAR* leaf )
    {
        return numEdgeTable[ getPrimalEdgesMask2( leaf ) ] ;
    }

    void setEdgeOffset( UCHAR* leaf, float pt, int count )
    {
        float * pts = ( float * ) ( leaf + leaf_node_bytes ) ;
#ifdef USE_HERMIT
        pts[ EDGE_FLOATS * count ] = pt ;
        pts[ EDGE_FLOATS * count + 1 ] = 0 ;
        pts[ EDGE_FLOATS * count + 2 ] = 0 ;
        pts[ EDGE_FLOATS * count + 3 ] = 0 ;
#else
        pts[ count ] = pt ;
#endif
    }

    void setEdgeOffsets( UCHAR* leaf, float pt[3], int len )
    {
        float * pts = ( float * ) ( leaf + leaf_node_bytes ) ;
        for ( int i = 0 ; i < len ; i ++ )
        {
            pts[i] = pt[i] ;
        }
    }

    float getEdgeOffset( UCHAR* leaf, int count )
    {
#ifdef USE_HERMIT
        return (( float * ) ( leaf + leaf_node_bytes ))[ 4 * count ] ;
#else
        return (( float * ) ( leaf + leaf_node_bytes ))[ count ] ;
#endif
    }

    UCHAR* updateEdgeOffsets( UCHAR* leaf, int oldlen, int newlen, float offs[3] )
    {
        // First, create a new leaf node
        UCHAR* nleaf = createLeaf( newlen ) ;
        for ( int i = 0 ; i < leaf_node_bytes ; i ++ )
        {
            nleaf[i] = leaf[i] ;
        }

        // Next, fill in the offsets
        setEdgeOffsets( nleaf, offs, newlen ) ;

        // Finally, delete the old leaf
        removeLeaf( oldlen, leaf ) ;

        return nleaf ;
    }

    void setOriginalIndex( UCHAR* leaf, int index )
    {
        ((int *)( leaf + leaf_node_bytes ))[ 0 ] = index ;
    }
    int getOriginalIndex( UCHAR* leaf )
    {
        return  ((int *)( leaf + leaf_node_bytes ))[ 0 ] ;
    }
#ifdef USE_HERMIT

    void setMinimizerIndex( UCHAR* leaf, int index )
    {
        ((int *)( leaf + leaf_node_bytes - leaf_extra_bytes - 4 ))[ 0 ] = index ;
    }

    int getMinimizerIndex( UCHAR* leaf )
    {
        return ((int *)( leaf + leaf_node_bytes - leaf_extra_bytes - 4 ))[ 0 ] ;
    }
    
    int getMinimizerIndex( UCHAR* leaf, int eind )
    {
        int add = manifold_table[ getSignMask( leaf ) ].pairs[ eind ][ 0 ] - 1 ;
        if ( add < 0 )
        {
            printf("Manifold components wrong!\n") ;
        }
        return ((int *)( leaf + leaf_node_bytes - leaf_extra_bytes - 4 ))[ 0 ] + add ;
    }

    void getMinimizerIndices( UCHAR* leaf, int eind, int inds[2] )
    {
        const int* add = manifold_table[ getSignMask( leaf ) ].pairs[ eind ] ;
        inds[0] = ((int *)( leaf + leaf_node_bytes - leaf_extra_bytes - 4 ))[ 0 ] + add[0] - 1 ;
        if ( add[0] == add[1] )
        {
            inds[1] = -1 ;
        }
        else
        {
            inds[1] = ((int *)( leaf + leaf_node_bytes - leaf_extra_bytes - 4 ))[ 0 ] + add[1] - 1 ;
        }
    }


    void setEdgeOffsetNormal( UCHAR* leaf, float pt, float a, float b, float c, int count )
    {
        float * pts = ( float * ) ( leaf + leaf_node_bytes ) ;
        pts[ 4 * count ] = pt ;
        pts[ 4 * count + 1 ] = a ;
        pts[ 4 * count + 2 ] = b ;
        pts[ 4 * count + 3 ] = c ;
    }

    float getEdgeOffsetNormal( UCHAR* leaf, int count, float& a, float& b, float& c )
    {
        float * pts = ( float * ) ( leaf + leaf_node_bytes ) ;
        a = pts[ 4 * count + 1 ] ;
        b = pts[ 4 * count + 2 ] ;
        c = pts[ 4 * count + 3 ] ;
        return pts[ 4 * count ] ;
    }

    void setEdgeOffsetsNormals( UCHAR* leaf, float pt[], float a[], float b[], float c[], int len )
    {
        float * pts = ( float * ) ( leaf + leaf_node_bytes ) ;
        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( UCHAR* leaf, int index, int st[3], int len, float pt[3], float nm[3] )
    {
        int count = getEdgeCount( leaf, index ) ;
        float * pts = ( float * ) ( leaf + leaf_node_bytes ) ;
        
        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( UCHAR* leaf, int index, float nm[3] )
    {
        int count = getEdgeCount( leaf, index ) ;
        float * pts = ( float * ) ( leaf + leaf_node_bytes ) ;
        
        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( UCHAR* leaf, int st[3], int len, float pts[12][3], float norms[12][3] )
    {
        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 ;
                UCHAR* 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 ;
                UCHAR* 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( UCHAR* leaf, int st[3], int len, float pts[12][3], 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 ) )
            {
                // 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 ;
                UCHAR* 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 ;
                UCHAR* 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( UCHAR* 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 ;
                UCHAR* 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 ;
                UCHAR* 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 ;
                }
            }
        }
    }


    UCHAR* updateEdgeOffsetsNormals( UCHAR* leaf, int oldlen, int newlen, float offs[3], float a[3], float b[3], float c[3] )
    {
        // First, create a new leaf node
        UCHAR* nleaf = createLeaf( newlen ) ;
        for ( int i = 0 ; i < leaf_node_bytes ; i ++ )
        {
            nleaf[i] = leaf[i] ;
        }

        // Next, fill in the offsets
        setEdgeOffsetsNormals( nleaf, offs, a, b, c, newlen ) ;

        // Finally, delete the old leaf
        removeLeaf( oldlen, leaf ) ;

        return nleaf ;
    }
#endif

    
    UCHAR* locateLeaf( int st[3] )
    {
        UCHAR* 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 = getChild( node, getChildCount( node, index ) ) ;
        }

        return node ;
    }
    
    UCHAR* locateLeaf( UCHAR* node, int len, int st[3] )
    {
        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 = getChild( node, getChildCount( node, index ) ) ;
        }

        return node ;
    }
    UCHAR* locateLeafCheck( int st[3] )
    {
        UCHAR* 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 ( ! hasChild( node, index ) )
            {
                return NULL ;
            }
            node = getChild( node, getChildCount( node, index ) ) ;
        }

        return node ;
    }
    UCHAR* locateParent( int len, int st[3], int& count )
    {
        UCHAR* node = root ;
        UCHAR* 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 = getChild( node, getChildCount( node, index ) ) ;
        }

        count = getChildCount( pre, index ) ;
        return pre ;
    }
    UCHAR* locateParent( UCHAR* papa, int len, int st[3], int& count )
    {
        UCHAR* node = papa ;
        UCHAR* 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 = getChild( node, getChildCount( node, index ) ) ;
        }

        count = getChildCount( pre, index ) ;
        return pre ;
    }
    /************ Operators for internal nodes ************/

    void printNode( UCHAR* node )
    {
        printf("Address: %p ", node ) ;
        printf("Leaf Mask: ") ;
        for ( int i = 0 ; i < 8 ; i ++ )
        {
            printf( "%d ", isLeaf( node, i ) ) ;
        }
        printf("Child Mask: ") ;
        for ( int i = 0 ; i < 8 ; i ++ )
        {
            printf( "%d ", hasChild( node, i ) ) ;
        }
        printf("\n") ;
    }

    int getSize ( int length )
    {
        return INTERNAL_NODE_BYTES + length * 4 ;   
    };

    int hasChild( UCHAR* node, int index )
    {
        return ( node[0] >> index ) & 1 ;
    };

    int isLeaf ( UCHAR* node, int index )
    {
        return ( node[1] >> index ) & 1 ;
    };

    UCHAR* getChild ( UCHAR* node, int count )
    {
        return (( UCHAR ** ) ( node + INTERNAL_NODE_BYTES )) [ count ] ;    
    };

    int getNumChildren( UCHAR* node )
    {
        return numChildrenTable[ node[0] ] ;
    };

    int getChildCount( UCHAR* node, int index )
    {
        return childrenCountTable[ node[0] ][ index ] ;
    }
    int getChildIndex( UCHAR* node, int count )
    {
        return childrenIndexTable[ node[0] ][ count ] ;
    }
    int* getChildCounts( UCHAR* node )
    {
        return childrenCountTable[ node[0] ] ;
    }

    void fillChildren( UCHAR* node, UCHAR* chd[ 8 ], int leaf[ 8 ] )
    {
        int count = 0 ;
        for ( int i = 0 ; i < 8 ; i ++ )
        {   
            leaf[ i ] = isLeaf( node, i ) ;
            if ( hasChild( node, i ) )
            {
                chd[ i ] = getChild( node, count ) ;
                count ++ ;
            }
            else
            {
                chd[ i ] = NULL ;
                leaf[ i ] = 0 ;
            }
        }
    }

    void setChild ( UCHAR* node, int count, UCHAR* chd )
    {
        (( UCHAR ** ) ( node + INTERNAL_NODE_BYTES )) [ count ] = chd ;
    }
    void setInternalChild ( UCHAR* node, int index, int count, UCHAR* chd )
    {
        setChild( node, count, chd ) ;
        node[0] |= ( 1 << index ) ;
    };
    void setLeafChild ( UCHAR* node, int index, int count, UCHAR* chd )
    {
        setChild( node, count, chd ) ;
        node[0] |= ( 1 << index ) ;
        node[1] |= ( 1 << index ) ;
    };

    UCHAR* addChild( UCHAR* node, int index, UCHAR* chd, int aLeaf )
    {
        // Create new internal node
        int num = getNumChildren( node ) ;
        UCHAR* rnode = createInternal( num + 1 ) ;

        // Establish children
        int i ;
        int count1 = 0, count2 = 0 ;
        for ( i = 0 ; i < 8 ; i ++ )
        {
            if ( i == index )
            {
                if ( aLeaf )
                {
                    setLeafChild( rnode, i, count2, chd ) ;
                }
                else
                {
                    setInternalChild( rnode, i, count2, chd ) ;
                }
                count2 ++ ;
            }
            else if ( hasChild( node, i ) )
            {
                if ( isLeaf( node, i ) )
                {
                    setLeafChild( rnode, i, count2, getChild( node, count1 ) ) ;
                }
                else
                {
                    setInternalChild( rnode, i, count2, getChild( node, count1 ) ) ;
                }
                count1 ++ ;
                count2 ++ ;
            }
        }

        removeInternal( num, node ) ;
        return rnode ;
    }

    UCHAR* createInternal( int length )
    {
        UCHAR* inode = alloc[ length ]->allocate( ) ;
        inode[0] = inode[1] = 0 ;
        return inode ;
    };
    UCHAR* createLeaf( int length )
    {
        if ( length > 3 )
        {
            printf("wierd");
        }
        UCHAR* lnode = leafalloc[ length ]->allocate( ) ;
        lnode[0] = lnode[1] = lnode[2] = 0 ;

        return lnode ;
    };

    void removeInternal ( int num, UCHAR* node )
    {
        alloc[ num ]->deallocate( node ) ;
    }

    void removeLeaf ( int num, UCHAR* leaf )
    {
        if ( num > 3 || num < 0 )
        {
            printf("wierd");
        }
        leafalloc[ num ]->deallocate( leaf ) ;
    }

    UCHAR* addLeafChild ( UCHAR* par, int index, int count, UCHAR* leaf )
    {
        int num = getNumChildren( par ) + 1 ;
        UCHAR* npar = createInternal( num ) ;
        npar[0] = par[0] ;
        npar[1] = par[1] ;
        
        if ( num == 1 )
        {
            setLeafChild( npar, index, 0, leaf ) ;
        }
        else
        {
            int i ;
            for ( i = 0 ; i < count ; i ++ )
            {
                setChild( npar, i, getChild( par, i ) ) ;
            }
            setLeafChild( npar, index, count, leaf ) ;
            for ( i = count + 1 ; i < num ; i ++ )
            {
                setChild( npar, i, getChild( par, i - 1 ) ) ;
            }
        }
        
        removeInternal( num-1, par ) ;
        return npar ;
    };

    UCHAR* addInternalChild ( UCHAR* par, int index, int count, UCHAR* node )
    {
        int num = getNumChildren( par ) + 1 ;
        UCHAR* npar = createInternal( num ) ;
        npar[0] = par[0] ;
        npar[1] = par[1] ;
        
        if ( num == 1 )
        {
            setInternalChild( npar, index, 0, node ) ;
        }
        else
        {
            int i ;
            for ( i = 0 ; i < count ; i ++ )
            {
                setChild( npar, i, getChild( par, i ) ) ;
            }
            setInternalChild( npar, index, count, node ) ;
            for ( i = count + 1 ; i < num ; i ++ )
            {
                setChild( npar, i, getChild( par, i - 1 ) ) ;
            }
        }
        
        removeInternal( num-1, par ) ;
        return npar ;
    };
};



#endif