forked from bartvdbraak/blender
2858 lines
64 KiB
C++
2858 lines
64 KiB
C++
/*
|
|
* ***** BEGIN GPL LICENSE BLOCK *****
|
|
*
|
|
* This program is free software; you can redistribute it and/or
|
|
* modify it under the terms of the GNU General Public License
|
|
* as published by the Free Software Foundation; either version 2
|
|
* of the License, or (at your option) any later version.
|
|
*
|
|
* This program is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
* GNU General Public License for more details.
|
|
*
|
|
* You should have received a copy of the GNU General Public License
|
|
* along with this program; if not, write to the Free Software Foundation,
|
|
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
*
|
|
* Contributor(s): Tao Ju, Nicholas Bishop
|
|
*
|
|
* ***** END GPL LICENSE BLOCK *****
|
|
*/
|
|
|
|
#ifdef WITH_CXX_GUARDEDALLOC
|
|
# include "MEM_guardedalloc.h"
|
|
#endif
|
|
|
|
#include "octree.h"
|
|
#include <Eigen/Dense>
|
|
#include <limits>
|
|
#include <time.h>
|
|
|
|
/**
|
|
* Implementations of Octree member functions.
|
|
*
|
|
* @author Tao Ju
|
|
*/
|
|
|
|
/* set to non-zero value to enable debugging output */
|
|
#define DC_DEBUG 0
|
|
|
|
#if DC_DEBUG
|
|
/* enable debug printfs */
|
|
#define dc_printf printf
|
|
#else
|
|
/* disable debug printfs */
|
|
#define dc_printf(...) do {} while (0)
|
|
#endif
|
|
|
|
Octree::Octree(ModelReader *mr,
|
|
DualConAllocOutput alloc_output_func,
|
|
DualConAddVert add_vert_func,
|
|
DualConAddQuad add_quad_func,
|
|
DualConFlags flags, DualConMode dualcon_mode, int depth,
|
|
float threshold, float sharpness)
|
|
: use_flood_fill(flags & DUALCON_FLOOD_FILL),
|
|
/* note on `use_manifold':
|
|
|
|
After playing around with this option, the only case I could
|
|
find where this option gives different results is on
|
|
relatively thin corners. Sometimes along these corners two
|
|
vertices from separate sides will be placed in the same
|
|
position, so hole gets filled with a 5-sided face, where two
|
|
of those vertices are in the same 3D location. If
|
|
`use_manifold' is disabled, then the modifier doesn't
|
|
generate two separate vertices so the results end up as all
|
|
quads.
|
|
|
|
Since the results are just as good with all quads, there
|
|
doesn't seem any reason to allow this to be toggled in
|
|
Blender. -nicholasbishop
|
|
*/
|
|
use_manifold(false),
|
|
hermite_num(sharpness),
|
|
mode(dualcon_mode),
|
|
alloc_output(alloc_output_func),
|
|
add_vert(add_vert_func),
|
|
add_quad(add_quad_func)
|
|
{
|
|
thresh = threshold;
|
|
reader = mr;
|
|
dimen = 1 << GRID_DIMENSION;
|
|
range = reader->getBoundingBox(origin);
|
|
nodeCount = nodeSpace = 0;
|
|
maxDepth = depth;
|
|
mindimen = (dimen >> maxDepth);
|
|
minshift = (GRID_DIMENSION - maxDepth);
|
|
buildTable();
|
|
|
|
maxTrianglePerCell = 0;
|
|
|
|
// Initialize memory
|
|
#ifdef IN_VERBOSE_MODE
|
|
dc_printf("Range: %f origin: %f, %f,%f \n", range, origin[0], origin[1], origin[2]);
|
|
dc_printf("Initialize memory...\n");
|
|
#endif
|
|
initMemory();
|
|
root = (Node *)createInternal(0);
|
|
|
|
// Read MC table
|
|
#ifdef IN_VERBOSE_MODE
|
|
dc_printf("Reading contour table...\n");
|
|
#endif
|
|
cubes = new Cubes();
|
|
|
|
}
|
|
|
|
Octree::~Octree()
|
|
{
|
|
delete cubes;
|
|
freeMemory();
|
|
}
|
|
|
|
void Octree::scanConvert()
|
|
{
|
|
// Scan triangles
|
|
#if DC_DEBUG
|
|
clock_t start, finish;
|
|
start = clock();
|
|
#endif
|
|
|
|
addAllTriangles();
|
|
resetMinimalEdges();
|
|
preparePrimalEdgesMask(&root->internal);
|
|
|
|
#if DC_DEBUG
|
|
finish = clock();
|
|
dc_printf("Time taken: %f seconds \n",
|
|
(double)(finish - start) / CLOCKS_PER_SEC);
|
|
#endif
|
|
|
|
// Generate signs
|
|
// Find holes
|
|
#if DC_DEBUG
|
|
dc_printf("Patching...\n");
|
|
start = clock();
|
|
#endif
|
|
trace();
|
|
#if DC_DEBUG
|
|
finish = clock();
|
|
dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC);
|
|
#ifdef IN_VERBOSE_MODE
|
|
dc_printf("Holes: %d Average Length: %f Max Length: %d \n", numRings, (float)totRingLengths / (float) numRings, maxRingLength);
|
|
#endif
|
|
#endif
|
|
|
|
// Check again
|
|
int tnumRings = numRings;
|
|
trace();
|
|
#ifdef IN_VERBOSE_MODE
|
|
dc_printf("Holes after patching: %d \n", numRings);
|
|
#endif
|
|
numRings = tnumRings;
|
|
|
|
#if DC_DEBUG
|
|
dc_printf("Building signs...\n");
|
|
start = clock();
|
|
#endif
|
|
buildSigns();
|
|
#if DC_DEBUG
|
|
finish = clock();
|
|
dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC);
|
|
#endif
|
|
|
|
if (use_flood_fill) {
|
|
/*
|
|
start = clock();
|
|
floodFill();
|
|
// Check again
|
|
tnumRings = numRings;
|
|
trace();
|
|
dc_printf("Holes after filling: %d \n", numRings);
|
|
numRings = tnumRings;
|
|
buildSigns();
|
|
finish = clock();
|
|
dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC);
|
|
*/
|
|
#if DC_DEBUG
|
|
start = clock();
|
|
dc_printf("Removing components...\n");
|
|
#endif
|
|
floodFill();
|
|
buildSigns();
|
|
// dc_printf("Checking...\n");
|
|
// floodFill();
|
|
#if DC_DEBUG
|
|
finish = clock();
|
|
dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC);
|
|
#endif
|
|
}
|
|
|
|
// Output
|
|
#if DC_DEBUG
|
|
start = clock();
|
|
#endif
|
|
writeOut();
|
|
#if DC_DEBUG
|
|
finish = clock();
|
|
#endif
|
|
// dc_printf("Time taken: %f seconds \n", (double)(finish - start) / CLOCKS_PER_SEC);
|
|
|
|
// Print info
|
|
#ifdef IN_VERBOSE_MODE
|
|
printMemUsage();
|
|
#endif
|
|
}
|
|
|
|
void Octree::initMemory()
|
|
{
|
|
leafalloc[0] = new MemoryAllocator<sizeof(LeafNode)>();
|
|
leafalloc[1] = new MemoryAllocator<sizeof(LeafNode) + sizeof(float) *EDGE_FLOATS>();
|
|
leafalloc[2] = new MemoryAllocator<sizeof(LeafNode) + sizeof(float) *EDGE_FLOATS * 2>();
|
|
leafalloc[3] = new MemoryAllocator<sizeof(LeafNode) + sizeof(float) *EDGE_FLOATS * 3>();
|
|
|
|
alloc[0] = new MemoryAllocator<sizeof(InternalNode)>();
|
|
alloc[1] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *)>();
|
|
alloc[2] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *) * 2>();
|
|
alloc[3] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *) * 3>();
|
|
alloc[4] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *) * 4>();
|
|
alloc[5] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *) * 5>();
|
|
alloc[6] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *) * 6>();
|
|
alloc[7] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *) * 7>();
|
|
alloc[8] = new MemoryAllocator<sizeof(InternalNode) + sizeof(Node *) * 8>();
|
|
}
|
|
|
|
void Octree::freeMemory()
|
|
{
|
|
for (int i = 0; i < 9; i++) {
|
|
alloc[i]->destroy();
|
|
delete alloc[i];
|
|
}
|
|
|
|
for (int i = 0; i < 4; i++) {
|
|
leafalloc[i]->destroy();
|
|
delete leafalloc[i];
|
|
}
|
|
}
|
|
|
|
void Octree::printMemUsage()
|
|
{
|
|
int totalbytes = 0;
|
|
dc_printf("********* Internal nodes: \n");
|
|
for (int i = 0; i < 9; i++) {
|
|
alloc[i]->printInfo();
|
|
|
|
totalbytes += alloc[i]->getAll() * alloc[i]->getBytes();
|
|
}
|
|
dc_printf("********* Leaf nodes: \n");
|
|
int totalLeafs = 0;
|
|
for (int i = 0; i < 4; i++) {
|
|
leafalloc[i]->printInfo();
|
|
|
|
totalbytes += leafalloc[i]->getAll() * leafalloc[i]->getBytes();
|
|
totalLeafs += leafalloc[i]->getAllocated();
|
|
}
|
|
|
|
dc_printf("Total allocated bytes on disk: %d \n", totalbytes);
|
|
dc_printf("Total leaf nodes: %d\n", totalLeafs);
|
|
}
|
|
|
|
void Octree::resetMinimalEdges()
|
|
{
|
|
cellProcParity(root, 0, maxDepth);
|
|
}
|
|
|
|
void Octree::addAllTriangles()
|
|
{
|
|
Triangle *trian;
|
|
int count = 0;
|
|
|
|
#if DC_DEBUG
|
|
int total = reader->getNumTriangles();
|
|
int unitcount = 1000;
|
|
dc_printf("\nScan converting to depth %d...\n", maxDepth);
|
|
#endif
|
|
|
|
srand(0);
|
|
|
|
while ((trian = reader->getNextTriangle()) != NULL) {
|
|
// Drop triangles
|
|
{
|
|
addTriangle(trian, count);
|
|
}
|
|
delete trian;
|
|
|
|
count++;
|
|
|
|
#if DC_DEBUG
|
|
if (count % unitcount == 0) {
|
|
putchar(13);
|
|
|
|
switch ((count / unitcount) % 4) {
|
|
case 0: dc_printf("-");
|
|
break;
|
|
case 1: dc_printf("/");
|
|
break;
|
|
case 2: dc_printf("|");
|
|
break;
|
|
case 3: dc_printf("\\");
|
|
break;
|
|
}
|
|
|
|
float percent = (float) count / total;
|
|
|
|
/*
|
|
int totbars = 50;
|
|
int bars =(int)(percent * totbars);
|
|
for(int i = 0; i < bars; i ++) {
|
|
putchar(219);
|
|
}
|
|
for(i = bars; i < totbars; i ++) {
|
|
putchar(176);
|
|
}
|
|
*/
|
|
|
|
dc_printf(" %d triangles: ", count);
|
|
dc_printf(" %f%% complete.", 100 * percent);
|
|
}
|
|
#endif
|
|
|
|
}
|
|
putchar(13);
|
|
}
|
|
|
|
/* 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;
|
|
|
|
/* 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;
|
|
}
|
|
|
|
/* 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 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;
|
|
}
|
|
|
|
#if 0
|
|
static void print_depth(int height, int maxDepth)
|
|
{
|
|
for (int i = 0; i < maxDepth - height; i++)
|
|
printf(" ");
|
|
}
|
|
#endif
|
|
|
|
InternalNode *Octree::addTriangle(InternalNode *node, CubeTriangleIsect *p, int height)
|
|
{
|
|
int i;
|
|
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];
|
|
tempdiff[2] += vertdiff[i][2];
|
|
|
|
/* Quick pruning using bounding box */
|
|
if (boxmask & (1 << i)) {
|
|
subp->shift(tempdiff);
|
|
tempdiff[0] = tempdiff[1] = tempdiff[2] = 0;
|
|
|
|
/* Pruning using intersection test */
|
|
if (subp->isIntersecting()) {
|
|
if (!node->has_child(i)) {
|
|
if (height == 1)
|
|
node = addLeafChild(node, i, count, createLeaf(0));
|
|
else
|
|
node = addInternalChild(node, i, count, createInternal(0));
|
|
}
|
|
Node *chd = node->get_child(count);
|
|
|
|
if (node->is_child_leaf(i))
|
|
node->set_child(count, (Node *)updateCell(&chd->leaf, subp));
|
|
else
|
|
node->set_child(count, (Node *)addTriangle(&chd->internal, subp, height - 1));
|
|
}
|
|
}
|
|
|
|
if (node->has_child(i))
|
|
count++;
|
|
}
|
|
|
|
delete subp;
|
|
|
|
return node;
|
|
}
|
|
|
|
LeafNode *Octree::updateCell(LeafNode *node, CubeTriangleIsect *p)
|
|
{
|
|
int i;
|
|
|
|
// Edge connectivity
|
|
int mask[3] = {0, 4, 8 };
|
|
int oldc = 0, newc = 0;
|
|
float offs[3];
|
|
float a[3], b[3], c[3];
|
|
|
|
for (i = 0; i < 3; i++) {
|
|
if (!getEdgeParity(node, mask[i])) {
|
|
if (p->isIntersectingPrimary(i)) {
|
|
// actualQuads ++;
|
|
setEdge(node, mask[i]);
|
|
offs[newc] = p->getIntersectionPrimary(i);
|
|
a[newc] = (float) p->inherit->norm[0];
|
|
b[newc] = (float) p->inherit->norm[1];
|
|
c[newc] = (float) p->inherit->norm[2];
|
|
newc++;
|
|
}
|
|
}
|
|
else {
|
|
offs[newc] = getEdgeOffsetNormal(node, oldc, a[newc], b[newc], c[newc]);
|
|
|
|
oldc++;
|
|
newc++;
|
|
}
|
|
}
|
|
|
|
if (newc > oldc) {
|
|
// New offsets added, update this node
|
|
node = updateEdgeOffsetsNormals(node, oldc, newc, offs, a, b, c);
|
|
}
|
|
|
|
return node;
|
|
}
|
|
|
|
void Octree::preparePrimalEdgesMask(InternalNode *node)
|
|
{
|
|
int count = 0;
|
|
for (int i = 0; i < 8; i++) {
|
|
if (node->has_child(i)) {
|
|
if (node->is_child_leaf(i))
|
|
createPrimalEdgesMask(&node->get_child(count)->leaf);
|
|
else
|
|
preparePrimalEdgesMask(&node->get_child(count)->internal);
|
|
|
|
count++;
|
|
}
|
|
}
|
|
}
|
|
|
|
void Octree::trace()
|
|
{
|
|
int st[3] = {0, 0, 0, };
|
|
numRings = 0;
|
|
totRingLengths = 0;
|
|
maxRingLength = 0;
|
|
|
|
PathList *chdpath = NULL;
|
|
root = trace(root, st, dimen, maxDepth, chdpath);
|
|
|
|
if (chdpath != NULL) {
|
|
dc_printf("there are incomplete rings.\n");
|
|
printPaths(chdpath);
|
|
};
|
|
}
|
|
|
|
Node *Octree::trace(Node *newnode, int *st, int len, int depth, PathList *& paths)
|
|
{
|
|
len >>= 1;
|
|
PathList *chdpaths[8];
|
|
Node *chd[8];
|
|
int nst[8][3];
|
|
int i, j;
|
|
|
|
// Get children paths
|
|
int chdleaf[8];
|
|
newnode->internal.fill_children(chd, chdleaf);
|
|
|
|
// int count = 0;
|
|
for (i = 0; i < 8; i++) {
|
|
for (j = 0; j < 3; j++) {
|
|
nst[i][j] = st[j] + len * vertmap[i][j];
|
|
}
|
|
|
|
if (chd[i] == NULL || newnode->internal.is_child_leaf(i)) {
|
|
chdpaths[i] = NULL;
|
|
}
|
|
else {
|
|
trace(chd[i], nst[i], len, depth - 1, chdpaths[i]);
|
|
}
|
|
}
|
|
|
|
// Get connectors on the faces
|
|
PathList *conn[12];
|
|
Node *nf[2];
|
|
int lf[2];
|
|
int df[2] = {depth - 1, depth - 1};
|
|
int *nstf[2];
|
|
|
|
newnode->internal.fill_children(chd, chdleaf);
|
|
|
|
for (i = 0; i < 12; i++) {
|
|
int c[2] = {cellProcFaceMask[i][0], cellProcFaceMask[i][1]};
|
|
|
|
for (int j = 0; j < 2; j++) {
|
|
lf[j] = chdleaf[c[j]];
|
|
nf[j] = chd[c[j]];
|
|
nstf[j] = nst[c[j]];
|
|
}
|
|
|
|
conn[i] = NULL;
|
|
|
|
findPaths((Node **)nf, lf, df, nstf, depth - 1, cellProcFaceMask[i][2], conn[i]);
|
|
|
|
//if(conn[i]) {
|
|
// printPath(conn[i]);
|
|
//}
|
|
}
|
|
|
|
// Connect paths
|
|
PathList *rings = NULL;
|
|
combinePaths(chdpaths[0], chdpaths[1], conn[8], rings);
|
|
combinePaths(chdpaths[2], chdpaths[3], conn[9], rings);
|
|
combinePaths(chdpaths[4], chdpaths[5], conn[10], rings);
|
|
combinePaths(chdpaths[6], chdpaths[7], conn[11], rings);
|
|
|
|
combinePaths(chdpaths[0], chdpaths[2], conn[4], rings);
|
|
combinePaths(chdpaths[4], chdpaths[6], conn[5], rings);
|
|
combinePaths(chdpaths[0], NULL, conn[6], rings);
|
|
combinePaths(chdpaths[4], NULL, conn[7], rings);
|
|
|
|
combinePaths(chdpaths[0], chdpaths[4], conn[0], rings);
|
|
combinePaths(chdpaths[0], NULL, conn[1], rings);
|
|
combinePaths(chdpaths[0], NULL, conn[2], rings);
|
|
combinePaths(chdpaths[0], NULL, conn[3], rings);
|
|
|
|
// By now, only chdpaths[0] and rings have contents
|
|
|
|
// Process rings
|
|
if (rings) {
|
|
// printPath(rings);
|
|
|
|
/* Let's count first */
|
|
PathList *trings = rings;
|
|
while (trings) {
|
|
numRings++;
|
|
totRingLengths += trings->length;
|
|
if (trings->length > maxRingLength) {
|
|
maxRingLength = trings->length;
|
|
}
|
|
trings = trings->next;
|
|
}
|
|
|
|
// printPath(rings);
|
|
newnode = patch(newnode, st, (len << 1), rings);
|
|
}
|
|
|
|
// Return incomplete paths
|
|
paths = chdpaths[0];
|
|
return newnode;
|
|
}
|
|
|
|
void Octree::findPaths(Node *node[2], int leaf[2], int depth[2], int *st[2], int maxdep, int dir, PathList *& paths)
|
|
{
|
|
if (!(node[0] && node[1])) {
|
|
return;
|
|
}
|
|
|
|
if (!(leaf[0] && leaf[1])) {
|
|
// Not at the bottom, recur
|
|
|
|
// Fill children nodes
|
|
int i, j;
|
|
Node *chd[2][8];
|
|
int chdleaf[2][8];
|
|
int nst[2][8][3];
|
|
|
|
for (j = 0; j < 2; j++) {
|
|
if (!leaf[j]) {
|
|
node[j]->internal.fill_children(chd[j], chdleaf[j]);
|
|
|
|
int len = (dimen >> (maxDepth - depth[j] + 1));
|
|
for (i = 0; i < 8; i++) {
|
|
for (int k = 0; k < 3; k++) {
|
|
nst[j][i][k] = st[j][k] + len * vertmap[i][k];
|
|
}
|
|
}
|
|
|
|
}
|
|
}
|
|
|
|
// 4 face calls
|
|
Node *nf[2];
|
|
int df[2];
|
|
int lf[2];
|
|
int *nstf[2];
|
|
for (i = 0; i < 4; i++) {
|
|
int c[2] = {faceProcFaceMask[dir][i][0], faceProcFaceMask[dir][i][1]};
|
|
for (int j = 0; j < 2; j++) {
|
|
if (leaf[j]) {
|
|
lf[j] = leaf[j];
|
|
nf[j] = node[j];
|
|
df[j] = depth[j];
|
|
nstf[j] = st[j];
|
|
}
|
|
else {
|
|
lf[j] = chdleaf[j][c[j]];
|
|
nf[j] = chd[j][c[j]];
|
|
df[j] = depth[j] - 1;
|
|
nstf[j] = nst[j][c[j]];
|
|
}
|
|
}
|
|
findPaths(nf, lf, df, nstf, maxdep - 1, faceProcFaceMask[dir][i][2], paths);
|
|
}
|
|
|
|
}
|
|
else {
|
|
// At the bottom, check this face
|
|
int ind = (depth[0] == maxdep ? 0 : 1);
|
|
int fcind = 2 * dir + (1 - ind);
|
|
if (getFaceParity((LeafNode *)node[ind], fcind)) {
|
|
// Add into path
|
|
PathElement *ele1 = new PathElement;
|
|
PathElement *ele2 = new PathElement;
|
|
|
|
ele1->pos[0] = st[0][0];
|
|
ele1->pos[1] = st[0][1];
|
|
ele1->pos[2] = st[0][2];
|
|
|
|
ele2->pos[0] = st[1][0];
|
|
ele2->pos[1] = st[1][1];
|
|
ele2->pos[2] = st[1][2];
|
|
|
|
ele1->next = ele2;
|
|
ele2->next = NULL;
|
|
|
|
PathList *lst = new PathList;
|
|
lst->head = ele1;
|
|
lst->tail = ele2;
|
|
lst->length = 2;
|
|
lst->next = paths;
|
|
paths = lst;
|
|
|
|
// int l =(dimen >> maxDepth);
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
void Octree::combinePaths(PathList *& list1, PathList *list2, PathList *paths, PathList *& rings)
|
|
{
|
|
// Make new list of paths
|
|
PathList *nlist = NULL;
|
|
|
|
// Search for each connectors in paths
|
|
PathList *tpaths = paths;
|
|
PathList *tlist, *pre;
|
|
while (tpaths) {
|
|
PathList *singlist = tpaths;
|
|
PathList *templist;
|
|
tpaths = tpaths->next;
|
|
singlist->next = NULL;
|
|
|
|
// Look for hookup in list1
|
|
tlist = list1;
|
|
pre = NULL;
|
|
while (tlist) {
|
|
if ((templist = combineSinglePath(list1, pre, tlist, singlist, NULL, singlist)) != NULL) {
|
|
singlist = templist;
|
|
continue;
|
|
}
|
|
pre = tlist;
|
|
tlist = tlist->next;
|
|
}
|
|
|
|
// Look for hookup in list2
|
|
tlist = list2;
|
|
pre = NULL;
|
|
while (tlist) {
|
|
if ((templist = combineSinglePath(list2, pre, tlist, singlist, NULL, singlist)) != NULL) {
|
|
singlist = templist;
|
|
continue;
|
|
}
|
|
pre = tlist;
|
|
tlist = tlist->next;
|
|
}
|
|
|
|
// Look for hookup in nlist
|
|
tlist = nlist;
|
|
pre = NULL;
|
|
while (tlist) {
|
|
if ((templist = combineSinglePath(nlist, pre, tlist, singlist, NULL, singlist)) != NULL) {
|
|
singlist = templist;
|
|
continue;
|
|
}
|
|
pre = tlist;
|
|
tlist = tlist->next;
|
|
}
|
|
|
|
// Add to nlist or rings
|
|
if (isEqual(singlist->head, singlist->tail)) {
|
|
PathElement *temp = singlist->head;
|
|
singlist->head = temp->next;
|
|
delete temp;
|
|
singlist->length--;
|
|
singlist->tail->next = singlist->head;
|
|
|
|
singlist->next = rings;
|
|
rings = singlist;
|
|
}
|
|
else {
|
|
singlist->next = nlist;
|
|
nlist = singlist;
|
|
}
|
|
|
|
}
|
|
|
|
// Append list2 and nlist to the end of list1
|
|
tlist = list1;
|
|
if (tlist != NULL) {
|
|
while (tlist->next != NULL) {
|
|
tlist = tlist->next;
|
|
}
|
|
tlist->next = list2;
|
|
}
|
|
else {
|
|
tlist = list2;
|
|
list1 = list2;
|
|
}
|
|
|
|
if (tlist != NULL) {
|
|
while (tlist->next != NULL) {
|
|
tlist = tlist->next;
|
|
}
|
|
tlist->next = nlist;
|
|
}
|
|
else {
|
|
tlist = nlist;
|
|
list1 = nlist;
|
|
}
|
|
|
|
}
|
|
|
|
PathList *Octree::combineSinglePath(PathList *& head1, PathList *pre1, PathList *& list1, PathList *& head2, PathList *pre2, PathList *& list2)
|
|
{
|
|
if (isEqual(list1->head, list2->head) || isEqual(list1->tail, list2->tail)) {
|
|
// Reverse the list
|
|
if (list1->length < list2->length) {
|
|
// Reverse list1
|
|
PathElement *prev = list1->head;
|
|
PathElement *next = prev->next;
|
|
prev->next = NULL;
|
|
while (next != NULL) {
|
|
PathElement *tnext = next->next;
|
|
next->next = prev;
|
|
|
|
prev = next;
|
|
next = tnext;
|
|
}
|
|
|
|
list1->tail = list1->head;
|
|
list1->head = prev;
|
|
}
|
|
else {
|
|
// Reverse list2
|
|
PathElement *prev = list2->head;
|
|
PathElement *next = prev->next;
|
|
prev->next = NULL;
|
|
while (next != NULL) {
|
|
PathElement *tnext = next->next;
|
|
next->next = prev;
|
|
|
|
prev = next;
|
|
next = tnext;
|
|
}
|
|
|
|
list2->tail = list2->head;
|
|
list2->head = prev;
|
|
}
|
|
}
|
|
|
|
if (isEqual(list1->head, list2->tail)) {
|
|
|
|
// Easy case
|
|
PathElement *temp = list1->head->next;
|
|
delete list1->head;
|
|
list2->tail->next = temp;
|
|
|
|
PathList *nlist = new PathList;
|
|
nlist->length = list1->length + list2->length - 1;
|
|
nlist->head = list2->head;
|
|
nlist->tail = list1->tail;
|
|
nlist->next = NULL;
|
|
|
|
deletePath(head1, pre1, list1);
|
|
deletePath(head2, pre2, list2);
|
|
|
|
return nlist;
|
|
}
|
|
else if (isEqual(list1->tail, list2->head)) {
|
|
// Easy case
|
|
PathElement *temp = list2->head->next;
|
|
delete list2->head;
|
|
list1->tail->next = temp;
|
|
|
|
PathList *nlist = new PathList;
|
|
nlist->length = list1->length + list2->length - 1;
|
|
nlist->head = list1->head;
|
|
nlist->tail = list2->tail;
|
|
nlist->next = NULL;
|
|
|
|
deletePath(head1, pre1, list1);
|
|
deletePath(head2, pre2, list2);
|
|
|
|
return nlist;
|
|
}
|
|
|
|
return NULL;
|
|
}
|
|
|
|
void Octree::deletePath(PathList *& head, PathList *pre, PathList *& curr)
|
|
{
|
|
PathList *temp = curr;
|
|
curr = temp->next;
|
|
delete temp;
|
|
|
|
if (pre == NULL) {
|
|
head = curr;
|
|
}
|
|
else {
|
|
pre->next = curr;
|
|
}
|
|
}
|
|
|
|
void Octree::printElement(PathElement *ele)
|
|
{
|
|
if (ele != NULL) {
|
|
dc_printf("(%d %d %d)", ele->pos[0], ele->pos[1], ele->pos[2]);
|
|
}
|
|
}
|
|
|
|
void Octree::printPath(PathList *path)
|
|
{
|
|
PathElement *n = path->head;
|
|
int same = 0;
|
|
|
|
#if DC_DEBUG
|
|
int len = (dimen >> maxDepth);
|
|
#endif
|
|
while (n && (same == 0 || n != path->head)) {
|
|
same++;
|
|
dc_printf("(%d %d %d)", n->pos[0] / len, n->pos[1] / len, n->pos[2] / len);
|
|
n = n->next;
|
|
}
|
|
|
|
if (n == path->head) {
|
|
dc_printf(" Ring!\n");
|
|
}
|
|
else {
|
|
dc_printf(" %p end!\n", n);
|
|
}
|
|
}
|
|
|
|
void Octree::printPath(PathElement *path)
|
|
{
|
|
PathElement *n = path;
|
|
int same = 0;
|
|
#if DC_DEBUG
|
|
int len = (dimen >> maxDepth);
|
|
#endif
|
|
while (n && (same == 0 || n != path)) {
|
|
same++;
|
|
dc_printf("(%d %d %d)", n->pos[0] / len, n->pos[1] / len, n->pos[2] / len);
|
|
n = n->next;
|
|
}
|
|
|
|
if (n == path) {
|
|
dc_printf(" Ring!\n");
|
|
}
|
|
else {
|
|
dc_printf(" %p end!\n", n);
|
|
}
|
|
|
|
}
|
|
|
|
|
|
void Octree::printPaths(PathList *path)
|
|
{
|
|
PathList *iter = path;
|
|
int i = 0;
|
|
while (iter != NULL) {
|
|
dc_printf("Path %d:\n", i);
|
|
printPath(iter);
|
|
iter = iter->next;
|
|
i++;
|
|
}
|
|
}
|
|
|
|
Node *Octree::patch(Node *newnode, int st[3], int len, PathList *rings)
|
|
{
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("Call to PATCH with rings: \n");
|
|
printPaths(rings);
|
|
#endif
|
|
|
|
/* Do nothing but couting
|
|
PathList* tlist = rings;
|
|
PathList* ttlist;
|
|
PathElement* telem, * ttelem;
|
|
while(tlist!= NULL) {
|
|
// printPath(tlist);
|
|
numRings ++;
|
|
totRingLengths += tlist->length;
|
|
if(tlist->length > maxRingLength) {
|
|
maxRingLength = tlist->length;
|
|
}
|
|
ttlist = tlist;
|
|
tlist = tlist->next;
|
|
}
|
|
return node;
|
|
*/
|
|
|
|
|
|
/* Pass onto separate calls in each direction */
|
|
if (len == mindimen) {
|
|
dc_printf("Error! should have no list by now.\n");
|
|
exit(0);
|
|
}
|
|
|
|
// YZ plane
|
|
PathList *xlists[2];
|
|
newnode = patchSplit(newnode, st, len, rings, 0, xlists[0], xlists[1]);
|
|
|
|
// XZ plane
|
|
PathList *ylists[4];
|
|
newnode = patchSplit(newnode, st, len, xlists[0], 1, ylists[0], ylists[1]);
|
|
newnode = patchSplit(newnode, st, len, xlists[1], 1, ylists[2], ylists[3]);
|
|
|
|
// XY plane
|
|
PathList *zlists[8];
|
|
newnode = patchSplit(newnode, st, len, ylists[0], 2, zlists[0], zlists[1]);
|
|
newnode = patchSplit(newnode, st, len, ylists[1], 2, zlists[2], zlists[3]);
|
|
newnode = patchSplit(newnode, st, len, ylists[2], 2, zlists[4], zlists[5]);
|
|
newnode = patchSplit(newnode, st, len, ylists[3], 2, zlists[6], zlists[7]);
|
|
|
|
// Recur
|
|
len >>= 1;
|
|
int count = 0;
|
|
for (int i = 0; i < 8; i++) {
|
|
if (zlists[i] != NULL) {
|
|
int nori[3] = {
|
|
st[0] + len * vertmap[i][0],
|
|
st[1] + len * vertmap[i][1],
|
|
st[2] + len * vertmap[i][2]
|
|
};
|
|
patch(newnode->internal.get_child(count), nori, len, zlists[i]);
|
|
}
|
|
|
|
if (newnode->internal.has_child(i)) {
|
|
count++;
|
|
}
|
|
}
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("Return from PATCH\n");
|
|
#endif
|
|
return newnode;
|
|
|
|
}
|
|
|
|
|
|
Node *Octree::patchSplit(Node *newnode, int st[3], int len, PathList *rings,
|
|
int dir, PathList *& nrings1, PathList *& nrings2)
|
|
{
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("Call to PATCHSPLIT with direction %d and rings: \n", dir);
|
|
printPaths(rings);
|
|
#endif
|
|
|
|
nrings1 = NULL;
|
|
nrings2 = NULL;
|
|
PathList *tmp;
|
|
while (rings != NULL) {
|
|
// Process this ring
|
|
newnode = patchSplitSingle(newnode, st, len, rings->head, dir, nrings1, nrings2);
|
|
|
|
// Delete this ring from the group
|
|
tmp = rings;
|
|
rings = rings->next;
|
|
delete tmp;
|
|
}
|
|
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("Return from PATCHSPLIT with \n");
|
|
dc_printf("Rings gourp 1:\n");
|
|
printPaths(nrings1);
|
|
dc_printf("Rings group 2:\n");
|
|
printPaths(nrings2);
|
|
#endif
|
|
|
|
return newnode;
|
|
}
|
|
|
|
Node *Octree::patchSplitSingle(Node *newnode, int st[3], int len, PathElement *head, int dir, PathList *& nrings1, PathList *& nrings2)
|
|
{
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("Call to PATCHSPLITSINGLE with direction %d and path: \n", dir);
|
|
printPath(head);
|
|
#endif
|
|
|
|
if (head == NULL) {
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("Return from PATCHSPLITSINGLE with head==NULL.\n");
|
|
#endif
|
|
return newnode;
|
|
}
|
|
else {
|
|
// printPath(head);
|
|
}
|
|
|
|
// Walk along the ring to find pair of intersections
|
|
PathElement *pre1 = NULL;
|
|
PathElement *pre2 = NULL;
|
|
int side = findPair(head, st[dir] + len / 2, dir, pre1, pre2);
|
|
|
|
/*
|
|
if(pre1 == pre2) {
|
|
int edgelen =(dimen >> maxDepth);
|
|
dc_printf("Location: %d %d %d Direction: %d Reso: %d\n", st[0]/edgelen, st[1]/edgelen, st[2]/edgelen, dir, len/edgelen);
|
|
printPath(head);
|
|
exit(0);
|
|
}
|
|
*/
|
|
|
|
if (side) {
|
|
// Entirely on one side
|
|
PathList *nring = new PathList();
|
|
nring->head = head;
|
|
|
|
if (side == -1) {
|
|
nring->next = nrings1;
|
|
nrings1 = nring;
|
|
}
|
|
else {
|
|
nring->next = nrings2;
|
|
nrings2 = nring;
|
|
}
|
|
}
|
|
else {
|
|
// Break into two parts
|
|
PathElement *nxt1 = pre1->next;
|
|
PathElement *nxt2 = pre2->next;
|
|
pre1->next = nxt2;
|
|
pre2->next = nxt1;
|
|
|
|
newnode = connectFace(newnode, st, len, dir, pre1, pre2);
|
|
|
|
if (isEqual(pre1, pre1->next)) {
|
|
if (pre1 == pre1->next) {
|
|
delete pre1;
|
|
pre1 = NULL;
|
|
}
|
|
else {
|
|
PathElement *temp = pre1->next;
|
|
pre1->next = temp->next;
|
|
delete temp;
|
|
}
|
|
}
|
|
if (isEqual(pre2, pre2->next)) {
|
|
if (pre2 == pre2->next) {
|
|
delete pre2;
|
|
pre2 = NULL;
|
|
}
|
|
else {
|
|
PathElement *temp = pre2->next;
|
|
pre2->next = temp->next;
|
|
delete temp;
|
|
}
|
|
}
|
|
|
|
compressRing(pre1);
|
|
compressRing(pre2);
|
|
|
|
// Recur
|
|
newnode = patchSplitSingle(newnode, st, len, pre1, dir, nrings1, nrings2);
|
|
newnode = patchSplitSingle(newnode, st, len, pre2, dir, nrings1, nrings2);
|
|
|
|
}
|
|
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("Return from PATCHSPLITSINGLE with \n");
|
|
dc_printf("Rings gourp 1:\n");
|
|
printPaths(nrings1);
|
|
dc_printf("Rings group 2:\n");
|
|
printPaths(nrings2);
|
|
#endif
|
|
|
|
return newnode;
|
|
}
|
|
|
|
Node *Octree::connectFace(Node *newnode, int st[3], int len, int dir,
|
|
PathElement *f1, PathElement *f2)
|
|
{
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("Call to CONNECTFACE with direction %d and length %d path: \n", dir, len);
|
|
dc_printf("Path(low side): \n");
|
|
printPath(f1);
|
|
// checkPath(f1);
|
|
dc_printf("Path(high side): \n");
|
|
printPath(f2);
|
|
// checkPath(f2);
|
|
#endif
|
|
|
|
// Setup 2D
|
|
int pos = st[dir] + len / 2;
|
|
int xdir = (dir + 1) % 3;
|
|
int ydir = (dir + 2) % 3;
|
|
|
|
// Use existing intersections on f1 and f2
|
|
int x1, y1, x2, y2;
|
|
float p1, q1, p2, q2;
|
|
|
|
getFacePoint(f2->next, dir, x1, y1, p1, q1);
|
|
getFacePoint(f2, dir, x2, y2, p2, q2);
|
|
|
|
float dx = x2 + p2 - x1 - p1;
|
|
float dy = y2 + q2 - y1 - q1;
|
|
|
|
// Do adapted Bresenham line drawing
|
|
float rx = p1, ry = q1;
|
|
int incx = 1, incy = 1;
|
|
int lx = x1, ly = y1;
|
|
int hx = x2, hy = y2;
|
|
int choice;
|
|
if (x2 < x1) {
|
|
incx = -1;
|
|
rx = 1 - rx;
|
|
lx = x2;
|
|
hx = x1;
|
|
}
|
|
if (y2 < y1) {
|
|
incy = -1;
|
|
ry = 1 - ry;
|
|
ly = y2;
|
|
hy = y1;
|
|
}
|
|
|
|
float sx = dx * incx;
|
|
float sy = dy * incy;
|
|
|
|
int ori[3];
|
|
ori[dir] = pos / mindimen;
|
|
ori[xdir] = x1;
|
|
ori[ydir] = y1;
|
|
int walkdir;
|
|
int inc;
|
|
float alpha;
|
|
|
|
PathElement *curEleN = f1;
|
|
PathElement *curEleP = f2->next;
|
|
Node *nodeN = NULL, *nodeP = NULL;
|
|
LeafNode *curN = locateLeaf(&newnode->internal, len, f1->pos);
|
|
LeafNode *curP = locateLeaf(&newnode->internal, len, f2->next->pos);
|
|
if (curN == NULL || curP == NULL) {
|
|
exit(0);
|
|
}
|
|
int stN[3], stP[3];
|
|
int lenN, lenP;
|
|
|
|
/* Unused code, leaving for posterity
|
|
|
|
float stpt[3], edpt[3];
|
|
stpt[dir] = edpt[dir] =(float) pos;
|
|
stpt[xdir] =(x1 + p1) * mindimen;
|
|
stpt[ydir] =(y1 + q1) * mindimen;
|
|
edpt[xdir] =(x2 + p2) * mindimen;
|
|
edpt[ydir] =(y2 + q2) * mindimen;
|
|
*/
|
|
while (ori[xdir] != x2 || ori[ydir] != y2) {
|
|
int next;
|
|
if (sy * (1 - rx) > sx * (1 - ry)) {
|
|
choice = 1;
|
|
next = ori[ydir] + incy;
|
|
if (next < ly || next > hy) {
|
|
choice = 4;
|
|
next = ori[xdir] + incx;
|
|
}
|
|
}
|
|
else {
|
|
choice = 2;
|
|
next = ori[xdir] + incx;
|
|
if (next < lx || next > hx) {
|
|
choice = 3;
|
|
next = ori[ydir] + incy;
|
|
}
|
|
}
|
|
|
|
if (choice & 1) {
|
|
ori[ydir] = next;
|
|
if (choice == 1) {
|
|
rx += (sy == 0 ? 0 : (1 - ry) * sx / sy);
|
|
ry = 0;
|
|
}
|
|
|
|
walkdir = 2;
|
|
inc = incy;
|
|
alpha = x2 < x1 ? 1 - rx : rx;
|
|
}
|
|
else {
|
|
ori[xdir] = next;
|
|
if (choice == 2) {
|
|
ry += (sx == 0 ? 0 : (1 - rx) * sy / sx);
|
|
rx = 0;
|
|
}
|
|
|
|
walkdir = 1;
|
|
inc = incx;
|
|
alpha = y2 < y1 ? 1 - ry : ry;
|
|
}
|
|
|
|
|
|
|
|
// Get the exact location of the marcher
|
|
int nori[3] = {ori[0] * mindimen, ori[1] * mindimen, ori[2] * mindimen};
|
|
float spt[3] = {(float) nori[0], (float) nori[1], (float) nori[2]};
|
|
spt[(dir + (3 - walkdir)) % 3] += alpha * mindimen;
|
|
if (inc < 0) {
|
|
spt[(dir + walkdir) % 3] += mindimen;
|
|
}
|
|
|
|
// dc_printf("new x,y: %d %d\n", ori[xdir] / edgelen, ori[ydir] / edgelen);
|
|
// dc_printf("nori: %d %d %d alpha: %f walkdir: %d\n", nori[0], nori[1], nori[2], alpha, walkdir);
|
|
// dc_printf("%f %f %f\n", spt[0], spt[1], spt[2]);
|
|
|
|
// Locate the current cells on both sides
|
|
newnode = locateCell(&newnode->internal, st, len, nori, dir, 1, nodeN, stN, lenN);
|
|
newnode = locateCell(&newnode->internal, st, len, nori, dir, 0, nodeP, stP, lenP);
|
|
|
|
updateParent(&newnode->internal, len, st);
|
|
|
|
int flag = 0;
|
|
// Add the cells to the rings and fill in the patch
|
|
PathElement *newEleN;
|
|
if (curEleN->pos[0] != stN[0] || curEleN->pos[1] != stN[1] || curEleN->pos[2] != stN[2]) {
|
|
if (curEleN->next->pos[0] != stN[0] || curEleN->next->pos[1] != stN[1] || curEleN->next->pos[2] != stN[2]) {
|
|
newEleN = new PathElement;
|
|
newEleN->next = curEleN->next;
|
|
newEleN->pos[0] = stN[0];
|
|
newEleN->pos[1] = stN[1];
|
|
newEleN->pos[2] = stN[2];
|
|
|
|
curEleN->next = newEleN;
|
|
}
|
|
else {
|
|
newEleN = curEleN->next;
|
|
}
|
|
curN = patchAdjacent(&newnode->internal, len, curEleN->pos, curN,
|
|
newEleN->pos, (LeafNode *)nodeN, walkdir,
|
|
inc, dir, 1, alpha);
|
|
|
|
curEleN = newEleN;
|
|
flag++;
|
|
}
|
|
|
|
PathElement *newEleP;
|
|
if (curEleP->pos[0] != stP[0] || curEleP->pos[1] != stP[1] || curEleP->pos[2] != stP[2]) {
|
|
if (f2->pos[0] != stP[0] || f2->pos[1] != stP[1] || f2->pos[2] != stP[2]) {
|
|
newEleP = new PathElement;
|
|
newEleP->next = curEleP;
|
|
newEleP->pos[0] = stP[0];
|
|
newEleP->pos[1] = stP[1];
|
|
newEleP->pos[2] = stP[2];
|
|
|
|
f2->next = newEleP;
|
|
}
|
|
else {
|
|
newEleP = f2;
|
|
}
|
|
curP = patchAdjacent(&newnode->internal, len, curEleP->pos, curP,
|
|
newEleP->pos, (LeafNode *)nodeP, walkdir,
|
|
inc, dir, 0, alpha);
|
|
|
|
|
|
|
|
curEleP = newEleP;
|
|
flag++;
|
|
}
|
|
|
|
|
|
/*
|
|
if(flag == 0) {
|
|
dc_printf("error: non-synchronized patching! at \n");
|
|
}
|
|
*/
|
|
}
|
|
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("Return from CONNECTFACE with \n");
|
|
dc_printf("Path(low side):\n");
|
|
printPath(f1);
|
|
checkPath(f1);
|
|
dc_printf("Path(high side):\n");
|
|
printPath(f2);
|
|
checkPath(f2);
|
|
#endif
|
|
|
|
|
|
return newnode;
|
|
}
|
|
|
|
LeafNode *Octree::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)
|
|
{
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("Before patching.\n");
|
|
printInfo(st1);
|
|
printInfo(st2);
|
|
dc_printf("-----------------%d %d %d; %d %d %d\n", st1[0], st2[1], st1[2], st2[0], st2[1], st2[2]);
|
|
#endif
|
|
|
|
// Get edge index on each leaf
|
|
int edgedir = (dir + (3 - walkdir)) % 3;
|
|
int incdir = (dir + walkdir) % 3;
|
|
int ind1 = (edgedir == 1 ? (dir + 3 - edgedir) % 3 - 1 : 2 - (dir + 3 - edgedir) % 3);
|
|
int ind2 = (edgedir == 1 ? (incdir + 3 - edgedir) % 3 - 1 : 2 - (incdir + 3 - edgedir) % 3);
|
|
|
|
int eind1 = ((edgedir << 2) | (side << ind1) | ((inc > 0 ? 1 : 0) << ind2));
|
|
int eind2 = ((edgedir << 2) | (side << ind1) | ((inc > 0 ? 0 : 1) << ind2));
|
|
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("Index 1: %d Alpha 1: %f Index 2: %d Alpha 2: %f\n", eind1, alpha, eind2, alpha);
|
|
/*
|
|
if(alpha < 0 || alpha > 1) {
|
|
dc_printf("Index 1: %d Alpha 1: %f Index 2: %d Alpha 2: %f\n", eind1, alpha, eind2, alpha);
|
|
printInfo(st1);
|
|
printInfo(st2);
|
|
}
|
|
*/
|
|
#endif
|
|
|
|
// Flip edge parity
|
|
LeafNode *nleaf1 = flipEdge(leaf1, eind1, alpha);
|
|
LeafNode *nleaf2 = flipEdge(leaf2, eind2, alpha);
|
|
|
|
// Update parent link
|
|
updateParent(node, len, st1, nleaf1);
|
|
updateParent(node, len, st2, nleaf2);
|
|
// updateParent(nleaf1, mindimen, st1);
|
|
// updateParent(nleaf2, mindimen, st2);
|
|
|
|
/*
|
|
float m[3];
|
|
dc_printf("Adding new point: %f %f %f\n", spt[0], spt[1], spt[2]);
|
|
getMinimizer(leaf1, m);
|
|
dc_printf("Cell %d now has minimizer %f %f %f\n", leaf1, m[0], m[1], m[2]);
|
|
getMinimizer(leaf2, m);
|
|
dc_printf("Cell %d now has minimizer %f %f %f\n", leaf2, m[0], m[1], m[2]);
|
|
*/
|
|
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("After patching.\n");
|
|
printInfo(st1);
|
|
printInfo(st2);
|
|
#endif
|
|
return nleaf2;
|
|
}
|
|
|
|
Node *Octree::locateCell(InternalNode *node, int st[3], int len, int ori[3], int dir, int side, Node *& rleaf, int rst[3], int& rlen)
|
|
{
|
|
#ifdef IN_DEBUG_MODE
|
|
// dc_printf("Call to LOCATECELL with node ");
|
|
// printNode(node);
|
|
#endif
|
|
|
|
int i;
|
|
len >>= 1;
|
|
int ind = 0;
|
|
for (i = 0; i < 3; i++) {
|
|
ind <<= 1;
|
|
if (i == dir && side == 1) {
|
|
ind |= (ori[i] <= (st[i] + len) ? 0 : 1);
|
|
}
|
|
else {
|
|
ind |= (ori[i] < (st[i] + len) ? 0 : 1);
|
|
}
|
|
}
|
|
|
|
#ifdef IN_DEBUG_MODE
|
|
// dc_printf("In LOCATECELL index of ori(%d %d %d) with dir %d side %d in st(%d %d %d, %d) is: %d\n",
|
|
// ori[0], ori[1], ori[2], dir, side, st[0], st[1], st[2], len, ind);
|
|
#endif
|
|
|
|
rst[0] = st[0] + vertmap[ind][0] * len;
|
|
rst[1] = st[1] + vertmap[ind][1] * len;
|
|
rst[2] = st[2] + vertmap[ind][2] * len;
|
|
|
|
if (node->has_child(ind)) {
|
|
int count = node->get_child_count(ind);
|
|
Node *chd = node->get_child(count);
|
|
if (node->is_child_leaf(ind)) {
|
|
rleaf = chd;
|
|
rlen = len;
|
|
}
|
|
else {
|
|
// Recur
|
|
node->set_child(count, locateCell(&chd->internal, rst, len, ori, dir, side, rleaf, rst, rlen));
|
|
}
|
|
}
|
|
else {
|
|
// Create a new child here
|
|
if (len == mindimen) {
|
|
LeafNode *chd = createLeaf(0);
|
|
node = addChild(node, ind, (Node *)chd, 1);
|
|
rleaf = (Node *)chd;
|
|
rlen = len;
|
|
}
|
|
else {
|
|
// Subdivide the empty cube
|
|
InternalNode *chd = createInternal(0);
|
|
node = addChild(node, ind,
|
|
locateCell(chd, rst, len, ori, dir, side, rleaf, rst, rlen), 0);
|
|
}
|
|
}
|
|
|
|
#ifdef IN_DEBUG_MODE
|
|
// dc_printf("Return from LOCATECELL with node ");
|
|
// printNode(newnode);
|
|
#endif
|
|
return (Node *)node;
|
|
}
|
|
|
|
void Octree::checkElement(PathElement * /*ele*/)
|
|
{
|
|
/*
|
|
if(ele != NULL && locateLeafCheck(ele->pos) != ele->node) {
|
|
dc_printf("Screwed! at pos: %d %d %d\n", ele->pos[0]>>minshift, ele->pos[1]>>minshift, ele->pos[2]>>minshift);
|
|
exit(0);
|
|
}
|
|
*/
|
|
}
|
|
|
|
void Octree::checkPath(PathElement *path)
|
|
{
|
|
PathElement *n = path;
|
|
int same = 0;
|
|
while (n && (same == 0 || n != path)) {
|
|
same++;
|
|
checkElement(n);
|
|
n = n->next;
|
|
}
|
|
|
|
}
|
|
|
|
void Octree::testFacePoint(PathElement *e1, PathElement *e2)
|
|
{
|
|
int i;
|
|
PathElement *e = NULL;
|
|
for (i = 0; i < 3; i++) {
|
|
if (e1->pos[i] != e2->pos[i]) {
|
|
if (e1->pos[i] < e2->pos[i]) {
|
|
e = e2;
|
|
}
|
|
else {
|
|
e = e1;
|
|
}
|
|
break;
|
|
}
|
|
}
|
|
|
|
int x, y;
|
|
float p, q;
|
|
dc_printf("Test.");
|
|
getFacePoint(e, i, x, y, p, q);
|
|
}
|
|
|
|
void Octree::getFacePoint(PathElement *leaf, int dir, int& x, int& y, float& p, float& q)
|
|
{
|
|
// Find average intersections
|
|
float avg[3] = {0, 0, 0};
|
|
float off[3];
|
|
int num = 0, num2 = 0;
|
|
|
|
LeafNode *leafnode = locateLeaf(leaf->pos);
|
|
for (int i = 0; i < 4; i++) {
|
|
int edgeind = faceMap[dir * 2][i];
|
|
int nst[3];
|
|
for (int j = 0; j < 3; j++) {
|
|
nst[j] = leaf->pos[j] + mindimen * vertmap[edgemap[edgeind][0]][j];
|
|
}
|
|
|
|
if (getEdgeIntersectionByIndex(nst, edgeind / 4, off, 1)) {
|
|
avg[0] += off[0];
|
|
avg[1] += off[1];
|
|
avg[2] += off[2];
|
|
num++;
|
|
}
|
|
if (getEdgeParity(leafnode, edgeind)) {
|
|
num2++;
|
|
}
|
|
}
|
|
if (num == 0) {
|
|
dc_printf("Wrong! dir: %d pos: %d %d %d num: %d\n", dir, leaf->pos[0] >> minshift, leaf->pos[1] >> minshift, leaf->pos[2] >> minshift, num2);
|
|
avg[0] = (float) leaf->pos[0];
|
|
avg[1] = (float) leaf->pos[1];
|
|
avg[2] = (float) leaf->pos[2];
|
|
}
|
|
else {
|
|
|
|
avg[0] /= num;
|
|
avg[1] /= num;
|
|
avg[2] /= num;
|
|
|
|
//avg[0] =(float) leaf->pos[0];
|
|
//avg[1] =(float) leaf->pos[1];
|
|
//avg[2] =(float) leaf->pos[2];
|
|
}
|
|
|
|
int xdir = (dir + 1) % 3;
|
|
int ydir = (dir + 2) % 3;
|
|
|
|
float xf = avg[xdir];
|
|
float yf = avg[ydir];
|
|
|
|
#ifdef IN_DEBUG_MODE
|
|
// Is it outside?
|
|
// PathElement* leaf = leaf1->len < leaf2->len ? leaf1 : leaf2;
|
|
/*
|
|
float* m =(leaf == leaf1 ? m1 : m2);
|
|
if(xf < leaf->pos[xdir] ||
|
|
yf < leaf->pos[ydir] ||
|
|
xf > leaf->pos[xdir] + leaf->len ||
|
|
yf > leaf->pos[ydir] + leaf->len) {
|
|
dc_printf("Outside cube(%d %d %d), %d : %d %d %f %f\n", leaf->pos[0], leaf->pos[1], leaf->pos[2], leaf->len,
|
|
pos, dir, xf, yf);
|
|
|
|
// For now, snap to cell
|
|
xf = m[xdir];
|
|
yf = m[ydir];
|
|
}
|
|
*/
|
|
|
|
/*
|
|
if(alpha < 0 || alpha > 1 ||
|
|
xf < leaf->pos[xdir] || xf > leaf->pos[xdir] + leaf->len ||
|
|
yf < leaf->pos[ydir] || yf > leaf->pos[ydir] + leaf->len) {
|
|
dc_printf("Alpha: %f Address: %d and %d\n", alpha, leaf1->node, leaf2->node);
|
|
dc_printf("GETFACEPOINT result:(%d %d %d) %d min:(%f %f %f);(%d %d %d) %d min:(%f %f %f).\n",
|
|
leaf1->pos[0], leaf1->pos[1], leaf1->pos[2], leaf1->len, m1[0], m1[1], m1[2],
|
|
leaf2->pos[0], leaf2->pos[1], leaf2->pos[2], leaf2->len, m2[0], m2[1], m2[2]);
|
|
dc_printf("Face point at dir %d pos %d: %f %f\n", dir, pos, xf, yf);
|
|
}
|
|
*/
|
|
#endif
|
|
|
|
|
|
// Get the integer and float part
|
|
x = ((leaf->pos[xdir]) >> minshift);
|
|
y = ((leaf->pos[ydir]) >> minshift);
|
|
|
|
p = (xf - leaf->pos[xdir]) / mindimen;
|
|
q = (yf - leaf->pos[ydir]) / mindimen;
|
|
|
|
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("Face point at dir %d : %f %f\n", dir, xf, yf);
|
|
#endif
|
|
}
|
|
|
|
int Octree::findPair(PathElement *head, int pos, int dir, PathElement *& pre1, PathElement *& pre2)
|
|
{
|
|
int side = getSide(head, pos, dir);
|
|
PathElement *cur = head;
|
|
PathElement *anchor;
|
|
PathElement *ppre1, *ppre2;
|
|
|
|
// Start from this face, find a pair
|
|
anchor = cur;
|
|
ppre1 = cur;
|
|
cur = cur->next;
|
|
while (cur != anchor && (getSide(cur, pos, dir) == side)) {
|
|
ppre1 = cur;
|
|
cur = cur->next;
|
|
}
|
|
if (cur == anchor) {
|
|
// No pair found
|
|
return side;
|
|
}
|
|
|
|
side = getSide(cur, pos, dir);
|
|
ppre2 = cur;
|
|
cur = cur->next;
|
|
while (getSide(cur, pos, dir) == side) {
|
|
ppre2 = cur;
|
|
cur = cur->next;
|
|
}
|
|
|
|
|
|
// Switch pre1 and pre2 if we start from the higher side
|
|
if (side == -1) {
|
|
cur = ppre1;
|
|
ppre1 = ppre2;
|
|
ppre2 = cur;
|
|
}
|
|
|
|
pre1 = ppre1;
|
|
pre2 = ppre2;
|
|
|
|
return 0;
|
|
}
|
|
|
|
int Octree::getSide(PathElement *e, int pos, int dir)
|
|
{
|
|
return (e->pos[dir] < pos ? -1 : 1);
|
|
}
|
|
|
|
int Octree::isEqual(PathElement *e1, PathElement *e2)
|
|
{
|
|
return (e1->pos[0] == e2->pos[0] && e1->pos[1] == e2->pos[1] && e1->pos[2] == e2->pos[2]);
|
|
}
|
|
|
|
void Octree::compressRing(PathElement *& ring)
|
|
{
|
|
if (ring == NULL) {
|
|
return;
|
|
}
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("Call to COMPRESSRING with path: \n");
|
|
printPath(ring);
|
|
#endif
|
|
|
|
PathElement *cur = ring->next->next;
|
|
PathElement *pre = ring->next;
|
|
PathElement *prepre = ring;
|
|
PathElement *anchor = prepre;
|
|
|
|
do {
|
|
while (isEqual(cur, prepre)) {
|
|
// Delete
|
|
if (cur == prepre) {
|
|
// The ring has shrinked to a point
|
|
delete pre;
|
|
delete cur;
|
|
anchor = NULL;
|
|
break;
|
|
}
|
|
else {
|
|
prepre->next = cur->next;
|
|
delete pre;
|
|
delete cur;
|
|
pre = prepre->next;
|
|
cur = pre->next;
|
|
anchor = prepre;
|
|
}
|
|
}
|
|
|
|
if (anchor == NULL) {
|
|
break;
|
|
}
|
|
|
|
prepre = pre;
|
|
pre = cur;
|
|
cur = cur->next;
|
|
} while (prepre != anchor);
|
|
|
|
ring = anchor;
|
|
|
|
#ifdef IN_DEBUG_MODE
|
|
dc_printf("Return from COMPRESSRING with path: \n");
|
|
printPath(ring);
|
|
#endif
|
|
}
|
|
|
|
void Octree::buildSigns()
|
|
{
|
|
// First build a lookup table
|
|
// dc_printf("Building up look up table...\n");
|
|
int size = 1 << 12;
|
|
unsigned char table[1 << 12];
|
|
for (int i = 0; i < size; i++) {
|
|
table[i] = 0;
|
|
}
|
|
for (int i = 0; i < 256; i++) {
|
|
int ind = 0;
|
|
for (int j = 11; j >= 0; j--) {
|
|
ind <<= 1;
|
|
if (((i >> edgemap[j][0]) & 1) ^ ((i >> edgemap[j][1]) & 1)) {
|
|
ind |= 1;
|
|
}
|
|
}
|
|
|
|
table[ind] = i;
|
|
}
|
|
|
|
// Next, traverse the grid
|
|
int sg = 1;
|
|
int cube[8];
|
|
buildSigns(table, root, 0, sg, cube);
|
|
}
|
|
|
|
void Octree::buildSigns(unsigned char table[], Node *node, int isLeaf, int sg, int rvalue[8])
|
|
{
|
|
if (node == NULL) {
|
|
for (int i = 0; i < 8; i++) {
|
|
rvalue[i] = sg;
|
|
}
|
|
return;
|
|
}
|
|
|
|
if (isLeaf == 0) {
|
|
// Internal node
|
|
Node *chd[8];
|
|
int leaf[8];
|
|
node->internal.fill_children(chd, leaf);
|
|
|
|
// Get the signs at the corners of the first cube
|
|
rvalue[0] = sg;
|
|
int oris[8];
|
|
buildSigns(table, chd[0], leaf[0], sg, oris);
|
|
|
|
// Get the rest
|
|
int cube[8];
|
|
for (int i = 1; i < 8; i++) {
|
|
buildSigns(table, chd[i], leaf[i], oris[i], cube);
|
|
rvalue[i] = cube[i];
|
|
}
|
|
|
|
}
|
|
else {
|
|
// Leaf node
|
|
generateSigns(&node->leaf, table, sg);
|
|
|
|
for (int i = 0; i < 8; i++) {
|
|
rvalue[i] = getSign(&node->leaf, i);
|
|
}
|
|
}
|
|
}
|
|
|
|
void Octree::floodFill()
|
|
{
|
|
// int threshold =(int)((dimen/mindimen) *(dimen/mindimen) * 0.5f);
|
|
int st[3] = {0, 0, 0};
|
|
|
|
// First, check for largest component
|
|
// size stored in -threshold
|
|
clearProcessBits(root, maxDepth);
|
|
int threshold = floodFill(root, st, dimen, maxDepth, 0);
|
|
|
|
// Next remove
|
|
dc_printf("Largest component: %d\n", threshold);
|
|
threshold *= thresh;
|
|
dc_printf("Removing all components smaller than %d\n", threshold);
|
|
|
|
int st2[3] = {0, 0, 0};
|
|
clearProcessBits(root, maxDepth);
|
|
floodFill(root, st2, dimen, maxDepth, threshold);
|
|
|
|
}
|
|
|
|
void Octree::clearProcessBits(Node *node, int height)
|
|
{
|
|
int i;
|
|
|
|
if (height == 0) {
|
|
// Leaf cell,
|
|
for (i = 0; i < 12; i++) {
|
|
setOutProcess(&node->leaf, i);
|
|
}
|
|
}
|
|
else {
|
|
// Internal cell, recur
|
|
int count = 0;
|
|
for (i = 0; i < 8; i++) {
|
|
if (node->internal.has_child(i)) {
|
|
clearProcessBits(node->internal.get_child(count), height - 1);
|
|
count++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
int Octree::floodFill(LeafNode *leaf, int st[3], int len, int /*height*/, int threshold)
|
|
{
|
|
int i, j;
|
|
int maxtotal = 0;
|
|
|
|
// Leaf cell,
|
|
int par, inp;
|
|
|
|
// Test if the leaf has intersection edges
|
|
for (i = 0; i < 12; i++) {
|
|
par = getEdgeParity(leaf, i);
|
|
inp = isInProcess(leaf, i);
|
|
|
|
if (par == 1 && inp == 0) {
|
|
// Intersection edge, hasn't been processed
|
|
// Let's start filling
|
|
GridQueue *queue = new GridQueue();
|
|
int total = 1;
|
|
|
|
// Set to in process
|
|
int mst[3];
|
|
mst[0] = st[0] + vertmap[edgemap[i][0]][0] * len;
|
|
mst[1] = st[1] + vertmap[edgemap[i][0]][1] * len;
|
|
mst[2] = st[2] + vertmap[edgemap[i][0]][2] * len;
|
|
int mdir = i / 4;
|
|
setInProcessAll(mst, mdir);
|
|
|
|
// Put this edge into queue
|
|
queue->pushQueue(mst, mdir);
|
|
|
|
// Queue processing
|
|
int nst[3], dir;
|
|
while (queue->popQueue(nst, dir) == 1) {
|
|
// dc_printf("nst: %d %d %d, dir: %d\n", nst[0]/mindimen, nst[1]/mindimen, nst[2]/mindimen, dir);
|
|
// locations
|
|
int stMask[3][3] = {
|
|
{0, 0 - len, 0 - len},
|
|
{0 - len, 0, 0 - len},
|
|
{0 - len, 0 - len, 0}
|
|
};
|
|
int cst[2][3];
|
|
for (j = 0; j < 3; j++) {
|
|
cst[0][j] = nst[j];
|
|
cst[1][j] = nst[j] + stMask[dir][j];
|
|
}
|
|
|
|
// cells
|
|
LeafNode *cs[2];
|
|
for (j = 0; j < 2; j++) {
|
|
cs[j] = locateLeaf(cst[j]);
|
|
}
|
|
|
|
// Middle sign
|
|
int s = getSign(cs[0], 0);
|
|
|
|
// Masks
|
|
int fcCells[4] = {1, 0, 1, 0};
|
|
int fcEdges[3][4][3] = {
|
|
{{9, 2, 11}, {8, 1, 10}, {5, 1, 7}, {4, 2, 6}},
|
|
{{10, 6, 11}, {8, 5, 9}, {1, 5, 3}, {0, 6, 2}},
|
|
{{6, 10, 7}, {4, 9, 5}, {2, 9, 3}, {0, 10, 1}}
|
|
};
|
|
|
|
// Search for neighboring connected intersection edges
|
|
for (int find = 0; find < 4; find++) {
|
|
int cind = fcCells[find];
|
|
int eind, edge;
|
|
if (s == 0) {
|
|
// Original order
|
|
for (eind = 0; eind < 3; eind++) {
|
|
edge = fcEdges[dir][find][eind];
|
|
if (getEdgeParity(cs[cind], edge) == 1) {
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
else {
|
|
// Inverse order
|
|
for (eind = 2; eind >= 0; eind--) {
|
|
edge = fcEdges[dir][find][eind];
|
|
if (getEdgeParity(cs[cind], edge) == 1) {
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (eind == 3 || eind == -1) {
|
|
dc_printf("Wrong! this is not a consistent sign. %d\n", eind);
|
|
}
|
|
else {
|
|
int est[3];
|
|
est[0] = cst[cind][0] + vertmap[edgemap[edge][0]][0] * len;
|
|
est[1] = cst[cind][1] + vertmap[edgemap[edge][0]][1] * len;
|
|
est[2] = cst[cind][2] + vertmap[edgemap[edge][0]][2] * len;
|
|
int edir = edge / 4;
|
|
|
|
if (isInProcess(cs[cind], edge) == 0) {
|
|
setInProcessAll(est, edir);
|
|
queue->pushQueue(est, edir);
|
|
// dc_printf("Pushed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir);
|
|
total++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
dc_printf("Size of component: %d ", total);
|
|
|
|
if (threshold == 0) {
|
|
// Measuring stage
|
|
if (total > maxtotal) {
|
|
maxtotal = total;
|
|
}
|
|
dc_printf(".\n");
|
|
delete queue;
|
|
continue;
|
|
}
|
|
|
|
if (total >= threshold) {
|
|
dc_printf("Maintained.\n");
|
|
delete queue;
|
|
continue;
|
|
}
|
|
dc_printf("Less then %d, removing...\n", threshold);
|
|
|
|
// We have to remove this noise
|
|
|
|
// Flip parity
|
|
// setOutProcessAll(mst, mdir);
|
|
flipParityAll(mst, mdir);
|
|
|
|
// Put this edge into queue
|
|
queue->pushQueue(mst, mdir);
|
|
|
|
// Queue processing
|
|
while (queue->popQueue(nst, dir) == 1) {
|
|
// dc_printf("nst: %d %d %d, dir: %d\n", nst[0]/mindimen, nst[1]/mindimen, nst[2]/mindimen, dir);
|
|
// locations
|
|
int stMask[3][3] = {
|
|
{0, 0 - len, 0 - len},
|
|
{0 - len, 0, 0 - len},
|
|
{0 - len, 0 - len, 0}
|
|
};
|
|
int cst[2][3];
|
|
for (j = 0; j < 3; j++) {
|
|
cst[0][j] = nst[j];
|
|
cst[1][j] = nst[j] + stMask[dir][j];
|
|
}
|
|
|
|
// cells
|
|
LeafNode *cs[2];
|
|
for (j = 0; j < 2; j++)
|
|
cs[j] = locateLeaf(cst[j]);
|
|
|
|
// Middle sign
|
|
int s = getSign(cs[0], 0);
|
|
|
|
// Masks
|
|
int fcCells[4] = {1, 0, 1, 0};
|
|
int fcEdges[3][4][3] = {
|
|
{{9, 2, 11}, {8, 1, 10}, {5, 1, 7}, {4, 2, 6}},
|
|
{{10, 6, 11}, {8, 5, 9}, {1, 5, 3}, {0, 6, 2}},
|
|
{{6, 10, 7}, {4, 9, 5}, {2, 9, 3}, {0, 10, 1}}
|
|
};
|
|
|
|
// Search for neighboring connected intersection edges
|
|
for (int find = 0; find < 4; find++) {
|
|
int cind = fcCells[find];
|
|
int eind, edge;
|
|
if (s == 0) {
|
|
// Original order
|
|
for (eind = 0; eind < 3; eind++) {
|
|
edge = fcEdges[dir][find][eind];
|
|
if (isInProcess(cs[cind], edge) == 1) {
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
else {
|
|
// Inverse order
|
|
for (eind = 2; eind >= 0; eind--) {
|
|
edge = fcEdges[dir][find][eind];
|
|
if (isInProcess(cs[cind], edge) == 1) {
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (eind == 3 || eind == -1) {
|
|
dc_printf("Wrong! this is not a consistent sign. %d\n", eind);
|
|
}
|
|
else {
|
|
int est[3];
|
|
est[0] = cst[cind][0] + vertmap[edgemap[edge][0]][0] * len;
|
|
est[1] = cst[cind][1] + vertmap[edgemap[edge][0]][1] * len;
|
|
est[2] = cst[cind][2] + vertmap[edgemap[edge][0]][2] * len;
|
|
int edir = edge / 4;
|
|
|
|
if (getEdgeParity(cs[cind], edge) == 1) {
|
|
flipParityAll(est, edir);
|
|
queue->pushQueue(est, edir);
|
|
// dc_printf("Pushed: est: %d %d %d, edir: %d\n", est[0]/len, est[1]/len, est[2]/len, edir);
|
|
total++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
delete queue;
|
|
}
|
|
}
|
|
|
|
return maxtotal;
|
|
}
|
|
|
|
int Octree::floodFill(Node *node, int st[3], int len, int height, int threshold)
|
|
{
|
|
int i;
|
|
int maxtotal = 0;
|
|
|
|
if (height == 0) {
|
|
maxtotal = floodFill(&node->leaf, st, len, height, threshold);
|
|
}
|
|
else {
|
|
// Internal cell, recur
|
|
int count = 0;
|
|
len >>= 1;
|
|
for (i = 0; i < 8; i++) {
|
|
if (node->internal.has_child(i)) {
|
|
int nst[3];
|
|
nst[0] = st[0] + vertmap[i][0] * len;
|
|
nst[1] = st[1] + vertmap[i][1] * len;
|
|
nst[2] = st[2] + vertmap[i][2] * len;
|
|
|
|
int d = floodFill(node->internal.get_child(count), nst, len, height - 1, threshold);
|
|
if (d > maxtotal) {
|
|
maxtotal = d;
|
|
}
|
|
count++;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
return maxtotal;
|
|
|
|
}
|
|
|
|
void Octree::writeOut()
|
|
{
|
|
int numQuads = 0;
|
|
int numVertices = 0;
|
|
int numEdges = 0;
|
|
|
|
countIntersection(root, maxDepth, numQuads, numVertices, numEdges);
|
|
|
|
dc_printf("Vertices counted: %d Polys counted: %d \n", numVertices, numQuads);
|
|
output_mesh = alloc_output(numVertices, numQuads);
|
|
int offset = 0;
|
|
int st[3] = {0, 0, 0};
|
|
|
|
// First, output vertices
|
|
offset = 0;
|
|
actualVerts = 0;
|
|
actualQuads = 0;
|
|
|
|
generateMinimizer(root, st, dimen, maxDepth, offset);
|
|
cellProcContour(root, 0, maxDepth);
|
|
dc_printf("Vertices written: %d Quads written: %d \n", offset, actualQuads);
|
|
}
|
|
|
|
void Octree::countIntersection(Node *node, int height, int& nedge, int& ncell, int& nface)
|
|
{
|
|
if (height > 0) {
|
|
int total = node->internal.get_num_children();
|
|
for (int i = 0; i < total; i++) {
|
|
countIntersection(node->internal.get_child(i), height - 1, nedge, ncell, nface);
|
|
}
|
|
}
|
|
else {
|
|
nedge += getNumEdges2(&node->leaf);
|
|
|
|
int smask = getSignMask(&node->leaf);
|
|
|
|
if (use_manifold) {
|
|
int comps = manifold_table[smask].comps;
|
|
ncell += comps;
|
|
}
|
|
else {
|
|
if (smask > 0 && smask < 255) {
|
|
ncell++;
|
|
}
|
|
}
|
|
|
|
for (int i = 0; i < 3; i++) {
|
|
if (getFaceEdgeNum(&node->leaf, i * 2)) {
|
|
nface++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/* from http://eigen.tuxfamily.org/bz/show_bug.cgi?id=257 */
|
|
template<typename _Matrix_Type_>
|
|
void pseudoInverse(const _Matrix_Type_ &a,
|
|
_Matrix_Type_ &result,
|
|
double epsilon = std::numeric_limits<typename _Matrix_Type_::Scalar>::epsilon())
|
|
{
|
|
Eigen::JacobiSVD< _Matrix_Type_ > svd = a.jacobiSvd(Eigen::ComputeFullU |
|
|
Eigen::ComputeFullV);
|
|
|
|
typename _Matrix_Type_::Scalar tolerance = epsilon * std::max(a.cols(),
|
|
a.rows()) *
|
|
svd.singularValues().array().abs().maxCoeff();
|
|
|
|
result = svd.matrixV() *
|
|
_Matrix_Type_((svd.singularValues().array().abs() >
|
|
tolerance).select(svd.singularValues().
|
|
array().inverse(), 0)).asDiagonal() *
|
|
svd.matrixU().adjoint();
|
|
}
|
|
|
|
static void solve_least_squares(const float halfA[], const float b[],
|
|
const float midpoint[], float rvalue[])
|
|
{
|
|
/* calculate pseudo-inverse */
|
|
Eigen::MatrixXf A(3, 3), pinv(3, 3);
|
|
A << halfA[0], halfA[1], halfA[2],
|
|
halfA[1], halfA[3], halfA[4],
|
|
halfA[2], halfA[4], halfA[5];
|
|
pseudoInverse(A, pinv);
|
|
|
|
Eigen::Vector3f b2(b), mp(midpoint), result;
|
|
b2 = b2 + A * -mp;
|
|
result = pinv * b2 + mp;
|
|
|
|
for (int i = 0; i < 3; i++)
|
|
rvalue[i] = result(i);
|
|
}
|
|
|
|
static void mass_point(float mp[3], const float pts[12][3], const int parity[12])
|
|
{
|
|
int ec = 0;
|
|
|
|
for (int i = 0; i < 12; i++) {
|
|
if (parity[i]) {
|
|
const float *p = pts[i];
|
|
|
|
mp[0] += p[0];
|
|
mp[1] += p[1];
|
|
mp[2] += p[2];
|
|
|
|
ec++;
|
|
}
|
|
}
|
|
|
|
if (ec == 0) {
|
|
return;
|
|
}
|
|
mp[0] /= ec;
|
|
mp[1] /= ec;
|
|
mp[2] /= ec;
|
|
}
|
|
|
|
static void minimize(float rvalue[3], float mp[3], const float pts[12][3],
|
|
const float norms[12][3], const int parity[12])
|
|
{
|
|
float ata[6] = {0, 0, 0, 0, 0, 0};
|
|
float atb[3] = {0, 0, 0};
|
|
int ec = 0;
|
|
|
|
for (int i = 0; i < 12; i++) {
|
|
// if(getEdgeParity(leaf, i))
|
|
if (parity[i]) {
|
|
const float *norm = norms[i];
|
|
const float *p = pts[i];
|
|
|
|
// QEF
|
|
ata[0] += (float)(norm[0] * norm[0]);
|
|
ata[1] += (float)(norm[0] * norm[1]);
|
|
ata[2] += (float)(norm[0] * norm[2]);
|
|
ata[3] += (float)(norm[1] * norm[1]);
|
|
ata[4] += (float)(norm[1] * norm[2]);
|
|
ata[5] += (float)(norm[2] * norm[2]);
|
|
|
|
const float pn = p[0] * norm[0] + p[1] * norm[1] + p[2] * norm[2];
|
|
|
|
atb[0] += (float)(norm[0] * pn);
|
|
atb[1] += (float)(norm[1] * pn);
|
|
atb[2] += (float)(norm[2] * pn);
|
|
|
|
// Minimizer
|
|
mp[0] += p[0];
|
|
mp[1] += p[1];
|
|
mp[2] += p[2];
|
|
|
|
ec++;
|
|
}
|
|
}
|
|
|
|
if (ec == 0) {
|
|
return;
|
|
}
|
|
mp[0] /= ec;
|
|
mp[1] /= ec;
|
|
mp[2] /= ec;
|
|
|
|
// Solve least squares
|
|
solve_least_squares(ata, atb, mp, rvalue);
|
|
}
|
|
|
|
void Octree::computeMinimizer(const LeafNode *leaf, int st[3], int len,
|
|
float rvalue[3]) const
|
|
{
|
|
// First, gather all edge intersections
|
|
float pts[12][3], norms[12][3];
|
|
int parity[12];
|
|
fillEdgeIntersections(leaf, st, len, pts, norms, parity);
|
|
|
|
switch (mode) {
|
|
case DUALCON_CENTROID:
|
|
rvalue[0] = (float) st[0] + len / 2;
|
|
rvalue[1] = (float) st[1] + len / 2;
|
|
rvalue[2] = (float) st[2] + len / 2;
|
|
break;
|
|
|
|
case DUALCON_MASS_POINT:
|
|
rvalue[0] = rvalue[1] = rvalue[2] = 0;
|
|
mass_point(rvalue, pts, parity);
|
|
break;
|
|
|
|
default: {
|
|
// Sharp features */
|
|
|
|
// construct QEF and minimizer
|
|
float mp[3] = {0, 0, 0};
|
|
minimize(rvalue, mp, pts, norms, parity);
|
|
|
|
/* Restraining the location of the minimizer */
|
|
float nh1 = hermite_num * len;
|
|
float nh2 = (1 + hermite_num) * len;
|
|
|
|
if (rvalue[0] < st[0] - nh1 ||
|
|
rvalue[1] < st[1] - nh1 ||
|
|
rvalue[2] < st[2] - nh1 ||
|
|
|
|
rvalue[0] > st[0] + nh2 ||
|
|
rvalue[1] > st[1] + nh2 ||
|
|
rvalue[2] > st[2] + nh2)
|
|
{
|
|
// Use mass point instead
|
|
rvalue[0] = mp[0];
|
|
rvalue[1] = mp[1];
|
|
rvalue[2] = mp[2];
|
|
}
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
void Octree::generateMinimizer(Node *node, int st[3], int len, int height, int& offset)
|
|
{
|
|
int i, j;
|
|
|
|
if (height == 0) {
|
|
// Leaf cell, generate
|
|
|
|
// First, find minimizer
|
|
float rvalue[3];
|
|
rvalue[0] = (float) st[0] + len / 2;
|
|
rvalue[1] = (float) st[1] + len / 2;
|
|
rvalue[2] = (float) st[2] + len / 2;
|
|
computeMinimizer(&node->leaf, st, len, rvalue);
|
|
|
|
// Update
|
|
//float fnst[3];
|
|
for (j = 0; j < 3; j++) {
|
|
rvalue[j] = rvalue[j] * range / dimen + origin[j];
|
|
//fnst[j] = st[j] * range / dimen + origin[j];
|
|
}
|
|
|
|
int mult = 0, smask = getSignMask(&node->leaf);
|
|
|
|
if (use_manifold) {
|
|
mult = manifold_table[smask].comps;
|
|
}
|
|
else {
|
|
if (smask > 0 && smask < 255) {
|
|
mult = 1;
|
|
}
|
|
}
|
|
|
|
for (j = 0; j < mult; j++) {
|
|
add_vert(output_mesh, rvalue);
|
|
}
|
|
|
|
// Store the index
|
|
setMinimizerIndex(&node->leaf, offset);
|
|
|
|
offset += mult;
|
|
}
|
|
else {
|
|
// Internal cell, recur
|
|
int count = 0;
|
|
len >>= 1;
|
|
for (i = 0; i < 8; i++) {
|
|
if (node->internal.has_child(i)) {
|
|
int nst[3];
|
|
nst[0] = st[0] + vertmap[i][0] * len;
|
|
nst[1] = st[1] + vertmap[i][1] * len;
|
|
nst[2] = st[2] + vertmap[i][2] * len;
|
|
|
|
generateMinimizer(node->internal.get_child(count),
|
|
nst, len, height - 1, offset);
|
|
count++;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void Octree::processEdgeWrite(Node *node[4], int /*depth*/[4], int /*maxdep*/, int dir)
|
|
{
|
|
//int color = 0;
|
|
|
|
int i = 3;
|
|
{
|
|
if (getEdgeParity((LeafNode *)(node[i]), processEdgeMask[dir][i])) {
|
|
int flip = 0;
|
|
int edgeind = processEdgeMask[dir][i];
|
|
if (getSign((LeafNode *)node[i], edgemap[edgeind][1]) > 0) {
|
|
flip = 1;
|
|
}
|
|
|
|
int num = 0;
|
|
{
|
|
int ind[8];
|
|
if (use_manifold) {
|
|
int vind[2];
|
|
int seq[4] = {0, 1, 3, 2};
|
|
for (int k = 0; k < 4; k++) {
|
|
getMinimizerIndices((LeafNode *)(node[seq[k]]), processEdgeMask[dir][seq[k]], vind);
|
|
ind[num] = vind[0];
|
|
num++;
|
|
|
|
if (vind[1] != -1) {
|
|
ind[num] = vind[1];
|
|
num++;
|
|
if (flip == 0) {
|
|
ind[num - 1] = vind[0];
|
|
ind[num - 2] = vind[1];
|
|
}
|
|
}
|
|
}
|
|
|
|
/* we don't use the manifold option, but if it is
|
|
ever enabled again note that it can output
|
|
non-quads */
|
|
}
|
|
else {
|
|
if (flip) {
|
|
ind[0] = getMinimizerIndex((LeafNode *)(node[2]));
|
|
ind[1] = getMinimizerIndex((LeafNode *)(node[3]));
|
|
ind[2] = getMinimizerIndex((LeafNode *)(node[1]));
|
|
ind[3] = getMinimizerIndex((LeafNode *)(node[0]));
|
|
}
|
|
else {
|
|
ind[0] = getMinimizerIndex((LeafNode *)(node[0]));
|
|
ind[1] = getMinimizerIndex((LeafNode *)(node[1]));
|
|
ind[2] = getMinimizerIndex((LeafNode *)(node[3]));
|
|
ind[3] = getMinimizerIndex((LeafNode *)(node[2]));
|
|
}
|
|
|
|
add_quad(output_mesh, ind);
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
else {
|
|
return;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Octree::edgeProcContour(Node *node[4], int leaf[4], int depth[4], int maxdep, int dir)
|
|
{
|
|
if (!(node[0] && node[1] && node[2] && node[3])) {
|
|
return;
|
|
}
|
|
if (leaf[0] && leaf[1] && leaf[2] && leaf[3]) {
|
|
processEdgeWrite(node, depth, maxdep, dir);
|
|
}
|
|
else {
|
|
int i, j;
|
|
Node *chd[4][8];
|
|
for (j = 0; j < 4; j++) {
|
|
for (i = 0; i < 8; i++) {
|
|
chd[j][i] = ((!leaf[j]) && node[j]->internal.has_child(i)) ?
|
|
node[j]->internal.get_child(
|
|
node[j]->internal.get_child_count(i)) : NULL;
|
|
}
|
|
}
|
|
|
|
// 2 edge calls
|
|
Node *ne[4];
|
|
int le[4];
|
|
int de[4];
|
|
for (i = 0; i < 2; i++) {
|
|
int c[4] = {edgeProcEdgeMask[dir][i][0],
|
|
edgeProcEdgeMask[dir][i][1],
|
|
edgeProcEdgeMask[dir][i][2],
|
|
edgeProcEdgeMask[dir][i][3]};
|
|
|
|
for (int j = 0; j < 4; j++) {
|
|
if (leaf[j]) {
|
|
le[j] = leaf[j];
|
|
ne[j] = node[j];
|
|
de[j] = depth[j];
|
|
}
|
|
else {
|
|
le[j] = node[j]->internal.is_child_leaf(c[j]);
|
|
ne[j] = chd[j][c[j]];
|
|
de[j] = depth[j] - 1;
|
|
}
|
|
}
|
|
|
|
edgeProcContour(ne, le, de, maxdep - 1, edgeProcEdgeMask[dir][i][4]);
|
|
}
|
|
|
|
}
|
|
}
|
|
|
|
void Octree::faceProcContour(Node *node[2], int leaf[2], int depth[2], int maxdep, int dir)
|
|
{
|
|
if (!(node[0] && node[1])) {
|
|
return;
|
|
}
|
|
|
|
if (!(leaf[0] && leaf[1])) {
|
|
int i, j;
|
|
// Fill children nodes
|
|
Node *chd[2][8];
|
|
for (j = 0; j < 2; j++) {
|
|
for (i = 0; i < 8; i++) {
|
|
chd[j][i] = ((!leaf[j]) && node[j]->internal.has_child(i)) ?
|
|
node[j]->internal.get_child(
|
|
node[j]->internal.get_child_count(i)) : NULL;
|
|
}
|
|
}
|
|
|
|
// 4 face calls
|
|
Node *nf[2];
|
|
int df[2];
|
|
int lf[2];
|
|
for (i = 0; i < 4; i++) {
|
|
int c[2] = {faceProcFaceMask[dir][i][0], faceProcFaceMask[dir][i][1]};
|
|
for (int j = 0; j < 2; j++) {
|
|
if (leaf[j]) {
|
|
lf[j] = leaf[j];
|
|
nf[j] = node[j];
|
|
df[j] = depth[j];
|
|
}
|
|
else {
|
|
lf[j] = node[j]->internal.is_child_leaf(c[j]);
|
|
nf[j] = chd[j][c[j]];
|
|
df[j] = depth[j] - 1;
|
|
}
|
|
}
|
|
faceProcContour(nf, lf, df, maxdep - 1, faceProcFaceMask[dir][i][2]);
|
|
}
|
|
|
|
// 4 edge calls
|
|
int orders[2][4] = {{0, 0, 1, 1}, {0, 1, 0, 1}};
|
|
Node *ne[4];
|
|
int le[4];
|
|
int de[4];
|
|
|
|
for (i = 0; i < 4; i++) {
|
|
int c[4] = {faceProcEdgeMask[dir][i][1], faceProcEdgeMask[dir][i][2],
|
|
faceProcEdgeMask[dir][i][3], faceProcEdgeMask[dir][i][4]};
|
|
int *order = orders[faceProcEdgeMask[dir][i][0]];
|
|
|
|
for (int j = 0; j < 4; j++) {
|
|
if (leaf[order[j]]) {
|
|
le[j] = leaf[order[j]];
|
|
ne[j] = node[order[j]];
|
|
de[j] = depth[order[j]];
|
|
}
|
|
else {
|
|
le[j] = node[order[j]]->internal.is_child_leaf(c[j]);
|
|
ne[j] = chd[order[j]][c[j]];
|
|
de[j] = depth[order[j]] - 1;
|
|
}
|
|
}
|
|
|
|
edgeProcContour(ne, le, de, maxdep - 1, faceProcEdgeMask[dir][i][5]);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Octree::cellProcContour(Node *node, int leaf, int depth)
|
|
{
|
|
if (node == NULL) {
|
|
return;
|
|
}
|
|
|
|
if (!leaf) {
|
|
int i;
|
|
|
|
// Fill children nodes
|
|
Node *chd[8];
|
|
for (i = 0; i < 8; i++) {
|
|
chd[i] = ((!leaf) && node->internal.has_child(i)) ?
|
|
node->internal.get_child(node->internal.get_child_count(i)) : NULL;
|
|
}
|
|
|
|
// 8 Cell calls
|
|
for (i = 0; i < 8; i++) {
|
|
cellProcContour(chd[i], node->internal.is_child_leaf(i), depth - 1);
|
|
}
|
|
|
|
// 12 face calls
|
|
Node *nf[2];
|
|
int lf[2];
|
|
int df[2] = {depth - 1, depth - 1};
|
|
for (i = 0; i < 12; i++) {
|
|
int c[2] = {cellProcFaceMask[i][0], cellProcFaceMask[i][1]};
|
|
|
|
lf[0] = 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]];
|
|
|
|
faceProcContour(nf, lf, df, depth - 1, cellProcFaceMask[i][2]);
|
|
}
|
|
|
|
// 6 edge calls
|
|
Node *ne[4];
|
|
int le[4];
|
|
int de[4] = {depth - 1, depth - 1, depth - 1, depth - 1};
|
|
for (i = 0; i < 6; i++) {
|
|
int c[4] = {cellProcEdgeMask[i][0], cellProcEdgeMask[i][1], cellProcEdgeMask[i][2], cellProcEdgeMask[i][3]};
|
|
|
|
for (int j = 0; j < 4; j++) {
|
|
le[j] = node->internal.is_child_leaf(c[j]);
|
|
ne[j] = chd[c[j]];
|
|
}
|
|
|
|
edgeProcContour(ne, le, de, depth - 1, cellProcEdgeMask[i][4]);
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
void Octree::processEdgeParity(LeafNode *node[4], int /*depth*/[4], int /*maxdep*/, int dir)
|
|
{
|
|
int con = 0;
|
|
for (int i = 0; i < 4; i++) {
|
|
// Minimal cell
|
|
// if(op == 0)
|
|
{
|
|
if (getEdgeParity(node[i], processEdgeMask[dir][i])) {
|
|
con = 1;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (con == 1) {
|
|
for (int i = 0; i < 4; i++) {
|
|
setEdge(node[i], processEdgeMask[dir][i]);
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
void Octree::edgeProcParity(Node *node[4], int leaf[4], int depth[4], int maxdep, int dir)
|
|
{
|
|
if (!(node[0] && node[1] && node[2] && node[3])) {
|
|
return;
|
|
}
|
|
if (leaf[0] && leaf[1] && leaf[2] && leaf[3]) {
|
|
processEdgeParity((LeafNode **)node, depth, maxdep, dir);
|
|
}
|
|
else {
|
|
int i, j;
|
|
Node *chd[4][8];
|
|
for (j = 0; j < 4; j++) {
|
|
for (i = 0; i < 8; i++) {
|
|
chd[j][i] = ((!leaf[j]) && node[j]->internal.has_child(i)) ?
|
|
node[j]->internal.get_child( node[j]->internal.get_child_count(i)) : NULL;
|
|
}
|
|
}
|
|
|
|
// 2 edge calls
|
|
Node *ne[4];
|
|
int le[4];
|
|
int de[4];
|
|
for (i = 0; i < 2; i++) {
|
|
int c[4] = {edgeProcEdgeMask[dir][i][0],
|
|
edgeProcEdgeMask[dir][i][1],
|
|
edgeProcEdgeMask[dir][i][2],
|
|
edgeProcEdgeMask[dir][i][3]};
|
|
|
|
// int allleaf = 1;
|
|
for (int j = 0; j < 4; j++) {
|
|
|
|
if (leaf[j]) {
|
|
le[j] = leaf[j];
|
|
ne[j] = node[j];
|
|
de[j] = depth[j];
|
|
}
|
|
else {
|
|
le[j] = node[j]->internal.is_child_leaf(c[j]);
|
|
ne[j] = chd[j][c[j]];
|
|
de[j] = depth[j] - 1;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
edgeProcParity(ne, le, de, maxdep - 1, edgeProcEdgeMask[dir][i][4]);
|
|
}
|
|
|
|
}
|
|
}
|
|
|
|
void Octree::faceProcParity(Node *node[2], int leaf[2], int depth[2], int maxdep, int dir)
|
|
{
|
|
if (!(node[0] && node[1])) {
|
|
return;
|
|
}
|
|
|
|
if (!(leaf[0] && leaf[1])) {
|
|
int i, j;
|
|
// Fill children nodes
|
|
Node *chd[2][8];
|
|
for (j = 0; j < 2; j++) {
|
|
for (i = 0; i < 8; i++) {
|
|
chd[j][i] = ((!leaf[j]) && node[j]->internal.has_child(i)) ?
|
|
node[j]->internal.get_child(
|
|
node[j]->internal.get_child_count(i)) : NULL;
|
|
}
|
|
}
|
|
|
|
// 4 face calls
|
|
Node *nf[2];
|
|
int df[2];
|
|
int lf[2];
|
|
for (i = 0; i < 4; i++) {
|
|
int c[2] = {faceProcFaceMask[dir][i][0], faceProcFaceMask[dir][i][1]};
|
|
for (int j = 0; j < 2; j++) {
|
|
if (leaf[j]) {
|
|
lf[j] = leaf[j];
|
|
nf[j] = node[j];
|
|
df[j] = depth[j];
|
|
}
|
|
else {
|
|
lf[j] = node[j]->internal.is_child_leaf(c[j]);
|
|
nf[j] = chd[j][c[j]];
|
|
df[j] = depth[j] - 1;
|
|
}
|
|
}
|
|
faceProcParity(nf, lf, df, maxdep - 1, faceProcFaceMask[dir][i][2]);
|
|
}
|
|
|
|
// 4 edge calls
|
|
int orders[2][4] = {{0, 0, 1, 1}, {0, 1, 0, 1}};
|
|
Node *ne[4];
|
|
int le[4];
|
|
int de[4];
|
|
|
|
for (i = 0; i < 4; i++) {
|
|
int c[4] = {faceProcEdgeMask[dir][i][1], faceProcEdgeMask[dir][i][2],
|
|
faceProcEdgeMask[dir][i][3], faceProcEdgeMask[dir][i][4]};
|
|
int *order = orders[faceProcEdgeMask[dir][i][0]];
|
|
|
|
for (int j = 0; j < 4; j++) {
|
|
if (leaf[order[j]]) {
|
|
le[j] = leaf[order[j]];
|
|
ne[j] = node[order[j]];
|
|
de[j] = depth[order[j]];
|
|
}
|
|
else {
|
|
le[j] = node[order[j]]->internal.is_child_leaf(c[j]);
|
|
ne[j] = chd[order[j]][c[j]];
|
|
de[j] = depth[order[j]] - 1;
|
|
}
|
|
}
|
|
|
|
edgeProcParity(ne, le, de, maxdep - 1, faceProcEdgeMask[dir][i][5]);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Octree::cellProcParity(Node *node, int leaf, int depth)
|
|
{
|
|
if (node == NULL) {
|
|
return;
|
|
}
|
|
|
|
if (!leaf) {
|
|
int i;
|
|
|
|
// Fill children nodes
|
|
Node *chd[8];
|
|
for (i = 0; i < 8; i++) {
|
|
chd[i] = ((!leaf) && node->internal.has_child(i)) ?
|
|
node->internal.get_child(node->internal.get_child_count(i)) : NULL;
|
|
}
|
|
|
|
// 8 Cell calls
|
|
for (i = 0; i < 8; i++) {
|
|
cellProcParity(chd[i], node->internal.is_child_leaf(i), depth - 1);
|
|
}
|
|
|
|
// 12 face calls
|
|
Node *nf[2];
|
|
int lf[2];
|
|
int df[2] = {depth - 1, depth - 1};
|
|
for (i = 0; i < 12; i++) {
|
|
int c[2] = {cellProcFaceMask[i][0], cellProcFaceMask[i][1]};
|
|
|
|
lf[0] = 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]];
|
|
|
|
faceProcParity(nf, lf, df, depth - 1, cellProcFaceMask[i][2]);
|
|
}
|
|
|
|
// 6 edge calls
|
|
Node *ne[4];
|
|
int le[4];
|
|
int de[4] = {depth - 1, depth - 1, depth - 1, depth - 1};
|
|
for (i = 0; i < 6; i++) {
|
|
int c[4] = {cellProcEdgeMask[i][0], cellProcEdgeMask[i][1], cellProcEdgeMask[i][2], cellProcEdgeMask[i][3]};
|
|
|
|
for (int j = 0; j < 4; j++) {
|
|
le[j] = node->internal.is_child_leaf(c[j]);
|
|
ne[j] = chd[c[j]];
|
|
}
|
|
|
|
edgeProcParity(ne, le, de, depth - 1, cellProcEdgeMask[i][4]);
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
/* definitions for global arrays */
|
|
const int edgemask[3] = {5, 3, 6};
|
|
|
|
const int faceMap[6][4] = {
|
|
{4, 8, 5, 9},
|
|
{6, 10, 7, 11},
|
|
{0, 8, 1, 10},
|
|
{2, 9, 3, 11},
|
|
{0, 4, 2, 6},
|
|
{1, 5, 3, 7}
|
|
};
|
|
|
|
const int cellProcFaceMask[12][3] = {
|
|
{0, 4, 0},
|
|
{1, 5, 0},
|
|
{2, 6, 0},
|
|
{3, 7, 0},
|
|
{0, 2, 1},
|
|
{4, 6, 1},
|
|
{1, 3, 1},
|
|
{5, 7, 1},
|
|
{0, 1, 2},
|
|
{2, 3, 2},
|
|
{4, 5, 2},
|
|
{6, 7, 2}
|
|
};
|
|
|
|
const int cellProcEdgeMask[6][5] = {
|
|
{0, 1, 2, 3, 0},
|
|
{4, 5, 6, 7, 0},
|
|
{0, 4, 1, 5, 1},
|
|
{2, 6, 3, 7, 1},
|
|
{0, 2, 4, 6, 2},
|
|
{1, 3, 5, 7, 2}
|
|
};
|
|
|
|
const int faceProcFaceMask[3][4][3] = {
|
|
{{4, 0, 0},
|
|
{5, 1, 0},
|
|
{6, 2, 0},
|
|
{7, 3, 0}},
|
|
{{2, 0, 1},
|
|
{6, 4, 1},
|
|
{3, 1, 1},
|
|
{7, 5, 1}},
|
|
{{1, 0, 2},
|
|
{3, 2, 2},
|
|
{5, 4, 2},
|
|
{7, 6, 2}}
|
|
};
|
|
|
|
const int faceProcEdgeMask[3][4][6] = {
|
|
{{1, 4, 0, 5, 1, 1},
|
|
{1, 6, 2, 7, 3, 1},
|
|
{0, 4, 6, 0, 2, 2},
|
|
{0, 5, 7, 1, 3, 2}},
|
|
{{0, 2, 3, 0, 1, 0},
|
|
{0, 6, 7, 4, 5, 0},
|
|
{1, 2, 0, 6, 4, 2},
|
|
{1, 3, 1, 7, 5, 2}},
|
|
{{1, 1, 0, 3, 2, 0},
|
|
{1, 5, 4, 7, 6, 0},
|
|
{0, 1, 5, 0, 4, 1},
|
|
{0, 3, 7, 2, 6, 1}}
|
|
};
|
|
|
|
const int edgeProcEdgeMask[3][2][5] = {
|
|
{{3, 2, 1, 0, 0},
|
|
{7, 6, 5, 4, 0}},
|
|
{{5, 1, 4, 0, 1},
|
|
{7, 3, 6, 2, 1}},
|
|
{{6, 4, 2, 0, 2},
|
|
{7, 5, 3, 1, 2}},
|
|
};
|
|
|
|
const int processEdgeMask[3][4] = {
|
|
{3, 2, 1, 0},
|
|
{7, 5, 6, 4},
|
|
{11, 10, 9, 8}
|
|
};
|
|
|
|
const int dirCell[3][4][3] = {
|
|
{{0, -1, -1},
|
|
{0, -1, 0},
|
|
{0, 0, -1},
|
|
{0, 0, 0}},
|
|
{{-1, 0, -1},
|
|
{-1, 0, 0},
|
|
{0, 0, -1},
|
|
{0, 0, 0}},
|
|
{{-1, -1, 0},
|
|
{-1, 0, 0},
|
|
{0, -1, 0},
|
|
{0, 0, 0}}
|
|
};
|
|
|
|
const int dirEdge[3][4] = {
|
|
{3, 2, 1, 0},
|
|
{7, 6, 5, 4},
|
|
{11, 10, 9, 8}
|
|
};
|
|
|
|
int InternalNode::numChildrenTable[256];
|
|
int InternalNode::childrenCountTable[256][8];
|
|
int InternalNode::childrenIndexTable[256][8];
|