diff options
Diffstat (limited to 'intern/dualcon/intern/octree.cpp')
| -rw-r--r-- | intern/dualcon/intern/octree.cpp | 4311 |
1 files changed, 4311 insertions, 0 deletions
diff --git a/intern/dualcon/intern/octree.cpp b/intern/dualcon/intern/octree.cpp new file mode 100644 index 00000000000..90dbf5376a2 --- /dev/null +++ b/intern/dualcon/intern/octree.cpp @@ -0,0 +1,4311 @@ +/* + * ***** 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 ***** + */ + +#include "octree.h" +#include <Eigen/Dense> +#include <limits> +#include <time.h> + +/** + * Implementations of Octree member functions. + * + * @author Tao Ju + */ + +/* set to non-zero value to enable debugging output */ +#define DC_DEBUG 0 + +#if DC_DEBUG +/* enable debug printfs */ +#define dc_printf printf +#else +/* disable debug printfs */ +#define dc_printf(...) do {} while(0) +#endif + +Octree::Octree( ModelReader* mr, + DualConAllocOutput alloc_output_func, + DualConAddVert add_vert_func, + DualConAddQuad add_quad_func, + DualConFlags flags, DualConMode dualcon_mode, int depth, + float threshold, float sharpness ) + : use_flood_fill(flags & DUALCON_FLOOD_FILL), + /* note on `use_manifold': + + After playing around with this option, the only case I could + find where this option gives different results is on + relatively thin corners. Sometimes along these corners two + vertices from seperate sides will be placed in the same + position, so hole gets filled with a 5-sided face, where two + of those vertices are in the same 3D location. If + `use_manifold' is disabled, then the modifier doesn't + generate two separate vertices so the results end up as all + quads. + + Since the results are just as good with all quads, there + doesn't seem any reason to allow this to be toggled in + Blender. -nicholasbishop + */ + use_manifold(false), + hermite_num(sharpness), + mode(dualcon_mode), + alloc_output(alloc_output_func), + add_vert(add_vert_func), + add_quad(add_quad_func) +{ + this->thresh = threshold ; + this->reader = mr ; + this->dimen = 1 << GRID_DIMENSION ; + this->range = reader->getBoundingBox( this->origin ) ; + this->nodeCount = this->nodeSpace = 0; + this->maxDepth = depth ; + this->mindimen = ( dimen >> maxDepth ) ; + this->minshift = ( GRID_DIMENSION - maxDepth ) ; + this->buildTable( ) ; + + flood_bytes = use_flood_fill ? FLOOD_FILL_BYTES : 0; + leaf_extra_bytes = flood_bytes + CINDY_BYTES; + +#ifdef USE_HERMIT + leaf_node_bytes = 7 + leaf_extra_bytes; +#else + leaf_node_bytes = 3 + leaf_extra_bytes; +#endif + +#ifdef QIANYI + dc_printf("Origin: (%f %f %f), Dimension: %f\n", origin[0], origin[1], origin[2], range) ; +#endif + + this->maxTrianglePerCell = 0 ; + + // Initialize memory +#ifdef IN_VERBOSE_MODE + dc_printf("Range: %f origin: %f, %f,%f \n", range, origin[0], origin[1], origin[2] ) ; + dc_printf("Initialize memory...\n") ; +#endif + initMemory( ) ; + this->root = createInternal( 0 ) ; + + // Read MC table +#ifdef IN_VERBOSE_MODE + dc_printf("Reading contour table...\n") ; +#endif + this->cubes = new Cubes(); + +} + +Octree::~Octree( ) +{ + freeMemory( ) ; +} + +void Octree::scanConvert() +{ + // Scan triangles +#if DC_DEBUG + clock_t start, finish ; + start = clock( ) ; +#endif + + this->addTrian( ) ; + this->resetMinimalEdges( ) ; + this->preparePrimalEdgesMask( this->root ) ; + +#if DC_DEBUG + finish = clock( ) ; + dc_printf("Time taken: %f seconds \n", + (double)(finish - start) / CLOCKS_PER_SEC ) ; +#endif + + // Generate signs + // Find holes +#if DC_DEBUG + dc_printf("Patching...\n") ; + start = clock( ) ; +#endif + this->trace( ) ; +#if DC_DEBUG + finish = clock( ) ; + dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC ) ; +#ifdef IN_VERBOSE_MODE + dc_printf("Holes: %d Average Length: %f Max Length: %d \n", numRings, (float)totRingLengths / (float) numRings, maxRingLength ) ; +#endif +#endif + + // Check again + int tnumRings = numRings ; + this->trace( ) ; +#ifdef IN_VERBOSE_MODE + dc_printf("Holes after patching: %d \n", numRings) ; +#endif + numRings = tnumRings ; + +#if DC_DEBUG + dc_printf("Building signs...\n") ; + start = clock( ) ; +#endif + this->buildSigns( ) ; +#if DC_DEBUG + finish = clock( ) ; + dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC ) ; +#endif + + if(use_flood_fill) { + /* + start = clock( ) ; + this->floodFill( ) ; + // Check again + tnumRings = numRings ; + this->trace( ) ; + dc_printf("Holes after filling: %d \n", numRings) ; + numRings = tnumRings ; + this->buildSigns( ) ; + finish = clock( ) ; + dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC ) ; + */ +#if DC_DEBUG + start = clock( ) ; + dc_printf("Removing components...\n"); +#endif + this->floodFill( ) ; + this->buildSigns( ) ; + // dc_printf("Checking...\n"); + // this->floodFill( ) ; +#if DC_DEBUG + finish = clock( ) ; + dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC ) ; +#endif + } + + // Output +#ifdef OUTPUT_REPAIRED +#if DC_DEBUG + start = clock( ) ; +#endif + writeOut(); +#if DC_DEBUG + finish = clock( ) ; +#endif + // dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC ) ; +#ifdef CINDY + this->writeTags( "tags.txt" ) ; + dc_printf("Tags output to tags.txt\n") ; +#endif + +#endif + + // Print info +#ifdef IN_VERBOSE_MODE + printMemUsage( ) ; +#endif +} + +#if 0 +void Octree::writeOut( char* fname ) +{ + dc_printf( "\n" ) ; + if ( strstr( fname, ".ply" ) != NULL ) + { + dc_printf("Writing PLY file format.\n") ; + this->outType = 1 ; + writePLY( fname ) ; + } + else if ( strstr( fname, ".off" ) != NULL ) + { + dc_printf("Writing OFF file format.\n") ; + this->outType = 0 ; + writeOFF( fname ) ; + } + else if ( strstr( fname, ".sof" ) != NULL ) + { + dc_printf("Writing Signed Octree File format.\n") ; + this->outType = 2 ; + writeOctree( fname ) ; + } + else if ( strstr( fname, ".dcf" ) != NULL ) + { +#ifdef USE_HERMIT + dc_printf("Writing Dual Contouring File format.\n") ; + this->outType = 3 ; + writeDCF( fname ) ; +#else + dc_printf("Can not write Dual Contouring File format in non-DC mode.\n") ; +#endif + } +#ifdef USE_HERMIT + else if ( strstr( fname, ".sog" ) != NULL ) + { + dc_printf("Writing signed octree with geometry.\n") ; + this->outType = 4 ; + writeOctreeGeom( fname ) ; + } +#endif + /* + else if ( strstr( fname, ".sof" ) != NULL ) + { + dc_printf("Writing SOF file format.\n") ; + this->outType = 2 ; + writeOctree( fname ) ; + } + */ + else + { + dc_printf("Unknown output format.\n") ; + } + +} +#endif + +void Octree::initMemory( ) +{ +#ifdef USE_HERMIT + const int leaf_node_bytes = 7; +#else + const int leaf_node_bytes = 3; +#endif + + if(use_flood_fill) { + const int bytes = leaf_node_bytes + CINDY_BYTES + FLOOD_FILL_BYTES; + this->leafalloc[ 0 ] = new MemoryAllocator< bytes > ( ) ; + this->leafalloc[ 1 ] = new MemoryAllocator< bytes + EDGE_BYTES > ( ) ; + this->leafalloc[ 2 ] = new MemoryAllocator< bytes + EDGE_BYTES * 2 > ( ) ; + this->leafalloc[ 3 ] = new MemoryAllocator< bytes + EDGE_BYTES * 3 > ( ) ; + } + else { + const int bytes = leaf_node_bytes + CINDY_BYTES; + this->leafalloc[ 0 ] = new MemoryAllocator< bytes > ( ) ; + this->leafalloc[ 1 ] = new MemoryAllocator< bytes + EDGE_BYTES > ( ) ; + this->leafalloc[ 2 ] = new MemoryAllocator< bytes + EDGE_BYTES * 2 > ( ) ; + this->leafalloc[ 3 ] = new MemoryAllocator< bytes + EDGE_BYTES * 3 > ( ) ; + } + + this->alloc[ 0 ] = new MemoryAllocator< INTERNAL_NODE_BYTES > ( ) ; + this->alloc[ 1 ] = new MemoryAllocator< INTERNAL_NODE_BYTES + POINTER_BYTES > ( ) ; + this->alloc[ 2 ] = new MemoryAllocator< INTERNAL_NODE_BYTES + POINTER_BYTES*2 > ( ) ; + this->alloc[ 3 ] = new MemoryAllocator< INTERNAL_NODE_BYTES + POINTER_BYTES*3 > ( ) ; + this->alloc[ 4 ] = new MemoryAllocator< INTERNAL_NODE_BYTES + POINTER_BYTES*4 > ( ) ; + this->alloc[ 5 ] = new MemoryAllocator< INTERNAL_NODE_BYTES + POINTER_BYTES*5 > ( ) ; + this->alloc[ 6 ] = new MemoryAllocator< INTERNAL_NODE_BYTES + POINTER_BYTES*6 > ( ) ; + this->alloc[ 7 ] = new MemoryAllocator< INTERNAL_NODE_BYTES + POINTER_BYTES*7 > ( ) ; + this->alloc[ 8 ] = new MemoryAllocator< INTERNAL_NODE_BYTES + POINTER_BYTES*8 > ( ) ; +} + +void Octree::freeMemory( ) +{ + for ( int i = 0 ; i < 9 ; i ++ ) + { + alloc[i]->destroy() ; + delete alloc[i] ; + } + + for ( int i = 0 ; i < 4 ; i ++ ) + { + leafalloc[i]->destroy() ; + delete leafalloc[i] ; + } +} + +void Octree::printMemUsage( ) +{ + int totalbytes = 0 ; + dc_printf("********* Internal nodes: \n") ; + for ( int i = 0 ; i < 9 ; i ++ ) + { + this->alloc[ i ]->printInfo() ; + + totalbytes += alloc[i]->getAll( ) * alloc[i]->getBytes() ; + } + dc_printf("********* Leaf nodes: \n") ; + int totalLeafs = 0 ; + for ( int i = 0 ; i < 4 ; i ++ ) + { + this->leafalloc[ i ]->printInfo() ; + + totalbytes += leafalloc[i]->getAll( ) * leafalloc[i]->getBytes() ; + totalLeafs += leafalloc[i]->getAllocated() ; + } + + dc_printf("Total allocated bytes on disk: %d \n", totalbytes) ; + dc_printf("Total leaf nodes: %d\n", totalLeafs ) ; +} + +void Octree::resetMinimalEdges( ) +{ + this->cellProcParity( this->root, 0, maxDepth ) ; +} + +void Octree::writeModel( char* fname ) +{ + reader->reset() ; + + int nFace = reader->getNumTriangles() ; + Triangle* trian ; + // int unitcount = 10000; + int count = 0 ; + int nVert = nFace * 3 ; + FILE* modelfout = fopen( "model.off", "w" ) ; + fprintf( modelfout, "OFF\n" ) ; + fprintf( modelfout, "%d %d 0\n", nVert, nFace ) ; + + //int total = this->reader->getNumTriangles() ; + dc_printf( "Start writing model to OFF...\n" ) ; + srand(0) ; + while ( ( trian = reader->getNextTriangle() ) != NULL ) + { + // Drop polygons + { + int i, j ; + + // Blow up the triangle + float mid[3] = {0, 0, 0} ; + for ( i = 0 ; i < 3 ; i ++ ) + for ( j = 0 ; j < 3 ; j ++ ) + { + trian->vt[i][j] = dimen * ( trian->vt[i][j] - origin[j] ) / range ; + + mid[j] += trian->vt[i][j] / 3 ; + } + + // Generate projections + // LONG cube[2][3] = { { 0, 0, 0 }, { dimen, dimen, dimen } } ; + int trig[3][3] ; + + // Blowing up the triangle to the grid + for ( i = 0 ; i < 3 ; i ++ ) + for ( j = 0 ; j < 3 ; j ++ ) + { + trig[i][j] = (int) (trian->vt[i][j]) ; + // Perturb end points, if set so + } + + + for ( i = 0 ; i < 3 ; i ++ ) + { + fprintf( modelfout, "%f %f %f\n", + (float)(((double) trig[i][0] / dimen) * range + origin[0]) , + (float)(((double) trig[i][1] / dimen) * range + origin[1]) , + (float)(((double) trig[i][2] / dimen) * range + origin[2]) ) ; + } + } + delete trian ; + + count ++ ; + + } + + for ( int i = 0 ; i < nFace ; i ++ ) + { + fprintf( modelfout, "3 %d %d %d\n", 3 * i + 2, 3 * i + 1, 3 * i ) ; + } + + fclose( modelfout ) ; + +} + +#ifdef CINDY +void Octree::writeTags( char* fname ) +{ + FILE* fout = fopen( fname, "w" ) ; + + clearCindyBits( root, maxDepth ) ; + readVertices() ; + outputTags( root, maxDepth, fout ) ; + + fclose ( fout ) ; +} + +void Octree::readVertices( ) +{ + int total = this->reader->getNumVertices() ; + reader->reset() ; + float v[3] ; + int st[3] = {0,0,0}; + int unitcount = 1000 ; + dc_printf( "\nRead in original %d vertices...\n", total ) ; + + for ( int i = 0 ; i < total ; i ++ ) + { + reader->getNextVertex( v ) ; + // Blowing up the triangle to the grid + float mid[3] = {0, 0, 0} ; + for ( int j = 0 ; j < 3 ; j ++ ) + { + v[j] = dimen * ( v[j] - origin[j] ) / range ; + } + +// dc_printf("vertex: %f %f %f, dimen: %d\n", v[0], v[1], v[2], dimen ) ; + readVertex ( root, st, dimen, maxDepth, v, i ) ; + + + if ( i % unitcount == 0 ) + { + putchar ( 13 ) ; + + switch ( ( i / unitcount ) % 4 ) + { + case 0 : dc_printf("-"); + break ; + case 1 : dc_printf("/") ; + break ; + case 2 : dc_printf("|"); + break ; + case 3 : dc_printf("\\") ; + break ; + } + + float percent = (float) i / total ; + /* + int totbars = 50 ; + int bars = (int)( percent * totbars ) ; + for ( int i = 0 ; i < bars ; i ++ ) + { + putchar( 219 ) ; + } + for ( i = bars ; i < totbars ; i ++ ) + { + putchar( 176 ) ; + } + */ + + dc_printf(" %d vertices: ", i ) ; + dc_printf( " %f%% complete.", 100 * percent ) ; + } + + } + putchar ( 13 ) ; + dc_printf(" \n"); +} + +void Octree::readVertex( UCHAR* node, int st[3], int len, int height, float v[3], int index ) +{ + int nst[3] ; + nst[0] = ( (int) v[0] / mindimen ) * mindimen ; + nst[1] = ( (int) v[1] / mindimen ) * mindimen ; + nst[2] = ( (int) v[2] / mindimen ) * mindimen ; + + UCHAR* cell = this->locateLeafCheck( nst ) ; + if ( cell == NULL ) + { + dc_printf("Cell %d %d %d is not found!\n", nst[0]/ mindimen, nst[1]/ mindimen, nst[2]/ mindimen) ; + return ; + } + + setOriginalIndex( cell, index ) ; + + + /* + int i ; + if ( height == 0 ) + { + // Leaf cell, assign index + dc_printf("Setting: %d\n", index ) ; + setOriginalIndex( node, index ) ; + } + else + { + len >>= 1 ; + // Internal cell, check and recur + int x = ( v[0] > st[0] + len ) ? 1 : 0 ; + int y = ( v[1] > st[1] + len ) ? 1 : 0 ; + int z = ( v[2] > st[2] + len ) ? 1 : 0 ; + int child = x * 4 + y * 2 + z ; + + int count = 0 ; + for ( i = 0 ; i < 8 ; i ++ ) + { + if ( i == child && hasChild( node, i ) ) + { + int nst[3] ; + nst[0] = st[0] + vertmap[i][0] * len ; + nst[1] = st[1] + vertmap[i][1] * len ; + nst[2] = st[2] + vertmap[i][2] * len ; + + dc_printf("Depth: %d -- child %d vertex: %f %f %f in %f %f %f\n", height - 1, child, v[0]/mindimen, v[1]/mindimen, v[2]/mindimen, + nst[0]/mindimen, nst[1]/mindimen, nst[2]/mindimen, len/mindimen ) ; + + readVertex( getChild( node, count ), nst, len, height - 1, v, index ) ; + count ++ ; + } + } + } + */ +} + +void Octree::outputTags( UCHAR* node, int height, FILE* fout ) +{ + int i ; + + if ( height == 0 ) + { + // Leaf cell, generate + int smask = getSignMask( node ) ; + + if(use_manifold) { + int comps = manifold_table[ smask ].comps ; + if ( comps != 1 ) + { + return ; + } + } + else + { + if ( smask == 0 || smask == 255 ) + { + return ; + } + } + + int rindex = getMinimizerIndex( node ) ; + int oindex = getOriginalIndex( node ) ; + + if ( oindex >= 0 ) + { + fprintf( fout, "%d: %d\n", rindex, oindex ) ; + } + + } + else + { + // Internal cell, recur + int count = 0 ; + for ( i = 0 ; i < 8 ; i ++ ) + { + if ( hasChild( node, i ) ) + { + outputTags( getChild( node, count ), height - 1, fout ) ; + count ++ ; + } + } + } +} + +void Octree::clearCindyBits( UCHAR* node, int height ) +{ + int i; + + if ( height == 0 ) + { + // Leaf cell, + { + setOriginalIndex( node, - 1 ) ; + } + } + else + { + // Internal cell, recur + int count = 0 ; + for ( i = 0 ; i < 8 ; i ++ ) + { + if ( hasChild( node, i ) ) + { + clearCindyBits( getChild( node, count ), height - 1 ) ; + count ++ ; + } + } + } +} +#endif + +void Octree::addTrian( ) +{ + Triangle* trian ; + int count = 0 ; + +#if DC_DEBUG + int total = this->reader->getNumTriangles() ; + int unitcount = 1000 ; + dc_printf( "\nScan converting to depth %d...\n", maxDepth ) ; +#endif + + srand(0) ; + + while ( ( trian = reader->getNextTriangle() ) != NULL ) + { + // Drop triangles + { + addTrian ( trian, count ) ; + } + delete trian ; + + count ++ ; + +#if DC_DEBUG + if ( count % unitcount == 0 ) + { + putchar ( 13 ) ; + + switch ( ( count / unitcount ) % 4 ) + { + case 0 : dc_printf("-"); + break ; + case 1 : dc_printf("/") ; + break ; + case 2 : dc_printf("|"); + break ; + case 3 : dc_printf("\\") ; + break ; + } + + float percent = (float) count / total ; + + /* + int totbars = 50 ; + int bars = (int)( percent * totbars ) ; + for ( int i = 0 ; i < bars ; i ++ ) + { + putchar( 219 ) ; + } + for ( i = bars ; i < totbars ; i ++ ) + { + putchar( 176 ) ; + } + */ + + dc_printf(" %d triangles: ", count ) ; + dc_printf( " %f%% complete.", 100 * percent ) ; + } +#endif + + } + putchar ( 13 ) ; +} + +void Octree::addTrian( Triangle* trian, int triind ) +{ + int i, j ; + + // Blowing up the triangle to the grid + float mid[3] = {0, 0, 0} ; + for ( i = 0 ; i < 3 ; i ++ ) + for ( j = 0 ; j < 3 ; j ++ ) + { + trian->vt[i][j] = dimen * ( trian->vt[i][j] - origin[j] ) / range ; + mid[j] += trian->vt[i][j] / 3 ; + } + + // Generate projections + LONG cube[2][3] = { { 0, 0, 0 }, { dimen, dimen, dimen } } ; + LONG trig[3][3] ; + + for ( i = 0 ; i < 3 ; i ++ ) + for ( j = 0 ; j < 3 ; j ++ ) + { + trig[i][j] = (LONG) (trian->vt[i][j]) ; + // Perturb end points, if set so + } + + // Add to the octree + // int start[3] = { 0, 0, 0 } ; + LONG errorvec = (LONG) ( 0 ) ; + Projections* proj = new Projections( cube, trig, errorvec, triind ) ; + root = addTrian( root, proj, maxDepth ) ; + + delete proj->inherit ; + delete proj ; +} + + +UCHAR* Octree::addTrian( UCHAR* node, Projections* p, int height ) +{ + int i ; + int vertdiff[8][3] = {{0,0,0},{0,0,1},{0,1,-1},{0,0,1},{1,-1,-1},{0,0,1},{0,1,-1},{0,0,1}} ; + UCHAR boxmask = p->getBoxMask( ) ; + Projections* subp = new Projections( p ) ; + + int count = 0 ; + int tempdiff[3] = {0,0,0} ; + for ( i = 0 ; i < 8 ; i ++ ) + { + tempdiff[0] += vertdiff[i][0] ; + tempdiff[1] += vertdiff[i][1] ; + tempdiff[2] += vertdiff[i][2] ; + + /* Quick pruning using bounding box */ + if ( boxmask & ( 1 << i ) ) + { + subp->shift( tempdiff ) ; + tempdiff[0] = tempdiff[1] = tempdiff[2] = 0 ; + + /* Pruning using intersection test */ + if ( subp->isIntersecting() ) + // if ( subp->getIntersectionMasks( cedgemask, edgemask ) ) + { + if ( ! hasChild( node, i ) ) + { + if ( height == 1 ) + { + node = addLeafChild( node, i, count, createLeaf(0) ) ; + } + else + { + node = addInternalChild( node, i, count, createInternal(0) ) ; + } + } + UCHAR* chd = getChild( node, count ) ; + + if ( ! isLeaf( node, i ) ) + { + // setChild( node, count, addTrian ( chd, subp, height - 1, vertmask[i], edgemask ) ) ; + setChild( node, count, addTrian ( chd, subp, height - 1 ) ) ; + } + else + { + setChild( node, count, updateCell( chd, subp ) ) ; + } + } + } + + if ( hasChild( node, i ) ) + { + count ++ ; + } + } + + delete subp ; + return node ; +} + +UCHAR* Octree::updateCell( UCHAR* node, Projections* p ) +{ + int i ; + + // Edge connectivity + int mask[3] = { 0, 4, 8 } ; + int oldc = 0, newc = 0 ; + float offs[3] ; +#ifdef USE_HERMIT + float a[3], b[3], c[3] ; +#endif + + for ( i = 0 ; i < 3 ; i ++ ) + { + if ( ! getEdgeParity( node, mask[i] ) ) + { + if ( p->isIntersectingPrimary( i ) ) + { + // this->actualQuads ++ ; + setEdge( node, mask[i] ) ; + offs[ newc ] = p->getIntersectionPrimary( i ) ; +#ifdef USE_HERMIT + a[ newc ] = (float) p->inherit->norm[0] ; + b[ newc ] = (float) p->inherit->norm[1] ; + c[ newc ] = (float) p->inherit->norm[2] ; +#endif + newc ++ ; + } + } + else + { +#ifndef USE_HERMIT + offs[ newc ] = getEdgeOffset( node, oldc ) ; +#else + offs[ newc ] = getEdgeOffsetNormal( node, oldc, a[ newc ], b[ newc ], c[ newc ] ) ; +#endif + +// if ( p->isIntersectingPrimary( i ) ) + { + // dc_printf("Multiple intersections!\n") ; + +// setPatchEdge( node, i ) ; + } + + oldc ++ ; + newc ++ ; + } + } + + if ( newc > oldc ) + { + // New offsets added, update this node +#ifndef USE_HERMIT + node = updateEdgeOffsets( node, oldc, newc, offs ) ; +#else + node = updateEdgeOffsetsNormals( node, oldc, newc, offs, a, b, c ) ; +#endif + } + + + + return node ; +} + +void Octree::preparePrimalEdgesMask( UCHAR* node ) +{ + int count = 0 ; + for ( int i = 0 ; i < 8 ; i ++ ) + { + if ( hasChild( node, i ) ) + { + if ( isLeaf( node, i ) ) + { + createPrimalEdgesMask( getChild( node, count ) ) ; + } + else + { + preparePrimalEdgesMask( getChild( node, count ) ) ; + } + + count ++ ; + } + } +} + +void Octree::trace( ) +{ + int st[3] = { 0, 0, 0, } ; + this->numRings = 0 ; + this->totRingLengths = 0 ; + this->maxRingLength = 0 ; + + PathList* chdpath = NULL ; + this->root = trace( this->root, st, dimen, maxDepth, chdpath ) ; + + if ( chdpath != NULL ) + { + dc_printf("there are incomplete rings.\n") ; + printPaths( chdpath ) ; + }; +} + +UCHAR* Octree::trace( UCHAR* node, int* st, int len, int depth, PathList*& paths) +{ + UCHAR* newnode = node ; + len >>= 1 ; + PathList* chdpaths[ 8 ] ; + UCHAR* chd[ 8 ] ; + int nst[ 8 ][ 3 ] ; + int i, j ; + + // Get children paths + int chdleaf[ 8 ] ; + fillChildren( newnode, chd, chdleaf ) ; + + // int count = 0 ; + for ( i = 0 ; i < 8 ; i ++ ) + { + for ( j = 0 ; j < 3 ; j ++ ) + { + nst[ i ][ j ] = st[ j ] + len * vertmap[ i ][ j ] ; + } + + if ( chd[ i ] == NULL || isLeaf( node, i ) ) + { + chdpaths[ i ] = NULL ; + } + else + { + trace( chd[ i ], nst[i], len, depth - 1, chdpaths[ i ] ) ; + } + } + + // Get connectors on the faces + PathList* conn[ 12 ] ; + UCHAR* nf[2] ; + int lf[2] ; + int df[2] = { depth - 1, depth - 1 } ; + int* nstf[ 2 ]; + + fillChildren( newnode, chd, chdleaf ) ; + + for ( i = 0 ; i < 12 ; i ++ ) + { + int c[ 2 ] = { cellProcFaceMask[ i ][ 0 ], cellProcFaceMask[ i ][ 1 ] }; + + for ( int j = 0 ; j < 2 ; j ++ ) + { + lf[j] = chdleaf[ c[j] ] ; + nf[j] = chd[ c[j] ] ; + nstf[j] = nst[ c[j] ] ; + } + + conn[ i ] = NULL ; + + findPaths( nf, lf, df, nstf, depth - 1, cellProcFaceMask[ i ][ 2 ], conn[ i ] ) ; + + //if ( conn[i] ) + //{ + // printPath( conn[i] ) ; + //} + } + + // Connect paths + PathList* rings = NULL ; + combinePaths( chdpaths[0], chdpaths[1], conn[8], rings ) ; + combinePaths( chdpaths[2], chdpaths[3], conn[9], rings ) ; + combinePaths( chdpaths[4], chdpaths[5], conn[10], rings ) ; + combinePaths( chdpaths[6], chdpaths[7], conn[11], rings ) ; + + combinePaths( chdpaths[0], chdpaths[2], conn[4], rings ) ; + combinePaths( chdpaths[4], chdpaths[6], conn[5], rings ) ; + combinePaths( chdpaths[0], NULL, conn[6], rings ) ; + combinePaths( chdpaths[4], NULL, conn[7], rings ) ; + + combinePaths( chdpaths[0], chdpaths[4], conn[0], rings ) ; + combinePaths( chdpaths[0], NULL, conn[1], rings ) ; + combinePaths( chdpaths[0], NULL, conn[2], rings ) ; + combinePaths( chdpaths[0], NULL, conn[3], rings ) ; + + // By now, only chdpaths[0] and rings have contents + + // Process rings + if ( rings ) + { + // printPath( rings ) ; + + /* Let's count first */ + PathList* trings = rings ; + while ( trings ) + { + this->numRings ++ ; + this->totRingLengths += trings->length ; + if ( trings->length > this->maxRingLength ) + { + this->maxRingLength = trings->length ; + } + trings = trings->next ; + } + + // printPath( rings ) ; + newnode = patch( newnode, st, ( len << 1 ), rings ) ; + } + + // Return incomplete paths + paths = chdpaths[0] ; + return newnode ; +} + +void Octree::findPaths( UCHAR* node[2], int leaf[2], int depth[2], int* st[2], int maxdep, int dir, PathList*& paths ) +{ + if ( ! ( node[0] && node[1] ) ) + { + return ; + } + + if ( ! ( leaf[0] && leaf[1] ) ) + { + // Not at the bottom, recur + + // Fill children nodes + int i, j ; + UCHAR* chd[ 2 ][ 8 ] ; + int chdleaf[ 2 ][ 8 ] ; + int nst[ 2 ][ 8 ][ 3 ] ; + + for ( j = 0 ; j < 2 ; j ++ ) + { + if ( ! leaf[j] ) + { + fillChildren( node[j], chd[j], chdleaf[j] ) ; + + int len = ( dimen >> ( maxDepth - depth[j] + 1 ) ) ; + for ( i = 0 ; i < 8 ; i ++ ) + { + for ( int k = 0 ; k < 3 ; k ++ ) + { + nst[ j ][ i ][ k ] = st[ j ][ k ] + len * vertmap[ i ][ k ] ; + } + } + + } + } + + // 4 face calls + UCHAR* nf[2] ; + int df[2] ; + int lf[2] ; + int* nstf[2] ; + for ( i = 0 ; i < 4 ; i ++ ) + { + int c[2] = { faceProcFaceMask[ dir ][ i ][ 0 ], faceProcFaceMask[ dir ][ i ][ 1 ] }; + for ( int j = 0 ; j < 2 ; j ++ ) + { + if ( leaf[j] ) + { + lf[j] = leaf[j] ; + nf[j] = node[j] ; + df[j] = depth[j] ; + nstf[j] = st[j] ; + } + else + { + lf[j] = chdleaf[ j ][ c[j] ] ; + nf[j] = chd[ j ][ c[j] ] ; + df[j] = depth[j] - 1 ; + nstf[j] = nst[ j ][ c[j] ] ; + } + } + findPaths( nf, lf, df, nstf, maxdep - 1, faceProcFaceMask[ dir ][ i ][ 2 ], paths ) ; + } + + } + else + { + // At the bottom, check this face + int ind = ( depth[0] == maxdep ? 0 : 1 ) ; + int fcind = 2 * dir + ( 1 - ind ) ; + if ( getFaceParity( node[ ind ], fcind ) ) + { + // Add into path + PathElement* ele1 = new PathElement ; + PathElement* ele2 = new PathElement ; + + ele1->pos[0] = st[0][0] ; + ele1->pos[1] = st[0][1] ; + ele1->pos[2] = st[0][2] ; + + ele2->pos[0] = st[1][0] ; + ele2->pos[1] = st[1][1] ; + ele2->pos[2] = st[1][2] ; + + ele1->next = ele2 ; + ele2->next = NULL ; + + PathList* lst = new PathList ; + lst->head = ele1 ; + lst->tail = ele2 ; + lst->length = 2 ; + lst->next = paths ; + paths = lst ; + + // int l = ( dimen >> maxDepth ) ; + } + } + +} + +void Octree::combinePaths( PathList*& list1, PathList* list2, PathList* paths, PathList*& rings ) +{ + // Make new list of paths + PathList* nlist = NULL ; + + // Search for each connectors in paths + PathList* tpaths = paths ; + PathList* tlist, * pre ; + while ( tpaths ) + { + PathList* singlist = tpaths ; + PathList* templist ; + tpaths = tpaths->next ; + singlist->next = NULL ; + + // Look for hookup in list1 + tlist = list1 ; + pre = NULL ; + while ( tlist ) + { + if ( (templist = combineSinglePath( list1, pre, tlist, singlist, NULL, singlist )) != NULL ) + { + singlist = templist ; + continue ; + } + pre = tlist ; + tlist = tlist->next ; + } + + // Look for hookup in list2 + tlist = list2 ; + pre = NULL ; + while ( tlist ) + { + if ( (templist = combineSinglePath( list2, pre, tlist, singlist, NULL, singlist )) != NULL ) + { + singlist = templist ; + continue ; + } + pre = tlist ; + tlist = tlist->next ; + } + + // Look for hookup in nlist + tlist = nlist ; + pre = NULL ; + while ( tlist ) + { + if ( (templist = combineSinglePath( nlist, pre, tlist, singlist, NULL, singlist )) != NULL ) + { + singlist = templist ; + continue ; + } + pre = tlist ; + tlist = tlist->next ; + } + + // Add to nlist or rings + if ( isEqual( singlist->head, singlist->tail ) ) + { + PathElement* temp = singlist->head ; + singlist->head = temp->next ; + delete temp ; + singlist->length -- ; + singlist->tail->next = singlist->head ; + + singlist->next = rings ; + rings = singlist ; + } + else + { + singlist->next = nlist ; + nlist = singlist ; + } + + } + + // Append list2 and nlist to the end of list1 + tlist = list1 ; + if ( tlist != NULL ) + { + while ( tlist->next != NULL ) + { + tlist = tlist->next ; + } + tlist->next = list2 ; + } + else + { + tlist = list2 ; + list1 = list2 ; + } + + if ( tlist != NULL ) + { + while ( tlist->next != NULL ) + { + tlist = tlist->next ; + } + tlist->next = nlist ; + } + else + { + tlist = nlist ; + list1 = nlist ; + } + +} + +PathList* Octree::combineSinglePath( PathList*& head1, PathList* pre1, PathList*& list1, PathList*& head2, PathList* pre2, PathList*& list2 ) +{ + if ( isEqual( list1->head, list2->head ) || isEqual( list1->tail, list2->tail ) ) + { + // Reverse the list + if ( list1->length < list2->length ) + { + // Reverse list1 + PathElement* prev = list1->head ; + PathElement* next = prev->next ; + prev->next = NULL ; + while ( next != NULL ) + { + PathElement* tnext = next->next ; + next->next = prev ; + + prev = next ; + next = tnext ; + } + + list1->tail = list1->head ; + list1->head = prev ; + } + else + { + // Reverse list2 + PathElement* prev = list2->head ; + PathElement* next = prev->next ; + prev->next = NULL ; + while ( next != NULL ) + { + PathElement* tnext = next->next ; + next->next = prev ; + + prev = next ; + next = tnext ; + } + + list2->tail = list2->head ; + list2->head = prev ; + } + } + + if ( isEqual( list1->head, list2->tail ) ) + { + + // Easy case + PathElement* temp = list1->head->next ; + delete list1->head ; + list2->tail->next = temp ; + + PathList* nlist = new PathList ; + nlist->length = list1->length + list2->length - 1 ; + nlist->head = list2->head ; + nlist->tail = list1->tail ; + nlist->next = NULL ; + + deletePath( head1, pre1, list1 ) ; + deletePath( head2, pre2, list2 ) ; + + return nlist ; + } + else if ( isEqual( list1->tail, list2->head ) ) + { + // Easy case + PathElement* temp = list2->head->next ; + delete list2->head ; + list1->tail->next = temp ; + + PathList* nlist = new PathList ; + nlist->length = list1->length + list2->length - 1 ; + nlist->head = list1->head ; + nlist->tail = list2->tail ; + nlist->next = NULL ; + + deletePath( head1, pre1, list1 ) ; + deletePath( head2, pre2, list2 ) ; + + return nlist ; + } + + return NULL ; +} + +void Octree::deletePath( PathList*& head, PathList* pre, PathList*& curr ) +{ + PathList* temp = curr ; + curr = temp->next ; + delete temp ; + + if ( pre == NULL ) + { + head = curr ; + } + else + { + pre->next = curr ; + } +} + +void Octree::printElement( PathElement* ele ) +{ + if ( ele != NULL ) + { + dc_printf( " (%d %d %d)", ele->pos[0], ele->pos[1], ele->pos[2] ) ; + } +} + +void Octree::printPath( PathList* path ) +{ + PathElement* n = path->head; + int same = 0 ; + +#if DC_DEBUG + int len = ( dimen >> maxDepth ) ; +#endif + while ( n && ( same == 0 || n != path->head ) ) + { + same ++ ; + dc_printf( " (%d %d %d)", n->pos[0] / len, n->pos[1] / len, n->pos[2] / len ) ; + n = n->next ; + } + + if ( n == path->head ) + { + dc_printf(" Ring!\n") ; + } + else + { + dc_printf(" %p end!\n", n) ; + } +} + +void Octree::printPath( PathElement* path ) +{ + PathElement *n = path; + int same = 0 ; +#if DC_DEBUG + int len = ( dimen >> maxDepth ) ; +#endif + while ( n && ( same == 0 || n != path ) ) + { + same ++ ; + dc_printf( " (%d %d %d)", n->pos[0] / len, n->pos[1] / len, n->pos[2] / len ) ; + n = n->next ; + } + + if ( n == path ) + { + dc_printf(" Ring!\n") ; + } + else + { + dc_printf(" %p end!\n", n) ; + } + +} + + +void Octree::printPaths( PathList* path ) +{ + PathList* iter = path ; + int i = 0 ; + while ( iter != NULL ) + { + dc_printf("Path %d:\n", i) ; + printPath( iter ) ; + iter = iter->next ; + i ++ ; + } +} + +UCHAR* Octree::patch( UCHAR* node, int st[3], int len, PathList* rings ) +{ +#ifdef IN_DEBUG_MODE + dc_printf("Call to PATCH with rings: \n"); + printPaths( rings ) ; +#endif + + /* Do nothing but couting + PathList* tlist = rings ; + PathList* ttlist ; + PathElement* telem, * ttelem ; + while ( tlist!= NULL ) + { + // printPath( tlist ) ; + this->numRings ++ ; + this->totRingLengths += tlist->length ; + if ( tlist->length > this->maxRingLength ) + { + this->maxRingLength = tlist->length ; + } + ttlist = tlist ; + tlist = tlist->next ; + } + return node ; + */ + + + /* Pass onto separate calls in each direction */ + UCHAR* newnode = node ; + if ( len == mindimen ) + { + dc_printf("Error! should have no list by now.\n") ; + exit(0) ; + } + + // YZ plane + PathList* xlists[2] ; + newnode = patchSplit( newnode, st, len, rings, 0, xlists[0], xlists[1] ) ; + + // XZ plane + PathList* ylists[4] ; + newnode = patchSplit( newnode, st, len, xlists[0], 1, ylists[0], ylists[1] ) ; + newnode = patchSplit( newnode, st, len, xlists[1], 1, ylists[2], ylists[3] ) ; + + // XY plane + PathList* zlists[8] ; + newnode = patchSplit( newnode, st, len, ylists[0], 2, zlists[0], zlists[1] ) ; + newnode = patchSplit( newnode, st, len, ylists[1], 2, zlists[2], zlists[3] ) ; + newnode = patchSplit( newnode, st, len, ylists[2], 2, zlists[4], zlists[5] ) ; + newnode = patchSplit( newnode, st, len, ylists[3], 2, zlists[6], zlists[7] ) ; + + // Recur + len >>= 1 ; + int count = 0 ; + for ( int i = 0 ; i < 8 ; i ++ ) + { + if ( zlists[i] != NULL ) + { + int nori[3] = { + st[0] + len * vertmap[i][0] , + st[1] + len * vertmap[i][1] , + st[2] + len * vertmap[i][2] } ; + patch( getChild( newnode , count ), nori, len, zlists[i] ) ; + } + + if ( hasChild( newnode, i ) ) + { + count ++ ; + } + } +#ifdef IN_DEBUG_MODE + dc_printf("Return from PATCH\n") ; +#endif + return newnode ; + +} + + +UCHAR* Octree::patchSplit( UCHAR* node, int st[3], int len, PathList* rings, int dir, PathList*& nrings1, PathList*& nrings2 ) +{ +#ifdef IN_DEBUG_MODE + dc_printf("Call to PATCHSPLIT with direction %d and rings: \n", dir); + printPaths( rings ) ; +#endif + + UCHAR* newnode = node ; + nrings1 = NULL ; + nrings2 = NULL ; + PathList* tmp ; + while ( rings != NULL ) + { + // Process this ring + newnode = patchSplitSingle( newnode, st, len, rings->head, dir, nrings1, nrings2 ) ; + + // Delete this ring from the group + tmp = rings ; + rings = rings->next ; + delete tmp ; + } + +#ifdef IN_DEBUG_MODE + dc_printf("Return from PATCHSPLIT with \n"); + dc_printf("Rings gourp 1:\n") ; + printPaths( nrings1 ) ; + dc_printf("Rings group 2:\n") ; + printPaths( nrings2 ) ; +#endif + + return newnode ; +} + +UCHAR* Octree::patchSplitSingle( UCHAR* node, int st[3], int len, PathElement* head, int dir, PathList*& nrings1, PathList*& nrings2 ) +{ +#ifdef IN_DEBUG_MODE + dc_printf("Call to PATCHSPLITSINGLE with direction %d and path: \n", dir ); + printPath( head ) ; +#endif + + UCHAR* newnode = node ; + + if ( head == NULL ) + { +#ifdef IN_DEBUG_MODE + dc_printf("Return from PATCHSPLITSINGLE with head==NULL.\n") ; +#endif + return newnode; + } + else + { + // printPath( head ) ; + } + + // Walk along the ring to find pair of intersections + PathElement* pre1 = NULL ; + PathElement* pre2 = NULL ; + int side = findPair ( head, st[ dir ] + len / 2 , dir, pre1, pre2 ) ; + + /* + if ( pre1 == pre2 ) + { + int edgelen = ( dimen >> maxDepth ) ; + dc_printf("Location: %d %d %d Direction: %d Reso: %d\n", st[0]/edgelen, st[1]/edgelen, st[2]/edgelen, dir, len/edgelen) ; + printPath( head ) ; + exit( 0 ) ; + } + */ + + if ( side ) + { + // Entirely on one side + PathList* nring = new PathList( ) ; + nring->head = head ; + + if ( side == -1 ) + { + nring->next = nrings1 ; + nrings1 = nring ; + } + else + { + nring->next = nrings2 ; + nrings2 = nring ; + } + } + else + { + // Break into two parts + PathElement* nxt1 = pre1->next ; + PathElement* nxt2 = pre2->next ; + pre1->next = nxt2 ; + pre2->next = nxt1 ; + + newnode = connectFace( newnode, st, len, dir, pre1, pre2 ) ; + + if ( isEqual( pre1, pre1->next ) ) + { + if ( pre1 == pre1->next ) + { + delete pre1 ; + pre1 = NULL ; + } + else + { + PathElement* temp = pre1->next ; + pre1->next = temp->next ; + delete temp ; + } + } + if ( isEqual( pre2, pre2->next ) ) + { + if ( pre2 == pre2->next ) + { + delete pre2 ; + pre2 = NULL ; + } + else + { + PathElement* temp = pre2->next ; + pre2->next = temp->next ; + delete temp ; + } + } + + compressRing ( pre1 ) ; + compressRing ( pre2 ) ; + + // Recur + newnode = patchSplitSingle( newnode, st, len, pre1, dir, nrings1, nrings2 ) ; + newnode = patchSplitSingle( newnode, st, len, pre2, dir, nrings1, nrings2 ) ; + + } + +#ifdef IN_DEBUG_MODE + dc_printf("Return from PATCHSPLITSINGLE with \n"); + dc_printf("Rings gourp 1:\n") ; + printPaths( nrings1 ) ; + dc_printf("Rings group 2:\n") ; + printPaths( nrings2 ) ; +#endif + + return newnode ; +} + +UCHAR* Octree::connectFace( UCHAR* node, int st[3], int len, int dir, PathElement* f1, PathElement* f2 ) +{ +#ifdef IN_DEBUG_MODE + dc_printf("Call to CONNECTFACE with direction %d and length %d path: \n", dir, len ); + dc_printf("Path (low side): \n" ) ; + printPath( f1 ) ; +// checkPath( f1 ) ; + dc_printf("Path (high side): \n" ) ; + printPath( f2 ) ; +// checkPath( f2 ) ; +#endif + + UCHAR* newnode = node ; + + // Setup 2D + int pos = st[ dir ] + len / 2 ; + int xdir = ( dir + 1 ) % 3 ; + int ydir = ( dir + 2 ) % 3 ; + + // Use existing intersections on f1 and f2 + int x1, y1, x2, y2 ; + float p1, q1, p2, q2 ; + + getFacePoint( f2->next, dir, x1, y1, p1, q1 ) ; + getFacePoint( f2, dir, x2, y2, p2, q2 ) ; + + float dx = x2 + p2 - x1 - p1 ; + float dy = y2 + q2 - y1 - q1 ; + + // Do adapted Bresenham line drawing + float rx = p1, ry = q1 ; + int incx = 1, incy = 1 ; + int lx = x1, ly = y1 ; + int hx = x2, hy = y2 ; + int choice ; + if ( x2 < x1 ) + { + incx = -1 ; + rx = 1 - rx ; + lx = x2 ; + hx = x1 ; + } + if ( y2 < y1 ) + { + incy = -1 ; + ry = 1 - ry ; + ly = y2 ; + hy = y1 ; + } + + float sx = dx * incx ; + float sy = dy * incy ; + + int ori[3] ; + ori[ dir ] = pos / mindimen ; + ori[ xdir ] = x1 ; + ori[ ydir ] = y1 ; + int walkdir ; + int inc ; + float alpha ; + + PathElement* curEleN = f1 ; + PathElement* curEleP = f2->next ; + UCHAR *nodeN = NULL, *nodeP = NULL ; + UCHAR *curN = locateLeaf( newnode, len, f1->pos ) ; + UCHAR *curP = locateLeaf( newnode, len, f2->next->pos ) ; + if ( curN == NULL || curP == NULL ) + { + exit(0) ; + } + int stN[3], stP[3] ; + int lenN, lenP ; + + /* Unused code, leaving for posterity + + float stpt[3], edpt[3] ; + stpt[ dir ] = edpt[ dir ] = (float) pos ; + stpt[ xdir ] = ( x1 + p1 ) * mindimen ; + stpt[ ydir ] = ( y1 + q1 ) * mindimen ; + edpt[ xdir ] = ( x2 + p2 ) * mindimen ; + edpt[ ydir ] = ( y2 + q2 ) * mindimen ; + */ + while( ori[ xdir ] != x2 || ori[ ydir ] != y2 ) + { + int next ; + if ( sy * (1 - rx) > sx * (1 - ry) ) + { + choice = 1 ; + next = ori[ ydir ] + incy ; + if ( next < ly || next > hy ) + { + choice = 4 ; + next = ori[ xdir ] + incx ; + } + } + else + { + choice = 2 ; + next = ori[ xdir ] + incx ; + if ( next < lx || next > hx ) + { + choice = 3 ; + next = ori[ ydir ] + incy ; + } + } + + if ( choice & 1 ) + { + ori[ ydir ] = next ; + if ( choice == 1 ) + { + rx += ( sy == 0 ? 0 : (1 - ry) * sx / sy ) ; + ry = 0 ; + } + + walkdir = 2 ; + inc = incy ; + alpha = x2 < x1 ? 1 - rx : rx ; + } + else + { + ori[ xdir ] = next ; + if ( choice == 2 ) + { + ry += ( sx == 0 ? 0 : (1 - rx) * sy / sx ) ; + rx = 0 ; + } + + walkdir = 1 ; + inc = incx ; + alpha = y2 < y1 ? 1 - ry : ry ; + } + + + + // Get the exact location of the marcher + int nori[3] = { ori[0] * mindimen, ori[1] * mindimen, ori[2] * mindimen } ; + float spt[3] = { (float) nori[0], (float) nori[1], (float) nori[2] } ; + spt[ ( dir + ( 3 - walkdir ) ) % 3 ] += alpha * mindimen ; + if ( inc < 0 ) + { + spt[ ( dir + walkdir ) % 3 ] += mindimen ; + } + + // dc_printf("new x,y: %d %d\n", ori[ xdir ] / edgelen, ori[ ydir ] / edgelen ) ; + // dc_printf("nori: %d %d %d alpha: %f walkdir: %d\n", nori[0], nori[1], nori[2], alpha, walkdir ) ; + // dc_printf("%f %f %f\n", spt[0], spt[1], spt[2] ) ; + + // Locate the current cells on both sides + newnode = locateCell( newnode, st, len, nori, dir, 1, nodeN, stN, lenN ) ; + newnode = locateCell( newnode, st, len, nori, dir, 0, nodeP, stP, lenP ) ; + + updateParent( newnode, len, st ) ; + + int flag = 0 ; + // Add the cells to the rings and fill in the patch + PathElement* newEleN ; + if ( curEleN->pos[0] != stN[0] || curEleN->pos[1] != stN[1] || curEleN->pos[2] != stN[2] ) + { + if ( curEleN->next->pos[0] != stN[0] || curEleN->next->pos[1] != stN[1] || curEleN->next->pos[2] != stN[2] ) + { + newEleN = new PathElement ; + newEleN->next = curEleN->next ; + newEleN->pos[0] = stN[0] ; + newEleN->pos[1] = stN[1] ; + newEleN->pos[2] = stN[2] ; + + curEleN->next = newEleN ; + } + else + { + newEleN = curEleN->next ; + } + curN = patchAdjacent( newnode, len, curEleN->pos, curN, newEleN->pos, nodeN, walkdir, inc, dir, 1, alpha ) ; + + curEleN = newEleN ; + flag ++ ; + } + + PathElement* newEleP ; + if ( curEleP->pos[0] != stP[0] || curEleP->pos[1] != stP[1] || curEleP->pos[2] != stP[2] ) + { + if ( f2->pos[0] != stP[0] || f2->pos[1] != stP[1] || f2->pos[2] != stP[2] ) + { + newEleP = new PathElement ; + newEleP->next = curEleP ; + newEleP->pos[0] = stP[0] ; + newEleP->pos[1] = stP[1] ; + newEleP->pos[2] = stP[2] ; + + f2->next = newEleP ; + } + else + { + newEleP = f2 ; + } + curP = patchAdjacent( newnode, len, curEleP->pos, curP, newEleP->pos, nodeP, walkdir, inc, dir, 0, alpha ) ; + + + + curEleP = newEleP ; + flag ++ ; + } + + + /* + if ( flag == 0 ) + { + dc_printf("error: non-synchronized patching! at \n") ; + } + */ + } + +#ifdef IN_DEBUG_MODE + dc_printf("Return from CONNECTFACE with \n"); + dc_printf("Path (low side):\n") ; + printPath( f1 ) ; + checkPath( f1 ) ; + dc_printf("Path (high side):\n") ; + printPath( f2 ) ; + checkPath( f2 ) ; +#endif + + + return newnode ; +} + +UCHAR* Octree::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 ) +{ +#ifdef IN_DEBUG_MODE + dc_printf("Before patching.\n") ; + printInfo( st1 ) ; + printInfo( st2 ) ; + dc_printf("-----------------%d %d %d ; %d %d %d\n", st1[0], st2[1], st1[2], st2[0], st2[1], st2[2] ) ; +#endif + + // Get edge index on each leaf + int edgedir = ( dir + ( 3 - walkdir ) ) % 3 ; + int incdir = ( dir + walkdir ) % 3 ; + int ind1 = ( edgedir == 1 ? ( dir + 3 - edgedir ) % 3 - 1 : 2 - ( dir + 3 - edgedir ) % 3 ) ; + int ind2 = ( edgedir == 1 ? ( incdir + 3 - edgedir ) % 3 - 1 : 2 - ( incdir + 3 - edgedir ) % 3 ) ; + + int eind1 = ( ( edgedir << 2 ) | ( side << ind1 ) | ( ( inc > 0 ? 1 : 0 ) << ind2 ) ) ; + int eind2 = ( ( edgedir << 2 ) | ( side << ind1 ) | ( ( inc > 0 ? 0 : 1 ) << ind2 ) ) ; + +#ifdef IN_DEBUG_MODE + dc_printf("Index 1: %d Alpha 1: %f Index 2: %d Alpha 2: %f\n", eind1, alpha, eind2, alpha ) ; + /* + if ( alpha < 0 || alpha > 1 ) + { + dc_printf("Index 1: %d Alpha 1: %f Index 2: %d Alpha 2: %f\n", eind1, alpha, eind2, alpha ) ; + printInfo( st1 ) ; + printInfo( st2 ) ; + } + */ +#endif + + // Flip edge parity + UCHAR* nleaf1 = flipEdge( leaf1, eind1, alpha ) ; + UCHAR* nleaf2 = flipEdge( leaf2, eind2, alpha ) ; + + // Update parent link + updateParent( node, len, st1, nleaf1 ) ; + updateParent( node, len, st2, nleaf2 ) ; + // updateParent( nleaf1, mindimen, st1 ) ; + // updateParent( nleaf2, mindimen, st2 ) ; + + /* + float m[3] ; + dc_printf("Adding new point: %f %f %f\n", spt[0], spt[1], spt[2] ) ; + getMinimizer( leaf1, m ) ; + dc_printf("Cell %d now has minimizer %f %f %f\n", leaf1, m[0], m[1], m[2] ) ; + getMinimizer( leaf2, m ) ; + dc_printf("Cell %d now has minimizer %f %f %f\n", leaf2, m[0], m[1], m[2] ) ; + */ + +#ifdef IN_DEBUG_MODE + dc_printf("After patching.\n") ; + printInfo( st1 ) ; + printInfo( st2 ) ; +#endif + return nleaf2 ; +} + +UCHAR* Octree::locateCell( UCHAR* node, int st[3], int len, int ori[3], int dir, int side, UCHAR*& rleaf, int rst[3], int& rlen ) +{ +#ifdef IN_DEBUG_MODE + // dc_printf("Call to LOCATECELL with node ") ; + // printNode( node ) ; +#endif + UCHAR* newnode = node ; + int i ; + len >>= 1 ; + int ind = 0 ; + for ( i = 0 ; i < 3 ; i ++ ) + { + ind <<= 1 ; + if ( i == dir && side == 1 ) + { + ind |= ( ori[ i ] <= ( st[ i ] + len ) ? 0 : 1 ) ; + } + else + { + ind |= ( ori[ i ] < ( st[ i ] + len ) ? 0 : 1 ) ; + } + } + +#ifdef IN_DEBUG_MODE + // dc_printf("In LOCATECELL index of ori (%d %d %d) with dir %d side %d in st (%d %d %d, %d) is: %d\n", + // ori[0], ori[1], ori[2], dir, side, st[0], st[1], st[2], len, ind ) ; +#endif + + rst[0] = st[0] + vertmap[ ind ][ 0 ] * len ; + rst[1] = st[1] + vertmap[ ind ][ 1 ] * len ; + rst[2] = st[2] + vertmap[ ind ][ 2 ] * len ; + + if ( hasChild( newnode, ind ) ) + { + int count = getChildCount( newnode, ind ) ; + UCHAR* chd = getChild( newnode, count ) ; + if ( isLeaf( newnode, ind ) ) + { + rleaf = chd ; + rlen = len ; + } + else + { + // Recur + setChild( newnode, count, locateCell( chd, rst, len, ori, dir, side, rleaf, rst, rlen ) ) ; + } + } + else + { + // Create a new child here + if ( len == this->mindimen ) + { + UCHAR* chd = createLeaf( 0 ) ; + newnode = addChild( newnode, ind, chd, 1 ) ; + rleaf = chd ; + rlen = len ; + } + else + { + // Subdivide the empty cube + UCHAR* chd = createInternal( 0 ) ; + newnode = addChild( newnode, ind, locateCell( chd, rst, len, ori, dir, side, rleaf, rst, rlen ), 0 ) ; + } + } + +#ifdef IN_DEBUG_MODE + // dc_printf("Return from LOCATECELL with node ") ; + // printNode( newnode ) ; +#endif + return newnode ; +} + +void Octree::checkElement( PathElement* ele ) +{ + /* + if ( ele != NULL && locateLeafCheck( ele->pos ) != ele->node ) + { + dc_printf("Screwed! at pos: %d %d %d\n", ele->pos[0]>>minshift, ele->pos[1]>>minshift, ele->pos[2]>>minshift); + exit( 0 ) ; + } + */ +} + +void Octree::checkPath( PathElement* path ) +{ + PathElement *n = path ; + int same = 0 ; + while ( n && ( same == 0 || n != path ) ) + { + same ++ ; + checkElement( n ) ; + n = n->next ; + } + +} + +void Octree::testFacePoint( PathElement* e1, PathElement* e2 ) +{ + int i ; + PathElement * e = NULL ; + for ( i = 0 ; i < 3 ; i ++ ) + { + if ( e1->pos[i] != e2->pos[i] ) + { + if ( e1->pos[i] < e2->pos[i] ) + { + e = e2 ; + } + else + { + e = e1 ; + } + break ; + } + } + + int x, y ; + float p, q ; + dc_printf("Test.") ; + getFacePoint( e, i, x, y, p, q ) ; +} + +void Octree::getFacePoint( PathElement* leaf, int dir, int& x, int& y, float& p, float& q ) +{ + // Find average intersections + float avg[3] = {0, 0, 0} ; + float off[3] ; + int num = 0, num2 = 0 ; + + UCHAR* leafnode = locateLeaf( leaf->pos ) ; + for ( int i = 0 ; i < 4 ; i ++ ) + { + int edgeind = faceMap[ dir * 2 ][ i ] ; + int nst[3] ; + for ( int j = 0 ; j < 3 ; j ++ ) + { + nst[j] = leaf->pos[j] + mindimen * vertmap[ edgemap[ edgeind][ 0 ] ][ j ] ; + } + + if ( getEdgeIntersectionByIndex( nst, edgeind / 4, off, 1 ) ) + { + avg[0] += off[0] ; + avg[1] += off[1] ; + avg[2] += off[2] ; + num ++ ; + } + if ( getEdgeParity( leafnode, edgeind ) ) + { + num2 ++ ; + } + } + if ( num == 0 ) + { + dc_printf("Wrong! dir: %d pos: %d %d %d num: %d\n", dir, leaf->pos[0]>>minshift, leaf->pos[1]>>minshift, leaf->pos[2]>>minshift, num2); + avg[0] = (float) leaf->pos[0] ; + avg[1] = (float) leaf->pos[1] ; + avg[2] = (float) leaf->pos[2] ; + } + else + { + + avg[0] /= num ; + avg[1] /= num ; + avg[2] /= num ; + + //avg[0] = (float) leaf->pos[0]; + //avg[1] = (float) leaf->pos[1]; + //avg[2] = (float) leaf->pos[2]; + } + + int xdir = ( dir + 1 ) % 3 ; + int ydir = ( dir + 2 ) % 3 ; + + float xf = avg[ xdir ] ; + float yf = avg[ ydir ] ; + +#ifdef IN_DEBUG_MODE + // Is it outside? + // PathElement* leaf = leaf1->len < leaf2->len ? leaf1 : leaf2 ; + /* + float* m = ( leaf == leaf1 ? m1 : m2 ) ; + if ( xf < leaf->pos[ xdir ] || + yf < leaf->pos[ ydir ] || + xf > leaf->pos[ xdir ] + leaf->len || + yf > leaf->pos[ ydir ] + leaf->len) + { + dc_printf("Outside cube (%d %d %d), %d : %d %d %f %f\n", leaf->pos[ 0 ], leaf->pos[1], leaf->pos[2], leaf->len, + pos, dir, xf, yf) ; + + // For now, snap to cell + xf = m[ xdir ] ; + yf = m[ ydir ] ; + } + */ + + /* + if ( alpha < 0 || alpha > 1 || + xf < leaf->pos[xdir] || xf > leaf->pos[xdir] + leaf->len || + yf < leaf->pos[ydir] || yf > leaf->pos[ydir] + leaf->len ) + { + dc_printf("Alpha: %f Address: %d and %d\n", alpha, leaf1->node, leaf2->node ) ; + dc_printf("GETFACEPOINT result: (%d %d %d) %d min: (%f %f %f) ;(%d %d %d) %d min: (%f %f %f).\n", + leaf1->pos[0], leaf1->pos[1], leaf1->pos[2], leaf1->len, m1[0], m1[1], m1[2], + leaf2->pos[0], leaf2->pos[1], leaf2->pos[2], leaf2->len, m2[0], m2[1], m2[2]); + dc_printf("Face point at dir %d pos %d: %f %f\n", dir, pos, xf, yf ) ; + } + */ +#endif + + + // Get the integer and float part + x = ( ( leaf->pos[ xdir ] ) >> minshift ) ; + y = ( ( leaf->pos[ ydir ] ) >> minshift ) ; + + p = ( xf - leaf->pos[ xdir ] ) / mindimen ; + q = ( yf - leaf->pos[ ydir ] ) / mindimen ; + + +#ifdef IN_DEBUG_MODE + dc_printf("Face point at dir %d : %f %f\n", dir, xf, yf ) ; +#endif +} + +int Octree::findPair( PathElement* head, int pos, int dir, PathElement*& pre1, PathElement*& pre2 ) +{ + int side = getSide ( head, pos, dir ) ; + PathElement* cur = head ; + PathElement* anchor ; + PathElement* ppre1, *ppre2 ; + + // Start from this face, find a pair + anchor = cur ; + ppre1 = cur ; + cur = cur->next ; + while ( cur != anchor && ( getSide( cur, pos, dir ) == side ) ) + { + ppre1 = cur ; + cur = cur->next ; + } + if ( cur == anchor ) + { + // No pair found + return side ; + } + + side = getSide( cur, pos, dir ) ; + ppre2 = cur ; + cur = cur->next ; + while ( getSide( cur, pos, dir ) == side ) + { + ppre2 = cur ; + cur = cur->next ; + } + + + // Switch pre1 and pre2 if we start from the higher side + if ( side == -1 ) + { + cur = ppre1 ; + ppre1 = ppre2 ; + ppre2 = cur ; + } + + pre1 = ppre1 ; + pre2 = ppre2 ; + + return 0 ; +} + +int Octree::getSide( PathElement* e, int pos, int dir ) +{ + return ( e->pos[ dir ] < pos ? -1 : 1 ) ; +} + +int Octree::isEqual( PathElement* e1, PathElement* e2 ) +{ + return ( e1->pos[0] == e2->pos[0] && e1->pos[1] == e2->pos[1] && e1->pos[2] == e2->pos[2] ) ; +} + +void Octree::compressRing( PathElement*& ring ) +{ + if ( ring == NULL ) + { + return ; + } +#ifdef IN_DEBUG_MODE + dc_printf("Call to COMPRESSRING with path: \n" ); + printPath( ring ) ; +#endif + + PathElement* cur = ring->next->next ; + PathElement* pre = ring->next ; + PathElement* prepre = ring ; + PathElement* anchor = prepre ; + + do + { + while ( isEqual( cur, prepre ) ) + { + // Delete + if ( cur == prepre ) + { + // The ring has shrinked to a point + delete pre ; + delete cur ; + anchor = NULL ; + break ; + } + else + { + prepre->next = cur->next ; + delete pre ; + delete cur ; + pre = prepre->next ; + cur = pre->next ; + anchor = prepre ; + } + } + + if ( anchor == NULL ) + { + break ; + } + + prepre = pre ; + pre = cur ; + cur = cur->next ; + } while ( prepre != anchor ) ; + + ring = anchor ; + +#ifdef IN_DEBUG_MODE + dc_printf("Return from COMPRESSRING with path: \n" ); + printPath( ring ) ; +#endif +} + +void Octree::buildSigns( ) +{ + // First build a lookup table + // dc_printf("Building up look up table...\n") ; + int size = 1 << 12 ; + unsigned char table[ 1 << 12 ] ; + for ( int i = 0 ; i < size ; i ++ ) + { + table[i] = 0 ; + } + for ( int i = 0 ; i < 256 ; i ++ ) + { + int ind = 0 ; + for ( int j = 11 ; j >= 0 ; j -- ) + { + ind <<= 1 ; + if ( ( ( i >> edgemap[j][0] ) & 1 ) ^ ( ( i >> edgemap[j][1] ) & 1 ) ) + { + ind |= 1 ; + } + } + + table[ ind ] = i ; + } + + // Next, traverse the grid + int sg = 1 ; + int cube[8] ; + buildSigns( table, this->root, 0, sg, cube ) ; +} + +void Octree::buildSigns( unsigned char table[], UCHAR* node, int isLeaf, int sg, int rvalue[8] ) +{ + if ( node == NULL ) + { + for ( int i = 0 ; i < 8 ; i ++ ) + { + rvalue[i] = sg ; + } + return ; + } + + if ( isLeaf == 0 ) + { + // Internal node + UCHAR* chd[8] ; + int leaf[8] ; + fillChildren( node, chd, leaf ) ; + + // Get the signs at the corners of the first cube + rvalue[0] = sg ; + int oris[8] ; + buildSigns( table, chd[0], leaf[0], sg, oris ) ; + + // Get the rest + int cube[8] ; + for ( int i = 1 ; i < 8 ; i ++ ) + { + buildSigns( table, chd[i], leaf[i], oris[i], cube ) ; + rvalue[i] = cube[i] ; + } + + } + else + { + // Leaf node + generateSigns( node, table, sg ) ; + + for ( int i = 0 ; i < 8 ; i ++ ) + { + rvalue[i] = getSign( node, i ) ; + } + } +} + +void Octree::floodFill( ) +{ + // int threshold = (int) ((dimen/mindimen) * (dimen/mindimen) * 0.5f) ; + int st[3] = { 0, 0, 0 } ; + + // First, check for largest component + // size stored in -threshold + this->clearProcessBits( root, maxDepth ) ; + int threshold = this->floodFill( root, st, dimen, maxDepth, 0 ) ; + + // Next remove + dc_printf("Largest component: %d\n", threshold); + threshold *= thresh ; + dc_printf("Removing all components smaller than %d\n", threshold) ; + + int st2[3] = { 0, 0, 0 } ; + this->clearProcessBits( root, maxDepth ) ; + this->floodFill( root, st2, dimen, maxDepth, threshold ) ; + +} + +void Octree::clearProcessBits( UCHAR* node, int height ) +{ + int i; + + if ( height == 0 ) + { + // Leaf cell, + for ( i = 0 ; i < 12 ; i ++ ) + { + setOutProcess( node, i ) ; + } + } + else + { + // Internal cell, recur + int count = 0 ; + for ( i = 0 ; i < 8 ; i ++ ) + { + if ( hasChild( node, i ) ) + { + clearProcessBits( getChild( node, count ), height - 1 ) ; + count ++ ; + } + } + } +} + +/* +void Octree::floodFill( UCHAR* node, int st[3], int len, int height, int threshold ) +{ + int i, j; + + if ( height == 0 ) + { + // Leaf cell, + int par, inp ; + + // Test if the leaf has intersection edges + for ( i = 0 ; i < 12 ; i ++ ) + { + par = getEdgeParity( node, i ) ; + inp = isInProcess( node, i ) ; + + if ( par == 1 && inp == 0 ) + { + // Intersection edge, hasn't been processed + // Let's start filling + GridQueue* queue = new GridQueue() ; + int total = 1 ; + + // Set to in process + int mst[3] ; + mst[0] = st[0] + vertmap[edgemap[i][0]][0] * len ; + mst[1] = st[1] + vertmap[edgemap[i][0]][1] * len ; + mst[2] = st[2] + vertmap[edgemap[i][0]][2] * len; + int mdir = i / 4 ; + setInProcessAll( mst, mdir ) ; + + // Put this edge into queue + queue->pushQueue( mst, mdir ) ; + + // Queue processing + int nst[3], dir ; + while ( queue->popQueue( nst, dir ) == 1 ) + { + // dc_printf("nst: %d %d %d, dir: %d\n", nst[0]/mindimen, nst[1]/mindimen, nst[2]/mindimen, dir) ; + // locations + int stMask[3][3] = { + { 0, 0 - len, 0 - len }, + { 0 - len, 0, 0 - len }, + { 0 - len, 0 - len, 0 } + }; + int cst[2][3] ; + for ( j = 0 ; j < 3 ; j ++ ) + { + cst[0][j] = nst[j] ; + cst[1][j] = nst[j] + stMask[ dir ][ j ] ; + } + + // cells + UCHAR* cs[2] ; + for ( j = 0 ; j < 2 ; j ++ ) + { + cs[ j ] = locateLeaf( cst[j] ) ; + } + + // Middle sign + int s = getSign( cs[0], 0 ) ; + + // Masks + int fcCells[4] = {1,0,1,0}; + int fcEdges[3][4][3] = { + {{9,2,11},{8,1,10},{5,1,7},{4,2,6}}, + {{10,6,11},{8,5,9},{1,5,3},{0,6,2}}, + {{6,10,7},{4,9,5},{2,9,3},{0,10,1}} + }; + + // Search for neighboring connected intersection edges + for ( int find = 0 ; find < 4 ; find ++ ) + { + int cind = fcCells[find] ; + int eind, edge ; + if ( s == 0 ) + { + // Original order + for ( eind = 0 ; eind < 3 ; eind ++ ) + { + edge = fcEdges[dir][find][eind] ; + if ( getEdgeParity( cs[cind], edge ) == 1 ) + { + break ; + } + } + } + else + { + // Inverse order + for ( eind = 2 ; eind >= 0 ; eind -- ) + { + edge = fcEdges[dir][find][eind] ; + if ( getEdgeParity( cs[cind], edge ) == 1 ) + { + break ; + } + } + } + + if ( eind == 3 || eind == -1 ) + { + dc_printf("Wrong! this is not a consistent sign. %d\n", eind ); + } + else + { + int est[3] ; + est[0] = cst[cind][0] + vertmap[edgemap[edge][0]][0] * len ; + est[1] = cst[cind][1] + vertmap[edgemap[edge][0]][1] * len ; + est[2] = cst[cind][2] + vertmap[edgemap[edge][0]][2] * len ; + int edir = edge / 4 ; + + if ( isInProcess( cs[cind], edge ) == 0 ) + { + setInProcessAll( est, edir ) ; + queue->pushQueue( est, edir ) ; + // dc_printf("Pushed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; + total ++ ; + } + else + { + // dc_printf("Processed, not pushed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; + } + } + + } + + } + + dc_printf("Size of component: %d ", total) ; + + if ( total > threshold ) + { + dc_printf("Maintained.\n") ; + continue ; + } + dc_printf("Less then %d, removing...\n", threshold) ; + + // We have to remove this noise + + // Flip parity + // setOutProcessAll( mst, mdir ) ; + flipParityAll( mst, mdir ) ; + + // Put this edge into queue + queue->pushQueue( mst, mdir ) ; + + // Queue processing + while ( queue->popQueue( nst, dir ) == 1 ) + { + // dc_printf("nst: %d %d %d, dir: %d\n", nst[0]/mindimen, nst[1]/mindimen, nst[2]/mindimen, dir) ; + // locations + int stMask[3][3] = { + { 0, 0 - len, 0 - len }, + { 0 - len, 0, 0 - len }, + { 0 - len, 0 - len, 0 } + }; + int cst[2][3] ; + for ( j = 0 ; j < 3 ; j ++ ) + { + cst[0][j] = nst[j] ; + cst[1][j] = nst[j] + stMask[ dir ][ j ] ; + } + + // cells + UCHAR* cs[2] ; + for ( j = 0 ; j < 2 ; j ++ ) + { + cs[ j ] = locateLeaf( cst[j] ) ; + } + + // Middle sign + int s = getSign( cs[0], 0 ) ; + + // Masks + int fcCells[4] = {1,0,1,0}; + int fcEdges[3][4][3] = { + {{9,2,11},{8,1,10},{5,1,7},{4,2,6}}, + {{10,6,11},{8,5,9},{1,5,3},{0,6,2}}, + {{6,10,7},{4,9,5},{2,9,3},{0,10,1}} + }; + + // Search for neighboring connected intersection edges + for ( int find = 0 ; find < 4 ; find ++ ) + { + int cind = fcCells[find] ; + int eind, edge ; + if ( s == 0 ) + { + // Original order + for ( eind = 0 ; eind < 3 ; eind ++ ) + { + edge = fcEdges[dir][find][eind] ; + if ( isInProcess( cs[cind], edge ) == 1 ) + { + break ; + } + } + } + else + { + // Inverse order + for ( eind = 2 ; eind >= 0 ; eind -- ) + { + edge = fcEdges[dir][find][eind] ; + if ( isInProcess( cs[cind], edge ) == 1 ) + { + break ; + } + } + } + + if ( eind == 3 || eind == -1 ) + { + dc_printf("Wrong! this is not a consistent sign. %d\n", eind ); + } + else + { + int est[3] ; + est[0] = cst[cind][0] + vertmap[edgemap[edge][0]][0] * len ; + est[1] = cst[cind][1] + vertmap[edgemap[edge][0]][1] * len ; + est[2] = cst[cind][2] + vertmap[edgemap[edge][0]][2] * len ; + int edir = edge / 4 ; + + if ( getEdgeParity( cs[cind], edge ) == 1 ) + { + flipParityAll( est, edir ) ; + queue->pushQueue( est, edir ) ; + // dc_printf("Pushed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; + total ++ ; + } + else + { + // dc_printf("Processed, not pushed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; + } + } + + } + + } + + } + } + } + else + { + // Internal cell, recur + int count = 0 ; + len >>= 1 ; + for ( i = 0 ; i < 8 ; i ++ ) + { + if ( hasChild( node, i ) ) + { + int nst[3] ; + nst[0] = st[0] + vertmap[i][0] * len ; + nst[1] = st[1] + vertmap[i][1] * len ; + nst[2] = st[2] + vertmap[i][2] * len ; + + floodFill( getChild( node, count ), nst, len, height - 1, threshold ) ; + count ++ ; + } + } + } +} +*/ + +int Octree::floodFill( UCHAR* node, int st[3], int len, int height, int threshold ) +{ + int i, j; + int maxtotal = 0 ; + + if ( height == 0 ) + { + // Leaf cell, + int par, inp ; + + // Test if the leaf has intersection edges + for ( i = 0 ; i < 12 ; i ++ ) + { + par = getEdgeParity( node, i ) ; + inp = isInProcess( node, i ) ; + + if ( par == 1 && inp == 0 ) + { + // Intersection edge, hasn't been processed + // Let's start filling + GridQueue* queue = new GridQueue() ; + int total = 1 ; + + // Set to in process + int mst[3] ; + mst[0] = st[0] + vertmap[edgemap[i][0]][0] * len ; + mst[1] = st[1] + vertmap[edgemap[i][0]][1] * len ; + mst[2] = st[2] + vertmap[edgemap[i][0]][2] * len; + int mdir = i / 4 ; + setInProcessAll( mst, mdir ) ; + + // Put this edge into queue + queue->pushQueue( mst, mdir ) ; + + // Queue processing + int nst[3], dir ; + while ( queue->popQueue( nst, dir ) == 1 ) + { + // dc_printf("nst: %d %d %d, dir: %d\n", nst[0]/mindimen, nst[1]/mindimen, nst[2]/mindimen, dir) ; + // locations + int stMask[3][3] = { + { 0, 0 - len, 0 - len }, + { 0 - len, 0, 0 - len }, + { 0 - len, 0 - len, 0 } + }; + int cst[2][3] ; + for ( j = 0 ; j < 3 ; j ++ ) + { + cst[0][j] = nst[j] ; + cst[1][j] = nst[j] + stMask[ dir ][ j ] ; + } + + // cells + UCHAR* cs[2] ; + for ( j = 0 ; j < 2 ; j ++ ) + { + cs[ j ] = locateLeaf( cst[j] ) ; + } + + // Middle sign + int s = getSign( cs[0], 0 ) ; + + // Masks + int fcCells[4] = {1,0,1,0}; + int fcEdges[3][4][3] = { + {{9,2,11},{8,1,10},{5,1,7},{4,2,6}}, + {{10,6,11},{8,5,9},{1,5,3},{0,6,2}}, + {{6,10,7},{4,9,5},{2,9,3},{0,10,1}} + }; + + // Search for neighboring connected intersection edges + for ( int find = 0 ; find < 4 ; find ++ ) + { + int cind = fcCells[find] ; + int eind, edge ; + if ( s == 0 ) + { + // Original order + for ( eind = 0 ; eind < 3 ; eind ++ ) + { + edge = fcEdges[dir][find][eind] ; + if ( getEdgeParity( cs[cind], edge ) == 1 ) + { + break ; + } + } + } + else + { + // Inverse order + for ( eind = 2 ; eind >= 0 ; eind -- ) + { + edge = fcEdges[dir][find][eind] ; + if ( getEdgeParity( cs[cind], edge ) == 1 ) + { + break ; + } + } + } + + if ( eind == 3 || eind == -1 ) + { + dc_printf("Wrong! this is not a consistent sign. %d\n", eind ); + } + else + { + int est[3] ; + est[0] = cst[cind][0] + vertmap[edgemap[edge][0]][0] * len ; + est[1] = cst[cind][1] + vertmap[edgemap[edge][0]][1] * len ; + est[2] = cst[cind][2] + vertmap[edgemap[edge][0]][2] * len ; + int edir = edge / 4 ; + + if ( isInProcess( cs[cind], edge ) == 0 ) + { + setInProcessAll( est, edir ) ; + queue->pushQueue( est, edir ) ; + // dc_printf("Pushed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; + total ++ ; + } + else + { + // dc_printf("Processed, not pushed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; + } + } + + } + + } + + dc_printf("Size of component: %d ", total) ; + + if ( threshold == 0 ) + { + // Measuring stage + if ( total > maxtotal ) + { + maxtotal = total ; + } + dc_printf(".\n") ; + continue ; + } + + if ( total >= threshold ) + { + dc_printf("Maintained.\n") ; + continue ; + } + dc_printf("Less then %d, removing...\n", threshold) ; + + // We have to remove this noise + + // Flip parity + // setOutProcessAll( mst, mdir ) ; + flipParityAll( mst, mdir ) ; + + // Put this edge into queue + queue->pushQueue( mst, mdir ) ; + + // Queue processing + while ( queue->popQueue( nst, dir ) == 1 ) + { + // dc_printf("nst: %d %d %d, dir: %d\n", nst[0]/mindimen, nst[1]/mindimen, nst[2]/mindimen, dir) ; + // locations + int stMask[3][3] = { + { 0, 0 - len, 0 - len }, + { 0 - len, 0, 0 - len }, + { 0 - len, 0 - len, 0 } + }; + int cst[2][3] ; + for ( j = 0 ; j < 3 ; j ++ ) + { + cst[0][j] = nst[j] ; + cst[1][j] = nst[j] + stMask[ dir ][ j ] ; + } + + // cells + UCHAR* cs[2] ; + for ( j = 0 ; j < 2 ; j ++ ) + { + cs[ j ] = locateLeaf( cst[j] ) ; + } + + // Middle sign + int s = getSign( cs[0], 0 ) ; + + // Masks + int fcCells[4] = {1,0,1,0}; + int fcEdges[3][4][3] = { + {{9,2,11},{8,1,10},{5,1,7},{4,2,6}}, + {{10,6,11},{8,5,9},{1,5,3},{0,6,2}}, + {{6,10,7},{4,9,5},{2,9,3},{0,10,1}} + }; + + // Search for neighboring connected intersection edges + for ( int find = 0 ; find < 4 ; find ++ ) + { + int cind = fcCells[find] ; + int eind, edge ; + if ( s == 0 ) + { + // Original order + for ( eind = 0 ; eind < 3 ; eind ++ ) + { + edge = fcEdges[dir][find][eind] ; + if ( isInProcess( cs[cind], edge ) == 1 ) + { + break ; + } + } + } + else + { + // Inverse order + for ( eind = 2 ; eind >= 0 ; eind -- ) + { + edge = fcEdges[dir][find][eind] ; + if ( isInProcess( cs[cind], edge ) == 1 ) + { + break ; + } + } + } + + if ( eind == 3 || eind == -1 ) + { + dc_printf("Wrong! this is not a consistent sign. %d\n", eind ); + } + else + { + int est[3] ; + est[0] = cst[cind][0] + vertmap[edgemap[edge][0]][0] * len ; + est[1] = cst[cind][1] + vertmap[edgemap[edge][0]][1] * len ; + est[2] = cst[cind][2] + vertmap[edgemap[edge][0]][2] * len ; + int edir = edge / 4 ; + + if ( getEdgeParity( cs[cind], edge ) == 1 ) + { + flipParityAll( est, edir ) ; + queue->pushQueue( est, edir ) ; + // dc_printf("Pushed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; + total ++ ; + } + else + { + // dc_printf("Processed, not pushed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir) ; + } + } + + } + + } + + } + } + + } + else + { + // Internal cell, recur + int count = 0 ; + len >>= 1 ; + for ( i = 0 ; i < 8 ; i ++ ) + { + if ( hasChild( node, i ) ) + { + int nst[3] ; + nst[0] = st[0] + vertmap[i][0] * len ; + nst[1] = st[1] + vertmap[i][1] * len ; + nst[2] = st[2] + vertmap[i][2] * len ; + + int d = floodFill( getChild( node, count ), nst, len, height - 1, threshold ) ; + if ( d > maxtotal) + { + maxtotal = d ; + } + count ++ ; + } + } + } + + + return maxtotal ; + +} + +void Octree::writeOut() +{ + int numQuads = 0 ; + int numVertices = 0 ; + int numEdges = 0 ; +#ifdef USE_HERMIT + countIntersection( root, maxDepth, numQuads, numVertices, numEdges ) ; +#else + countIntersection( root, maxDepth, numQuads, numVertices ) ; + numEdges = numQuads * 3 / 2 ; +#endif + dc_printf("Vertices counted: %d Polys counted: %d \n", numVertices, numQuads ) ; + output_mesh = alloc_output(numVertices, numQuads); + int offset = 0 ; + int st[3] = { 0, 0, 0 } ; + + // First, output vertices + offset = 0 ; + actualVerts = 0 ; + actualQuads = 0 ; +#ifdef USE_HERMIT + generateMinimizer( root, st, dimen, maxDepth, offset ) ; + cellProcContour( this->root, 0, maxDepth ) ; + dc_printf("Vertices written: %d Quads written: %d \n", offset, actualQuads ) ; +#else + writeVertex( root, st, dimen, maxDepth, offset, out ) ; + writeQuad( root, st, dimen, maxDepth, out ) ; + dc_printf("Vertices written: %d Triangles written: %d \n", offset, actualQuads ) ; +#endif +} + +#if 0 +void Octree::writePLY( char* fname ) +{ + int numQuads = 0 ; + int numVertices = 0 ; + int numEdges = 0 ; +#ifdef USE_HERMIT + countIntersection( root, maxDepth, numQuads, numVertices, numEdges ) ; +#else + countIntersection( root, maxDepth, numQuads, numVertices ) ; + numEdges = numQuads * 3 / 2 ; +#endif + // int euler = numVertices + numQuads - numEdges ; + // int genus = ( 2 - euler ) / 2 ; + // dc_printf("%d vertices %d quads %d edges\n", numVertices, numQuads, numEdges ) ; + // dc_printf("Genus: %d Euler: %d\n", genus, euler ) ; + + FILE* fout = fopen ( fname, "wb" ) ; + dc_printf("Vertices counted: %d Polys counted: %d \n", numVertices, numQuads ) ; + PLYWriter::writeHeader( fout, numVertices, numQuads ) ; + int offset = 0 ; + int st[3] = { 0, 0, 0 } ; + + // First, output vertices + offset = 0 ; + actualVerts = 0 ; + actualQuads = 0 ; +#ifdef USE_HERMIT + generateMinimizer( root, st, dimen, maxDepth, offset, fout ) ; +#ifdef TESTMANIFOLD + testfout = fopen("test.txt", "w"); + fprintf(testfout, "{"); +#endif + cellProcContour( this->root, 0, maxDepth, fout ) ; +#ifdef TESTMANIFOLD + fprintf(testfout, "}"); + fclose( testfout ) ; +#endif + dc_printf("Vertices written: %d Quads written: %d \n", offset, actualQuads ) ; +#else + writeVertex( root, st, dimen, maxDepth, offset, fout ) ; + writeQuad( root, st, dimen, maxDepth, fout ) ; + dc_printf("Vertices written: %d Triangles written: %d \n", offset, actualQuads ) ; +#endif + + + fclose( fout ) ; +} +#endif + +void Octree::writeOctree( char* fname ) +{ + FILE* fout = fopen ( fname, "wb" ) ; + + int sized = ( 1 << maxDepth ) ; + fwrite( &sized, sizeof( int ), 1, fout ) ; + writeOctree( fout, root, maxDepth ) ; + dc_printf("Grid dimension: %d\n", sized ) ; + + + fclose( fout ) ; +} +void Octree::writeOctree( FILE* fout, UCHAR* node, int depth ) +{ + char type ; + if ( depth > 0 ) + { + type = 0 ; + fwrite( &type, sizeof( char ), 1, fout ) ; + + // Get sign at the center + char sg = (char) getSign( getChild( node, 0 ), depth - 1, 7 - getChildIndex( node, 0 ) ) ; + + int t = 0 ; + for ( int i = 0 ; i < 8 ; i ++ ) + { + if ( hasChild( node, i ) ) + { + writeOctree( fout, getChild( node, t ), depth - 1 ) ; + t ++ ; + } + else + { + type = 1 ; + fwrite( &type, sizeof( char ), 1, fout ) ; + fwrite( &sg, sizeof( char ), 1, fout ) ; + } + } + } + else + { + type = 2 ; + fwrite( &type, sizeof( char ), 1, fout ) ; + fwrite( &(node[2]), sizeof ( UCHAR ), 1, fout ); + } +} + +#ifdef USE_HERMIT +#if 0 +void Octree::writeOctreeGeom( char* fname ) +{ + FILE* fout = fopen ( fname, "wb" ) ; + + // Write header + char header[]="SOG.Format 1.0"; + int nlen = 128 - 4 * 4 - strlen(header) - 1 ; + char* header2 = new char[ nlen ]; + for ( int i = 0 ; i < nlen ; i ++ ) + { + header2[i] = '\0'; + } + fwrite( header, sizeof( char ), strlen(header) + 1, fout ) ; + fwrite( origin, sizeof( float ), 3, fout ) ; + fwrite( &range, sizeof( float ), 1, fout ) ; + fwrite( header2, sizeof( char ), nlen, fout ) ; + + + int sized = ( 1 << maxDepth ) ; + int st[3] = {0,0,0}; + fwrite( &sized, sizeof( int ), 1, fout ) ; + + writeOctreeGeom( fout, root, st, dimen, maxDepth ) ; + dc_printf("Grid dimension: %d\n", sized ) ; + + + fclose( fout ) ; +} +#endif +void Octree::writeOctreeGeom( FILE* fout, UCHAR* node, int st[3], int len, int depth ) +{ + char type ; + if ( depth > 0 ) + { + type = 0 ; + fwrite( &type, sizeof( char ), 1, fout ) ; + + // Get sign at the center + char sg = (char) getSign( getChild( node, 0 ), depth - 1, 7 - getChildIndex( node, 0 ) ) ; + + int t = 0 ; + len >>= 1 ; + for ( int i = 0 ; i < 8 ; i ++ ) + { + if ( hasChild( node, i ) ) + { + int nst[3] ; + nst[0] = st[0] + vertmap[i][0] * len ; + nst[1] = st[1] + vertmap[i][1] * len ; + nst[2] = st[2] + vertmap[i][2] * len ; + writeOctreeGeom( fout, getChild( node, t ), nst, len, depth - 1 ) ; + t ++ ; + } + else + { + type = 1 ; + fwrite( &type, sizeof( char ), 1, fout ) ; + fwrite( &sg, sizeof( char ), 1, fout ) ; + } + } + } + else + { + type = 2 ; + fwrite( &type, sizeof( char ), 1, fout ) ; + fwrite( &(node[2]), sizeof ( UCHAR ), 1, fout ); + + // Compute minimizer + // First, find minimizer + float rvalue[3] ; + rvalue[0] = (float) st[0] + len / 2 ; + rvalue[1] = (float) st[1] + len / 2 ; + rvalue[2] = (float) st[2] + len / 2 ; + computeMinimizer( node, st, len, rvalue ) ; + + // Update + // float flen = len * range / dimen ; + for ( int j = 0 ; j < 3 ; j ++ ) + { + rvalue[ j ] = rvalue[ j ] * range / dimen + origin[ j ] ; + } + + fwrite( rvalue, sizeof ( float ), 3, fout ); + } +} +#endif + +#ifdef USE_HERMIT +void Octree::writeDCF( char* fname ) +{ + FILE* fout = fopen ( fname, "wb" ) ; + + // Writing out version + char version[10] = "multisign"; + fwrite ( &version, sizeof ( char ), 10, fout ); + + // Writing out size + int sized = ( 1 << maxDepth ) ; + fwrite( &sized, sizeof( int ), 1, fout ) ; + fwrite( &sized, sizeof( int ), 1, fout ) ; + fwrite( &sized, sizeof( int ), 1, fout ) ; + + int st[3] = {0, 0, 0} ; + writeDCF( fout, root, maxDepth, st, dimen ) ; + + dc_printf("Grid dimension: %d\n", sized ) ; + fclose( fout ) ; +} + +void Octree::writeDCF( FILE* fout, UCHAR* node, int height, int st[3], int len ) +{ + nodetype type ; + if ( height > 0 ) + { + type = 0 ; + len >>= 1 ; + fwrite( &type, sizeof( nodetype ), 1, fout ) ; + + // Get sign at the center + signtype sg = 1 - (signtype) getSign( getChild( node, 0 ), height - 1, 7 - getChildIndex( node, 0 ) ) ; + + int t = 0 ; + for ( int i = 0 ; i < 8 ; i ++ ) + { + if ( hasChild( node, i ) ) + { + int nst[3] ; + nst[0] = st[0] + vertmap[i][0] * len ; + nst[1] = st[1] + vertmap[i][1] * len ; + nst[2] = st[2] + vertmap[i][2] * len ; + + + writeDCF( fout, getChild( node, t ), height - 1, nst, len ) ; + t ++ ; + } + else + { + type = 1 ; + fwrite( &type, sizeof( nodetype ), 1, fout ) ; + fwrite ( &(sg), sizeof ( signtype ), 1, fout ); + } + } + } + else + { + type = 2 ; + fwrite( &type, sizeof( nodetype ), 1, fout ) ; + + // Write signs + signtype sgn[8] ; + for ( int i = 0 ; i < 8 ; i ++ ) + { + sgn[ i ] = 1 - (signtype) getSign( node, i ) ; + } + fwrite (sgn, sizeof (signtype), 8, fout ); + + // Write edge data + float pts[12], norms[12][3] ; + int parity[12] ; + fillEdgeOffsetsNormals( node, st, len, pts, norms, parity ) ; + + numtype zero = 0, one = 1 ; + for ( int i = 0 ; i < 12 ; i ++ ) + { + int par = getEdgeParity( node, i ) ; + // Let's check first + if ( par ) + { + if ( sgn[ edgemap[i][0] ] == sgn[ edgemap[i][1] ] ) + { + dc_printf("Wrong! Parity: %d Sign: %d %d\n", parity[i], sgn[ edgemap[i][0] ], sgn[ edgemap[i][1] ]); + exit(0) ; + } + if ( parity[ i ] == 0 ) + { + dc_printf("Wrong! No intersection found.\n"); + exit(0) ; + } + fwrite( &one, sizeof ( numtype ) , 1, fout ) ; + fwrite( &(pts[i]), sizeof( float ), 1, fout ) ; + fwrite( norms[i], sizeof( float ), 3, fout ) ; + + } + else + { + if ( sgn[ edgemap[i][0] ] != sgn[ edgemap[i][1] ] ) + { + dc_printf("Wrong! Parity: %d Sign: %d %d\n", parity[i], sgn[ edgemap[i][0] ], sgn[ edgemap[i][1] ]); + exit(0) ; + } + fwrite ( &zero, sizeof ( numtype ) , 1, fout ); + } + } + } +} +#endif + + +void Octree::writeOpenEdges( FILE* fout ) +{ + // Total number of rings + fprintf( fout, "%d\n", numRings ) ; + dc_printf("Number of rings to write: %d\n", numRings) ; + + // Write each ring + PathList* tlist = ringList ; + for ( int i = 0 ; i < numRings ; i ++ ) + { + fprintf(fout, "%d\n", tlist->length) ; + // dc_printf("Ring length: %d\n", tlist->length ) ; + PathElement* cur = tlist->head ; + for ( int j = 0 ; j < tlist->length ; j ++ ) + { + float cent[3] ; + float flen = mindimen * range / dimen ; + for ( int k = 0 ; k < 3 ; k ++ ) + { + cent[ k ] = cur->pos[ k ] * range / dimen + origin[ k ] + flen / 2 ; + } + fprintf(fout, "%f %f %f\n", cent[0], cent[1], cent[2]) ; + cur = cur->next ; + } + + tlist = tlist->next ; + } +} + +#ifndef USE_HERMIT +void Octree::countIntersection( UCHAR* node, int height, int& nquad, int& nvert ) +{ + if ( height > 0 ) + { + int total = getNumChildren( node ) ; + for ( int i = 0 ; i < total ; i ++ ) + { + countIntersection( getChild( node, i ), height - 1, nquad, nvert ) ; + } + } + else + { + int mask = getSignMask( node ) ; + nvert += getNumEdges2( node ) ; + nquad += cubes->getNumTriangle( mask ) ; + + } +} + +void Octree::writeVertex( UCHAR* node, int st[3], int len, int height, int& offset, FILE* fout ) +{ + int i ; + + if ( height == 0 ) + { + // Leaf cell, generate + int emap[] = { 0, 4, 8 } ; + for ( int i = 0 ; i < 3 ; i ++ ) + { + if ( getEdgeParity( node, emap[i] ) ) + { + // Get intersection location + int count = getEdgeCount( node, i ) ; + float off = getEdgeOffset( node, count ) ; + + float rvalue[3] ; + rvalue[0] = (float) st[0] ; + rvalue[1] = (float) st[1] ; + rvalue[2] = (float) st[2] ; + rvalue[i] += off * mindimen ; + + // Update + float fnst[3] ; + float flen = len * range / dimen ; + for ( int j = 0 ; j < 3 ; j ++ ) + { + rvalue[ j ] = rvalue[ j ] * range / dimen + origin[ j ] ; + fnst[ j ] = st[ j ] * range / dimen + origin[ j ] ; + } + + if ( this->outType == 0 ) + { + fprintf( fout, "%f %f %f\n", rvalue[0], rvalue[1], rvalue[2] ) ; + } + else if ( this->outType == 1 ) + { + PLYWriter::writeVertex( fout, rvalue ) ; + } + + // Store the index + setEdgeIntersectionIndex( node, count, offset ) ; + offset ++ ; + } + } + + } + else + { + // Internal cell, recur + int count = 0 ; + len >>= 1 ; + for ( i = 0 ; i < 8 ; i ++ ) + { + if ( hasChild( node, i ) ) + { + int nst[3] ; + nst[0] = st[0] + vertmap[i][0] * len ; + nst[1] = st[1] + vertmap[i][1] * len ; + nst[2] = st[2] + vertmap[i][2] * len ; + + writeVertex( getChild( node, count ), nst, len, height - 1, offset, fout ) ; + count ++ ; + } + } + } +} + +void Octree::writeQuad( UCHAR* node, int st[3], int len, int height, FILE* fout ) +{ + int i ; + if ( height == 0 ) + { + int mask = getSignMask( node ) ; + int num = cubes->getNumTriangle( mask ) ; + int indices[12] ; + fillEdgeIntersectionIndices( node, st, len, indices ) ; + int einds[3], ind[3] ; + + //int flag1 = 0 ; + //int flag2 = 0 ; + for ( i = 0 ; i < num ; i ++ ) + { + int color = 0 ; + cubes->getTriangle( mask, i, einds ) ; + // dc_printf("(%d %d %d) ", einds[0], einds[1], einds[2] ) ; + + for ( int j = 0 ; j < 3 ; j ++ ) + { + ind[j] = indices[ einds[j] ] ; + /* + if ( ind[j] == 78381 ) + { + flag1 = 1 ; + } + if ( ind[j] == 78384 ) + { + flag2 = 1 ; + } + */ + } + + if ( this->outType == 0 ) + { + // OFF + int numpoly = ( color ? -3 : 3 ) ; + fprintf(fout, "%d %d %d %d\n", numpoly, ind[0], ind[1], ind[2] ) ; + actualQuads ++ ; + } + else if ( this->outType == 1 ) + { + // PLY + PLYWriter::writeFace( fout, 3, ind ) ; + actualQuads ++ ; + } + } + + /* + if (flag1 && flag2) + { + dc_printf("%d\n", mask); + cubes->printTriangles( mask ) ; + } + */ + } + else + { + // Internal cell, recur + int count = 0 ; + len >>= 1 ; + for ( i = 0 ; i < 8 ; i ++ ) + { + if ( hasChild( node, i ) ) + { + int nst[3] ; + nst[0] = st[0] + vertmap[i][0] * len ; + nst[1] = st[1] + vertmap[i][1] * len ; + nst[2] = st[2] + vertmap[i][2] * len ; + + writeQuad( getChild( node, count ), nst, len, height - 1, fout ) ; + count ++ ; + } + } + } +} + +#endif + + +#ifdef USE_HERMIT +void Octree::countIntersection( UCHAR* node, int height, int& nedge, int& ncell, int& nface ) +{ + if ( height > 0 ) + { + int total = getNumChildren( node ) ; + for ( int i = 0 ; i < total ; i ++ ) + { + countIntersection( getChild( node, i ), height - 1, nedge, ncell, nface ) ; + } + } + else + { + nedge += getNumEdges2( node ) ; + + int smask = getSignMask( node ) ; + + if(use_manifold) + { + int comps = manifold_table[ smask ].comps ; + ncell += comps ; + } + else { + if ( smask > 0 && smask < 255 ) + { + ncell ++ ; + } + } + + for ( int i = 0 ; i < 3 ; i ++ ) + { + if ( getFaceEdgeNum( node, i * 2 ) ) + { + nface ++ ; + } + } + } +} + +/* from http://eigen.tuxfamily.org/bz/show_bug.cgi?id=257 */ +template<typename _Matrix_Type_> +void pseudoInverse(const _Matrix_Type_ &a, + _Matrix_Type_ &result, + double epsilon = std::numeric_limits<typename _Matrix_Type_::Scalar>::epsilon()) +{ + Eigen::JacobiSVD< _Matrix_Type_ > svd = a.jacobiSvd(Eigen::ComputeFullU | + Eigen::ComputeFullV); + + typename _Matrix_Type_::Scalar tolerance = epsilon * std::max(a.cols(), + a.rows()) * + svd.singularValues().array().abs().maxCoeff(); + + result = svd.matrixV() * + _Matrix_Type_((svd.singularValues().array().abs() > + tolerance).select(svd.singularValues(). + array().inverse(), 0)).asDiagonal() * + svd.matrixU().adjoint(); +} + +void solve_least_squares(const float halfA[], const float b[], + const float midpoint[], float rvalue[]) +{ + /* calculate pseudo-inverse */ + Eigen::MatrixXf A(3, 3), pinv(3, 3); + A << halfA[0], halfA[1], halfA[2], + halfA[1], halfA[3], halfA[4], + halfA[2], halfA[4], halfA[5]; + pseudoInverse(A, pinv); + + Eigen::Vector3f b2(b), mp(midpoint), result; + b2 = b2 + A * -mp; + result = pinv * b2 + mp; + + for(int i = 0; i < 3; i++) + rvalue[i] = result(i); +} + +void minimize(float rvalue[3], float mp[3], const float pts[12][3], + const float norms[12][3], const int parity[12]) +{ + float ata[6] = { 0, 0, 0, 0, 0, 0 }; + float atb[3] = { 0, 0, 0 } ; + int ec = 0 ; + + for ( int i = 0 ; i < 12 ; i ++ ) + { + // if ( getEdgeParity( leaf, i) ) + if ( parity[ i ] ) + { + const float* norm = norms[i] ; + const float* p = pts[i] ; + + // QEF + ata[ 0 ] += (float) ( norm[ 0 ] * norm[ 0 ] ); + ata[ 1 ] += (float) ( norm[ 0 ] * norm[ 1 ] ); + ata[ 2 ] += (float) ( norm[ 0 ] * norm[ 2 ] ); + ata[ 3 ] += (float) ( norm[ 1 ] * norm[ 1 ] ); + ata[ 4 ] += (float) ( norm[ 1 ] * norm[ 2 ] ); + ata[ 5 ] += (float) ( norm[ 2 ] * norm[ 2 ] ); + + double pn = p[0] * norm[0] + p[1] * norm[1] + p[2] * norm[2] ; + + atb[ 0 ] += (float) ( norm[ 0 ] * pn ) ; + atb[ 1 ] += (float) ( norm[ 1 ] * pn ) ; + atb[ 2 ] += (float) ( norm[ 2 ] * pn ) ; + + // Minimizer + mp[0] += p[0] ; + mp[1] += p[1] ; + mp[2] += p[2] ; + + ec ++ ; + } + } + + if ( ec == 0 ) + { + return ; + } + mp[0] /= ec ; + mp[1] /= ec ; + mp[2] /= ec ; + + // Solve least squares + solve_least_squares(ata, atb, mp, rvalue); +} + +void Octree::computeMinimizer( UCHAR* leaf, int st[3], int len, float rvalue[3] ) +{ + // First, gather all edge intersections + float pts[12][3], norms[12][3] ; + // fillEdgeIntersections( leaf, st, len, pts, norms ) ; + int parity[12] ; + fillEdgeIntersections( leaf, st, len, pts, norms, parity ) ; + + // Next, construct QEF and minimizer + float mp[3] = {0, 0, 0}; + minimize(rvalue, mp, pts, norms, parity); + + /* Restraining the location of the minimizer */ + float nh1 = hermite_num * len ; + float nh2 = ( 1 + hermite_num ) * len ; + if((mode == DUALCON_MASS_POINT || mode == DUALCON_CENTROID) || + ( rvalue[0] < st[0] - nh1 || rvalue[1] < st[1] - nh1 || rvalue[2] < st[2] - nh1 || + rvalue[0] > st[0] + nh2 || rvalue[1] > st[1] + nh2 || rvalue[2] > st[2] + nh2 )) + { + if(mode == DUALCON_CENTROID) { + // Use centroids + rvalue[0] = (float) st[0] + len / 2 ; + rvalue[1] = (float) st[1] + len / 2 ; + rvalue[2] = (float) st[2] + len / 2 ; + } + else { + // Use mass point instead + rvalue[0] = mp[0] ; + rvalue[1] = mp[1] ; + rvalue[2] = mp[2] ; + } + } +} + +void Octree::generateMinimizer( UCHAR* node, int st[3], int len, int height, int& offset ) +{ + int i, j ; + + if ( height == 0 ) + { + // Leaf cell, generate + + // First, find minimizer + float rvalue[3] ; + rvalue[0] = (float) st[0] + len / 2 ; + rvalue[1] = (float) st[1] + len / 2 ; + rvalue[2] = (float) st[2] + len / 2 ; + computeMinimizer( node, st, len, rvalue ) ; + + // Update + //float fnst[3] ; + for ( j = 0 ; j < 3 ; j ++ ) + { + rvalue[ j ] = rvalue[ j ] * range / dimen + origin[ j ] ; + //fnst[ j ] = st[ j ] * range / dimen + origin[ j ] ; + } + + int mult = 0, smask = getSignMask( node ) ; + + if(use_manifold) + { + mult = manifold_table[ smask ].comps ; + } + else + { + if ( smask > 0 && smask < 255 ) + { + mult = 1 ; + } + } + + for ( j = 0 ; j < mult ; j ++ ) + { + add_vert(output_mesh, rvalue); + } + + // Store the index + setMinimizerIndex( node, offset ) ; + + offset += mult ; + } + else + { + // Internal cell, recur + int count = 0 ; + len >>= 1 ; + for ( i = 0 ; i < 8 ; i ++ ) + { + if ( hasChild( node, i ) ) + { + int nst[3] ; + nst[0] = st[0] + vertmap[i][0] * len ; + nst[1] = st[1] + vertmap[i][1] * len ; + nst[2] = st[2] + vertmap[i][2] * len ; + + generateMinimizer( getChild( node, count ), nst, len, height - 1, offset ) ; + count ++ ; + } + } + } +} + +void Octree::processEdgeWrite( UCHAR* node[4], int depth[4], int maxdep, int dir ) +{ + //int color = 0 ; + + int i = 3 ; + { + if ( getEdgeParity( node[i], processEdgeMask[dir][i] ) ) + { + int flip = 0 ; + int edgeind = processEdgeMask[dir][i] ; + if ( getSign( node[i], edgemap[ edgeind ][ 1 ] ) > 0 ) + { + flip = 1 ; + } + + int num = 0 ; + { + int ind[8]; + if(use_manifold) + { + /* Deprecated + int ind[4] = { + getMinimizerIndex( node[0], processEdgeMask[dir][0] ), + getMinimizerIndex( node[1], processEdgeMask[dir][1] ), + getMinimizerIndex( node[3], processEdgeMask[dir][3] ), + getMinimizerIndex( node[2], processEdgeMask[dir][2] ) + } ; + num = 4 ; + */ + int vind[2] ; + int seq[4] = {0,1,3,2}; + for ( int k = 0 ; k < 4 ; k ++ ) + { + getMinimizerIndices( node[seq[k]], processEdgeMask[dir][seq[k]], vind ) ; + ind[num] = vind[0] ; + num ++ ; + + if ( vind[1] != -1 ) + { + ind[num] = vind[1] ; + num ++ ; + if ( flip == 0 ) + { + ind[num-1] = vind[0] ; + ind[num-2] = vind[1] ; + } + } + } +#ifdef TESTMANIFOLD + if ( num != 4 ) + { + dc_printf("Polygon: %d\n", num); + } + for ( k = 0 ; k < num ; k ++ ) + { + fprintf(testfout, "{%d,%d},", ind[k], ind[(k+1)%num] ); + } +#endif + + /* we don't use the manifold option, but if it is + ever enabled again note that it can output + non-quads */ + } + else { + if(flip) { + ind[0] = getMinimizerIndex( node[2] ); + ind[1] = getMinimizerIndex( node[3] ); + ind[2] = getMinimizerIndex( node[1] ); + ind[3] = getMinimizerIndex( node[0] ); + } + else { + ind[0] = getMinimizerIndex( node[0] ); + ind[1] = getMinimizerIndex( node[1] ); + ind[2] = getMinimizerIndex( node[3] ); + ind[3] = getMinimizerIndex( node[2] ); + } + + add_quad(output_mesh, ind); + } + /* + if ( this->outType == 0 ) + { + // OFF + + num = ( color ? -num : num ) ; + + fprintf(fout, "%d ", num); + + if ( flip ) + { + for ( int k = num - 1 ; k >= 0 ; k -- ) + { + fprintf(fout, "%d ", ind[k] ) ; + } + } + else + { + for ( int k = 0 ; k < num ; k ++ ) + { + fprintf(fout, "%d ", ind[k] ) ; + } + } + + fprintf( fout, "\n") ; + + actualQuads ++ ; + } + else if ( this->outType == 1 ) + { + // PLY + + if ( flip ) + { + int tind[8]; + for ( int k = num - 1 ; k >= 0 ; k -- ) + { + tind[k] = ind[num-1-k] ; + } + // PLYWriter::writeFace( fout, num, tind ) ; + } + else + { + // PLYWriter::writeFace( fout, num, ind ) ; + } + + actualQuads ++ ; + }*/ + } + return ; + } + else + { + return ; + } + } +} + + +void Octree::edgeProcContour( UCHAR* node[4], int leaf[4], int depth[4], int maxdep, int dir ) +{ + if ( ! ( node[0] && node[1] && node[2] && node[3] ) ) + { + return ; + } + if ( leaf[0] && leaf[1] && leaf[2] && leaf[3] ) + { + processEdgeWrite( node, depth, maxdep, dir ) ; + } + else + { + int i, j ; + UCHAR* chd[ 4 ][ 8 ] ; + for ( j = 0 ; j < 4 ; j ++ ) + { + for ( i = 0 ; i < 8 ; i ++ ) + { + chd[ j ][ i ] = ((!leaf[j]) && hasChild( node[j], i ) )? getChild( node[j], getChildCount( node[j], i ) ) : NULL ; + } + } + + // 2 edge calls + UCHAR* ne[4] ; + int le[4] ; + int de[4] ; + for ( i = 0 ; i < 2 ; i ++ ) + { + int c[ 4 ] = { edgeProcEdgeMask[ dir ][ i ][ 0 ], + edgeProcEdgeMask[ dir ][ i ][ 1 ], + edgeProcEdgeMask[ dir ][ i ][ 2 ], + edgeProcEdgeMask[ dir ][ i ][ 3 ] } ; + + for ( int j = 0 ; j < 4 ; j ++ ) + { + if ( leaf[j] ) + { + le[j] = leaf[j] ; + ne[j] = node[j] ; + de[j] = depth[j] ; + } + else + { + le[j] = isLeaf( node[j], c[j] ) ; + ne[j] = chd[ j ][ c[j] ] ; + de[j] = depth[j] - 1 ; + } + } + + edgeProcContour( ne, le, de, maxdep - 1, edgeProcEdgeMask[ dir ][ i ][ 4 ] ) ; + } + + } +} + +void Octree::faceProcContour( UCHAR* node[2], int leaf[2], int depth[2], int maxdep, int dir ) +{ + if ( ! ( node[0] && node[1] ) ) + { + return ; + } + + if ( ! ( leaf[0] && leaf[1] ) ) + { + int i, j ; + // Fill children nodes + UCHAR* chd[ 2 ][ 8 ] ; + for ( j = 0 ; j < 2 ; j ++ ) + { + for ( i = 0 ; i < 8 ; i ++ ) + { + chd[ j ][ i ] = ((!leaf[j]) && hasChild( node[j], i )) ? getChild( node[j], getChildCount( node[j], i ) ) : NULL ; + } + } + + // 4 face calls + UCHAR* nf[2] ; + int df[2] ; + int lf[2] ; + for ( i = 0 ; i < 4 ; i ++ ) + { + int c[2] = { faceProcFaceMask[ dir ][ i ][ 0 ], faceProcFaceMask[ dir ][ i ][ 1 ] }; + for ( int j = 0 ; j < 2 ; j ++ ) + { + if ( leaf[j] ) + { + lf[j] = leaf[j] ; + nf[j] = node[j] ; + df[j] = depth[j] ; + } + else + { + lf[j] = isLeaf( node[j], c[j] ) ; + nf[j] = chd[ j ][ c[j] ] ; + df[j] = depth[j] - 1 ; + } + } + faceProcContour( nf, lf, df, maxdep - 1, faceProcFaceMask[ dir ][ i ][ 2 ] ) ; + } + + // 4 edge calls + int orders[2][4] = {{ 0, 0, 1, 1 }, { 0, 1, 0, 1 }} ; + UCHAR* ne[4] ; + int le[4] ; + int de[4] ; + + for ( i = 0 ; i < 4 ; i ++ ) + { + int c[4] = { faceProcEdgeMask[ dir ][ i ][ 1 ], faceProcEdgeMask[ dir ][ i ][ 2 ], + faceProcEdgeMask[ dir ][ i ][ 3 ], faceProcEdgeMask[ dir ][ i ][ 4 ] }; + int* order = orders[ faceProcEdgeMask[ dir ][ i ][ 0 ] ] ; + + for ( int j = 0 ; j < 4 ; j ++ ) + { + if ( leaf[order[j]] ) + { + le[j] = leaf[order[j]] ; + ne[j] = node[order[j]] ; + de[j] = depth[order[j]] ; + } + else + { + le[j] = isLeaf( node[order[j]], c[j] ) ; + ne[j] = chd[ order[ j ] ][ c[j] ] ; + de[j] = depth[order[j]] - 1 ; + } + } + + edgeProcContour( ne, le, de, maxdep - 1, faceProcEdgeMask[ dir ][ i ][ 5 ] ) ; + } + } +} + + +void Octree::cellProcContour( UCHAR* node, int leaf, int depth ) +{ + if ( node == NULL ) + { + return ; + } + + if ( ! leaf ) + { + int i ; + + // Fill children nodes + UCHAR* chd[ 8 ] ; + for ( i = 0 ; i < 8 ; i ++ ) + { + chd[ i ] = ((!leaf) && hasChild( node, i )) ? getChild( node, getChildCount( node, i ) ) : NULL ; + } + + // 8 Cell calls + for ( i = 0 ; i < 8 ; i ++ ) + { + cellProcContour( chd[ i ], isLeaf( node, i ), depth - 1 ) ; + } + + // 12 face calls + UCHAR* nf[2] ; + int lf[2] ; + int df[2] = { depth - 1, depth - 1 } ; + for ( i = 0 ; i < 12 ; i ++ ) + { + int c[ 2 ] = { cellProcFaceMask[ i ][ 0 ], cellProcFaceMask[ i ][ 1 ] }; + + lf[0] = isLeaf( node, c[0] ) ; + lf[1] = isLeaf( node, c[1] ) ; + + nf[0] = chd[ c[0] ] ; + nf[1] = chd[ c[1] ] ; + + faceProcContour( nf, lf, df, depth - 1, cellProcFaceMask[ i ][ 2 ] ) ; + } + + // 6 edge calls + UCHAR* ne[4] ; + int le[4] ; + int de[4] = { depth - 1, depth - 1, depth - 1, depth - 1 } ; + for ( i = 0 ; i < 6 ; i ++ ) + { + int c[ 4 ] = { cellProcEdgeMask[ i ][ 0 ], cellProcEdgeMask[ i ][ 1 ], cellProcEdgeMask[ i ][ 2 ], cellProcEdgeMask[ i ][ 3 ] }; + + for ( int j = 0 ; j < 4 ; j ++ ) + { + le[j] = isLeaf( node, c[j] ) ; + ne[j] = chd[ c[j] ] ; + } + + edgeProcContour( ne, le, de, depth - 1, cellProcEdgeMask[ i ][ 4 ] ) ; + } + } + +} + +#endif + + + +void Octree::processEdgeParity( UCHAR* node[4], int depth[4], int maxdep, int dir ) +{ + int con = 0 ; + for ( int i = 0 ; i < 4 ; i ++ ) + { + // Minimal cell + // if ( op == 0 ) + { + if ( getEdgeParity( node[i], processEdgeMask[dir][i] ) ) + { + con = 1 ; + break ; + } + } + } + + if ( con == 1 ) + { + for ( int i = 0 ; i < 4 ; i ++ ) + { + setEdge( node[ i ], processEdgeMask[dir][i] ) ; + } + } + +} + +void Octree::edgeProcParity( UCHAR* node[4], int leaf[4], int depth[4], int maxdep, int dir ) +{ + if ( ! ( node[0] && node[1] && node[2] && node[3] ) ) + { + return ; + } + if ( leaf[0] && leaf[1] && leaf[2] && leaf[3] ) + { + processEdgeParity( node, depth, maxdep, dir ) ; + } + else + { + int i, j ; + UCHAR* chd[ 4 ][ 8 ] ; + for ( j = 0 ; j < 4 ; j ++ ) + { + for ( i = 0 ; i < 8 ; i ++ ) + { + chd[ j ][ i ] = ((!leaf[j]) && hasChild( node[j], i ) )? getChild( node[j], getChildCount( node[j], i ) ) : NULL ; + } + } + + // 2 edge calls + UCHAR* ne[4] ; + int le[4] ; + int de[4] ; + for ( i = 0 ; i < 2 ; i ++ ) + { + int c[ 4 ] = { edgeProcEdgeMask[ dir ][ i ][ 0 ], + edgeProcEdgeMask[ dir ][ i ][ 1 ], + edgeProcEdgeMask[ dir ][ i ][ 2 ], + edgeProcEdgeMask[ dir ][ i ][ 3 ] } ; + + // int allleaf = 1 ; + for ( int j = 0 ; j < 4 ; j ++ ) + { + + if ( leaf[j] ) + { + le[j] = leaf[j] ; + ne[j] = node[j] ; + de[j] = depth[j] ; + } + else + { + le[j] = isLeaf( node[j], c[j] ) ; + ne[j] = chd[ j ][ c[j] ] ; + de[j] = depth[j] - 1 ; + + } + + } + + edgeProcParity( ne, le, de, maxdep - 1, edgeProcEdgeMask[ dir ][ i ][ 4 ] ) ; + } + + } +} + +void Octree::faceProcParity( UCHAR* node[2], int leaf[2], int depth[2], int maxdep, int dir ) +{ + if ( ! ( node[0] && node[1] ) ) + { + return ; + } + + if ( ! ( leaf[0] && leaf[1] ) ) + { + int i, j ; + // Fill children nodes + UCHAR* chd[ 2 ][ 8 ] ; + for ( j = 0 ; j < 2 ; j ++ ) + { + for ( i = 0 ; i < 8 ; i ++ ) + { + chd[ j ][ i ] = ((!leaf[j]) && hasChild( node[j], i )) ? getChild( node[j], getChildCount( node[j], i ) ) : NULL ; + } + } + + // 4 face calls + UCHAR* nf[2] ; + int df[2] ; + int lf[2] ; + for ( i = 0 ; i < 4 ; i ++ ) + { + int c[2] = { faceProcFaceMask[ dir ][ i ][ 0 ], faceProcFaceMask[ dir ][ i ][ 1 ] }; + for ( int j = 0 ; j < 2 ; j ++ ) + { + if ( leaf[j] ) + { + lf[j] = leaf[j] ; + nf[j] = node[j] ; + df[j] = depth[j] ; + } + else + { + lf[j] = isLeaf( node[j], c[j] ) ; + nf[j] = chd[ j ][ c[j] ] ; + df[j] = depth[j] - 1 ; + } + } + faceProcParity( nf, lf, df, maxdep - 1, faceProcFaceMask[ dir ][ i ][ 2 ] ) ; + } + + // 4 edge calls + int orders[2][4] = {{ 0, 0, 1, 1 }, { 0, 1, 0, 1 }} ; + UCHAR* ne[4] ; + int le[4] ; + int de[4] ; + + for ( i = 0 ; i < 4 ; i ++ ) + { + int c[4] = { faceProcEdgeMask[ dir ][ i ][ 1 ], faceProcEdgeMask[ dir ][ i ][ 2 ], + faceProcEdgeMask[ dir ][ i ][ 3 ], faceProcEdgeMask[ dir ][ i ][ 4 ] }; + int* order = orders[ faceProcEdgeMask[ dir ][ i ][ 0 ] ] ; + + for ( int j = 0 ; j < 4 ; j ++ ) + { + if ( leaf[order[j]] ) + { + le[j] = leaf[order[j]] ; + ne[j] = node[order[j]] ; + de[j] = depth[order[j]] ; + } + else + { + le[j] = isLeaf( node[order[j]], c[j] ) ; + ne[j] = chd[ order[ j ] ][ c[j] ] ; + de[j] = depth[order[j]] - 1 ; + } + } + + edgeProcParity( ne, le, de, maxdep - 1, faceProcEdgeMask[ dir ][ i ][ 5 ] ) ; + } + } +} + + +void Octree::cellProcParity( UCHAR* node, int leaf, int depth ) +{ + if ( node == NULL ) + { + return ; + } + + if ( ! leaf ) + { + int i ; + + // Fill children nodes + UCHAR* chd[ 8 ] ; + for ( i = 0 ; i < 8 ; i ++ ) + { + chd[ i ] = ((!leaf) && hasChild( node, i )) ? getChild( node, getChildCount( node, i ) ) : NULL ; + } + + // 8 Cell calls + for ( i = 0 ; i < 8 ; i ++ ) + { + cellProcParity( chd[ i ], isLeaf( node, i ), depth - 1 ) ; + } + + // 12 face calls + UCHAR* nf[2] ; + int lf[2] ; + int df[2] = { depth - 1, depth - 1 } ; + for ( i = 0 ; i < 12 ; i ++ ) + { + int c[ 2 ] = { cellProcFaceMask[ i ][ 0 ], cellProcFaceMask[ i ][ 1 ] }; + + lf[0] = isLeaf( node, c[0] ) ; + lf[1] = isLeaf( node, c[1] ) ; + + nf[0] = chd[ c[0] ] ; + nf[1] = chd[ c[1] ] ; + + faceProcParity( nf, lf, df, depth - 1, cellProcFaceMask[ i ][ 2 ] ) ; + } + + // 6 edge calls + UCHAR* ne[4] ; + int le[4] ; + int de[4] = { depth - 1, depth - 1, depth - 1, depth - 1 } ; + for ( i = 0 ; i < 6 ; i ++ ) + { + int c[ 4 ] = { cellProcEdgeMask[ i ][ 0 ], cellProcEdgeMask[ i ][ 1 ], cellProcEdgeMask[ i ][ 2 ], cellProcEdgeMask[ i ][ 3 ] }; + + for ( int j = 0 ; j < 4 ; j ++ ) + { + le[j] = isLeaf( node, c[j] ) ; + ne[j] = chd[ c[j] ] ; + } + + edgeProcParity( ne, le, de, depth - 1, cellProcEdgeMask[ i ][ 4 ] ) ; + } + } + +} + +/* definitions for global arrays */ +const int edgemask[3] = {5, 3, 6}; + +const int faceMap[6][4] = { + {4, 8, 5, 9}, + {6, 10, 7, 11}, + {0, 8, 1, 10}, + {2, 9, 3, 11}, + {0, 4, 2, 6}, + {1, 5, 3, 7} +}; + +const int cellProcFaceMask[12][3] = { + {0, 4, 0}, + {1, 5, 0}, + {2, 6, 0}, + {3, 7, 0}, + {0, 2, 1}, + {4, 6, 1}, + {1, 3, 1}, + {5, 7, 1}, + {0, 1, 2}, + {2, 3, 2}, + {4, 5, 2}, + {6, 7, 2} +}; + +const int cellProcEdgeMask[6][5] = { + {0, 1, 2, 3, 0}, + {4, 5, 6, 7, 0}, + {0, 4, 1, 5, 1}, + {2, 6, 3, 7, 1}, + {0, 2, 4, 6, 2}, + {1, 3, 5, 7, 2} +}; + +const int faceProcFaceMask[3][4][3] = { + {{4, 0, 0}, + {5, 1, 0}, + {6, 2, 0}, + {7, 3, 0}}, + {{2, 0, 1}, + {6, 4, 1}, + {3, 1, 1}, + {7, 5, 1}}, + {{1, 0, 2}, + {3, 2, 2}, + {5, 4, 2}, + {7, 6, 2}} +}; + +const int faceProcEdgeMask[3][4][6] = { + {{1, 4, 0, 5, 1, 1}, + {1, 6, 2, 7, 3, 1}, + {0, 4, 6, 0, 2, 2}, + {0, 5, 7, 1, 3, 2}}, + {{0, 2, 3, 0, 1, 0}, + {0, 6, 7, 4, 5, 0}, + {1, 2, 0, 6, 4, 2}, + {1, 3, 1, 7, 5, 2}}, + {{1, 1, 0, 3, 2, 0}, + {1, 5, 4, 7, 6, 0}, + {0, 1, 5, 0, 4, 1}, + {0, 3, 7, 2, 6, 1}} +}; + +const int edgeProcEdgeMask[3][2][5] = { + {{3, 2, 1, 0, 0}, + {7, 6, 5, 4, 0}}, + {{5, 1, 4, 0, 1}, + {7, 3, 6, 2, 1}}, + {{6, 4, 2, 0, 2}, + {7, 5, 3, 1, 2}}, +}; + +const int processEdgeMask[3][4] = { + {3, 2, 1, 0}, + {7, 5, 6, 4}, + {11, 10, 9, 8} +}; + +const int dirCell[3][4][3] = { + {{0, -1, -1}, + {0, -1, 0}, + {0, 0, -1}, + {0, 0, 0}}, + {{-1, 0, -1}, + {-1, 0, 0}, + {0, 0, -1}, + {0, 0, 0}}, + {{-1, -1, 0}, + {-1, 0, 0}, + {0, -1, 0}, + {0, 0, 0}} +}; + +const int dirEdge[3][4] = { + {3, 2, 1, 0}, + {7, 6, 5, 4}, + {11, 10, 9, 8} +}; |