diff options
Diffstat (limited to 'intern/dualcon')
-rw-r--r-- | intern/dualcon/intern/Projections.cpp | 305 | ||||
-rw-r--r-- | intern/dualcon/intern/Projections.h | 815 | ||||
-rw-r--r-- | intern/dualcon/intern/octree.cpp | 132 | ||||
-rw-r--r-- | intern/dualcon/intern/octree.h | 2001 |
4 files changed, 1389 insertions, 1864 deletions
diff --git a/intern/dualcon/intern/Projections.cpp b/intern/dualcon/intern/Projections.cpp index e50065c004d..7e7d5e0081c 100644 --- a/intern/dualcon/intern/Projections.cpp +++ b/intern/dualcon/intern/Projections.cpp @@ -71,3 +71,308 @@ const int facemap[6][4] = { {0, 2, 4, 6}, {1, 3, 5, 7} }; + +/** + * Method to perform cross-product + */ +static void crossProduct(int64_t res[3], const int64_t a[3], const int64_t b[3]) +{ + res[0] = a[1] * b[2] - a[2] * b[1]; + res[1] = a[2] * b[0] - a[0] * b[2]; + res[2] = a[0] * b[1] - a[1] * b[0]; +} + +static void crossProduct(double res[3], const double a[3], const double b[3]) +{ + res[0] = a[1] * b[2] - a[2] * b[1]; + res[1] = a[2] * b[0] - a[0] * b[2]; + res[2] = a[0] * b[1] - a[1] * b[0]; +} + +/** + * Method to perform dot product + */ +int64_t dotProduct(const int64_t a[3], const int64_t b[3]) +{ + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; +} + +void normalize(double a[3]) +{ + double mag = a[0] * a[0] + a[1] * a[1] + a[2] * a[2]; + if (mag > 0) { + mag = sqrt(mag); + a[0] /= mag; + a[1] /= mag; + a[2] /= mag; + } +} + +/* Create projection axes for cube+triangle intersection testing. + * 0, 1, 2: cube face normals + * + * 3: triangle normal + * + * 4, 5, 6, + * 7, 8, 9, + * 10, 11, 12: cross of each triangle edge vector with each cube + * face normal + */ +static void create_projection_axes(int64_t axes[NUM_AXES][3], const int64_t tri[3][3]) +{ + /* Cube face normals */ + axes[0][0] = 1; + axes[0][1] = 0; + axes[0][2] = 0; + axes[1][0] = 0; + axes[1][1] = 1; + axes[1][2] = 0; + axes[2][0] = 0; + axes[2][1] = 0; + axes[2][2] = 1; + + /* Get triangle edge vectors */ + int64_t tri_edges[3][3]; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) + tri_edges[i][j] = tri[(i + 1) % 3][j] - tri[i][j]; + } + + /* Triangle normal */ + crossProduct(axes[3], tri_edges[0], tri_edges[1]); + + // Face edges and triangle edges + int ct = 4; + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + crossProduct(axes[ct], axes[j], tri_edges[i]); + ct++; + } + } +} + +/** + * Construction from a cube (axes aligned) and triangle + */ +CubeTriangleIsect::CubeTriangleIsect(int64_t cube[2][3], int64_t tri[3][3], int64_t error, int triind) +{ + int i; + inherit = new TriangleProjection; + inherit->index = triind; + + int64_t axes[NUM_AXES][3]; + create_projection_axes(axes, tri); + + /* Normalize face normal and store */ + double dedge1[] = {(double)tri[1][0] - (double)tri[0][0], + (double)tri[1][1] - (double)tri[0][1], + (double)tri[1][2] - (double)tri[0][2]}; + double dedge2[] = {(double)tri[2][0] - (double)tri[1][0], + (double)tri[2][1] - (double)tri[1][1], + (double)tri[2][2] - (double)tri[1][2]}; + crossProduct(inherit->norm, dedge1, dedge2); + normalize(inherit->norm); + + int64_t cubeedge[3][3]; + for (i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + cubeedge[i][j] = 0; + } + cubeedge[i][i] = cube[1][i] - cube[0][i]; + } + + /* Project the cube on to each axis */ + for (int axis = 0; axis < NUM_AXES; axis++) { + CubeProjection &cube_proj = cubeProj[axis]; + + /* Origin */ + cube_proj.origin = dotProduct(axes[axis], cube[0]); + + /* 3 direction vectors */ + for (i = 0; i < 3; i++) + cube_proj.edges[i] = dotProduct(axes[axis], cubeedge[i]); + + /* Offsets of 2 ends of cube projection */ + int64_t max = 0; + int64_t min = 0; + for (i = 1; i < 8; i++) { + int64_t proj = (vertmap[i][0] * cube_proj.edges[0] + + vertmap[i][1] * cube_proj.edges[1] + + vertmap[i][2] * cube_proj.edges[2]); + if (proj > max) { + max = proj; + } + if (proj < min) { + min = proj; + } + } + cube_proj.min = min; + cube_proj.max = max; + + } + + /* Project the triangle on to each axis */ + for (int axis = 0; axis < NUM_AXES; axis++) { + const int64_t vts[3] = {dotProduct(axes[axis], tri[0]), + dotProduct(axes[axis], tri[1]), + dotProduct(axes[axis], tri[2])}; + + // Triangle + inherit->tri_proj[axis][0] = vts[0]; + inherit->tri_proj[axis][1] = vts[0]; + for (i = 1; i < 3; i++) { + if (vts[i] < inherit->tri_proj[axis][0]) + inherit->tri_proj[axis][0] = vts[i]; + + if (vts[i] > inherit->tri_proj[axis][1]) + inherit->tri_proj[axis][1] = vts[i]; + } + } +} + +/** + * Construction + * from a parent CubeTriangleIsect object and the index of the children + */ +CubeTriangleIsect::CubeTriangleIsect(CubeTriangleIsect *parent) +{ + // Copy inheritable projections + this->inherit = parent->inherit; + + // Shrink cube projections + for (int i = 0; i < NUM_AXES; i++) { + cubeProj[i].origin = parent->cubeProj[i].origin; + + for (int j = 0; j < 3; j++) + cubeProj[i].edges[j] = parent->cubeProj[i].edges[j] >> 1; + + cubeProj[i].min = parent->cubeProj[i].min >> 1; + cubeProj[i].max = parent->cubeProj[i].max >> 1; + } +} + +unsigned char CubeTriangleIsect::getBoxMask( ) +{ + int i, j, k; + int bmask[3][2] = {{0, 0}, {0, 0}, {0, 0}}; + unsigned char boxmask = 0; + int64_t child_len = cubeProj[0].edges[0] >> 1; + + for (i = 0; i < 3; i++) { + int64_t mid = cubeProj[i].origin + child_len; + + // Check bounding box + if (mid >= inherit->tri_proj[i][0]) { + bmask[i][0] = 1; + } + if (mid <= inherit->tri_proj[i][1]) { + bmask[i][1] = 1; + } + + } + + // Fill in masks + int ct = 0; + for (i = 0; i < 2; i++) { + for (j = 0; j < 2; j++) { + for (k = 0; k < 2; k++) { + boxmask |= ( (bmask[0][i] & bmask[1][j] & bmask[2][k]) << ct); + ct++; + } + } + } + + // Return bounding box masks + return boxmask; +} + + +/** + * Shifting a cube to a new origin + */ +void CubeTriangleIsect::shift(int off[3]) +{ + for (int i = 0; i < NUM_AXES; i++) { + cubeProj[i].origin += (off[0] * cubeProj[i].edges[0] + + off[1] * cubeProj[i].edges[1] + + off[2] * cubeProj[i].edges[2]); + } +} + +/** + * Method to test intersection of the triangle and the cube + */ +int CubeTriangleIsect::isIntersecting() const +{ + for (int i = 0; i < NUM_AXES; i++) { + /* + int64_t proj0 = cubeProj[i][0] + + vertmap[inherit->cubeEnds[i][0]][0] * cubeProj[i][1] + + vertmap[inherit->cubeEnds[i][0]][1] * cubeProj[i][2] + + vertmap[inherit->cubeEnds[i][0]][2] * cubeProj[i][3] ; + int64_t proj1 = cubeProj[i][0] + + vertmap[inherit->cubeEnds[i][1]][0] * cubeProj[i][1] + + vertmap[inherit->cubeEnds[i][1]][1] * cubeProj[i][2] + + vertmap[inherit->cubeEnds[i][1]][2] * cubeProj[i][3] ; + */ + + int64_t proj0 = cubeProj[i].origin + cubeProj[i].min; + int64_t proj1 = cubeProj[i].origin + cubeProj[i].max; + + if (proj0 > inherit->tri_proj[i][1] || + proj1 < inherit->tri_proj[i][0]) { + return 0; + } + } + + return 1; +} + +int CubeTriangleIsect::isIntersectingPrimary(int edgeInd) const +{ + for (int i = 0; i < NUM_AXES; i++) { + + int64_t proj0 = cubeProj[i].origin; + int64_t proj1 = cubeProj[i].origin + cubeProj[i].edges[edgeInd]; + + if (proj0 < proj1) { + if (proj0 > inherit->tri_proj[i][1] || + proj1 < inherit->tri_proj[i][0]) { + return 0; + } + } + else { + if (proj1 > inherit->tri_proj[i][1] || + proj0 < inherit->tri_proj[i][0]) { + return 0; + } + } + + } + + // printf( "Intersecting: %d %d\n", edgemap[edgeInd][0], edgemap[edgeInd][1] ) ; + return 1; +} + +float CubeTriangleIsect::getIntersectionPrimary(int edgeInd) const +{ + int i = 3; + + + int64_t proj0 = cubeProj[i].origin; + int64_t proj1 = cubeProj[i].origin + cubeProj[i].edges[edgeInd]; + int64_t proj2 = inherit->tri_proj[i][1]; + int64_t d = proj1 - proj0; + double alpha; + + if (d == 0) + alpha = 0.5; + else { + alpha = (double)((proj2 - proj0)) / (double)d; + + if (alpha < 0 || alpha > 1) + alpha = 0.5; + } + + return (float)alpha; +} diff --git a/intern/dualcon/intern/Projections.h b/intern/dualcon/intern/Projections.h index 7740b0d1634..2cc3320f8e4 100644 --- a/intern/dualcon/intern/Projections.h +++ b/intern/dualcon/intern/Projections.h @@ -32,808 +32,99 @@ #if defined(_WIN32) && !defined(__MINGW32__) #define isnan(n) _isnan(n) #define LONG __int64 +#define int64_t __int64 #else #include <stdint.h> -#define LONG int64_t #endif -#define UCHAR unsigned char /** - * Structures and classes for computing projections of triangles - * onto separating axes during scan conversion - * - * @author Tao Ju - */ - +* Structures and classes for computing projections of triangles onto +* separating axes during scan conversion +* +* @author Tao Ju +*/ extern const int vertmap[8][3]; extern const int centmap[3][3][3][2]; extern const int edgemap[12][2]; extern const int facemap[6][4]; +/* Axes: + * 0, 1, 2: cube face normals + * + * 3: triangle normal + * + * 4, 5, 6, + * 7, 8, 9, + * 10, 11, 12: cross of each triangle edge vector with each cube + * face normal + */ +#define NUM_AXES 13 + /** * Structure for the projections inheritable from parent */ -struct InheritableProjections { - /// Projections of triangle - LONG trigProj[13][2]; - - /// Projections of triangle vertices on primary axes - LONG trigVertProj[13][3]; - - /// Projections of triangle edges - LONG trigEdgeProj[13][3][2]; +struct TriangleProjection { + /// Projections of triangle (min and max) + int64_t tri_proj[NUM_AXES][2]; /// Normal of the triangle double norm[3]; - double normA, normB; - - /// End points along each axis - //int cubeEnds[13][2] ; - /// Error range on each axis - /// LONG errorProj[13]; - -#ifdef CONTAINS_INDEX /// Index of polygon int index; -#endif +}; + +/* This is a projection for the cube against a single projection + axis, see CubeTriangleIsect.cubeProj */ +struct CubeProjection { + int64_t origin; + int64_t edges[3]; + int64_t min, max; }; /** * Class for projections of cube / triangle vertices on the separating axes */ -class Projections +class CubeTriangleIsect { public: -/// Inheritable portion -InheritableProjections *inherit; + /// Inheritable portion + TriangleProjection *inherit; -/// Projections of the cube vertices -LONG cubeProj[13][6]; + /// Projections of the cube vertices + CubeProjection cubeProj[NUM_AXES]; public: + CubeTriangleIsect() {} -Projections( ) -{ -} - -/** - * Construction - * from a cube (axes aligned) and triangle - */ -Projections(LONG cube[2][3], LONG trig[3][3], LONG error, int triind) -{ - int i, j; - inherit = new InheritableProjections; -#ifdef CONTAINS_INDEX - inherit->index = triind; -#endif - /// Create axes - LONG axes[13][3]; - - // Cube faces - axes[0][0] = 1; - axes[0][1] = 0; - axes[0][2] = 0; - - axes[1][0] = 0; - axes[1][1] = 1; - axes[1][2] = 0; - - axes[2][0] = 0; - axes[2][1] = 0; - axes[2][2] = 1; - - // Triangle face - LONG trigedge[3][3]; - for (i = 0; i < 3; i++) - { - for (j = 0; j < 3; j++) - { - trigedge[i][j] = trig[(i + 1) % 3][j] - trig[i][j]; - } - } - crossProduct(trigedge[0], trigedge[1], axes[3]); - - /// Normalize face normal and store - double dedge1[] = { (double) trig[1][0] - (double) trig[0][0], - (double) trig[1][1] - (double) trig[0][1], - (double) trig[1][2] - (double) trig[0][2] }; - double dedge2[] = { (double) trig[2][0] - (double) trig[1][0], - (double) trig[2][1] - (double) trig[1][1], - (double) trig[2][2] - (double) trig[1][2] }; - crossProduct(dedge1, dedge2, inherit->norm); - normalize(inherit->norm); -// inherit->normA = norm[ 0 ] ; -// inherit->normB = norm[ 2 ] > 0 ? norm[ 1 ] : 2 + norm[ 1 ] ; - - // Face edges and triangle edges - int ct = 4; - for (i = 0; i < 3; i++) - for (j = 0; j < 3; j++) - { - crossProduct(axes[j], trigedge[i], axes[ct]); - ct++; - } - - /// Generate projections - LONG cubeedge[3][3]; - for (i = 0; i < 3; i++) - { - for (j = 0; j < 3; j++) - { - cubeedge[i][j] = 0; - } - cubeedge[i][i] = cube[1][i] - cube[0][i]; - } - - for (j = 0; j < 13; j++) - { - // Origin - cubeProj[j][0] = dotProduct(axes[j], cube[0]); - - // 3 direction vectors - for (i = 1; i < 4; i++) - { - cubeProj[j][i] = dotProduct(axes[j], cubeedge[i - 1]); - } - - // Offsets of 2 ends of cube projection - LONG max = 0; - LONG min = 0; - for (i = 1; i < 8; i++) - { - LONG proj = vertmap[i][0] * cubeProj[j][1] + vertmap[i][1] * cubeProj[j][2] + vertmap[i][2] * cubeProj[j][3]; - if (proj > max) - { - max = proj; - } - if (proj < min) - { - min = proj; - } - } - cubeProj[j][4] = min; - cubeProj[j][5] = max; - - } - - for (j = 0; j < 13; j++) - { - LONG vts[3] = { dotProduct(axes[j], trig[0]), - dotProduct(axes[j], trig[1]), - dotProduct(axes[j], trig[2]) }; - - // Vertex - inherit->trigVertProj[j][0] = vts[0]; - inherit->trigVertProj[j][1] = vts[1]; - inherit->trigVertProj[j][2] = vts[2]; - - // Edge - for (i = 0; i < 3; i++) - { - if (vts[i] < vts[(i + 1) % 3]) - { - inherit->trigEdgeProj[j][i][0] = vts[i]; - inherit->trigEdgeProj[j][i][1] = vts[(i + 1) % 3]; - } - else { - inherit->trigEdgeProj[j][i][1] = vts[i]; - inherit->trigEdgeProj[j][i][0] = vts[(i + 1) % 3]; - } - } - - // Triangle - inherit->trigProj[j][0] = vts[0]; - inherit->trigProj[j][1] = vts[0]; - for (i = 1; i < 3; i++) - { - if (vts[i] < inherit->trigProj[j][0]) - { - inherit->trigProj[j][0] = vts[i]; - } - if (vts[i] > inherit->trigProj[j][1]) - { - inherit->trigProj[j][1] = vts[i]; - } - } - } - -} - -/** - * Construction - * from a parent Projections object and the index of the children - */ -Projections (Projections *parent) -{ - // Copy inheritable projections - this->inherit = parent->inherit; - - // Shrink cube projections - for (int i = 0; i < 13; i++) - { - cubeProj[i][0] = parent->cubeProj[i][0]; - for (int j = 1; j < 6; j++) - { - cubeProj[i][j] = parent->cubeProj[i][j] >> 1; - } - } -}; - -Projections (Projections *parent, int box[3], int depth) -{ - int mask = (1 << depth) - 1; - int nbox[3] = { box[0] & mask, box[1] & mask, box[2] & mask }; - - // Copy inheritable projections - this->inherit = parent->inherit; - - // Shrink cube projections - for (int i = 0; i < 13; i++) - { - for (int j = 1; j < 6; j++) - { - cubeProj[i][j] = parent->cubeProj[i][j] >> depth; - } - - cubeProj[i][0] = parent->cubeProj[i][0] + nbox[0] * cubeProj[i][1] + nbox[1] * cubeProj[i][2] + nbox[2] * cubeProj[i][3]; - } -}; - -/** - * Testing intersection based on vertex/edge masks - */ -int getIntersectionMasks(UCHAR cedgemask, UCHAR& edgemask) -{ - int i, j; - edgemask = cedgemask; - - // Pre-processing - /* - if ( cvertmask & 1 ) - { - edgemask |= 5 ; - } - if ( cvertmask & 2 ) - { - edgemask |= 3 ; - } - if ( cvertmask & 4 ) - { - edgemask |= 6 ; - } - + /** + * Construction from a cube (axes aligned) and triangle */ - - // Test axes for edge intersection - UCHAR bit = 1; - for (j = 0; j < 3; j++) - { - if (edgemask & bit) - { - for (i = 0; i < 13; i++) - { - LONG proj0 = cubeProj[i][0] + cubeProj[i][4]; - LONG proj1 = cubeProj[i][0] + cubeProj[i][5]; - - if (proj0 > inherit->trigEdgeProj[i][j][1] || - proj1 < inherit->trigEdgeProj[i][j][0]) - { - edgemask &= (~bit); - break; - } - } - } - bit <<= 1; - } - - /* - if ( edgemask != 0 ) - { - printf("%d %d\n", cedgemask, edgemask) ; - } + CubeTriangleIsect(int64_t cube[2][3], int64_t trig[3][3], int64_t error, int triind); + + /** + * Construction from a parent CubeTriangleIsect object and the index of + * the children */ + CubeTriangleIsect(CubeTriangleIsect *parent); + + unsigned char getBoxMask( ); - // Test axes for triangle intersection - if (edgemask) - { - return 1; - } - - for (i = 3; i < 13; i++) - { - LONG proj0 = cubeProj[i][0] + cubeProj[i][4]; - LONG proj1 = cubeProj[i][0] + cubeProj[i][5]; - - if (proj0 > inherit->trigProj[i][1] || - proj1 < inherit->trigProj[i][0]) - { - return 0; - } - } - - return 1; -} - -/** - * Retrieving children masks using PRIMARY AXES - */ -UCHAR getChildrenMasks(UCHAR cvertmask, UCHAR vertmask[8]) -{ - int i, j, k; - int bmask[3][2] = {{0, 0}, {0, 0}, {0, 0}}; - int vmask[3][3][2] = {{{0, 0}, {0, 0}, {0, 0}}, {{0, 0}, {0, 0}, {0, 0}}, {{0, 0}, {0, 0}, {0, 0}}}; - UCHAR boxmask = 0; - LONG len = cubeProj[0][1] >> 1; - - for (i = 0; i < 3; i++) - { - LONG mid = cubeProj[i][0] + len; - - // Check bounding box - if (mid >= inherit->trigProj[i][0]) - { - bmask[i][0] = 1; - } - if (mid <= inherit->trigProj[i][1]) - { - bmask[i][1] = 1; - } - - // Check vertex mask - if (cvertmask) - { - for (j = 0; j < 3; j++) - { - if (cvertmask & (1 << j) ) - { - // Only check if it's contained this node - if (mid >= inherit->trigVertProj[i][j]) - { - vmask[i][j][0] = 1; - } - if (mid <= inherit->trigVertProj[i][j]) - { - vmask[i][j][1] = 1; - } - } - } - } - - /* - // Check edge mask - if ( cedgemask ) - { - for ( j = 0 ; j < 3 ; j ++ ) - { - if ( cedgemask & ( 1 << j ) ) - { - // Only check if it's contained this node - if ( mid >= inherit->trigEdgeProj[i][j][0] ) - { - emask[i][j][0] = 1 ; - } - if ( mid <= inherit->trigEdgeProj[i][j][1] ) - { - emask[i][j][1] = 1 ; - } - } - } - } - */ - - } - - // Fill in masks - int ct = 0; - for (i = 0; i < 2; i++) - for (j = 0; j < 2; j++) - for (k = 0; k < 2; k++) - { - boxmask |= ( (bmask[0][i] & bmask[1][j] & bmask[2][k]) << ct); - vertmask[ct] = ((vmask[0][0][i] & vmask[1][0][j] & vmask[2][0][k]) | - ((vmask[0][1][i] & vmask[1][1][j] & vmask[2][1][k]) << 1) | - ((vmask[0][2][i] & vmask[1][2][j] & vmask[2][2][k]) << 2) ); - /* - edgemask[ct] = (( emask[0][0][i] & emask[1][0][j] & emask[2][0][k] ) | - (( emask[0][1][i] & emask[1][1][j] & emask[2][1][k] ) << 1 ) | - (( emask[0][2][i] & emask[1][2][j] & emask[2][2][k] ) << 2 ) ) ; - edgemask[ct] = cedgemask ; - */ - ct++; - } - - // Return bounding box masks - return boxmask; -} - -UCHAR getBoxMask( ) -{ - int i, j, k; - int bmask[3][2] = {{0, 0}, {0, 0}, {0, 0}}; - UCHAR boxmask = 0; - LONG len = cubeProj[0][1] >> 1; - - for (i = 0; i < 3; i++) - { - LONG mid = cubeProj[i][0] + len; - - // Check bounding box - if (mid >= inherit->trigProj[i][0]) - { - bmask[i][0] = 1; - } - if (mid <= inherit->trigProj[i][1]) - { - bmask[i][1] = 1; - } - - } - - // Fill in masks - int ct = 0; - for (i = 0; i < 2; i++) - for (j = 0; j < 2; j++) - for (k = 0; k < 2; k++) - { - boxmask |= ( (bmask[0][i] & bmask[1][j] & bmask[2][k]) << ct); - ct++; - } - - // Return bounding box masks - return boxmask; -} - - -/** - * Get projections for sub-cubes (simple axes) - */ -void getSubProjectionsSimple(Projections *p[8]) -{ - // Process the axes cooresponding to the triangle's normal - int ind = 3; - LONG len = cubeProj[0][1] >> 1; - LONG trigproj[3] = { cubeProj[ind][1] >> 1, cubeProj[ind][2] >> 1, cubeProj[ind][3] >> 1 }; - - int ct = 0; - for (int i = 0; i < 2; i++) - for (int j = 0; j < 2; j++) - for (int k = 0; k < 2; k++) - { - p[ct] = new Projections( ); - p[ct]->inherit = inherit; - - p[ct]->cubeProj[0][0] = cubeProj[0][0] + i * len; - p[ct]->cubeProj[1][0] = cubeProj[1][0] + j * len; - p[ct]->cubeProj[2][0] = cubeProj[2][0] + k * len; - p[ct]->cubeProj[0][1] = len; - - for (int m = 1; m < 4; m++) - { - p[ct]->cubeProj[ind][m] = trigproj[m - 1]; - } - p[ct]->cubeProj[ind][0] = cubeProj[ind][0] + i * trigproj[0] + j * trigproj[1] + k * trigproj[2]; - - ct++; - } -} - -/** - * Shifting a cube to a new origin - */ -void shift(int off[3]) -{ - for (int i = 0; i < 13; i++) - { - cubeProj[i][0] += off[0] * cubeProj[i][1] + off[1] * cubeProj[i][2] + off[2] * cubeProj[i][3]; - } -} - -void shiftNoPrimary(int off[3]) -{ - for (int i = 3; i < 13; i++) - { - cubeProj[i][0] += off[0] * cubeProj[i][1] + off[1] * cubeProj[i][2] + off[2] * cubeProj[i][3]; - } -} - -/** - * Method to test intersection of the triangle and the cube - */ -int isIntersecting( ) -{ - for (int i = 0; i < 13; i++) - { - /* - LONG proj0 = cubeProj[i][0] + - vertmap[inherit->cubeEnds[i][0]][0] * cubeProj[i][1] + - vertmap[inherit->cubeEnds[i][0]][1] * cubeProj[i][2] + - vertmap[inherit->cubeEnds[i][0]][2] * cubeProj[i][3] ; - LONG proj1 = cubeProj[i][0] + - vertmap[inherit->cubeEnds[i][1]][0] * cubeProj[i][1] + - vertmap[inherit->cubeEnds[i][1]][1] * cubeProj[i][2] + - vertmap[inherit->cubeEnds[i][1]][2] * cubeProj[i][3] ; - */ - - LONG proj0 = cubeProj[i][0] + cubeProj[i][4]; - LONG proj1 = cubeProj[i][0] + cubeProj[i][5]; - - if (proj0 > inherit->trigProj[i][1] || - proj1 < inherit->trigProj[i][0]) - { - return 0; - } - } - - return 1; -}; - -int isIntersectingNoPrimary( ) -{ - for (int i = 3; i < 13; i++) - { - /* - LONG proj0 = cubeProj[i][0] + - vertmap[inherit->cubeEnds[i][0]][0] * cubeProj[i][1] + - vertmap[inherit->cubeEnds[i][0]][1] * cubeProj[i][2] + - vertmap[inherit->cubeEnds[i][0]][2] * cubeProj[i][3] ; - LONG proj1 = cubeProj[i][0] + - vertmap[inherit->cubeEnds[i][1]][0] * cubeProj[i][1] + - vertmap[inherit->cubeEnds[i][1]][1] * cubeProj[i][2] + - vertmap[inherit->cubeEnds[i][1]][2] * cubeProj[i][3] ; - */ - - LONG proj0 = cubeProj[i][0] + cubeProj[i][4]; - LONG proj1 = cubeProj[i][0] + cubeProj[i][5]; - - if (proj0 > inherit->trigProj[i][1] || - proj1 < inherit->trigProj[i][0]) - { - return 0; - } - } - - return 1; -}; - -/** - * Method to test intersection of the triangle and one edge - */ -int isIntersecting(int edgeInd) -{ - for (int i = 0; i < 13; i++) - { - - LONG proj0 = cubeProj[i][0] + - vertmap[edgemap[edgeInd][0]][0] * cubeProj[i][1] + - vertmap[edgemap[edgeInd][0]][1] * cubeProj[i][2] + - vertmap[edgemap[edgeInd][0]][2] * cubeProj[i][3]; - LONG proj1 = cubeProj[i][0] + - vertmap[edgemap[edgeInd][1]][0] * cubeProj[i][1] + - vertmap[edgemap[edgeInd][1]][1] * cubeProj[i][2] + - vertmap[edgemap[edgeInd][1]][2] * cubeProj[i][3]; - - - if (proj0 < proj1) - { - if (proj0 > inherit->trigProj[i][1] || - proj1 < inherit->trigProj[i][0]) - { - return 0; - } - } - else { - if (proj1 > inherit->trigProj[i][1] || - proj0 < inherit->trigProj[i][0]) - { - return 0; - } - } - } - - // printf( "Intersecting: %d %d\n", edgemap[edgeInd][0], edgemap[edgeInd][1] ) ; - return 1; -}; - -/** - * Method to test intersection of one triangle edge and one cube face - */ -int isIntersecting(int edgeInd, int faceInd) -{ - for (int i = 0; i < 13; i++) - { - LONG trigproj0 = inherit->trigVertProj[i][edgeInd]; - LONG trigproj1 = inherit->trigVertProj[i][(edgeInd + 1) % 3]; - - if (trigproj0 < trigproj1) - { - int t1 = 1, t2 = 1; - for (int j = 0; j < 4; j++) - { - LONG proj = cubeProj[i][0] + - vertmap[facemap[faceInd][j]][0] * cubeProj[i][1] + - vertmap[facemap[faceInd][j]][1] * cubeProj[i][2] + - vertmap[facemap[faceInd][j]][2] * cubeProj[i][3]; - if (proj >= trigproj0) - { - t1 = 0; - } - if (proj <= trigproj1) - { - t2 = 0; - } - } - if (t1 || t2) - { - return 0; - } - } - else { - int t1 = 1, t2 = 1; - for (int j = 0; j < 4; j++) - { - LONG proj = cubeProj[i][0] + - vertmap[facemap[faceInd][j]][0] * cubeProj[i][1] + - vertmap[facemap[faceInd][j]][1] * cubeProj[i][2] + - vertmap[facemap[faceInd][j]][2] * cubeProj[i][3]; - if (proj >= trigproj1) - { - t1 = 0; - } - if (proj <= trigproj0) - { - t2 = 0; - } - } - if (t1 || t2) - { - return 0; - } - } - } - - return 1; -}; - - -int isIntersectingPrimary(int edgeInd) -{ - for (int i = 0; i < 13; i++) - { - - LONG proj0 = cubeProj[i][0]; - LONG proj1 = cubeProj[i][0] + cubeProj[i][edgeInd + 1]; - - if (proj0 < proj1) - { - if (proj0 > inherit->trigProj[i][1] || - proj1 < inherit->trigProj[i][0]) - { - return 0; - } - } - else { - if (proj1 > inherit->trigProj[i][1] || - proj0 < inherit->trigProj[i][0]) - { - return 0; - } - } - - } - - // printf( "Intersecting: %d %d\n", edgemap[edgeInd][0], edgemap[edgeInd][1] ) ; - return 1; -}; - -double getIntersection(int edgeInd) -{ - int i = 3; - - LONG proj0 = cubeProj[i][0] + - vertmap[edgemap[edgeInd][0]][0] * cubeProj[i][1] + - vertmap[edgemap[edgeInd][0]][1] * cubeProj[i][2] + - vertmap[edgemap[edgeInd][0]][2] * cubeProj[i][3]; - LONG proj1 = cubeProj[i][0] + - vertmap[edgemap[edgeInd][1]][0] * cubeProj[i][1] + - vertmap[edgemap[edgeInd][1]][1] * cubeProj[i][2] + - vertmap[edgemap[edgeInd][1]][2] * cubeProj[i][3]; - LONG proj2 = inherit->trigProj[i][1]; - - /* - if ( proj0 < proj1 ) - { - if ( proj2 < proj0 || proj2 > proj1 ) - { - return -1 ; - } - } - else - { - if ( proj2 < proj1 || proj2 > proj0 ) - { - return -1 ; - } - } + /** + * Shifting a cube to a new origin */ + void shift(int off[3]); - double alpha = (double)(proj2 - proj0) / (double)(proj1 - proj0); - /* - if ( alpha < 0 ) - { - alpha = 0.5 ; - } - else if ( alpha > 1 ) - { - alpha = 0.5 ; - } + /** + * Method to test intersection of the triangle and the cube */ + int isIntersecting() const; - return alpha; -}; - -float getIntersectionPrimary(int edgeInd) -{ - int i = 3; - - - LONG proj0 = cubeProj[i][0]; - LONG proj1 = cubeProj[i][0] + cubeProj[i][edgeInd + 1]; - LONG proj2 = inherit->trigProj[i][1]; - LONG d = proj1 - proj0; - double alpha; - - if (d == 0) - alpha = 0.5; - else { - alpha = (double)((proj2 - proj0)) / (double)d; - - if (alpha < 0 || alpha > 1) - alpha = 0.5; - } - - return (float)alpha; -}; - -/** - * Method to perform cross-product - */ -void crossProduct(LONG a[3], LONG b[3], LONG res[3]) -{ - res[0] = a[1] * b[2] - a[2] * b[1]; - res[1] = a[2] * b[0] - a[0] * b[2]; - res[2] = a[0] * b[1] - a[1] * b[0]; -} -void crossProduct(double a[3], double b[3], double res[3]) -{ - res[0] = a[1] * b[2] - a[2] * b[1]; - res[1] = a[2] * b[0] - a[0] * b[2]; - res[2] = a[0] * b[1] - a[1] * b[0]; -} - -/** - * Method to perform dot product - */ -LONG dotProduct(LONG a[3], LONG b[3]) -{ - return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; -} - -void normalize(double a[3]) -{ - double mag = a[0] * a[0] + a[1] * a[1] + a[2] * a[2]; - if (mag > 0) - { - mag = sqrt(mag); - a[0] /= mag; - a[1] /= mag; - a[2] /= mag; - } -} + int isIntersectingPrimary(int edgeInd) const; + float getIntersectionPrimary(int edgeInd) const; }; #endif diff --git a/intern/dualcon/intern/octree.cpp b/intern/dualcon/intern/octree.cpp index c5f50d8af4e..1ad502b018d 100644 --- a/intern/dualcon/intern/octree.cpp +++ b/intern/dualcon/intern/octree.cpp @@ -113,7 +113,7 @@ void Octree::scanConvert() start = clock(); #endif - addTrian(); + addAllTriangles(); resetMinimalEdges(); preparePrimalEdgesMask(&root->internal); @@ -257,7 +257,7 @@ void Octree::resetMinimalEdges() cellProcParity(root, 0, maxDepth); } -void Octree::addTrian() +void Octree::addAllTriangles() { Triangle *trian; int count = 0; @@ -273,7 +273,7 @@ void Octree::addTrian() while ((trian = reader->getNextTriangle()) != NULL) { // Drop triangles { - addTrian(trian, count); + addTriangle(trian, count); } delete trian; @@ -316,48 +316,60 @@ void Octree::addTrian() putchar(13); } -void Octree::addTrian(Triangle *trian, int triind) +/* Prepare a triangle for insertion into the octree; call the other + addTriangle() to (recursively) build the octree */ +void Octree::addTriangle(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++) { + /* Project the triangle's coordinates into the grid */ + 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 - } + /* Generate projections */ + int64_t cube[2][3] = {{0, 0, 0}, {dimen, dimen, dimen}}; + int64_t trig[3][3]; + for (i = 0; i < 3; i++) { + for (j = 0; j < 3; j++) + trig[i][j] = (int64_t)(trian->vt[i][j]); + } - // Add to the octree - // int start[3] = {0, 0, 0}; - LONG errorvec = (LONG)(0); - Projections *proj = new Projections(cube, trig, errorvec, triind); - root = (Node *)addTrian(&root->internal, proj, maxDepth); + /* Add triangle to the octree */ + int64_t errorvec = (int64_t)(0); + CubeTriangleIsect *proj = new CubeTriangleIsect(cube, trig, errorvec, triind); + root = (Node *)addTriangle(&root->internal, proj, maxDepth); delete proj->inherit; delete proj; } +void print_depth(int height, int maxDepth) +{ + for (int i = 0; i < maxDepth - height; i++) + printf(" "); +} -InternalNode *Octree::addTrian(InternalNode *node, Projections *p, int height) +InternalNode *Octree::addTriangle(InternalNode *node, CubeTriangleIsect *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); - + const 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}}; + unsigned char boxmask = p->getBoxMask(); + CubeTriangleIsect *subp = new CubeTriangleIsect(p); + int count = 0; int tempdiff[3] = {0, 0, 0}; + + /* Check triangle against each of the input node's children */ for (i = 0; i < 8; i++) { tempdiff[0] += vertdiff[i][0]; tempdiff[1] += vertdiff[i][1]; @@ -370,30 +382,23 @@ InternalNode *Octree::addTrian(InternalNode *node, Projections *p, int height) /* Pruning using intersection test */ if (subp->isIntersecting()) { - // if(subp->getIntersectionMasks(cedgemask, edgemask)) if (!hasChild(node, i)) { - if (height == 1) { + if (height == 1) node = addLeafChild(node, i, count, createLeaf(0)); - } - else { + else node = addInternalChild(node, i, count, createInternal(0)); - } } Node *chd = getChild(node, count); - if (!isLeaf(node, i)) { - // setChild(node, count, addTrian(chd, subp, height - 1, vertmask[i], edgemask)); - setChild(node, count, (Node *)addTrian(&chd->internal, subp, height - 1)); - } - else { + if (node->is_child_leaf(i)) setChild(node, count, (Node *)updateCell(&chd->leaf, subp)); - } + else + setChild(node, count, (Node *)addTriangle(&chd->internal, subp, height - 1)); } } - if (hasChild(node, i)) { + if (hasChild(node, i)) count++; - } } delete subp; @@ -401,7 +406,7 @@ InternalNode *Octree::addTrian(InternalNode *node, Projections *p, int height) return node; } -LeafNode *Octree::updateCell(LeafNode *node, Projections *p) +LeafNode *Octree::updateCell(LeafNode *node, CubeTriangleIsect *p) { int i; @@ -426,13 +431,6 @@ LeafNode *Octree::updateCell(LeafNode *node, Projections *p) else { offs[newc] = getEdgeOffsetNormal(node, oldc, a[newc], b[newc], c[newc]); -// if(p->isIntersectingPrimary(i)) - { - // dc_printf("Multiple intersections!\n"); - -// setPatchEdge(node, i); - } - oldc++; newc++; } @@ -451,7 +449,7 @@ void Octree::preparePrimalEdgesMask(InternalNode *node) int count = 0; for (int i = 0; i < 8; i++) { if (hasChild(node, i)) { - if (isLeaf(node, i)) + if (node->is_child_leaf(i)) createPrimalEdgesMask(&getChild(node, count)->leaf); else preparePrimalEdgesMask(&getChild(node, count)->internal); @@ -495,7 +493,7 @@ Node *Octree::trace(Node *newnode, int *st, int len, int depth, PathList *& path nst[i][j] = st[j] + len * vertmap[i][j]; } - if (chd[i] == NULL || isLeaf(&newnode->internal, i)) { + if (chd[i] == NULL || newnode->internal.is_child_leaf(i)) { chdpaths[i] = NULL; } else { @@ -1411,7 +1409,7 @@ Node *Octree::locateCell(InternalNode *node, int st[3], int len, int ori[3], int if (hasChild(node, ind)) { int count = getChildCount(node, ind); Node *chd = getChild(node, count); - if (isLeaf(node, ind)) { + if (node->is_child_leaf(ind)) { rleaf = chd; rlen = len; } @@ -2367,7 +2365,7 @@ void Octree::edgeProcContour(Node *node[4], int leaf[4], int depth[4], int maxde de[j] = depth[j]; } else { - le[j] = isLeaf(&node[j]->internal, c[j]); + le[j] = node[j]->internal.is_child_leaf(c[j]); ne[j] = chd[j][c[j]]; de[j] = depth[j] - 1; } @@ -2410,7 +2408,7 @@ void Octree::faceProcContour(Node *node[2], int leaf[2], int depth[2], int maxde df[j] = depth[j]; } else { - lf[j] = isLeaf(&node[j]->internal, c[j]); + lf[j] = node[j]->internal.is_child_leaf(c[j]); nf[j] = chd[j][c[j]]; df[j] = depth[j] - 1; } @@ -2436,7 +2434,7 @@ void Octree::faceProcContour(Node *node[2], int leaf[2], int depth[2], int maxde de[j] = depth[order[j]]; } else { - le[j] = isLeaf(&node[order[j]]->internal, c[j]); + le[j] = node[order[j]]->internal.is_child_leaf(c[j]); ne[j] = chd[order[j]][c[j]]; de[j] = depth[order[j]] - 1; } @@ -2467,7 +2465,7 @@ void Octree::cellProcContour(Node *node, int leaf, int depth) // 8 Cell calls for (i = 0; i < 8; i++) { - cellProcContour(chd[i], isLeaf(&node->internal, i), depth - 1); + cellProcContour(chd[i], node->internal.is_child_leaf(i), depth - 1); } // 12 face calls @@ -2477,8 +2475,8 @@ void Octree::cellProcContour(Node *node, int leaf, int depth) for (i = 0; i < 12; i++) { int c[2] = {cellProcFaceMask[i][0], cellProcFaceMask[i][1]}; - lf[0] = isLeaf(&node->internal, c[0]); - lf[1] = isLeaf(&node->internal, c[1]); + lf[0] = node->internal.is_child_leaf(c[0]); + lf[1] = node->internal.is_child_leaf(c[1]); nf[0] = chd[c[0]]; nf[1] = chd[c[1]]; @@ -2494,7 +2492,7 @@ void Octree::cellProcContour(Node *node, int leaf, int depth) 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->internal, c[j]); + le[j] = node->internal.is_child_leaf(c[j]); ne[j] = chd[c[j]]; } @@ -2563,7 +2561,7 @@ void Octree::edgeProcParity(Node *node[4], int leaf[4], int depth[4], int maxdep de[j] = depth[j]; } else { - le[j] = isLeaf(&node[j]->internal, c[j]); + le[j] = node[j]->internal.is_child_leaf(c[j]); ne[j] = chd[j][c[j]]; de[j] = depth[j] - 1; @@ -2608,7 +2606,7 @@ void Octree::faceProcParity(Node *node[2], int leaf[2], int depth[2], int maxdep df[j] = depth[j]; } else { - lf[j] = isLeaf(&node[j]->internal, c[j]); + lf[j] = node[j]->internal.is_child_leaf(c[j]); nf[j] = chd[j][c[j]]; df[j] = depth[j] - 1; } @@ -2634,7 +2632,7 @@ void Octree::faceProcParity(Node *node[2], int leaf[2], int depth[2], int maxdep de[j] = depth[order[j]]; } else { - le[j] = isLeaf((InternalNode *)(node[order[j]]), c[j]); + le[j] = node[order[j]]->internal.is_child_leaf(c[j]); ne[j] = chd[order[j]][c[j]]; de[j] = depth[order[j]] - 1; } @@ -2665,7 +2663,7 @@ void Octree::cellProcParity(Node *node, int leaf, int depth) // 8 Cell calls for (i = 0; i < 8; i++) { - cellProcParity(chd[i], isLeaf((InternalNode *)node, i), depth - 1); + cellProcParity(chd[i], node->internal.is_child_leaf(i), depth - 1); } // 12 face calls @@ -2675,8 +2673,8 @@ void Octree::cellProcParity(Node *node, int leaf, int depth) for (i = 0; i < 12; i++) { int c[2] = {cellProcFaceMask[i][0], cellProcFaceMask[i][1]}; - lf[0] = isLeaf((InternalNode *)node, c[0]); - lf[1] = isLeaf((InternalNode *)node, c[1]); + lf[0] = node->internal.is_child_leaf(c[0]); + lf[1] = node->internal.is_child_leaf(c[1]); nf[0] = chd[c[0]]; nf[1] = chd[c[1]]; @@ -2692,7 +2690,7 @@ void Octree::cellProcParity(Node *node, int leaf, int depth) 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((InternalNode *)node, c[j]); + le[j] = node->internal.is_child_leaf(c[j]); ne[j] = chd[c[j]]; } diff --git a/intern/dualcon/intern/octree.h b/intern/dualcon/intern/octree.h index 35d24a074bb..550d584baa7 100644 --- a/intern/dualcon/intern/octree.h +++ b/intern/dualcon/intern/octree.h @@ -65,6 +65,12 @@ struct InternalNode { /* Can have up to eight children */ Node *children[0]; + + /// Test if child is leaf + int is_child_leaf(int index) const + { + return (child_is_leaf >> index) & 1; + } }; @@ -145,1277 +151,1202 @@ struct PathList { */ class Octree { -public: -/* Public members */ - -/// Memory allocators -VirtualMemoryAllocator *alloc[9]; -VirtualMemoryAllocator *leafalloc[4]; + public: + /* Public members */ -/// Root node -Node *root; + /// Memory allocators + VirtualMemoryAllocator *alloc[9]; + VirtualMemoryAllocator *leafalloc[4]; -/// Model reader -ModelReader *reader; + /// Root node + Node *root; -/// Marching cubes table -Cubes *cubes; + /// Model reader + ModelReader *reader; -/// Length of grid -int dimen; -int mindimen, minshift; + /// Marching cubes table + Cubes *cubes; -/// Maximum depth -int maxDepth; + /// Length of grid + int dimen; + int mindimen, minshift; -/// The lower corner of the bounding box and the size -float origin[3]; -float range; + /// Maximum depth + int maxDepth; -/// Counting information -int nodeCount; -int nodeSpace; -int nodeCounts[9]; + /// The lower corner of the bounding box and the size + float origin[3]; + float range; -int actualQuads, actualVerts; + /// Counting information + int nodeCount; + int nodeSpace; + int nodeCounts[9]; -PathList *ringList; + int actualQuads, actualVerts; -int maxTrianglePerCell; -int outType; // 0 for OFF, 1 for PLY, 2 for VOL + PathList *ringList; -// For flood filling -int use_flood_fill; -float thresh; + int maxTrianglePerCell; + int outType; // 0 for OFF, 1 for PLY, 2 for VOL -int use_manifold; + // For flood filling + int use_flood_fill; + float thresh; -float hermite_num; + int use_manifold; -DualConMode mode; - -public: -/** - * Construtor - */ -Octree(ModelReader *mr, - DualConAllocOutput alloc_output_func, - DualConAddVert add_vert_func, - DualConAddQuad add_quad_func, - DualConFlags flags, DualConMode mode, int depth, - float threshold, float hermite_num); - -/** - * Destructor - */ -~Octree(); - -/** - * Scan convert - */ -void scanConvert(); + float hermite_num; -void *getOutputMesh() { - return output_mesh; -} + DualConMode mode; -private: -/* Helper functions */ - -/** - * Initialize memory allocators - */ -void initMemory(); - -/** - * Release memory - */ -void freeMemory(); - -/** - * Print memory usage - */ -void printMemUsage(); + public: + /** + * Construtor + */ + Octree(ModelReader *mr, + DualConAllocOutput alloc_output_func, + DualConAddVert add_vert_func, + DualConAddQuad add_quad_func, + DualConFlags flags, DualConMode mode, int depth, + float threshold, float hermite_num); + + /** + * Destructor + */ + ~Octree(); + /** + * Scan convert + */ + void scanConvert(); -/** - * Methods to set / restore minimum edges - */ -void resetMinimalEdges(); + void *getOutputMesh() { + return output_mesh; + } -void cellProcParity(Node *node, int leaf, int depth); -void faceProcParity(Node * node[2], int leaf[2], int depth[2], int maxdep, int dir); -void edgeProcParity(Node * node[4], int leaf[4], int depth[4], int maxdep, int dir); + private: + /* Helper functions */ -void processEdgeParity(LeafNode * node[4], int depths[4], int maxdep, int dir); + /** + * Initialize memory allocators + */ + void initMemory(); -/** - * Add triangles to the tree - */ -void addTrian(); -void addTrian(Triangle *trian, int triind); -InternalNode *addTrian(InternalNode *node, Projections *p, int height); + /** + * Release memory + */ + void freeMemory(); -/** - * Method to update minimizer in a cell: update edge intersections instead - */ -LeafNode *updateCell(LeafNode *node, Projections *p); + /** + * Print memory usage + */ + void printMemUsage(); -/* Routines to detect and patch holes */ -int numRings; -int totRingLengths; -int maxRingLength; -/** - * Entry routine. - */ -void trace(); -/** - * Trace the given node, find patches and fill them in - */ -Node *trace(Node *node, int *st, int len, int depth, PathList *& paths); -/** - * Look for path on the face and add to paths - */ -void findPaths(Node * node[2], int leaf[2], int depth[2], int *st[2], int maxdep, int dir, PathList * &paths); -/** - * Combine two list1 and list2 into list1 using connecting paths list3, - * while closed paths are appended to rings - */ -void combinePaths(PathList *& list1, PathList *list2, PathList *paths, PathList *& rings); -/** - * Helper function: combine current paths in list1 and list2 to a single path and append to list3 - */ -PathList *combineSinglePath(PathList *& head1, PathList *pre1, PathList *& list1, PathList *& head2, PathList *pre2, PathList *& list2); + /** + * Methods to set / restore minimum edges + */ + void resetMinimalEdges(); -/** - * Functions to patch rings in a node - */ -Node *patch(Node * node, int st[3], int len, PathList * rings); -Node *patchSplit(Node * node, int st[3], int len, PathList * rings, int dir, PathList * &nrings1, PathList * &nrings2); -Node *patchSplitSingle(Node * node, int st[3], int len, PathElement * head, int dir, PathList * &nrings1, PathList * &nrings2); -Node *connectFace(Node * node, int st[3], int len, int dir, PathElement * f1, PathElement * f2); -Node *locateCell(InternalNode * node, int st[3], int len, int ori[3], int dir, int side, Node * &rleaf, int rst[3], int& rlen); -void compressRing(PathElement *& ring); -void getFacePoint(PathElement *leaf, int dir, int& x, int& y, float& p, float& q); -LeafNode *patchAdjacent(InternalNode * node, int len, int st1[3], LeafNode * leaf1, int st2[3], LeafNode * leaf2, int walkdir, int inc, int dir, int side, float alpha); -int findPair(PathElement *head, int pos, int dir, PathElement *& pre1, PathElement *& pre2); -int getSide(PathElement *e, int pos, int dir); -int isEqual(PathElement *e1, PathElement *e2); -void preparePrimalEdgesMask(InternalNode *node); -void testFacePoint(PathElement *e1, PathElement *e2); + void cellProcParity(Node *node, int leaf, int depth); + void faceProcParity(Node * node[2], int leaf[2], int depth[2], int maxdep, int dir); + void edgeProcParity(Node * node[4], int leaf[4], int depth[4], int maxdep, int dir); -/** - * Path-related functions - */ -void deletePath(PathList *& head, PathList *pre, PathList *& curr); -void printPath(PathList *path); -void printPath(PathElement *path); -void printElement(PathElement *ele); -void printPaths(PathList *path); -void checkElement(PathElement *ele); -void checkPath(PathElement *path); + void processEdgeParity(LeafNode * node[4], int depths[4], int maxdep, int dir); + /** + * Add triangles to the tree + */ + void addAllTriangles(); + void addTriangle(Triangle *trian, int triind); + InternalNode *addTriangle(InternalNode *node, CubeTriangleIsect *p, int height); -/** - * Routines to build signs to create a partitioned volume - *(after patching rings) - */ -void buildSigns(); -void buildSigns(unsigned char table[], Node * node, int isLeaf, int sg, int rvalue[8]); + /** + * Method to update minimizer in a cell: update edge intersections instead + */ + LeafNode *updateCell(LeafNode *node, CubeTriangleIsect *p); -/************************************************************************/ -/* To remove disconnected components */ -/************************************************************************/ -void floodFill(); -void clearProcessBits(Node *node, int height); -int floodFill(LeafNode * leaf, int st[3], int len, int height, int threshold); -int floodFill(Node * node, int st[3], int len, int height, int threshold); + /* Routines to detect and patch holes */ + int numRings; + int totRingLengths; + int maxRingLength; -/** - * Write out polygon file - */ -void writeOut(); + /** + * Entry routine. + */ + void trace(); + /** + * Trace the given node, find patches and fill them in + */ + Node *trace(Node *node, int *st, int len, int depth, PathList *& paths); + /** + * Look for path on the face and add to paths + */ + void findPaths(Node * node[2], int leaf[2], int depth[2], int *st[2], int maxdep, int dir, PathList * &paths); + /** + * Combine two list1 and list2 into list1 using connecting paths list3, + * while closed paths are appended to rings + */ + void combinePaths(PathList *& list1, PathList *list2, PathList *paths, PathList *& rings); + /** + * Helper function: combine current paths in list1 and list2 to a single path and append to list3 + */ + PathList *combineSinglePath(PathList *& head1, PathList *pre1, PathList *& list1, PathList *& head2, PathList *pre2, PathList *& list2); -void countIntersection(Node *node, int height, int& nedge, int& ncell, int& nface); -void generateMinimizer(Node * node, int st[3], int len, int height, int& offset); -void computeMinimizer(LeafNode * leaf, int st[3], int len, float rvalue[3]); -/** - * Traversal functions to generate polygon model - * op: 0 for counting, 1 for writing OBJ, 2 for writing OFF, 3 for writing PLY - */ -void cellProcContour(Node *node, int leaf, int depth); -void faceProcContour(Node * node[2], int leaf[2], int depth[2], int maxdep, int dir); -void edgeProcContour(Node * node[4], int leaf[4], int depth[4], int maxdep, int dir); -void processEdgeWrite(Node * node[4], int depths[4], int maxdep, int dir); - -/* output callbacks/data */ -DualConAllocOutput alloc_output; -DualConAddVert add_vert; -DualConAddQuad add_quad; -void *output_mesh; - -private: -/************ Operators for all nodes ************/ - -/// Lookup table -int numChildrenTable[256]; -int childrenCountTable[256][8]; -int childrenIndexTable[256][8]; -int numEdgeTable[8]; -int edgeCountTable[8][3]; - -/// Build up lookup table -void buildTable() -{ - for (int i = 0; i < 256; i++) + /** + * Functions to patch rings in a node + */ + Node *patch(Node * node, int st[3], int len, PathList * rings); + Node *patchSplit(Node * node, int st[3], int len, PathList * rings, int dir, PathList * &nrings1, PathList * &nrings2); + Node *patchSplitSingle(Node * node, int st[3], int len, PathElement * head, int dir, PathList * &nrings1, PathList * &nrings2); + Node *connectFace(Node * node, int st[3], int len, int dir, PathElement * f1, PathElement * f2); + Node *locateCell(InternalNode * node, int st[3], int len, int ori[3], int dir, int side, Node * &rleaf, int rst[3], int& rlen); + void compressRing(PathElement *& ring); + void getFacePoint(PathElement *leaf, int dir, int& x, int& y, float& p, float& q); + LeafNode *patchAdjacent(InternalNode * node, int len, int st1[3], LeafNode * leaf1, int st2[3], LeafNode * leaf2, int walkdir, int inc, int dir, int side, float alpha); + int findPair(PathElement *head, int pos, int dir, PathElement *& pre1, PathElement *& pre2); + int getSide(PathElement *e, int pos, int dir); + int isEqual(PathElement *e1, PathElement *e2); + void preparePrimalEdgesMask(InternalNode *node); + void testFacePoint(PathElement *e1, PathElement *e2); + + /** + * Path-related functions + */ + void deletePath(PathList *& head, PathList *pre, PathList *& curr); + void printPath(PathList *path); + void printPath(PathElement *path); + void printElement(PathElement *ele); + void printPaths(PathList *path); + void checkElement(PathElement *ele); + void checkPath(PathElement *path); + + + /** + * Routines to build signs to create a partitioned volume + *(after patching rings) + */ + void buildSigns(); + void buildSigns(unsigned char table[], Node * node, int isLeaf, int sg, int rvalue[8]); + + /************************************************************************/ + /* To remove disconnected components */ + /************************************************************************/ + void floodFill(); + void clearProcessBits(Node *node, int height); + int floodFill(LeafNode * leaf, int st[3], int len, int height, int threshold); + int floodFill(Node * node, int st[3], int len, int height, int threshold); + + /** + * Write out polygon file + */ + void writeOut(); + + void countIntersection(Node *node, int height, int& nedge, int& ncell, int& nface); + void generateMinimizer(Node * node, int st[3], int len, int height, int& offset); + void computeMinimizer(LeafNode * leaf, int st[3], int len, float rvalue[3]); + /** + * Traversal functions to generate polygon model + * op: 0 for counting, 1 for writing OBJ, 2 for writing OFF, 3 for writing PLY + */ + void cellProcContour(Node *node, int leaf, int depth); + void faceProcContour(Node * node[2], int leaf[2], int depth[2], int maxdep, int dir); + void edgeProcContour(Node * node[4], int leaf[4], int depth[4], int maxdep, int dir); + void processEdgeWrite(Node * node[4], int depths[4], int maxdep, int dir); + + /* output callbacks/data */ + DualConAllocOutput alloc_output; + DualConAddVert add_vert; + DualConAddQuad add_quad; + void *output_mesh; + + private: + /************ Operators for all nodes ************/ + + /// Lookup table + int numChildrenTable[256]; + int childrenCountTable[256][8]; + int childrenIndexTable[256][8]; + int numEdgeTable[8]; + int edgeCountTable[8][3]; + + /// Build up lookup table + void buildTable() { - numChildrenTable[i] = 0; - int count = 0; - for (int j = 0; j < 8; j++) - { - numChildrenTable[i] += ((i >> j) & 1); - childrenCountTable[i][j] = count; - childrenIndexTable[i][count] = j; - count += ((i >> j) & 1); + for (int i = 0; i < 256; i++) { + numChildrenTable[i] = 0; + int count = 0; + for (int j = 0; j < 8; j++) { + numChildrenTable[i] += ((i >> j) & 1); + childrenCountTable[i][j] = count; + childrenIndexTable[i][count] = j; + count += ((i >> j) & 1); + } } - } - for (int i = 0; i < 8; i++) - { - numEdgeTable[i] = 0; - int count = 0; - for (int j = 0; j < 3; j++) - { - numEdgeTable[i] += ((i >> j) & 1); - edgeCountTable[i][j] = count; - count += ((i >> j) & 1); + for (int i = 0; i < 8; i++) { + numEdgeTable[i] = 0; + int count = 0; + for (int j = 0; j < 3; j++) { + numEdgeTable[i] += ((i >> j) & 1); + edgeCountTable[i][j] = count; + count += ((i >> j) & 1); + } } } -} -int getSign(Node *node, int height, int index) -{ - if (height == 0) + int getSign(Node *node, int height, int index) { - return getSign(&node->leaf, index); - } - else { - if (hasChild(&node->internal, index)) - { - return getSign(getChild(&node->internal, getChildCount(&node->internal, index)), - height - 1, - index); + if (height == 0) { + return getSign(&node->leaf, index); } else { - return getSign(getChild(&node->internal, 0), - height - 1, - 7 - getChildIndex(&node->internal, 0)); + if (hasChild(&node->internal, index)) { + return getSign(getChild(&node->internal, getChildCount(&node->internal, index)), + height - 1, + index); + } + else { + return getSign(getChild(&node->internal, 0), + height - 1, + 7 - getChildIndex(&node->internal, 0)); + } } } -} -/************ Operators for leaf nodes ************/ + /************ Operators for leaf nodes ************/ -void printInfo(int st[3]) -{ - printf("INFO AT: %d %d %d\n", st[0] >> minshift, st[1] >> minshift, st[2] >> minshift); - LeafNode *leaf = (LeafNode *)locateLeafCheck(st); - if (leaf) - printInfo(leaf); - else - printf("Leaf not exists!\n"); -} - -void printInfo(const LeafNode *leaf) -{ - /* - printf("Edge mask: "); - for(int i = 0; i < 12; i ++) - { - printf("%d ", getEdgeParity(leaf, i)); - } - printf("\n"); - printf("Stored edge mask: "); - for(i = 0; i < 3; i ++) - { - printf("%d ", getStoredEdgesParity(leaf, i)); - } - printf("\n"); - */ - printf("Sign mask: "); - for (int i = 0; i < 8; i++) + void printInfo(int st[3]) { - printf("%d ", getSign(leaf, i)); + printf("INFO AT: %d %d %d\n", st[0] >> minshift, st[1] >> minshift, st[2] >> minshift); + LeafNode *leaf = (LeafNode *)locateLeafCheck(st); + if (leaf) + printInfo(leaf); + else + printf("Leaf not exists!\n"); } - printf("\n"); - -} -/// Retrieve signs -int getSign(const LeafNode *leaf, int index) -{ - return ((leaf->signs >> index) & 1); -} - -/// Set sign -void setSign(LeafNode *leaf, int index) -{ - leaf->signs |= (1 << index); -} + void printInfo(const LeafNode *leaf) + { + /* + printf("Edge mask: "); + for(int i = 0; i < 12; i ++) + { + printf("%d ", getEdgeParity(leaf, i)); + } + printf("\n"); + printf("Stored edge mask: "); + for(i = 0; i < 3; i ++) + { + printf("%d ", getStoredEdgesParity(leaf, i)); + } + printf("\n"); + */ + printf("Sign mask: "); + for (int i = 0; i < 8; i++) { + printf("%d ", getSign(leaf, i)); + } + printf("\n"); -void setSign(LeafNode *leaf, int index, int sign) -{ - leaf->signs &= (~(1 << index)); - leaf->signs |= ((sign & 1) << index); -} + } -int getSignMask(const LeafNode *leaf) -{ - return leaf->signs; -} + /// Retrieve signs + int getSign(const LeafNode *leaf, int index) + { + return ((leaf->signs >> index) & 1); + } -void setInProcessAll(int st[3], int dir) -{ - int nst[3], eind; - for (int i = 0; i < 4; i++) + /// Set sign + void setSign(LeafNode *leaf, int index) { - nst[0] = st[0] + dirCell[dir][i][0] * mindimen; - nst[1] = st[1] + dirCell[dir][i][1] * mindimen; - nst[2] = st[2] + dirCell[dir][i][2] * mindimen; - eind = dirEdge[dir][i]; + leaf->signs |= (1 << index); + } - LeafNode *cell = locateLeafCheck(nst); - assert(cell); + void setSign(LeafNode *leaf, int index, int sign) + { + leaf->signs &= (~(1 << index)); + leaf->signs |= ((sign & 1) << index); + } - setInProcess(cell, eind); + int getSignMask(const LeafNode *leaf) + { + return leaf->signs; } -} -void flipParityAll(int st[3], int dir) -{ - int nst[3], eind; - for (int i = 0; i < 4; i++) + void setInProcessAll(int st[3], int dir) { - nst[0] = st[0] + dirCell[dir][i][0] * mindimen; - nst[1] = st[1] + dirCell[dir][i][1] * mindimen; - nst[2] = st[2] + dirCell[dir][i][2] * mindimen; - eind = dirEdge[dir][i]; + int nst[3], eind; + for (int i = 0; i < 4; i++) { + nst[0] = st[0] + dirCell[dir][i][0] * mindimen; + nst[1] = st[1] + dirCell[dir][i][1] * mindimen; + nst[2] = st[2] + dirCell[dir][i][2] * mindimen; + eind = dirEdge[dir][i]; + + LeafNode *cell = locateLeafCheck(nst); + assert(cell); - LeafNode *cell = locateLeaf(nst); - flipEdge(cell, eind); + setInProcess(cell, eind); + } } -} -void setInProcess(LeafNode *leaf, int eind) -{ - assert(eind >= 0 && eind <= 11); + void flipParityAll(int st[3], int dir) + { + int nst[3], eind; + for (int i = 0; i < 4; i++) { + nst[0] = st[0] + dirCell[dir][i][0] * mindimen; + nst[1] = st[1] + dirCell[dir][i][1] * mindimen; + nst[2] = st[2] + dirCell[dir][i][2] * mindimen; + eind = dirEdge[dir][i]; + + LeafNode *cell = locateLeaf(nst); + flipEdge(cell, eind); + } + } - leaf->flood_fill |= (1 << eind); -} + void setInProcess(LeafNode *leaf, int eind) + { + assert(eind >= 0 && eind <= 11); -void setOutProcess(LeafNode *leaf, int eind) -{ - assert(eind >= 0 && eind <= 11); + leaf->flood_fill |= (1 << eind); + } - leaf->flood_fill &= ~(1 << eind); -} + void setOutProcess(LeafNode *leaf, int eind) + { + assert(eind >= 0 && eind <= 11); -int isInProcess(LeafNode *leaf, int eind) -{ - assert(eind >= 0 && eind <= 11); + leaf->flood_fill &= ~(1 << eind); + } - return (leaf->flood_fill >> eind) & 1; -} + int isInProcess(LeafNode *leaf, int eind) + { + assert(eind >= 0 && eind <= 11); -/// Generate signs at the corners from the edge parity -void generateSigns(LeafNode *leaf, unsigned char table[], int start) -{ - leaf->signs = table[leaf->edge_parity]; + return (leaf->flood_fill >> eind) & 1; + } - if ((start ^ leaf->signs) & 1) + /// Generate signs at the corners from the edge parity + void generateSigns(LeafNode *leaf, unsigned char table[], int start) { - leaf->signs = ~(leaf->signs); + leaf->signs = table[leaf->edge_parity]; + + if ((start ^ leaf->signs) & 1) { + leaf->signs = ~(leaf->signs); + } } -} -/// Get edge parity -int getEdgeParity(LeafNode *leaf, int index) -{ - assert(index >= 0 && index <= 11); + /// Get edge parity + int getEdgeParity(LeafNode *leaf, int index) + { + assert(index >= 0 && index <= 11); - return (leaf->edge_parity >> index) & 1; -} + return (leaf->edge_parity >> index) & 1; + } -/// Get edge parity on a face -int getFaceParity(LeafNode *leaf, int index) -{ - int a = getEdgeParity(leaf, faceMap[index][0]) + + /// Get edge parity on a face + int getFaceParity(LeafNode *leaf, int index) + { + int a = getEdgeParity(leaf, faceMap[index][0]) + getEdgeParity(leaf, faceMap[index][1]) + getEdgeParity(leaf, faceMap[index][2]) + getEdgeParity(leaf, faceMap[index][3]); - return (a & 1); -} -int getFaceEdgeNum(LeafNode *leaf, int index) -{ - int a = getEdgeParity(leaf, faceMap[index][0]) + + return (a & 1); + } + int getFaceEdgeNum(LeafNode *leaf, int index) + { + int a = getEdgeParity(leaf, faceMap[index][0]) + getEdgeParity(leaf, faceMap[index][1]) + getEdgeParity(leaf, faceMap[index][2]) + getEdgeParity(leaf, faceMap[index][3]); - return a; -} - -/// Set edge parity -void flipEdge(LeafNode *leaf, int index) -{ - assert(index >= 0 && index <= 11); - - leaf->edge_parity ^= (1 << index); -} - -/// Set 1 -void setEdge(LeafNode *leaf, int index) -{ - assert(index >= 0 && index <= 11); + return a; + } - leaf->edge_parity |= (1 << index); -} + /// Set edge parity + void flipEdge(LeafNode *leaf, int index) + { + assert(index >= 0 && index <= 11); -/// Set 0 -void resetEdge(LeafNode *leaf, int index) -{ - assert(index >= 0 && index <= 11); + leaf->edge_parity ^= (1 << index); + } - leaf->edge_parity &= ~(1 << index); -} + /// Set 1 + void setEdge(LeafNode *leaf, int index) + { + assert(index >= 0 && index <= 11); -/// Flipping with a new intersection offset -void createPrimalEdgesMask(LeafNode *leaf) -{ - leaf->primary_edge_intersections = getPrimalEdgesMask2(leaf); -} + leaf->edge_parity |= (1 << index); + } -void setStoredEdgesParity(LeafNode *leaf, int pindex) -{ - assert(pindex <= 2 && pindex >= 0); + /// Set 0 + void resetEdge(LeafNode *leaf, int index) + { + assert(index >= 0 && index <= 11); - leaf->primary_edge_intersections |= (1 << pindex); -} -int getStoredEdgesParity(LeafNode *leaf, int pindex) -{ - assert(pindex <= 2 && pindex >= 0); + leaf->edge_parity &= ~(1 << index); + } - return (leaf->primary_edge_intersections >> pindex) & 1; -} + /// Flipping with a new intersection offset + void createPrimalEdgesMask(LeafNode *leaf) + { + leaf->primary_edge_intersections = getPrimalEdgesMask2(leaf); + } -LeafNode *flipEdge(LeafNode *leaf, int index, float alpha) -{ - flipEdge(leaf, index); + void setStoredEdgesParity(LeafNode *leaf, int pindex) + { + assert(pindex <= 2 && pindex >= 0); - if ((index & 3) == 0) + leaf->primary_edge_intersections |= (1 << pindex); + } + int getStoredEdgesParity(LeafNode *leaf, int pindex) { - int ind = index / 4; - if (getEdgeParity(leaf, index) && !getStoredEdgesParity(leaf, ind)) - { - // Create a new node - int num = getNumEdges(leaf) + 1; - setStoredEdgesParity(leaf, ind); - int count = getEdgeCount(leaf, ind); - LeafNode *nleaf = createLeaf(num); - *nleaf = *leaf; + assert(pindex <= 2 && pindex >= 0); - setEdgeOffset(nleaf, alpha, count); + return (leaf->primary_edge_intersections >> pindex) & 1; + } - if (num > 1) - { - float *pts = leaf->edge_intersections; - float *npts = nleaf->edge_intersections; - for (int i = 0; i < count; i++) - { - for (int j = 0; j < EDGE_FLOATS; j++) - { - npts[i * EDGE_FLOATS + j] = pts[i * EDGE_FLOATS + j]; + LeafNode *flipEdge(LeafNode *leaf, int index, float alpha) + { + flipEdge(leaf, index); + + if ((index & 3) == 0) { + int ind = index / 4; + if (getEdgeParity(leaf, index) && !getStoredEdgesParity(leaf, ind)) { + // Create a new node + int num = getNumEdges(leaf) + 1; + setStoredEdgesParity(leaf, ind); + int count = getEdgeCount(leaf, ind); + LeafNode *nleaf = createLeaf(num); + *nleaf = *leaf; + + setEdgeOffset(nleaf, alpha, count); + + if (num > 1) { + float *pts = leaf->edge_intersections; + float *npts = nleaf->edge_intersections; + for (int i = 0; i < count; i++) { + for (int j = 0; j < EDGE_FLOATS; j++) { + npts[i * EDGE_FLOATS + j] = pts[i * EDGE_FLOATS + j]; + } } - } - for (int i = count + 1; i < num; i++) - { - for (int j = 0; j < EDGE_FLOATS; j++) - { - npts[i * EDGE_FLOATS + j] = pts[(i - 1) * EDGE_FLOATS + j]; + for (int i = count + 1; i < num; i++) { + for (int j = 0; j < EDGE_FLOATS; j++) { + npts[i * EDGE_FLOATS + j] = pts[(i - 1) * EDGE_FLOATS + j]; + } } } - } - removeLeaf(num - 1, (LeafNode *)leaf); - leaf = nleaf; + removeLeaf(num - 1, (LeafNode *)leaf); + leaf = nleaf; + } } - } - return leaf; -} + return leaf; + } -/// Update parent link -void updateParent(InternalNode *node, int len, int st[3], LeafNode *leaf) -{ - // First, locate the parent - int count; - InternalNode *parent = locateParent(node, len, st, count); + /// Update parent link + void updateParent(InternalNode *node, int len, int st[3], LeafNode *leaf) + { + // First, locate the parent + int count; + InternalNode *parent = locateParent(node, len, st, count); - // Update - setChild(parent, count, (Node *)leaf); -} + // Update + setChild(parent, count, (Node *)leaf); + } -void updateParent(InternalNode *node, int len, int st[3]) -{ - if (len == dimen) + void updateParent(InternalNode *node, int len, int st[3]) { - root = (Node *)node; - return; - } + if (len == dimen) { + root = (Node *)node; + return; + } - // First, locate the parent - int count; - InternalNode *parent = locateParent(len, st, count); + // First, locate the parent + int count; + InternalNode *parent = locateParent(len, st, count); - // UPdate - setChild(parent, count, (Node *)node); -} + // UPdate + setChild(parent, count, (Node *)node); + } -/// Find edge intersection on a given edge -int getEdgeIntersectionByIndex(int st[3], int index, float pt[3], int check) -{ - // First, locat the leaf - LeafNode *leaf; - if (check) + /// Find edge intersection on a given edge + int getEdgeIntersectionByIndex(int st[3], int index, float pt[3], int check) { - leaf = locateLeafCheck(st); + // First, locat the leaf + LeafNode *leaf; + if (check) { + leaf = locateLeafCheck(st); + } + else { + leaf = locateLeaf(st); + } + + if (leaf && getStoredEdgesParity(leaf, index)) { + float off = getEdgeOffset(leaf, getEdgeCount(leaf, index)); + pt[0] = (float) st[0]; + pt[1] = (float) st[1]; + pt[2] = (float) st[2]; + pt[index] += off * mindimen; + + return 1; + } + else { + return 0; + } } - else { - leaf = locateLeaf(st); + + /// Retrieve number of edges intersected + int getPrimalEdgesMask(LeafNode *leaf) + { + return leaf->primary_edge_intersections; } - if (leaf && getStoredEdgesParity(leaf, index)) + int getPrimalEdgesMask2(LeafNode *leaf) { - float off = getEdgeOffset(leaf, getEdgeCount(leaf, index)); - pt[0] = (float) st[0]; - pt[1] = (float) st[1]; - pt[2] = (float) st[2]; - pt[index] += off * mindimen; + return (((leaf->edge_parity & 0x1) >> 0) | + ((leaf->edge_parity & 0x10) >> 3) | + ((leaf->edge_parity & 0x100) >> 6)); + } - return 1; + /// Get the count for a primary edge + int getEdgeCount(LeafNode *leaf, int index) + { + return edgeCountTable[getPrimalEdgesMask(leaf)][index]; } - else { - return 0; + int getNumEdges(LeafNode *leaf) + { + return numEdgeTable[getPrimalEdgesMask(leaf)]; } -} - -/// Retrieve number of edges intersected -int getPrimalEdgesMask(LeafNode *leaf) -{ - return leaf->primary_edge_intersections; -} - -int getPrimalEdgesMask2(LeafNode *leaf) -{ - return (((leaf->edge_parity & 0x1) >> 0) | - ((leaf->edge_parity & 0x10) >> 3) | - ((leaf->edge_parity & 0x100) >> 6)); -} - -/// Get the count for a primary edge -int getEdgeCount(LeafNode *leaf, int index) -{ - return edgeCountTable[getPrimalEdgesMask(leaf)][index]; -} -int getNumEdges(LeafNode *leaf) -{ - return numEdgeTable[getPrimalEdgesMask(leaf)]; -} - -int getNumEdges2(LeafNode *leaf) -{ - return numEdgeTable[getPrimalEdgesMask2(leaf)]; -} -/// Set edge intersection -void setEdgeOffset(LeafNode *leaf, float pt, int count) -{ - float *pts = leaf->edge_intersections; - pts[EDGE_FLOATS * count] = pt; - pts[EDGE_FLOATS * count + 1] = 0; - pts[EDGE_FLOATS * count + 2] = 0; - pts[EDGE_FLOATS * count + 3] = 0; -} - -/// Set multiple edge intersections -void setEdgeOffsets(LeafNode *leaf, float pt[3], int len) -{ - float *pts = leaf->edge_intersections; - for (int i = 0; i < len; i++) + int getNumEdges2(LeafNode *leaf) { - pts[i] = pt[i]; + return numEdgeTable[getPrimalEdgesMask2(leaf)]; } -} -/// Retrieve edge intersection -float getEdgeOffset(LeafNode *leaf, int count) -{ - return leaf->edge_intersections[4 * count]; -} - -/// Update method -LeafNode *updateEdgeOffsets(LeafNode *leaf, int oldlen, int newlen, float offs[3]) -{ - // First, create a new leaf node - LeafNode *nleaf = createLeaf(newlen); - *nleaf = *leaf; + /// Set edge intersection + void setEdgeOffset(LeafNode *leaf, float pt, int count) + { + float *pts = leaf->edge_intersections; + pts[EDGE_FLOATS * count] = pt; + pts[EDGE_FLOATS * count + 1] = 0; + pts[EDGE_FLOATS * count + 2] = 0; + pts[EDGE_FLOATS * count + 3] = 0; + } - // Next, fill in the offsets - setEdgeOffsets(nleaf, offs, newlen); + /// Set multiple edge intersections + void setEdgeOffsets(LeafNode *leaf, float pt[3], int len) + { + float *pts = leaf->edge_intersections; + for (int i = 0; i < len; i++) { + pts[i] = pt[i]; + } + } - // Finally, delete the old leaf - removeLeaf(oldlen, leaf); + /// Retrieve edge intersection + float getEdgeOffset(LeafNode *leaf, int count) + { + return leaf->edge_intersections[4 * count]; + } - return nleaf; -} + /// Update method + LeafNode *updateEdgeOffsets(LeafNode *leaf, int oldlen, int newlen, float offs[3]) + { + // First, create a new leaf node + LeafNode *nleaf = createLeaf(newlen); + *nleaf = *leaf; -/// Set minimizer index -void setMinimizerIndex(LeafNode *leaf, int index) -{ - leaf->minimizer_index = index; -} + // Next, fill in the offsets + setEdgeOffsets(nleaf, offs, newlen); -/// Get minimizer index -int getMinimizerIndex(LeafNode *leaf) -{ - return leaf->minimizer_index; -} + // Finally, delete the old leaf + removeLeaf(oldlen, leaf); -int getMinimizerIndex(LeafNode *leaf, int eind) -{ - int add = manifold_table[getSignMask(leaf)].pairs[eind][0] - 1; - assert(add >= 0); - return leaf->minimizer_index + add; -} + return nleaf; + } -void getMinimizerIndices(LeafNode *leaf, int eind, int inds[2]) -{ - const int *add = manifold_table[getSignMask(leaf)].pairs[eind]; - inds[0] = leaf->minimizer_index + add[0] - 1; - if (add[0] == add[1]) + /// Set minimizer index + void setMinimizerIndex(LeafNode *leaf, int index) { - inds[1] = -1; + leaf->minimizer_index = index; } - else { - inds[1] = leaf->minimizer_index + add[1] - 1; + + /// Get minimizer index + int getMinimizerIndex(LeafNode *leaf) + { + return leaf->minimizer_index; } -} + int getMinimizerIndex(LeafNode *leaf, int eind) + { + int add = manifold_table[getSignMask(leaf)].pairs[eind][0] - 1; + assert(add >= 0); + return leaf->minimizer_index + add; + } -/// Set edge intersection -void setEdgeOffsetNormal(LeafNode *leaf, float pt, float a, float b, float c, int count) -{ - float *pts = leaf->edge_intersections; - pts[4 * count] = pt; - pts[4 * count + 1] = a; - pts[4 * count + 2] = b; - pts[4 * count + 3] = c; -} - -float getEdgeOffsetNormal(LeafNode *leaf, int count, float& a, float& b, float& c) -{ - float *pts = leaf->edge_intersections; - a = pts[4 * count + 1]; - b = pts[4 * count + 2]; - c = pts[4 * count + 3]; - return pts[4 * count]; -} - -/// Set multiple edge intersections -void setEdgeOffsetsNormals(LeafNode *leaf, float pt[], float a[], float b[], float c[], int len) -{ - float *pts = leaf->edge_intersections; - for (int i = 0; i < len; i++) + void getMinimizerIndices(LeafNode *leaf, int eind, int inds[2]) { - if (pt[i] > 1 || pt[i] < 0) - { - printf("\noffset: %f\n", pt[i]); + const int *add = manifold_table[getSignMask(leaf)].pairs[eind]; + inds[0] = leaf->minimizer_index + add[0] - 1; + if (add[0] == add[1]) { + inds[1] = -1; + } + else { + inds[1] = leaf->minimizer_index + add[1] - 1; } - pts[i * 4] = pt[i]; - pts[i * 4 + 1] = a[i]; - pts[i * 4 + 2] = b[i]; - pts[i * 4 + 3] = c[i]; } -} - -/// Retrieve complete edge intersection -void getEdgeIntersectionByIndex(LeafNode *leaf, int index, int st[3], int len, float pt[3], float nm[3]) -{ - int count = getEdgeCount(leaf, index); - float *pts = leaf->edge_intersections; - - float off = pts[4 * count]; - - pt[0] = (float) st[0]; - pt[1] = (float) st[1]; - pt[2] = (float) st[2]; - pt[index] += (off * len); - - nm[0] = pts[4 * count + 1]; - nm[1] = pts[4 * count + 2]; - nm[2] = pts[4 * count + 3]; -} - -float getEdgeOffsetNormalByIndex(LeafNode *leaf, int index, float nm[3]) -{ - int count = getEdgeCount(leaf, index); - float *pts = leaf->edge_intersections; - - float off = pts[4 * count]; - - nm[0] = pts[4 * count + 1]; - nm[1] = pts[4 * count + 2]; - nm[2] = pts[4 * count + 3]; - return off; -} -void fillEdgeIntersections(LeafNode *leaf, int st[3], int len, float pts[12][3], float norms[12][3]) -{ - int i; - // int stt[3] = {0, 0, 0}; - - // The three primal edges are easy - int pmask[3] = {0, 4, 8}; - for (i = 0; i < 3; i++) + /// Set edge intersection + void setEdgeOffsetNormal(LeafNode *leaf, float pt, float a, float b, float c, int count) { - if (getEdgeParity(leaf, pmask[i])) - { - // getEdgeIntersectionByIndex(leaf, i, stt, 1, pts[pmask[i]], norms[pmask[i]]); - getEdgeIntersectionByIndex(leaf, i, st, len, pts[pmask[i]], norms[pmask[i]]); - } + float *pts = leaf->edge_intersections; + pts[4 * count] = pt; + pts[4 * count + 1] = a; + pts[4 * count + 2] = b; + pts[4 * count + 3] = c; } - // 3 face adjacent cubes - int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}}; - int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}}; - for (i = 0; i < 3; i++) + float getEdgeOffsetNormal(LeafNode *leaf, int count, float& a, float& b, float& c) { - int e1 = getEdgeParity(leaf, fmask[i][0]); - int e2 = getEdgeParity(leaf, fmask[i][1]); - if (e1 || e2) - { - int nst[3] = {st[0], st[1], st[2]}; - nst[i] += len; - // int nstt[3] = {0, 0, 0}; - // nstt[i] += 1; - LeafNode *node = locateLeaf(nst); + float *pts = leaf->edge_intersections; + a = pts[4 * count + 1]; + b = pts[4 * count + 2]; + c = pts[4 * count + 3]; + return pts[4 * count]; + } - if (e1) - { - // getEdgeIntersectionByIndex(node, femask[i][0], nstt, 1, pts[fmask[i][0]], norms[fmask[i][0]]); - getEdgeIntersectionByIndex(node, femask[i][0], nst, len, pts[fmask[i][0]], norms[fmask[i][0]]); - } - if (e2) - { - // getEdgeIntersectionByIndex(node, femask[i][1], nstt, 1, pts[fmask[i][1]], norms[fmask[i][1]]); - getEdgeIntersectionByIndex(node, femask[i][1], nst, len, pts[fmask[i][1]], norms[fmask[i][1]]); + /// Set multiple edge intersections + void setEdgeOffsetsNormals(LeafNode *leaf, const float pt[], + const float a[], const float b[], + const float c[], int len) + { + float *pts = leaf->edge_intersections; + for (int i = 0; i < len; i++) { + if (pt[i] > 1 || pt[i] < 0) { + printf("\noffset: %f\n", pt[i]); } + pts[i * 4] = pt[i]; + pts[i * 4 + 1] = a[i]; + pts[i * 4 + 2] = b[i]; + pts[i * 4 + 3] = c[i]; } } - // 3 edge adjacent cubes - int emask[3] = {3, 7, 11}; - int eemask[3] = {0, 1, 2}; - for (i = 0; i < 3; i++) + /// Retrieve complete edge intersection + void getEdgeIntersectionByIndex(LeafNode *leaf, int index, int st[3], int len, float pt[3], float nm[3]) { - if (getEdgeParity(leaf, emask[i])) - { - int nst[3] = {st[0] + len, st[1] + len, st[2] + len}; - nst[i] -= len; - // int nstt[3] = {1, 1, 1}; - // nstt[i] -= 1; - LeafNode *node = locateLeaf(nst); + int count = getEdgeCount(leaf, index); + float *pts = leaf->edge_intersections; - // getEdgeIntersectionByIndex(node, eemask[i], nstt, 1, pts[emask[i]], norms[emask[i]]); - getEdgeIntersectionByIndex(node, eemask[i], nst, len, pts[emask[i]], norms[emask[i]]); - } - } -} + float off = pts[4 * count]; + pt[0] = (float) st[0]; + pt[1] = (float) st[1]; + pt[2] = (float) st[2]; + pt[index] += (off * len); -void fillEdgeIntersections(LeafNode *leaf, int st[3], int len, float pts[12][3], float norms[12][3], int parity[12]) -{ - int i; - for (i = 0; i < 12; i++) - { - parity[i] = 0; + nm[0] = pts[4 * count + 1]; + nm[1] = pts[4 * count + 2]; + nm[2] = pts[4 * count + 3]; } - // int stt[3] = {0, 0, 0}; - // The three primal edges are easy - int pmask[3] = {0, 4, 8}; - for (i = 0; i < 3; i++) + float getEdgeOffsetNormalByIndex(LeafNode *leaf, int index, float nm[3]) { - if (getStoredEdgesParity(leaf, i)) - { - // getEdgeIntersectionByIndex(leaf, i, stt, 1, pts[pmask[i]], norms[pmask[i]]); - getEdgeIntersectionByIndex(leaf, i, st, len, pts[pmask[i]], norms[pmask[i]]); - parity[pmask[i]] = 1; - } + int count = getEdgeCount(leaf, index); + float *pts = leaf->edge_intersections; + + float off = pts[4 * count]; + + nm[0] = pts[4 * count + 1]; + nm[1] = pts[4 * count + 2]; + nm[2] = pts[4 * count + 3]; + + return off; } - // 3 face adjacent cubes - int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}}; - int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}}; - for (i = 0; i < 3; i++) + void fillEdgeIntersections(LeafNode *leaf, int st[3], int len, float pts[12][3], float norms[12][3]) { - { - int nst[3] = {st[0], st[1], st[2]}; - nst[i] += len; - // int nstt[3] = {0, 0, 0}; - // nstt[i] += 1; - LeafNode *node = locateLeafCheck(nst); - if (node == NULL) - { - continue; + int i; + // int stt[3] = {0, 0, 0}; + + // The three primal edges are easy + int pmask[3] = {0, 4, 8}; + for (i = 0; i < 3; i++) { + if (getEdgeParity(leaf, pmask[i])) { + // getEdgeIntersectionByIndex(leaf, i, stt, 1, pts[pmask[i]], norms[pmask[i]]); + getEdgeIntersectionByIndex(leaf, i, st, len, pts[pmask[i]], norms[pmask[i]]); } + } - int e1 = getStoredEdgesParity(node, femask[i][0]); - int e2 = getStoredEdgesParity(node, femask[i][1]); - - if (e1) - { - // getEdgeIntersectionByIndex(node, femask[i][0], nstt, 1, pts[fmask[i][0]], norms[fmask[i][0]]); - getEdgeIntersectionByIndex(node, femask[i][0], nst, len, pts[fmask[i][0]], norms[fmask[i][0]]); - parity[fmask[i][0]] = 1; - } - if (e2) - { - // getEdgeIntersectionByIndex(node, femask[i][1], nstt, 1, pts[fmask[i][1]], norms[fmask[i][1]]); - getEdgeIntersectionByIndex(node, femask[i][1], nst, len, pts[fmask[i][1]], norms[fmask[i][1]]); - parity[fmask[i][1]] = 1; + // 3 face adjacent cubes + int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}}; + int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}}; + for (i = 0; i < 3; i++) { + int e1 = getEdgeParity(leaf, fmask[i][0]); + int e2 = getEdgeParity(leaf, fmask[i][1]); + if (e1 || e2) { + int nst[3] = {st[0], st[1], st[2]}; + nst[i] += len; + // int nstt[3] = {0, 0, 0}; + // nstt[i] += 1; + LeafNode *node = locateLeaf(nst); + + if (e1) { + // getEdgeIntersectionByIndex(node, femask[i][0], nstt, 1, pts[fmask[i][0]], norms[fmask[i][0]]); + getEdgeIntersectionByIndex(node, femask[i][0], nst, len, pts[fmask[i][0]], norms[fmask[i][0]]); + } + if (e2) { + // getEdgeIntersectionByIndex(node, femask[i][1], nstt, 1, pts[fmask[i][1]], norms[fmask[i][1]]); + getEdgeIntersectionByIndex(node, femask[i][1], nst, len, pts[fmask[i][1]], norms[fmask[i][1]]); + } } } - } - // 3 edge adjacent cubes - int emask[3] = {3, 7, 11}; - int eemask[3] = {0, 1, 2}; - for (i = 0; i < 3; i++) - { -// if(getEdgeParity(leaf, emask[i])) - { - int nst[3] = {st[0] + len, st[1] + len, st[2] + len}; - nst[i] -= len; - // int nstt[3] = {1, 1, 1}; - // nstt[i] -= 1; - LeafNode *node = locateLeafCheck(nst); - if (node == NULL) - { - continue; - } + // 3 edge adjacent cubes + int emask[3] = {3, 7, 11}; + int eemask[3] = {0, 1, 2}; + for (i = 0; i < 3; i++) { + if (getEdgeParity(leaf, emask[i])) { + int nst[3] = {st[0] + len, st[1] + len, st[2] + len}; + nst[i] -= len; + // int nstt[3] = {1, 1, 1}; + // nstt[i] -= 1; + LeafNode *node = locateLeaf(nst); - if (getStoredEdgesParity(node, eemask[i])) - { // getEdgeIntersectionByIndex(node, eemask[i], nstt, 1, pts[emask[i]], norms[emask[i]]); getEdgeIntersectionByIndex(node, eemask[i], nst, len, pts[emask[i]], norms[emask[i]]); - parity[emask[i]] = 1; } } } -} -void fillEdgeOffsetsNormals(LeafNode *leaf, int st[3], int len, float pts[12], float norms[12][3], int parity[12]) -{ - int i; - for (i = 0; i < 12; i++) - { - parity[i] = 0; - } - // int stt[3] = {0, 0, 0}; - // The three primal edges are easy - int pmask[3] = {0, 4, 8}; - for (i = 0; i < 3; i++) + void fillEdgeIntersections(LeafNode *leaf, int st[3], int len, float pts[12][3], float norms[12][3], int parity[12]) { - if (getStoredEdgesParity(leaf, i)) - { - pts[pmask[i]] = getEdgeOffsetNormalByIndex(leaf, i, norms[pmask[i]]); - parity[pmask[i]] = 1; + int i; + for (i = 0; i < 12; i++) { + parity[i] = 0; + } + // int stt[3] = {0, 0, 0}; + + // The three primal edges are easy + int pmask[3] = {0, 4, 8}; + for (i = 0; i < 3; i++) { + if (getStoredEdgesParity(leaf, i)) { + // getEdgeIntersectionByIndex(leaf, i, stt, 1, pts[pmask[i]], norms[pmask[i]]); + getEdgeIntersectionByIndex(leaf, i, st, len, pts[pmask[i]], norms[pmask[i]]); + parity[pmask[i]] = 1; + } } - } - // 3 face adjacent cubes - int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}}; - int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}}; - for (i = 0; i < 3; i++) - { - { - int nst[3] = {st[0], st[1], st[2]}; - nst[i] += len; - // int nstt[3] = {0, 0, 0}; - // nstt[i] += 1; - LeafNode *node = locateLeafCheck(nst); - if (node == NULL) + // 3 face adjacent cubes + int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}}; + int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}}; + for (i = 0; i < 3; i++) { { - continue; - } + int nst[3] = {st[0], st[1], st[2]}; + nst[i] += len; + // int nstt[3] = {0, 0, 0}; + // nstt[i] += 1; + LeafNode *node = locateLeafCheck(nst); + if (node == NULL) { + continue; + } - int e1 = getStoredEdgesParity(node, femask[i][0]); - int e2 = getStoredEdgesParity(node, femask[i][1]); + int e1 = getStoredEdgesParity(node, femask[i][0]); + int e2 = getStoredEdgesParity(node, femask[i][1]); - if (e1) - { - pts[fmask[i][0]] = getEdgeOffsetNormalByIndex(node, femask[i][0], norms[fmask[i][0]]); - parity[fmask[i][0]] = 1; + if (e1) { + // getEdgeIntersectionByIndex(node, femask[i][0], nstt, 1, pts[fmask[i][0]], norms[fmask[i][0]]); + getEdgeIntersectionByIndex(node, femask[i][0], nst, len, pts[fmask[i][0]], norms[fmask[i][0]]); + parity[fmask[i][0]] = 1; + } + if (e2) { + // getEdgeIntersectionByIndex(node, femask[i][1], nstt, 1, pts[fmask[i][1]], norms[fmask[i][1]]); + getEdgeIntersectionByIndex(node, femask[i][1], nst, len, pts[fmask[i][1]], norms[fmask[i][1]]); + parity[fmask[i][1]] = 1; + } } - if (e2) + } + + // 3 edge adjacent cubes + int emask[3] = {3, 7, 11}; + int eemask[3] = {0, 1, 2}; + for (i = 0; i < 3; i++) { + // if(getEdgeParity(leaf, emask[i])) { - pts[fmask[i][1]] = getEdgeOffsetNormalByIndex(node, femask[i][1], norms[fmask[i][1]]); - parity[fmask[i][1]] = 1; + int nst[3] = {st[0] + len, st[1] + len, st[2] + len}; + nst[i] -= len; + // int nstt[3] = {1, 1, 1}; + // nstt[i] -= 1; + LeafNode *node = locateLeafCheck(nst); + if (node == NULL) { + continue; + } + + if (getStoredEdgesParity(node, eemask[i])) { + // getEdgeIntersectionByIndex(node, eemask[i], nstt, 1, pts[emask[i]], norms[emask[i]]); + getEdgeIntersectionByIndex(node, eemask[i], nst, len, pts[emask[i]], norms[emask[i]]); + parity[emask[i]] = 1; + } } } } - // 3 edge adjacent cubes - int emask[3] = {3, 7, 11}; - int eemask[3] = {0, 1, 2}; - for (i = 0; i < 3; i++) + void fillEdgeOffsetsNormals(LeafNode *leaf, int st[3], int len, float pts[12], float norms[12][3], int parity[12]) { -// if(getEdgeParity(leaf, emask[i])) - { - int nst[3] = {st[0] + len, st[1] + len, st[2] + len}; - nst[i] -= len; - // int nstt[3] = {1, 1, 1}; - // nstt[i] -= 1; - LeafNode *node = locateLeafCheck(nst); - if (node == NULL) - { - continue; + int i; + for (i = 0; i < 12; i++) { + parity[i] = 0; + } + // int stt[3] = {0, 0, 0}; + + // The three primal edges are easy + int pmask[3] = {0, 4, 8}; + for (i = 0; i < 3; i++) { + if (getStoredEdgesParity(leaf, i)) { + pts[pmask[i]] = getEdgeOffsetNormalByIndex(leaf, i, norms[pmask[i]]); + parity[pmask[i]] = 1; } + } - if (getStoredEdgesParity(node, eemask[i])) + // 3 face adjacent cubes + int fmask[3][2] = {{6, 10}, {2, 9}, {1, 5}}; + int femask[3][2] = {{1, 2}, {0, 2}, {0, 1}}; + for (i = 0; i < 3; i++) { { - pts[emask[i]] = getEdgeOffsetNormalByIndex(node, eemask[i], norms[emask[i]]); - parity[emask[i]] = 1; + int nst[3] = {st[0], st[1], st[2]}; + nst[i] += len; + // int nstt[3] = {0, 0, 0}; + // nstt[i] += 1; + LeafNode *node = locateLeafCheck(nst); + if (node == NULL) { + continue; + } + + int e1 = getStoredEdgesParity(node, femask[i][0]); + int e2 = getStoredEdgesParity(node, femask[i][1]); + + if (e1) { + pts[fmask[i][0]] = getEdgeOffsetNormalByIndex(node, femask[i][0], norms[fmask[i][0]]); + parity[fmask[i][0]] = 1; + } + if (e2) { + pts[fmask[i][1]] = getEdgeOffsetNormalByIndex(node, femask[i][1], norms[fmask[i][1]]); + parity[fmask[i][1]] = 1; + } } } - } -} + // 3 edge adjacent cubes + int emask[3] = {3, 7, 11}; + int eemask[3] = {0, 1, 2}; + for (i = 0; i < 3; i++) { + // if(getEdgeParity(leaf, emask[i])) + { + int nst[3] = {st[0] + len, st[1] + len, st[2] + len}; + nst[i] -= len; + // int nstt[3] = {1, 1, 1}; + // nstt[i] -= 1; + LeafNode *node = locateLeafCheck(nst); + if (node == NULL) { + continue; + } -/// Update method -LeafNode *updateEdgeOffsetsNormals(LeafNode *leaf, int oldlen, int newlen, float offs[3], float a[3], float b[3], float c[3]) -{ - // First, create a new leaf node - LeafNode *nleaf = createLeaf(newlen); - *nleaf = *leaf; + if (getStoredEdgesParity(node, eemask[i])) { + pts[emask[i]] = getEdgeOffsetNormalByIndex(node, eemask[i], norms[emask[i]]); + parity[emask[i]] = 1; + } + } + } + } - // Next, fill in the offsets - setEdgeOffsetsNormals(nleaf, offs, a, b, c, newlen); - // Finally, delete the old leaf - removeLeaf(oldlen, leaf); + /// Update method + LeafNode *updateEdgeOffsetsNormals(LeafNode *leaf, int oldlen, int newlen, float offs[3], float a[3], float b[3], float c[3]) + { + // First, create a new leaf node + LeafNode *nleaf = createLeaf(newlen); + *nleaf = *leaf; - return nleaf; -} + // Next, fill in the offsets + setEdgeOffsetsNormals(nleaf, offs, a, b, c, newlen); -/// Locate a leaf -/// WARNING: assuming this leaf already exists! + // Finally, delete the old leaf + removeLeaf(oldlen, leaf); -LeafNode *locateLeaf(int st[3]) -{ - Node *node = (Node *)root; - for (int i = GRID_DIMENSION - 1; i > GRID_DIMENSION - maxDepth - 1; i--) - { - int index = (((st[0] >> i) & 1) << 2) | - (((st[1] >> i) & 1) << 1) | - (((st[2] >> i) & 1)); - node = getChild(&node->internal, getChildCount(&node->internal, index)); + return nleaf; } - return &node->leaf; -} + /// Locate a leaf + /// WARNING: assuming this leaf already exists! -LeafNode *locateLeaf(InternalNode *parent, int len, int st[3]) -{ - Node *node = (Node *)parent; - int index; - for (int i = len / 2; i >= mindimen; i >>= 1) + LeafNode *locateLeaf(int st[3]) { - index = (((st[0] & i) ? 4 : 0) | - ((st[1] & i) ? 2 : 0) | - ((st[2] & i) ? 1 : 0)); - node = getChild(&node->internal, - getChildCount(&node->internal, index)); - } + Node *node = (Node *)root; + for (int i = GRID_DIMENSION - 1; i > GRID_DIMENSION - maxDepth - 1; i--) { + int index = (((st[0] >> i) & 1) << 2) | + (((st[1] >> i) & 1) << 1) | + (((st[2] >> i) & 1)); + node = getChild(&node->internal, getChildCount(&node->internal, index)); + } - return &node->leaf; -} + return &node->leaf; + } -LeafNode *locateLeafCheck(int st[3]) -{ - Node *node = (Node *)root; - for (int i = GRID_DIMENSION - 1; i > GRID_DIMENSION - maxDepth - 1; i--) - { - int index = (((st[0] >> i) & 1) << 2) | - (((st[1] >> i) & 1) << 1) | - (((st[2] >> i) & 1)); - if (!hasChild(&node->internal, index)) - { - return NULL; + LeafNode *locateLeaf(InternalNode *parent, int len, int st[3]) + { + Node *node = (Node *)parent; + int index; + for (int i = len / 2; i >= mindimen; i >>= 1) { + index = (((st[0] & i) ? 4 : 0) | + ((st[1] & i) ? 2 : 0) | + ((st[2] & i) ? 1 : 0)); + node = getChild(&node->internal, + getChildCount(&node->internal, index)); } - node = getChild(&node->internal, getChildCount(&node->internal, index)); - } - return &node->leaf; -} + return &node->leaf; + } -InternalNode *locateParent(int len, int st[3], int& count) -{ - InternalNode *node = (InternalNode *)root; - InternalNode *pre = NULL; - int index = 0; - for (int i = dimen / 2; i >= len; i >>= 1) + LeafNode *locateLeafCheck(int st[3]) { - index = (((st[0] & i) ? 4 : 0) | - ((st[1] & i) ? 2 : 0) | - ((st[2] & i) ? 1 : 0)); - pre = node; - node = &getChild(node, getChildCount(node, index))->internal; - } + Node *node = (Node *)root; + for (int i = GRID_DIMENSION - 1; i > GRID_DIMENSION - maxDepth - 1; i--) { + int index = (((st[0] >> i) & 1) << 2) | + (((st[1] >> i) & 1) << 1) | + (((st[2] >> i) & 1)); + if (!hasChild(&node->internal, index)) { + return NULL; + } + node = getChild(&node->internal, getChildCount(&node->internal, index)); + } - count = getChildCount(pre, index); - return pre; -} + return &node->leaf; + } -InternalNode *locateParent(InternalNode *parent, int len, int st[3], int& count) -{ - InternalNode *node = parent; - InternalNode *pre = NULL; - int index = 0; - for (int i = len / 2; i >= mindimen; i >>= 1) + InternalNode *locateParent(int len, int st[3], int& count) { - index = (((st[0] & i) ? 4 : 0) | - ((st[1] & i) ? 2 : 0) | - ((st[2] & i) ? 1 : 0)); - pre = node; - node = (InternalNode *)getChild(node, getChildCount(node, index)); - } + InternalNode *node = (InternalNode *)root; + InternalNode *pre = NULL; + int index = 0; + for (int i = dimen / 2; i >= len; i >>= 1) { + index = (((st[0] & i) ? 4 : 0) | + ((st[1] & i) ? 2 : 0) | + ((st[2] & i) ? 1 : 0)); + pre = node; + node = &getChild(node, getChildCount(node, index))->internal; + } - count = getChildCount(pre, index); - return pre; -} + count = getChildCount(pre, index); + return pre; + } -/************ Operators for internal nodes ************/ + InternalNode *locateParent(InternalNode *parent, int len, int st[3], int& count) + { + InternalNode *node = parent; + InternalNode *pre = NULL; + int index = 0; + for (int i = len / 2; i >= mindimen; i >>= 1) { + index = (((st[0] & i) ? 4 : 0) | + ((st[1] & i) ? 2 : 0) | + ((st[2] & i) ? 1 : 0)); + pre = node; + node = (InternalNode *)getChild(node, getChildCount(node, index)); + } -/// If child index exists -int hasChild(InternalNode *node, int index) -{ - return (node->has_child >> index) & 1; -} + count = getChildCount(pre, index); + return pre; + } -/// Test if child is leaf -int isLeaf(InternalNode *node, int index) -{ - return (node->child_is_leaf >> index) & 1; -} + /************ Operators for internal nodes ************/ -/// Get the pointer to child index -Node *getChild(InternalNode *node, int count) -{ - return node->children[count]; -}; + /// If child index exists + int hasChild(InternalNode *node, int index) + { + return (node->has_child >> index) & 1; + } -/// Get total number of children -int getNumChildren(InternalNode *node) -{ - return numChildrenTable[node->has_child]; -} + /// Get the pointer to child index + Node *getChild(InternalNode *node, int count) + { + return node->children[count]; + }; -/// Get the count of children -int getChildCount(InternalNode *node, int index) -{ - return childrenCountTable[node->has_child][index]; -} -int getChildIndex(InternalNode *node, int count) -{ - return childrenIndexTable[node->has_child][count]; -} -int *getChildCounts(InternalNode *node) -{ - return childrenCountTable[node->has_child]; -} + /// Get total number of children + int getNumChildren(InternalNode *node) + { + return numChildrenTable[node->has_child]; + } -/// Get all children -void fillChildren(InternalNode *node, Node *children[8], int leaf[8]) -{ - int count = 0; - for (int i = 0; i < 8; i++) - { - leaf[i] = isLeaf(node, i); - if (hasChild(node, i)) - { - children[i] = getChild(node, count); - count++; - } - else { - children[i] = NULL; - leaf[i] = 0; - } + /// Get the count of children + int getChildCount(InternalNode *node, int index) + { + return childrenCountTable[node->has_child][index]; + } + int getChildIndex(InternalNode *node, int count) + { + return childrenIndexTable[node->has_child][count]; + } + int *getChildCounts(InternalNode *node) + { + return childrenCountTable[node->has_child]; } -} -/// Sets the child pointer -void setChild(InternalNode *node, int count, Node *chd) -{ - node->children[count] = chd; -} -void setInternalChild(InternalNode *node, int index, int count, InternalNode *chd) -{ - setChild(node, count, (Node *)chd); - node->has_child |= (1 << index); -} -void setLeafChild(InternalNode *node, int index, int count, LeafNode *chd) -{ - setChild(node, count, (Node *)chd); - node->has_child |= (1 << index); - node->child_is_leaf |= (1 << index); -} - -/// Add a kid to an existing internal node -/// Fix me: can we do this without wasting memory ? -/// Fixed: using variable memory -InternalNode *addChild(InternalNode *node, int index, Node *child, int aLeaf) -{ - // Create new internal node - int num = getNumChildren(node); - InternalNode *rnode = createInternal(num + 1); - - // Establish children - int i; - int count1 = 0, count2 = 0; - for (i = 0; i < 8; i++) - { - if (i == index) - { - if (aLeaf) - { - setLeafChild(rnode, i, count2, &child->leaf); - } - else { - setInternalChild(rnode, i, count2, &child->internal); - } - count2++; - } - else if (hasChild(node, i)) - { - if (isLeaf(node, i)) - { - setLeafChild(rnode, i, count2, &getChild(node, count1)->leaf); + /// Get all children + void fillChildren(InternalNode *node, Node *children[8], int leaf[8]) + { + int count = 0; + for (int i = 0; i < 8; i++) { + leaf[i] = node->is_child_leaf(i); + if (hasChild(node, i)) { + children[i] = getChild(node, count); + count++; } else { - setInternalChild(rnode, i, count2, &getChild(node, count1)->internal); + children[i] = NULL; + leaf[i] = 0; } - count1++; - count2++; } } - removeInternal(num, node); - return rnode; -} + /// Sets the child pointer + void setChild(InternalNode *node, int count, Node *chd) + { + node->children[count] = chd; + } + void setInternalChild(InternalNode *node, int index, int count, InternalNode *chd) + { + setChild(node, count, (Node *)chd); + node->has_child |= (1 << index); + } + void setLeafChild(InternalNode *node, int index, int count, LeafNode *chd) + { + setChild(node, count, (Node *)chd); + node->has_child |= (1 << index); + node->child_is_leaf |= (1 << index); + } -/// Allocate a node -InternalNode *createInternal(int length) -{ - InternalNode *inode = (InternalNode *)alloc[length]->allocate(); - inode->has_child = 0; - inode->child_is_leaf = 0; - return inode; -} + /// Add a kid to an existing internal node + /// Fix me: can we do this without wasting memory ? + /// Fixed: using variable memory + InternalNode *addChild(InternalNode *node, int index, Node *child, int aLeaf) + { + // Create new internal node + int num = getNumChildren(node); + InternalNode *rnode = createInternal(num + 1); -LeafNode *createLeaf(int length) -{ - assert(length <= 3); + // Establish children + int i; + int count1 = 0, count2 = 0; + for (i = 0; i < 8; i++) { + if (i == index) { + if (aLeaf) { + setLeafChild(rnode, i, count2, &child->leaf); + } + else { + setInternalChild(rnode, i, count2, &child->internal); + } + count2++; + } + else if (hasChild(node, i)) { + if (node->is_child_leaf(i)) { + setLeafChild(rnode, i, count2, &getChild(node, count1)->leaf); + } + else { + setInternalChild(rnode, i, count2, &getChild(node, count1)->internal); + } + count1++; + count2++; + } + } - LeafNode *lnode = (LeafNode *)leafalloc[length]->allocate(); - lnode->edge_parity = 0; - lnode->primary_edge_intersections = 0; - lnode->signs = 0; + removeInternal(num, node); + return rnode; + } - return lnode; -} + /// Allocate a node + InternalNode *createInternal(int length) + { + InternalNode *inode = (InternalNode *)alloc[length]->allocate(); + inode->has_child = 0; + inode->child_is_leaf = 0; + return inode; + } -void removeInternal(int num, InternalNode *node) -{ - alloc[num]->deallocate(node); -} + LeafNode *createLeaf(int length) + { + assert(length <= 3); -void removeLeaf(int num, LeafNode *leaf) -{ - assert(num >= 0 && num <= 3); - leafalloc[num]->deallocate(leaf); -} + LeafNode *lnode = (LeafNode *)leafalloc[length]->allocate(); + lnode->edge_parity = 0; + lnode->primary_edge_intersections = 0; + lnode->signs = 0; -/// Add a leaf (by creating a new par node with the leaf added) -InternalNode *addLeafChild(InternalNode *par, int index, int count, - LeafNode *leaf) -{ - int num = getNumChildren(par) + 1; - InternalNode *npar = createInternal(num); - *npar = *par; + return lnode; + } - if (num == 1) + void removeInternal(int num, InternalNode *node) { - setLeafChild(npar, index, 0, leaf); + alloc[num]->deallocate(node); } - else { - int i; - for (i = 0; i < count; i++) - { - setChild(npar, i, getChild(par, i)); - } - setLeafChild(npar, index, count, leaf); - for (i = count + 1; i < num; i++) - { - setChild(npar, i, getChild(par, i - 1)); - } + + void removeLeaf(int num, LeafNode *leaf) + { + assert(num >= 0 && num <= 3); + leafalloc[num]->deallocate(leaf); } - removeInternal(num - 1, par); - return npar; -} + /// Add a leaf (by creating a new par node with the leaf added) + InternalNode *addLeafChild(InternalNode *par, int index, int count, + LeafNode *leaf) + { + int num = getNumChildren(par) + 1; + InternalNode *npar = createInternal(num); + *npar = *par; -InternalNode *addInternalChild(InternalNode *par, int index, int count, - InternalNode *node) -{ - int num = getNumChildren(par) + 1; - InternalNode *npar = createInternal(num); - *npar = *par; + if (num == 1) { + setLeafChild(npar, index, 0, leaf); + } + else { + int i; + for (i = 0; i < count; i++) { + setChild(npar, i, getChild(par, i)); + } + setLeafChild(npar, index, count, leaf); + for (i = count + 1; i < num; i++) { + setChild(npar, i, getChild(par, i - 1)); + } + } - if (num == 1) - { - setInternalChild(npar, index, 0, node); + removeInternal(num - 1, par); + return npar; } - else { - int i; - for (i = 0; i < count; i++) - { - setChild(npar, i, getChild(par, i)); + + InternalNode *addInternalChild(InternalNode *par, int index, int count, + InternalNode *node) + { + int num = getNumChildren(par) + 1; + InternalNode *npar = createInternal(num); + *npar = *par; + + if (num == 1) { + setInternalChild(npar, index, 0, node); } - setInternalChild(npar, index, count, node); - for (i = count + 1; i < num; i++) - { - setChild(npar, i, getChild(par, i - 1)); + else { + int i; + for (i = 0; i < count; i++) { + setChild(npar, i, getChild(par, i)); + } + setInternalChild(npar, index, count, node); + for (i = count + 1; i < num; i++) { + setChild(npar, i, getChild(par, i - 1)); + } } - } - removeInternal(num - 1, par); - return npar; -} + removeInternal(num - 1, par); + return npar; + } }; #endif |