-
sparker
thiago
Manta Interactive Ray Tracer Development Mailing List
Text archives Help
- From: knolla@sci.utah.edu
- To: manta@sci.utah.edu
- Subject: [MANTA] r1088 - in trunk: Model/Primitives scenes
- Date: Sat, 27 May 2006 15:09:59 -0600 (MDT)
Author: knolla
Date: Sat May 27 15:09:52 2006
New Revision: 1088
Modified:
trunk/Model/Primitives/IsosurfaceOctreeVolume.cc
trunk/Model/Primitives/OctreeVolume.cc
trunk/Model/Primitives/OctreeVolume.h
trunk/scenes/octisovol.cc
Log:
Major change to octree volume format. New format is 20% faster despite 1-2%
larger memory footprint. IsosurfaceOctreeVolume now allows multiresolution
traversals, which slow down the application considerably. This feature can be
enabled or disabled via OCTVOL_DYNAMIC_MULTIRES in OctreeVolume.
Modified: trunk/Model/Primitives/IsosurfaceOctreeVolume.cc
==============================================================================
--- trunk/Model/Primitives/IsosurfaceOctreeVolume.cc (original)
+++ trunk/Model/Primitives/IsosurfaceOctreeVolume.cc Sat May 27 15:09:52
2006
@@ -11,10 +11,73 @@
#define MIN4(a,b,c,d) min(min(a,b), min(c,d));
#define MAX4(a,b,c,d) max(max(a,b), max(c,d));
-
-using namespace Manta;
+
+//requires local variables:
+// NODE (either node or cap... anything that has .values[])
+// offset = Vec3i(0,0,1)
+// every other variable must be either locally or globall declared. See
caller examples below.
+#define octvol_fill_cell(NODE, cw) \
+ if (target_child & 1) \
+ this_rho = octdata->lookup_neighbor<0,0,1>(child_cell,
offset, stop_depth, prev_depth, index_trace); \
+ else \
+ this_rho = NODE.values[target_child | 1]; \
+ min_rho = MIN(min_rho, this_rho); \
+ max_rho = MAX(max_rho, this_rho); \
+ rho[0][0][1] = static_cast<float>(this_rho); \
+ offset.data[1] = cw; \
+ if (target_child & 3) \
+ this_rho = octdata->lookup_neighbor<0,1,1>(child_cell,
offset, stop_depth, prev_depth, index_trace); \
+ else \
+ this_rho = NODE.values[target_child | 3]; \
+ min_rho = MIN(min_rho, this_rho); \
+ max_rho = MAX(max_rho, this_rho); \
+ rho[0][1][1] = static_cast<float>(this_rho); \
+ offset.data[0] = cw; \
+ if (target_child & 7) \
+ this_rho = octdata->lookup_neighbor<1,1,1>(child_cell,
offset, stop_depth, prev_depth, index_trace); \
+ else \
+ this_rho = NODE.values[target_child | 7]; \
+ min_rho = MIN(min_rho, this_rho); \
+ max_rho = MAX(max_rho, this_rho); \
+ rho[1][1][1] = static_cast<float>(this_rho); \
+ offset.data[2] = 0; \
+ if (target_child & 6) \
+ this_rho = octdata->lookup_neighbor<1,1,0>(child_cell,
offset, stop_depth, prev_depth, index_trace); \
+ else \
+ this_rho = NODE.values[target_child | 6]; \
+ min_rho = MIN(min_rho, this_rho); \
+ max_rho = MAX(max_rho, this_rho); \
+ rho[1][1][0] = static_cast<float>(this_rho); \
+ offset.data[1] = 0; \
+ if (target_child & 4) \
+ this_rho = octdata->lookup_neighbor<1,0,0>(child_cell,
offset, stop_depth, prev_depth, index_trace); \
+ else \
+ this_rho = NODE.values[target_child | 4]; \
+ min_rho = MIN(min_rho, this_rho); \
+ max_rho = MAX(max_rho, this_rho); \
+ rho[1][0][0] = static_cast<float>(this_rho); \
+ offset.data[2] = cw; \
+ if (target_child & 5) \
+ this_rho = octdata->lookup_neighbor<1,0,1>(child_cell,
offset, stop_depth, prev_depth, index_trace); \
+ else \
+ this_rho = NODE.values[target_child | 5]; \
+ min_rho = MIN(min_rho, this_rho); \
+ max_rho = MAX(max_rho, this_rho); \
+ rho[1][0][1] = static_cast<float>(this_rho); \
+ offset.data[0] = 0; \
+ offset.data[1] = cw; \
+ offset.data[2] = 0; \
+ if (target_child & 2) \
+ this_rho = octdata->lookup_neighbor<0,1,0>(child_cell,
offset, stop_depth, prev_depth, index_trace); \
+ else \
+ this_rho = NODE.values[target_child | 2]; \
+ min_rho = MIN(min_rho, this_rho); \
+ max_rho = MAX(max_rho, this_rho); \
+ rho[0][1][0] = static_cast<float>(this_rho); \
static const int axis_table[] = {4, 2, 1};
+
+using namespace Manta;
IsosurfaceOctreeVolume::IsosurfaceOctreeVolume(OctreeVolume* _octdata,
Material* _matl)
: PrimitiveCommon(_matl), octdata(_octdata)
@@ -113,30 +176,23 @@
//tenter_padded = MAX(0.f, tenter_padded);
unsigned int index_trace[octdata->get_max_depth()+1];
+
+ int stop_depth = octdata->get_cap_depth() - 0;
- Vec3i cell(0,0,0);
-
- #if 0
- if (rays.hit(which_one, tenter_padded, PrimitiveCommon::getMaterial(),
this, 0))
- {
- rays.setNormal(which_one, Vector(-1,0,0));
- return;
- }
- #else
- single_traverse_node(rays, which_one, orig, dir, inv_dir, 0,
+ Vec3i cell(0,0,0);
+ single_traverse_node(rays, which_one, orig, dir, inv_dir, stop_depth,
0, 0, index_trace, cell, tenter_padded,
texit_padded);
- #endif
}
bool IsosurfaceOctreeVolume::single_traverse_node(RayPacket& rays, int
which_one,
const Vector& orig, const Vector& dir, const Vector&
inv_dir,
- int res, int depth, unsigned int node_index,
+ int stop_depth, int depth, unsigned int node_index,
unsigned int index_trace[], Vec3i& cell, const float
tenter,
const float texit) const
{
+ //cerr << "single_traverse_node, depth=" << depth << ", node_index=" <<
node_index << ", cell=" << cell.data[0] << "," << cell.data[1] << "," <<
cell.data[2] << endl;
OctNode& node = octdata->get_node(depth, node_index);
-
- //cerr << "single_traverse_node, depth=" << depth << ", node_index=" <<
node_index << ", cell=" << cell[0] << "," << cell[1] << "," << cell[2] <<
endl;
+
index_trace[depth] = node_index;
int child_bit = octdata->get_child_bit_depth(depth);
@@ -187,8 +243,6 @@
axis_isects.data[2] = 0;
}
}
-
- const int axis_table[] = {4, 2, 1};
float child_tenter = tenter;
float child_texit;
@@ -212,28 +266,63 @@
}
}
- if (octdata->get_isovalue() >= node.child_mins[target_child] &&
octdata->get_isovalue() <= node.child_maxs[target_child])
+ if (octdata->get_isovalue() >= node.mins[target_child] &&
octdata->get_isovalue() <= node.maxs[target_child])
{
- if (node.scalar_leaf & (1 << target_child))
+
+#ifdef OCTVOL_DYNAMIC_MULTIRES
+ if (depth == stop_depth)
+ {
+ Vector cmin(child_cell.data[0],
child_cell.data[1], child_cell.data[2]);
+ Vector cmax(child_cell.data[0] + child_bit,
child_cell.data[1] + child_bit, child_cell.data[2] + child_bit);
+
+ float rho[2][2][2];
+ ST min_rho, max_rho, this_rho;
+ min_rho = max_rho = this_rho =
node.values[target_child];
+ rho[0][0][0] = static_cast<float>(this_rho);
+ int prev_depth = depth-1;
+ Vec3i offset(0,0,child_bit);
+ octvol_fill_cell(node, child_bit);
+
+ if (octdata->get_isovalue() >= min_rho &&
octdata->get_isovalue() <= max_rho)
+ {
+ float hit_t;
+ if
(IsosurfaceImplicit::single_intersect(orig, dir, cmin, cmax, rho,
+ octdata->get_isovalue(),
child_tenter, child_texit, hit_t))
+ {
+ if (rays.hit(which_one,
hit_t, PrimitiveCommon::getMaterial(), this, 0))
+ {
+ Vector normal;
+ Vector phit = orig +
dir*hit_t;
+
IsosurfaceImplicit::single_normal(normal, cmin, cmax, phit, rho);
+ normal.normalize();
+
rays.setNormal(which_one, normal);
+ return true;
+ }
+ }
+ }
+ }
+ else
+#endif
+ if (node.offsets[target_child]==-1)
{
return false;
- if (single_traverse_leaf(rays, which_one, orig, dir,
inv_dir, res,
- next_depth, depth, node.children[target_child].scalar,
+ if (single_traverse_leaf(rays, which_one, orig, dir,
inv_dir, stop_depth,
+ next_depth, depth, node.values[target_child],
child_cell, index_trace, child_cell, child_tenter,
child_texit))
return true;
}
else
{
- unsigned int child_idx = node.children_start +
node.children[target_child].offset;
- if (next_depth+1 == octdata->get_max_depth()) //cap
+ unsigned int child_idx = node.children_start +
node.offsets[target_child];
+ if (depth == octdata->get_pre_cap_depth()) //cap
{
- if (single_traverse_cap(rays, which_one, orig, dir,
inv_dir, res, next_depth, child_idx,
+ if (single_traverse_cap(rays, which_one, orig, dir,
inv_dir, stop_depth, next_depth, child_idx,
index_trace, child_cell, child_tenter, child_texit))
return true;
}
else
{
- if (single_traverse_node(rays, which_one, orig, dir,
inv_dir, res, next_depth, child_idx,
+ if (single_traverse_node(rays, which_one, orig, dir,
inv_dir, stop_depth, next_depth, child_idx,
index_trace, child_cell, child_tenter, child_texit))
return true;
}
@@ -261,14 +350,14 @@
}
bool IsosurfaceOctreeVolume::single_traverse_leaf(RayPacket& rays, int
which_one,
- const Vector& orig, const Vector& dir, const Vector&
inv_dir, int res,
+ const Vector& orig, const Vector& dir, const Vector&
inv_dir, int stop_depth,
int depth, int leaf_depth, ST scalar, Vec3i&
leaf_base_cell,
unsigned int index_trace[], Vec3i& cell, const float
tenter, const float texit) const
{
//using pretty much the same algorithm, find which "implicit" cell
(-->voxel) we're in at octdata->get_max_depth()
- int child_bit = octdata->get_child_bit_depth(depth+res);
- int unsafe_zone = octdata->get_child_bit_depth(depth-1+res) -
octdata->get_child_bit_depth(octdata->get_max_depth()-1);
+ int child_bit = octdata->get_child_bit_depth(depth);
+ int unsafe_zone = octdata->get_child_bit_depth(depth-1) -
octdata->get_child_bit_depth(octdata->get_cap_depth());
Vector center(static_cast<float>(cell.data[0] | child_bit),
static_cast<float>(cell.data[1] | child_bit), static_cast<float>(cell.data[2]
| child_bit));
Vector tcenter = inv_dir * (center - orig);
@@ -316,8 +405,6 @@
}
}
- const int axis_table[] = {4, 2, 1};
-
float child_tenter = tenter;
float child_texit;
@@ -343,24 +430,24 @@
Vec3i local_child_cell = child_cell - leaf_base_cell;
if (local_child_cell.data[0] & unsafe_zone ||
local_child_cell.data[1] & unsafe_zone || local_child_cell.data[2] &
unsafe_zone)
{
- if (depth == octdata->get_max_depth()-1-res)
+ if (depth == stop_depth)
{
//try isosurface intersection
Vector cmin(child_cell.data[0], child_cell.data[1],
child_cell.data[2]);
- Vector cmax = cmin + Vector(1.f,1.f,1.f);
-
+ Vector cmax(child_cell.data[0] + child_bit,
child_cell.data[1] + child_bit, child_cell.data[2] + child_bit);
+
#ifdef USE_OCTREE_DATA
//use octree data
float rho[2][2][2];
ST min_rho, max_rho, this_rho;
min_rho = max_rho = this_rho = scalar;
rho[0][0][0] = static_cast<float>(this_rho);
-
+ Vec3i offset(0,0,child_bit);
+
//0,0,1
- Vec3i offset(0,0,1);
if (target_child & 1)
{
- this_rho = octdata->lookup_neighbor<0,0,1>(child_cell,
offset, res, leaf_depth, index_trace);
+ this_rho = octdata->lookup_neighbor<0,0,1>(child_cell,
offset, stop_depth, leaf_depth, index_trace);
min_rho = MIN(min_rho, this_rho);
max_rho = MAX(max_rho, this_rho);
}
@@ -369,10 +456,10 @@
rho[0][0][1] = static_cast<float>(this_rho);
//0,1,1
- offset.data[1] = 1;
+ offset.data[1] = child_bit;
if (target_child & 3)
{
- this_rho = octdata->lookup_neighbor<0,1,1>(child_cell,
offset, res, leaf_depth, index_trace);
+ this_rho = octdata->lookup_neighbor<0,1,1>(child_cell,
offset, stop_depth, leaf_depth, index_trace);
min_rho = MIN(min_rho, this_rho);
max_rho = MAX(max_rho, this_rho);
}
@@ -381,10 +468,10 @@
rho[0][1][1] = static_cast<float>(this_rho);
//1,1,1
- offset.data[0] = 1;
+ offset.data[0] = child_bit;
if (target_child & 7)
{
- this_rho = octdata->lookup_neighbor<1,1,1>(child_cell,
offset, res, leaf_depth, index_trace);
+ this_rho = octdata->lookup_neighbor<1,1,1>(child_cell,
offset, stop_depth, leaf_depth, index_trace);
min_rho = MIN(min_rho, this_rho);
max_rho = MAX(max_rho, this_rho);
}
@@ -396,7 +483,7 @@
offset.data[2] = 0;
if (target_child & 6)
{
- this_rho = octdata->lookup_neighbor<1,1,0>(child_cell,
offset, res, leaf_depth, index_trace);
+ this_rho = octdata->lookup_neighbor<1,1,0>(child_cell,
offset, stop_depth, leaf_depth, index_trace);
min_rho = MIN(min_rho, this_rho);
max_rho = MAX(max_rho, this_rho);
}
@@ -408,7 +495,7 @@
offset.data[1] = 0;
if (target_child & 4)
{
- this_rho = octdata->lookup_neighbor<1,0,0>(child_cell,
offset, res, leaf_depth, index_trace);
+ this_rho = octdata->lookup_neighbor<1,0,0>(child_cell,
offset, stop_depth, leaf_depth, index_trace);
min_rho = MIN(min_rho, this_rho);
max_rho = MAX(max_rho, this_rho);
}
@@ -417,10 +504,10 @@
rho[1][0][0] = static_cast<float>(this_rho);
//1,0,1
- offset.data[2] = 1;
+ offset.data[2] = child_bit;
if (target_child & 5)
{
- this_rho = octdata->lookup_neighbor<1,0,1>(child_cell,
offset, res, leaf_depth, index_trace);
+ this_rho = octdata->lookup_neighbor<1,0,1>(child_cell,
offset, stop_depth, leaf_depth, index_trace);
min_rho = MIN(min_rho, this_rho);
max_rho = MAX(max_rho, this_rho);
}
@@ -430,11 +517,11 @@
//0,1,0
offset.data[0] = 0;
- offset.data[1] = 1;
+ offset.data[1] = child_bit;
offset.data[2] = 0;
if (target_child & 2)
{
- this_rho = octdata->lookup_neighbor<0,1,0>(child_cell,
offset, res, leaf_depth, index_trace);
+ this_rho = octdata->lookup_neighbor<0,1,0>(child_cell,
offset, stop_depth, leaf_depth, index_trace);
min_rho = MIN(min_rho, this_rho);
max_rho = MAX(max_rho, this_rho);
}
@@ -460,13 +547,6 @@
#endif
if (octdata->get_isovalue() >= min_rho &&
octdata->get_isovalue() <= max_rho)
{
- #if 1
- if (rays.hit(which_one, child_tenter,
PrimitiveCommon::getMaterial(), this, 0))
- {
- rays.setNormal(which_one, Vector(-1,0,0));
- return true;
- }
- #else
float hit_t;
if (IsosurfaceImplicit::single_intersect(orig, dir,
cmin, cmax, rho,
octdata->get_isovalue(), child_tenter, child_texit,
hit_t))
@@ -481,12 +561,11 @@
return true;
}
}
- #endif
}
}
else //not at cap-level depth
{
- if (single_traverse_leaf(rays, which_one, orig, dir,
inv_dir, res,
+ if (single_traverse_leaf(rays, which_one, orig, dir,
inv_dir, stop_depth,
next_depth, leaf_depth, scalar, leaf_base_cell,
index_trace, child_cell, child_tenter, child_texit))
return true;
@@ -514,7 +593,7 @@
}
bool IsosurfaceOctreeVolume::single_traverse_cap(RayPacket& rays, int
which_one,
- const Vector& orig, const Vector& dir, const Vector&
inv_dir, int res,
+ const Vector& orig, const Vector& dir, const Vector&
inv_dir, int stop_depth,
int depth, unsigned int cap_index, unsigned int
index_trace[], Vec3i& cell,
const float tenter, const float texit) const
{
@@ -569,8 +648,6 @@
}
}
- const int axis_table[] = {4, 2, 1};
-
float child_tenter = tenter;
float child_texit;
@@ -593,87 +670,16 @@
//try isosurface intersection in this node
Vector cmin(child_cell.data[0], child_cell.data[1],
child_cell.data[2]);
- Vector cmax = cmin + Vector(1.f,1.f,1.f);
+ Vector cmax(child_cell.data[0] + 1, child_cell.data[1] + 1,
child_cell.data[2] + 1);
#ifdef USE_OCTREE_DATA
- //use octree data
float rho[2][2][2];
ST min_rho, max_rho, this_rho;
min_rho = max_rho = this_rho = cap.values[target_child];
rho[0][0][0] = static_cast<float>(this_rho);
int prev_depth = depth-1;
-
- //0,0,1
Vec3i offset(0,0,1);
- if (target_child & 1)
- this_rho = octdata->lookup_neighbor<0,0,1>(child_cell, offset,
res, prev_depth, index_trace);
- else
- this_rho = cap.values[target_child | 1];
- min_rho = MIN(min_rho, this_rho);
- max_rho = MAX(max_rho, this_rho);
- rho[0][0][1] = static_cast<float>(this_rho);
-
- //0,1,1
- offset.data[1] = 1;
- if (target_child & 3)
- this_rho = octdata->lookup_neighbor<0,1,1>(child_cell, offset,
res, prev_depth, index_trace);
- else
- this_rho = cap.values[target_child | 3];
- min_rho = MIN(min_rho, this_rho);
- max_rho = MAX(max_rho, this_rho);
- rho[0][1][1] = static_cast<float>(this_rho);
-
- //1,1,1
- offset.data[0] = 1;
- if (target_child & 7)
- this_rho = octdata->lookup_neighbor<1,1,1>(child_cell, offset,
res, prev_depth, index_trace);
- else
- this_rho = cap.values[target_child | 7];
- min_rho = MIN(min_rho, this_rho);
- max_rho = MAX(max_rho, this_rho);
- rho[1][1][1] = static_cast<float>(this_rho);
-
- //1,1,0
- offset.data[2] = 0;
- if (target_child & 6)
- this_rho = octdata->lookup_neighbor<1,1,0>(child_cell, offset,
res, prev_depth, index_trace);
- else
- this_rho = cap.values[target_child | 6];
- min_rho = MIN(min_rho, this_rho);
- max_rho = MAX(max_rho, this_rho);
- rho[1][1][0] = static_cast<float>(this_rho);
-
- //1,0,0
- offset.data[1] = 0;
- if (target_child & 4)
- this_rho = octdata->lookup_neighbor<1,0,0>(child_cell, offset,
res, prev_depth, index_trace);
- else
- this_rho = cap.values[target_child | 4];
- min_rho = MIN(min_rho, this_rho);
- max_rho = MAX(max_rho, this_rho);
- rho[1][0][0] = static_cast<float>(this_rho);
-
- //1,0,1
- offset.data[2] = 1;
- if (target_child & 5)
- this_rho = octdata->lookup_neighbor<1,0,1>(child_cell, offset,
res, prev_depth, index_trace);
- else
- this_rho = cap.values[target_child | 5];
- min_rho = MIN(min_rho, this_rho);
- max_rho = MAX(max_rho, this_rho);
- rho[1][0][1] = static_cast<float>(this_rho);
-
- //0,1,0
- offset.data[0] = 0;
- offset.data[1] = 1;
- offset.data[2] = 0;
- if (target_child & 2)
- this_rho = octdata->lookup_neighbor<0,1,0>(child_cell, offset,
res, prev_depth, index_trace);
- else
- this_rho = cap.values[target_child | 2];
- min_rho = MIN(min_rho, this_rho);
- max_rho = MAX(max_rho, this_rho);
- rho[0][1][0] = static_cast<float>(this_rho);
+ octvol_fill_cell(cap, 1);
#else
//use original grid data
float rho[2][2][2];
@@ -693,13 +699,6 @@
#endif
if (octdata->get_isovalue() >= min_rho && octdata->get_isovalue() <=
max_rho)
{
- #if 0
- if (rays.hit(which_one, child_tenter,
PrimitiveCommon::getMaterial(), this, 0))
- {
- rays.setNormal(which_one, Vector(-1,0,0));
- return true;
- }
- #else
float hit_t;
if (IsosurfaceImplicit::single_intersect(orig, dir, cmin, cmax,
rho,
octdata->get_isovalue(), child_tenter, child_texit, hit_t))
@@ -714,7 +713,6 @@
return true;
}
}
- #endif
}
if (axis==-1)
Modified: trunk/Model/Primitives/OctreeVolume.cc
==============================================================================
--- trunk/Model/Primitives/OctreeVolume.cc (original)
+++ trunk/Model/Primitives/OctreeVolume.cc Sat May 27 15:09:52 2006
@@ -56,7 +56,6 @@
int endTimestep,
double variance_threshold,
int kernelWidth,
- int mresLevels,
ST isomin,
ST isomax)
{
@@ -123,6 +122,8 @@
double ddepth_3 = log((double)nz) / log(2.0);
double ddepth = MAX(MAX(ddepth_1, ddepth_2), ddepth_3);
max_depth = (int)(ceil(ddepth));
+ cap_depth = max_depth - 1;
+ pre_cap_depth = max_depth - 2;
padded_dims.data[0] = padded_dims.data[1] = padded_dims.data[2] =
(float)max_depth_bit;
@@ -173,8 +174,8 @@
cerr << "indata(150,150,150)=" << indata(150,150,150) << endl;
//build the octree data from the single-res data.
- cerr << "Building multi-resolution octree data. (var thresh =
"<< variance_threshold << ", kernel width = " << kernelWidth << ", mres
levels = " << mresLevels << ")...";
- make_from_single_res_data(indata, ts, variance_threshold,
kernelWidth, mresLevels, isomin, isomax);
+ cerr << "Building multi-resolution octree data. (var thresh =
"<< variance_threshold << ", kernel width = " << kernelWidth << ")...";
+ make_from_single_res_data(indata, ts, variance_threshold,
kernelWidth, isomin, isomax);
cerr << "\nDone making from single res data" << endl;
}
@@ -192,19 +193,12 @@
{
for(int ts=0; ts<num_timesteps; ts++)
{
- for(int r=0; r<mres_levels+1; r++)
- {
- delete[] steps[ts].caps[r];
+ delete[] steps[ts].caps;
- for(int d=0; d<max_depth-1-r; d++)
- delete[] steps[ts].nodes[r][d];
- delete[] steps[ts].nodes[r];
- delete[] steps[ts].num_nodes[r];
- }
-
- delete[] steps[ts].caps;
- delete[] steps[ts].nodes;
- delete[] steps[ts].num_nodes;
+ for(int d=0; d<max_depth-1; d++)
+ delete[] steps[ts].nodes[d];
+ delete[] steps[ts].nodes;
+ delete[] steps[ts].num_nodes;
}
delete[] steps;
@@ -215,7 +209,7 @@
void OctreeVolume::make_from_single_res_data(Array3<ST>& originalData,
int ts, double
variance_threshold,
- int kernelWidth, int
mresLevels,
+ int kernelWidth,
ST isomin, ST isomax)
{
if (isomin >= isomax)
@@ -224,21 +218,12 @@
isomax = std::numeric_limits<ST>::max();
}
- cout << "Timestep " << ts << ": Building octree data of depth " <<
max_depth << ", variance threshold " << variance_threshold << ", kernel width
" << kernel_width << ", multiresolution levels " << mresLevels << ", isomin "
<< (int)isomin << ", isomax " << (int)isomax << ".\n";
-
- mres_levels = mresLevels;
-
+ cout << "Timestep " << ts << ": Building octree data of depth " <<
max_depth << ", variance threshold " << variance_threshold << ", kernel width
" << kernel_width << ", isomin " << (int)isomin << ", isomax " << (int)isomax
<< ".\n";
+
Array3<BuildNode>* buildnodes = new Array3<BuildNode>[max_depth];
- steps[ts].num_nodes = new unsigned int*[mres_levels+1];
- steps[ts].nodes = new OctNode**[mres_levels+1];
- for(int r=0; r<mres_levels+1; r++)
- {
- steps[ts].num_nodes[r] = new unsigned int[max_depth-1];
- steps[ts].nodes[r] = new OctNode*[max_depth-1];
- }
- steps[ts].num_caps = new unsigned int[mres_levels+1];
- steps[ts].caps = new OctCap*[mres_levels+1];
+ steps[ts].num_nodes = new unsigned int[max_depth-1];
+ steps[ts].nodes = new OctNode*[max_depth-1];
//1. start at the bottom. Create 2^(depth-1) nodes in each dim
//NOTE: d=0 is root node (coarsest depth)
@@ -253,326 +238,316 @@
unsigned int totalnodesalldepths = 0;
unsigned int totalmaxnodesalldepths = 0;
- for(int r=0; r<mres_levels+1; r++)
- {
- //build the buildnodes, from bottom up.
- for(int d=max_depth-1-r; d>=0; d--)
- {
- Array3<BuildNode>& bnodes = buildnodes[d];
- unsigned int ncaps = 0;
- unsigned int nnodes = 0;
-
- int nodes_dim = 1 << d;
- //cout << "\nd="<<d<<":";
- //cout << "nodes_dim = " << nodes_dim << endl;
-
- int dim_extents = 1 << (max_depth - d - r - 1);
- //cout << "dim_extents = " << dim_extents << endl;
-
- for(int i=0; i<nodes_dim; i++)
- for(int j=0; j<nodes_dim; j++)
- for(int k=0; k<nodes_dim; k++)
- {
- BuildNode& bnode = bnodes(i,j,k);
-
- if (r==0)
- {
- if (d==max_depth-1)
- {
- for(int pi=0; pi<2; pi++)
- for(int pj=0; pj<2; pj++)
- for(int pk=0; pk<2; pk++)
- {
- Vec3i dp((2*i + pi) *
dim_extents,
- (2*j + pj) *
dim_extents,
- (2*k + pk) *
dim_extents);
- ST odval =
lookup_safe(originalData, dp.data[0], dp.data[1], dp.data[2]);
- narrow_range(odval, isomin,
isomax);
- bnode.values[4*pi | 2*pj | pk] =
odval;
- }
- }
- else
- {
- for(int pi=0; pi<2; pi++)
- for(int pj=0; pj<2; pj++)
- for(int pk=0; pk<2; pk++)
- {
- bnode.values[4*pi | 2*pj |
pk] = buildnodes[d+1](2*i+pi, 2*j+pj, 2*k+pk).expected_value;
- }
- }
- }
-
- if (d == max_depth-1-r) //finest-level data
are automatically caps, *if* they are kept
- {
-
- double expectedValue = 0.0;
- for(int v=0; v<8; v++)
- expectedValue +=
static_cast<double>(bnode.values[v]);
- expectedValue *= 0.125;
-
- double variances[8];
- bnode.keepMe = false;
- bnode.isLeaf = true;
- for(int v=0; v<8; v++)
- {
- double diff =
static_cast<double>(bnode.values[v]) - expectedValue;
- variances[v] = diff * diff * 0.125;
-
- if (variances[v] >= variance_threshold)
- {
- bnode.keepMe = true;
- break;
- }
- }
- bnode.expected_value = (ST)expectedValue;
-
- //find min, max using a specified kernel
width
- // 0 = ignore this step
- // 1 = this node only
- // 2 = this node + 1 forward (forward
differences)
- // 3 = this node, 1 forward, 1 backward
(central differences)
- // 4 = this, 2 forward, 1 backward.
- // etc.
- Vec3i kmin;
- kmin.data[0] = 2*i;
- kmin.data[1] = 2*j;
- kmin.data[2] = 2*k;
- Vec3i kmax = kmin;
-
- kmin.data[0] -= (kernel_width-1)/2;
- kmin.data[1] -= (kernel_width-1)/2;
- kmin.data[2] -= (kernel_width-1)/2;
- kmax.data[0] += (kernel_width/2)+1;
- kmax.data[1] += (kernel_width/2)+1;
- kmax.data[2] += (kernel_width/2)+1;
- bnode.min = bnode.max = bnode.values[0];
- if (r==0)
- {
- for(int pi=kmin.data[0];
pi<=kmax.data[0]; pi++)
- for(int pj=kmin.data[1];
pj<=kmax.data[1]; pj++)
- for(int pk=kmin.data[2];
pk<=kmax.data[2]; pk++)
- {
- ST val =
lookup_safe(originalData, pi, pj, pk);
- narrow_range(val, isomin,
isomax);
- bnode.min = MIN(bnode.min,
val);
- bnode.max = MAX(bnode.max,
val);
- }
- }
- else
- {
- int maxbound = (1<<d+1) - 1;
- Array3<BuildNode>& children_bnodes =
buildnodes[d+1];
- kmin = Min(kmin, Vec3i(0,0,0));
- kmax = Max(kmax,
Vec3i(maxbound,maxbound,maxbound));
- for(int pi=kmin.data[0];
pi<=kmax.data[0]; pi++)
- for(int pj=kmin.data[1];
pj<=kmax.data[1]; pj++)
- for(int pk=kmin.data[2];
pk<=kmax.data[2]; pk++)
- {
- ST val = children_bnodes(pi,
pj, pk).expected_value;
- bnode.min = MIN(bnode.min,
val);
- bnode.max = MAX(bnode.max,
val);
- }
- }
-
- ncaps += bnode.keepMe;
- }
- else //if we're not at finest-level depth
- {
- //fill in the children
- Array3<BuildNode>& children_bnodes =
buildnodes[d+1];
- bnode.isLeaf = false;
- bnode.keepMe = true;
- bnode.scalar_leaf = 0;
-
- bnode.myIndex = 0;
- bnode.min = children_bnodes(2*i, 2*j,
2*k).min;
- bnode.max = children_bnodes(2*i, 2*j,
2*k).max;
- for(int pi=0; pi<2; pi++)
- for(int pj=0; pj<2; pj++)
- for(int pk=0; pk<2; pk++)
- {
- int cidx = 4*pi+2*pj+pk;
- BuildNode& child_bnode =
children_bnodes(2*i+pi, 2*j+pj, 2*k+pk);
-
- //find min, max using children
- bnode.min = MIN(bnode.min,
child_bnode.min);
- bnode.max = MAX(bnode.max,
child_bnode.max);
-
- if (child_bnode.keepMe == false)
//make it a scalar leaf
- {
- bnode.scalar_leaf |= (1 <<
cidx);
- bnode.children[cidx].scalar
= child_bnode.expected_value;
- }
- else
- {
- //update the pointer to this
child index
- bnode.children_start =
child_bnode.myStart;
- bnode.children[cidx].offset
= child_bnode.myIndex;
- }
- }
-
- double expectedValue = 0.0;
- for(int v=0; v<8; v++)
- expectedValue +=
static_cast<double>(bnode.values[v]);
- expectedValue *= 0.125;
- bnode.expected_value = (ST)expectedValue;
-
- //if all children are scalar leaves
- if (bnode.scalar_leaf & 255)
- {
- bnode.keepMe = ((bnode.max - bnode.min)
>= variance_threshold);
- }
- nnodes += bnode.keepMe;
- }
- }
-
- //this is just for stats...
- unsigned int maxnodesthisdepth = (1 << d);
- maxnodesthisdepth = maxnodesthisdepth * maxnodesthisdepth *
maxnodesthisdepth;
- totalmaxnodesalldepths += maxnodesthisdepth;
-
- if (d == max_depth-1-r)
- {
- cerr << "\nAllocating "<<ncaps<< " caps at depth " << d<<
"(" << 100.0f*ncaps/(float)maxnodesthisdepth <<"% full)";
- totalnodesalldepths += ncaps;
- steps[ts].num_caps[r] = ncaps;
- steps[ts].caps[r] = new OctCap[ncaps];
- }
- else
- {
- cerr << "\nAllocating "<<nnodes<< " nodes at depth " << d<<
"(" << 100.0f*nnodes/(float)maxnodesthisdepth <<"% full)";
- totalnodesalldepths += nnodes;
- steps[ts].num_nodes[r][d] = nnodes;
- steps[ts].nodes[r][d] = new OctNode[nnodes];
- }
-
- //fill in caps and nodes, storing blocks of 8 (sort of like
bricking)
- unsigned int cap = 0;
- unsigned int node = 0;
- if (d == max_depth-1-r) //bottom nodes
- {
- //NEW
- unsigned int actualnodesthislvl = 1 << d;
- actualnodesthislvl = actualnodesthislvl * actualnodesthislvl
* actualnodesthislvl;
- unsigned int n=0;
- unsigned int cap_start = n;
- while(n < actualnodesthislvl)
- {
- //find the x, y, z coordinates
- Vec3i coords(0,0,0);
- unsigned int remainder = n;
- for(int dc=d; dc>=0; dc--)
- {
- int divisor = (1 << 3*dc);
- int div = remainder / divisor;
- remainder = remainder % divisor;
- if (div)
- {
- coords.data[0] |= ((div&4)!=0) << dc;
- coords.data[1] |= ((div&2)!=0) << dc;
- coords.data[2] |= (div&1) << dc;
- }
- }
- BuildNode& bnode = bnodes(coords.data[0],
coords.data[1], coords.data[2]);
- if (bnode.keepMe)
- {
- OctCap& ocap = steps[ts].caps[r][cap];
- memcpy(ocap.values, bnode.values, 8*sizeof(SO));
- bnode.myIndex = (unsigned char)(cap - cap_start);
- bnode.myStart = cap_start;
- cap++;
- }
- n++;
- if (n%8==0)
- cap_start = cap;
- }
- }
- else if (d > 0) //inner nodes
- {
- unsigned int actualnodesthislvl = 1 << d;
- actualnodesthislvl = actualnodesthislvl * actualnodesthislvl
* actualnodesthislvl;
- unsigned int n=0;
- unsigned int node_start = n;
- while(n < actualnodesthislvl)
- {
- //find the x, y, z coordinates
- Vec3i coords(0,0,0);
- unsigned int remainder = n;
- for(int dc=d; dc>=0; dc--)
- {
- int divisor = (1 << 3*dc);
- int div = remainder / divisor;
- remainder = remainder % divisor;
- if (div)
- {
- coords.data[0] |= ((div&4)!=0) << dc;
- coords.data[1] |= ((div&2)!=0) << dc;
- coords.data[2] |= (div&1) << dc;
- }
- }
- BuildNode& bnode = bnodes(coords.data[0],
coords.data[1], coords.data[2]);
- if (bnode.keepMe)
- {
- OctNode& onode = steps[ts].nodes[r][d][node];
- Vec3i child_coords = coords * 2;
- #pragma unroll(8)
- for(int c=0; c<8; c++)
- {
- BuildNode& child_bnode =
buildnodes[d+1](child_coords.data[0] + ((c&4)!=0), child_coords.data[1] +
((c&2)!=0), child_coords.data[2] + (c&1));
- onode.child_mins[c] = child_bnode.min;
- onode.child_maxs[c] = child_bnode.max;
- }
-
- memcpy(onode.children, bnode.children, 8*sizeof(SO));
- onode.scalar_leaf = bnode.scalar_leaf;
- onode.children_start = bnode.children_start;
- bnode.myIndex = (unsigned char)(node - node_start);
- bnode.myStart = node_start;
- node++;
- }
- n++;
- if (n%8==0)
- node_start = node;
- }
- }
- else //root node; can't group nodes into 8 because we only
have 1!
- {
- BuildNode& bnode = bnodes(0,0,0);
- OctNode& onode = steps[ts].nodes[r][d][node];
- for(int c=0; c<8; c++)
- {
- BuildNode& child_bnode = buildnodes[1]((c&4)!=0,
(c&2)!=0, (c&1));
- onode.child_mins[c] = child_bnode.min;
- onode.child_maxs[c] = child_bnode.max;
- }
-
- memcpy(onode.children, bnode.children, 8*sizeof(SO));
- onode.scalar_leaf = 0;
- onode.children_start = bnode.children_start;
- bnode.myIndex = node;
- node++;
- }
- //cout<< "caps filled = " << cap << endl;
- //cout<< "nodes filled = " << node << endl;
- }
-
- //set the min and max
- OctNode& rootnode = steps[ts].nodes[r][0][0];
- datamin = rootnode.child_mins[0];
- datamax = rootnode.child_maxs[0];
- for(int c=1; c<8; c++)
- {
- datamin = MIN(datamin, rootnode.child_mins[c]);
- datamax = MAX(datamax, rootnode.child_maxs[c]);
- }
-
- datamin = MAX(datamin, isomin);
- datamax = MIN(datamax, isomax);
-
- //cout << "datamin = " << (int)datamin << ", datamax = " <<
(int)datamax << endl;
-
- cout << "\nMultires level " << r << ":\n Total octree nodes = " <<
totalnodesalldepths << "/" << totalmaxnodesalldepths << " ~= " <<
totalnodesalldepths*100/(float)totalmaxnodesalldepths << "% full.\n\n";
- }
+ //build the buildnodes, from bottom up.
+ for(int d=max_depth-1; d>=0; d--)
+ {
+ Array3<BuildNode>& bnodes = buildnodes[d];
+ unsigned int ncaps = 0;
+ unsigned int nnodes = 0;
+
+ int nodes_dim = 1 << d;
+ //cout << "\nd="<<d<<":";
+ //cout << "nodes_dim = " << nodes_dim << endl;
+
+ int dim_extents = 1 << (max_depth - d - 1);
+ //cout << "dim_extents = " << dim_extents << endl;
+
+ for(int i=0; i<nodes_dim; i++)
+ for(int j=0; j<nodes_dim; j++)
+ for(int k=0; k<nodes_dim; k++)
+ {
+ BuildNode& bnode = bnodes(i,j,k);
+
+ if (d==max_depth-1)
+ {
+ for(int pi=0; pi<2; pi++)
+ for(int pj=0; pj<2;
pj++)
+ for(int pk=0;
pk<2; pk++)
+ {
+ Vec3i
dp((2*i + pi) * dim_extents,
+
(2*j + pj) * dim_extents,
+
(2*k + pk) * dim_extents);
+ ST
odval = lookup_safe(originalData, dp.data[0], dp.data[1], dp.data[2]);
+
narrow_range(odval, isomin, isomax);
+
bnode.values[4*pi | 2*pj | pk] = odval;
+ }
+ }
+ else
+ {
+ for(int pi=0; pi<2; pi++)
+ for(int pj=0; pj<2;
pj++)
+ for(int pk=0;
pk<2; pk++)
+ {
+
bnode.values[4*pi | 2*pj | pk] = buildnodes[d+1](2*i+pi, 2*j+pj,
2*k+pk).expected_value;
+ }
+ }
+
+ if (d == max_depth-1)
//finest-level data are automatically caps, *if* they are kept
+ {
+
+ double expectedValue = 0.0;
+ for(int v=0; v<8; v++)
+ expectedValue +=
static_cast<double>(bnode.values[v]);
+ expectedValue *= 0.125;
+
+ double variances[8];
+ bnode.keepMe = false;
+ bnode.isLeaf = true;
+ for(int v=0; v<8; v++)
+ {
+ double diff =
static_cast<double>(bnode.values[v]) - expectedValue;
+ variances[v] = diff *
diff * 0.125;
+
+ if (variances[v] >=
variance_threshold)
+ {
+ bnode.keepMe
= true;
+ break;
+ }
+ }
+ bnode.expected_value =
(ST)expectedValue;
+
+ //find min, max using a
specified kernel width
+ // 0 = ignore this step
+ // 1 = this node only
+ // 2 = this node + 1 forward
(forward differences)
+ // 3 = this node, 1 forward,
1 backward (central differences)
+ // 4 = this, 2 forward, 1
backward.
+ // etc.
+ Vec3i kmin;
+ kmin.data[0] = 2*i;
+ kmin.data[1] = 2*j;
+ kmin.data[2] = 2*k;
+ Vec3i kmax = kmin;
+
+ kmin.data[0] -=
(kernel_width-1)/2;
+ kmin.data[1] -=
(kernel_width-1)/2;
+ kmin.data[2] -=
(kernel_width-1)/2;
+ kmax.data[0] +=
(kernel_width/2)+1;
+ kmax.data[1] +=
(kernel_width/2)+1;
+ kmax.data[2] +=
(kernel_width/2)+1;
+ bnode.min = bnode.max =
bnode.values[0];
+ for(int pi=kmin.data[0];
pi<=kmax.data[0]; pi++)
+ for(int
pj=kmin.data[1]; pj<=kmax.data[1]; pj++)
+ for(int
pk=kmin.data[2]; pk<=kmax.data[2]; pk++)
+ {
+ ST
val = lookup_safe(originalData, pi, pj, pk);
+
narrow_range(val, isomin, isomax);
+
bnode.min = MIN(bnode.min, val);
+
bnode.max = MAX(bnode.max, val);
+ }
+ ncaps += bnode.keepMe;
+ }
+ else //if we're not at
finest-level depth
+ {
+ //fill in the children
+ Array3<BuildNode>&
children_bnodes = buildnodes[d+1];
+ bnode.isLeaf = false;
+ bnode.keepMe = true;
+
+ bnode.myIndex = 0;
+ bnode.min =
children_bnodes(2*i, 2*j, 2*k).min;
+ bnode.max =
children_bnodes(2*i, 2*j, 2*k).max;
+ for(int pi=0; pi<2; pi++)
+ for(int pj=0; pj<2;
pj++)
+ for(int pk=0;
pk<2; pk++)
+ {
+ int
cidx = 4*pi+2*pj+pk;
+
BuildNode& child_bnode = children_bnodes(2*i+pi, 2*j+pj, 2*k+pk);
+
+
//find min, max using children
+
bnode.min = MIN(bnode.min, child_bnode.min);
+
bnode.max = MAX(bnode.max, child_bnode.max);
+
+ if
(child_bnode.keepMe == false) //make it a scalar leaf
+ {
+
bnode.offsets[cidx] = -1;
+ }
+ else
+ {
+
//update the pointer to this child index
+
bnode.children_start = child_bnode.myStart;
+
bnode.offsets[cidx] = child_bnode.myIndex;
+ }
+
+
bnode.values[cidx] = child_bnode.expected_value;
+
+ }
+
+ double expectedValue = 0.0;
+ for(int v=0; v<8; v++)
+ expectedValue +=
static_cast<double>(bnode.values[v]);
+ expectedValue *= 0.125;
+ bnode.expected_value =
(ST)expectedValue;
+
+ //if all children are scalar
leaves
+ char chsum = 0;
+ #pragma unroll(8)
+ for(int ch=0; ch<8; ch++)
+ chsum +=
bnode.offsets[ch];
+ if (chsum == -8)
+ {
+ bnode.keepMe =
((bnode.max - bnode.min) >= variance_threshold);
+ }
+
+ /*
+ if (bnode.scalar_leaf & 255)
+ {
+ bnode.keepMe =
((bnode.max - bnode.min) >= variance_threshold);
+ }
+ */
+
+ nnodes += bnode.keepMe;
+
+ } //end (not at finest-level
depth)
+ } //end (iterate over i,j,k)
+
+ //this is just for stats...
+ unsigned int maxnodesthisdepth = (1 << d);
+ maxnodesthisdepth = maxnodesthisdepth * maxnodesthisdepth *
maxnodesthisdepth;
+ totalmaxnodesalldepths += maxnodesthisdepth;
+
+ if (d == max_depth-1)
+ {
+ cerr << "\nAllocating "<<ncaps<< " caps at depth " <<
d<< "(" << 100.0f*ncaps/(float)maxnodesthisdepth <<"% full)";
+ totalnodesalldepths += ncaps;
+ steps[ts].num_caps = ncaps;
+ steps[ts].caps = new OctCap[ncaps];
+ }
+ else
+ {
+ cerr << "\nAllocating "<<nnodes<< " nodes at depth "
<< d<< "(" << 100.0f*nnodes/(float)maxnodesthisdepth <<"% full)";
+ totalnodesalldepths += nnodes;
+ steps[ts].num_nodes[d] = nnodes;
+ steps[ts].nodes[d] = new OctNode[nnodes];
+ }
+
+ //fill in caps and nodes, storing blocks of 8 (sort of like
bricking)
+ unsigned int cap = 0;
+ unsigned int node = 0;
+ if (d == max_depth-1) //bottom nodes
+ {
+ //NEW
+ unsigned int actualnodesthislvl = 1 << d;
+ actualnodesthislvl = actualnodesthislvl *
actualnodesthislvl * actualnodesthislvl;
+ unsigned int n=0;
+ unsigned int cap_start = n;
+ while(n < actualnodesthislvl)
+ {
+ //find the x, y, z coordinates
+ Vec3i coords(0,0,0);
+ unsigned int remainder = n;
+ for(int dc=d; dc>=0; dc--)
+ {
+ int divisor = (1 << 3*dc);
+ int div = remainder / divisor;
+ remainder = remainder % divisor;
+ if (div)
+ {
+ coords.data[0] |=
((div&4)!=0) << dc;
+ coords.data[1] |=
((div&2)!=0) << dc;
+ coords.data[2] |= (div&1) <<
dc;
+ }
+ }
+ BuildNode& bnode = bnodes(coords.data[0],
coords.data[1], coords.data[2]);
+ if (bnode.keepMe)
+ {
+ OctCap& ocap = steps[ts].caps[cap];
+ memcpy(ocap.values, bnode.values,
8*sizeof(ST));
+ bnode.myIndex = (unsigned char)(cap -
cap_start);
+ bnode.myStart = cap_start;
+ cap++;
+ }
+ n++;
+ if (n%8==0)
+ cap_start = cap;
+ }
+ }
+ else if (d > 0) //inner nodes
+ {
+ unsigned int actualnodesthislvl = 1 << d;
+ actualnodesthislvl = actualnodesthislvl *
actualnodesthislvl * actualnodesthislvl;
+ unsigned int n=0;
+ unsigned int node_start = n;
+ while(n < actualnodesthislvl)
+ {
+ //find the x, y, z coordinates
+ Vec3i coords(0,0,0);
+ unsigned int remainder = n;
+ for(int dc=d; dc>=0; dc--)
+ {
+ int divisor = (1 << 3*dc);
+ int div = remainder / divisor;
+ remainder = remainder % divisor;
+ if (div)
+ {
+ coords.data[0] |=
((div&4)!=0) << dc;
+ coords.data[1] |=
((div&2)!=0) << dc;
+ coords.data[2] |= (div&1) <<
dc;
+ }
+ }
+ BuildNode& bnode = bnodes(coords.data[0],
coords.data[1], coords.data[2]);
+ if (bnode.keepMe)
+ {
+ OctNode& onode =
steps[ts].nodes[d][node];
+ memcpy(onode.offsets, bnode.offsets,
8*sizeof(ST));
+
+ Vec3i child_coords = coords * 2;
+ #pragma unroll(8)
+ for(int c=0; c<8; c++)
+ {
+ BuildNode& child_bnode =
buildnodes[d+1](child_coords.data[0] + ((c&4)!=0), child_coords.data[1] +
((c&2)!=0), child_coords.data[2] + (c&1));
+ onode.mins[c] =
child_bnode.min;
+ onode.maxs[c] =
child_bnode.max;
+ }
+ onode.children_start =
bnode.children_start;
+ memcpy(onode.values, bnode.values,
8*sizeof(ST));
+
+ bnode.myIndex = (unsigned char)(node
- node_start);
+ bnode.myStart = node_start;
+ node++;
+ }
+ n++;
+ if (n%8==0)
+ node_start = node;
+ }
+ }
+ else //root node; can't group nodes into 8 because we only
have 1!
+ {
+ BuildNode& bnode = bnodes(0,0,0);
+ OctNode& onode = steps[ts].nodes[d][node];
+ memcpy(onode.offsets, bnode.offsets, 8);
+ for(int c=0; c<8; c++)
+ {
+ BuildNode& child_bnode =
buildnodes[1]((c&4)!=0, (c&2)!=0, (c&1));
+ onode.mins[c] = child_bnode.min;
+ onode.maxs[c] = child_bnode.max;
+ }
+
+ memcpy(onode.values, bnode.values, 8*sizeof(ST));
+ onode.children_start = bnode.children_start;
+ bnode.myIndex = node;
+ node++;
+ }
+ //cout<< "caps filled = " << cap << endl;
+ //cout<< "nodes filled = " << node << endl;
+ }
+
+ //set the min and max
+ OctNode& rootnode = steps[ts].nodes[0][0];
+ datamin = rootnode.mins[0];
+ datamax = rootnode.maxs[0];
+ for(int c=1; c<8; c++)
+ {
+ datamin = MIN(datamin, rootnode.mins[c]);
+ datamax = MAX(datamax, rootnode.maxs[c]);
+ }
+
+ datamin = MAX(datamin, isomin);
+ datamax = MIN(datamax, isomax);
+
+ //cout << "datamin = " << (int)datamin << ", datamax = " <<
(int)datamax << endl;
+
+ cout << "Total octree nodes = " << totalnodesalldepths << "/" <<
totalmaxnodesalldepths << " ~= " <<
totalnodesalldepths*100/(float)totalmaxnodesalldepths << "% full.\n\n";
delete[] buildnodes;
}
@@ -594,7 +569,9 @@
in >> datamin;
in >> datamax;
in >> num_timesteps;
- in >> mres_levels;
+
+ cap_depth = max_depth - 1;
+ pre_cap_depth = max_depth - 2;
child_bit_depth = new int[max_depth+1];
inv_child_bit_depth = new float[max_depth+1];
@@ -620,19 +597,12 @@
steps = new Timestep[num_timesteps];
for(int ts = 0; ts < num_timesteps; ts++)
{
- steps[ts].num_nodes = new unsigned int*[mres_levels+1];
- steps[ts].nodes = new OctNode**[mres_levels+1];
- steps[ts].caps = new OctCap*[mres_levels+1];
- steps[ts].num_caps = new unsigned int[mres_levels+1];
- for(int r=0; r<mres_levels+1; r++)
- {
- steps[ts].num_nodes[r] = new unsigned int[max_depth-1];
- for(int d=0; d<max_depth-1-r; d++)
- {
- in >> steps[ts].num_nodes[r][d];
- }
- in >> steps[ts].num_caps[r];
- }
+ steps[ts].num_nodes = new unsigned int[max_depth-1];
+ for(int d=0; d<max_depth-1; d++)
+ {
+ in >> steps[ts].num_nodes[d];
+ }
+ in >> steps[ts].num_caps;
}
in.close();
@@ -648,29 +618,26 @@
for(int ts = 0; ts < num_timesteps; ts++)
{
//read in the nodes
- for(int r=0; r<mres_levels+1; r++)
+ steps[ts].nodes = new OctNode*[max_depth-1];
+ for(int d=0; d<max_depth-1; d++)
{
- steps[ts].nodes[r] = new OctNode*[max_depth-1-r];
- for(int d=0; d<max_depth-1-r; d++)
- {
- steps[ts].nodes[r][d] = new
OctNode[steps[ts].num_nodes[r][d]];
- int num_read = fread((char*)steps[ts].nodes[r][d],
sizeof(OctNode), steps[ts].num_nodes[r][d],in2);
- cerr << "Read " << num_read << " nodes" << " at depth " << d
<< endl;
- if (num_read != steps[ts].num_nodes[r][d])
- {
- cerr << "OctreeVolume: error reading octnodes" << endl;
- exit(0);
- }
- }
- //read in the caps
- steps[ts].caps[r] = new OctCap[steps[ts].num_caps[r]];
- int num_read =
fread((char*)steps[ts].caps[r],sizeof(OctCap),steps[ts].num_caps[r],in2);
- cerr << "Read " << num_read << " caps" << endl;
- if (num_read != steps[ts].num_caps[r])
- {
- cerr << "OctreeVolume: error reading octcaps" << endl;
- exit(0);
- }
+ steps[ts].nodes[d] = new
OctNode[steps[ts].num_nodes[d]];
+ int num_read = fread((char*)steps[ts].nodes[d],
sizeof(OctNode), steps[ts].num_nodes[d],in2);
+ cerr << "Read " << num_read << " nodes" << " at depth
" << d << endl;
+ if (num_read != steps[ts].num_nodes[d])
+ {
+ cerr << "OctreeVolume: error reading
octnodes" << endl;
+ exit(0);
+ }
+ }
+ //read in the caps
+ steps[ts].caps = new OctCap[steps[ts].num_caps];
+ int num_read =
fread((char*)steps[ts].caps,sizeof(OctCap),steps[ts].num_caps,in2);
+ cerr << "Read " << num_read << " caps" << endl;
+ if (num_read != steps[ts].num_caps)
+ {
+ cerr << "OctreeVolume: error reading octcaps" << endl;
+ exit(0);
}
}
fclose(in2);
@@ -686,7 +653,7 @@
indata.resize(int_dims.data[0], int_dims.data[1], int_dims.data[2]);
cerr << "Reading indata at " << filename << "...";
cerr << "(" << indata.get_datasize() << " read/" << int_dims.data[0] *
int_dims.data[1] * int_dims.data[2] * sizeof(ST) << " expected)...";
- fread_big((char*)(&indata.get_dataptr()[0][0][0]), 1,
indata.get_datasize(), din);
+ fread_big((void*)(&indata.get_dataptr()[0][0][0]), 1,
indata.get_datasize(), din);
cerr << "done.\n";
fclose(din);
#endif
@@ -711,16 +678,12 @@
out << kernel_width << endl;
out << datamin << " " << datamax << endl;
out << num_timesteps << endl;
- out << mres_levels << endl;
for(int ts = 0; ts < num_timesteps; ts++)
{
- for(int r=0; r<mres_levels+1; r++)
- {
- for(int d=0; d<max_depth-1-r; d++)
- out << steps[ts].num_nodes[r][d] << " ";
- out << endl;
- out << steps[ts].num_caps[r] << endl;
- }
+ for(int d=0; d<max_depth-1; d++)
+ out << steps[ts].num_nodes[d] << " ";
+ out << endl;
+ out << steps[ts].num_caps << endl;
}
out.close();
@@ -731,25 +694,22 @@
for(int ts = 0; ts < num_timesteps; ts++)
{
- for(int r=0; r<mres_levels+1; r++)
+ //write out the nodes
+ for(int d=0; d<max_depth-1; d++)
{
- //write out the nodes
- for(int d=0; d<max_depth-1-r; d++)
- {
- int num_written =
fwrite(reinterpret_cast<void*>(steps[ts].nodes[r][d]), sizeof(OctNode),
steps[ts].num_nodes[r][d], out2);
- if (num_written != steps[ts].num_nodes[r][d])
- {
- cerr << "OctreeVolume: error writing octnodes" << endl;
- exit(0);
- }
- }
- //write out the caps
- int num_written =
fwrite(reinterpret_cast<void*>(steps[ts].caps[r]), sizeof(OctCap),
steps[ts].num_caps[r], out2);
- if (num_written != steps[ts].num_caps[r])
- {
- cerr << "OctreeVolume: error writing octcaps" << endl;
- exit(0);
- }
+ int num_written =
fwrite(reinterpret_cast<void*>(steps[ts].nodes[d]), sizeof(OctNode),
steps[ts].num_nodes[d], out2);
+ if (num_written != steps[ts].num_nodes[d])
+ {
+ cerr << "OctreeVolume: error writing
octnodes" << endl;
+ exit(0);
+ }
+ }
+ //write out the caps
+ int num_written =
fwrite(reinterpret_cast<void*>(steps[ts].caps), sizeof(OctCap),
steps[ts].num_caps, out2);
+ if (num_written != steps[ts].num_caps)
+ {
+ cerr << "OctreeVolume: error writing octcaps" << endl;
+ exit(0);
}
}
fclose(out2);
Modified: trunk/Model/Primitives/OctreeVolume.h
==============================================================================
--- trunk/Model/Primitives/OctreeVolume.h (original)
+++ trunk/Model/Primitives/OctreeVolume.h Sat May 27 15:09:52 2006
@@ -35,6 +35,9 @@
//enable this only when testing octree data and original data side by side
//#define TEST_INDATA
+//enable dynamic multires via the "stop_depth"
+//#define OCTVOL_DYNAMIC_MULTIRES
+
#include <Core/Geometry/BBox.h>
#include <SCIRun/Core/Containers/Array3.h>
#include <Core/Geometry/vecdefs.h>
@@ -62,31 +65,28 @@
{
ST values[8];
};
-
- struct OctNode
- {
- SO children[8];
- ST child_mins[8];
- ST child_maxs[8];
- unsigned int children_start;
- unsigned char scalar_leaf; //mask: tells us if a
child octant is a scalar leaf node
- char padding[3]; //this will
make an OctNode 32 bytes exactly
- };
+
+ struct OctNode
+ {
+ char offsets[8]; //-1 if leaf; 0-7 if child exists
+ ST mins[8];
+ ST maxs[8];
+ unsigned int children_start;
+ ST values[8];
+ };
//node used in build
struct BuildNode
{
- unsigned int myStart;
- unsigned char myIndex;
+ char offsets[8];
+ ST values[8];
+ unsigned int myStart;
unsigned int children_start;
- SO children[8];
+ ST expected_value;
+ ST min, max;
+ unsigned char myIndex;
bool keepMe;
bool isLeaf;
- ST expected_value;
- unsigned char scalar_leaf;
-
- ST values[8];
- ST min, max;
};
//class ST = the scalar value type (e.g. unsigned char)
@@ -95,11 +95,15 @@
{
public:
-#ifdef TEST_INDATA
- SCIRun::Array3<ST> indata;
-#endif
+ // depth of the caps minus 1; hence always == max_depth-2
+ // traversals often check this value
+ int pre_cap_depth;
+
//the maximum depth of the octree, finest resolution
int max_depth;
+
+ //the depth of caps; always == max_depth-2
+ int cap_depth;
//the width of the kernel that is allegedly reading the octree data.
// This determines the ghosting overlap, if any, of the min/max
values
@@ -122,10 +126,10 @@
struct Timestep
{
- OctNode*** nodes; //interior nodes, e.g. from depth = 0 to
max_depth-2
- OctCap** caps; //caps *would* be leaves in a full
octree.
- unsigned int** num_nodes;
- unsigned int* num_caps;
+ OctNode** nodes; //interior nodes, e.g. from depth = 0 to
max_depth-2
+ OctCap* caps; //caps *would* be leaves in a full
octree.
+ unsigned int* num_nodes;
+ unsigned int num_caps;
};
Timestep* steps;
@@ -141,6 +145,10 @@
ST isoval;
bool use_adaptive_res;
/* END VARS SHOULD BE IN OCTISOVOLUME */
+
+ #ifdef TEST_INDATA
+ SCIRun::Array3<ST> indata;
+ #endif
inline float getMaxTime() const
{
@@ -213,6 +221,17 @@
{
return max_depth;
}
+
+ inline int get_cap_depth() const
+ {
+ return cap_depth;
+ }
+
+
+ inline int get_pre_cap_depth() const
+ {
+ return pre_cap_depth;
+ }
OctreeVolume();
@@ -226,7 +245,6 @@
int endTimestep,
double variance_threshold,
int kernelWidth,
- int mresLevels,
ST isomin,
ST isomax);
@@ -234,83 +252,93 @@
inline OctNode& get_node(int depth, unsigned int index) const
{
- return
steps[current_timestep].nodes[multires_level][depth][index];
+ return steps[current_timestep].nodes[depth][index];
}
inline OctCap& get_cap(unsigned int index) const
{
- return steps[current_timestep].caps[multires_level][index];
+ return steps[current_timestep].caps[index];
}
inline ST operator()(Vec3i& v) const
{
- return lookup_node(v, 0, 0, 0);
+ return lookup_node(v, cap_depth, 0, 0);
}
inline ST operator()(int d1, int d2, int d3) const
{
//do a lookup in octree data, on [0, 1 << max_depth]
Vec3i p(d1, d2, d3);
- return lookup_node(p, 0, 0, 0);
+ return lookup_node(p, cap_depth, 0, 0);
}
inline ST locate_trace(int d1, int d2, int d3, unsigned int*
index_trace) const
{
//do a lookup in octree data
Vec3i p(d1, d2, d3);
- return lookup_node_trace(p, 0, 0, 0, index_trace);
+ return lookup_node_trace(p, cap_depth, 0, 0, index_trace);
}
- inline ST lookup_node_trace(Vec3i& p, int res, int depth, unsigned
int index, unsigned int* index_trace) const
+ inline ST lookup_node_trace(Vec3i& p, int stop_depth, int depth,
unsigned int index, unsigned int* index_trace) const
{
for(;;)
{
- OctNode& node =
steps[current_timestep].nodes[res][depth][index];
+ OctNode& node = steps[current_timestep].nodes[depth][index];
index_trace[depth] = index;
-
- int child_bit = child_bit_depth[depth+res];
+
+ int child_bit = child_bit_depth[depth];
int target_child = (((p.data[0] & child_bit)!=0) << 2) |
(((p.data[1] & child_bit)!=0) << 1) | ((p.data[2] & child_bit)!=0);
- if (node.scalar_leaf & (1 << target_child))
- return node.children[target_child].scalar;
- if (depth == max_depth-2-res)
+ #ifdef OCTVOL_DYNAMIC_MULTIRES
+ if (depth == stop_depth)
+ return node.values[target_child];
+ #endif
+
+ if (node.offsets[target_child] == -1)
+ return node.values[target_child];
+ if (depth == pre_cap_depth)
{
- index = node.children_start +
node.children[target_child].offset;
+ index = node.children_start + node.offsets[target_child];
index_trace[depth+1] = index;
- int target_child = ((p.data[0] & 1) << 2) | ((p.data[1]
& 1) << 1) | (p.data[2] & 1);
- return
steps[current_timestep].caps[res][index].values[target_child];
+ target_child = ((p.data[0] & 1) << 2) | ((p.data[1] & 1)
<< 1) | (p.data[2] & 1);
+ return
steps[current_timestep].caps[index].values[target_child];
}
-
- index = node.children_start +
node.children[target_child].offset;
+
+ index = node.children_start + node.offsets[target_child];
depth++;
}
}
- inline ST lookup_node(Vec3i& p, int res, int depth, unsigned int
index) const
+ inline ST lookup_node(Vec3i& p, int stop_depth, int depth, unsigned
int index) const
{
for(;;)
{
- OctNode& node =
steps[current_timestep].nodes[res][depth][index];
+ OctNode& node = steps[current_timestep].nodes[depth][index];
- int child_bit = child_bit_depth[depth+res];
+ int child_bit = child_bit_depth[depth];
int target_child = (((p.data[0] & child_bit)!=0) << 2) |
(((p.data[1] & child_bit)!=0) << 1) | ((p.data[2] & child_bit)!=0);
+
+ #ifdef OCTVOL_DYNAMIC_MULTIRES
+ if (depth == stop_depth)
+ return node.values[target_child];
+ #endif
- if (node.scalar_leaf & (1 << target_child))
- return node.children[target_child].scalar;
- if (depth == max_depth-2-res)
+ if (node.offsets[target_child] == -1)
+ return node.values[target_child];
+ if (depth == pre_cap_depth)
{
- index = node.children_start +
node.children[target_child].offset;
+ index = node.children_start + node.offsets[target_child];
int target_child = ((p.data[0] & 1) << 2) | ((p.data[1]
& 1) << 1) | (p.data[2] & 1);
- return
steps[current_timestep].caps[res][index].values[target_child];
+ return
steps[current_timestep].caps[index].values[target_child];
}
- index = node.children_start +
node.children[target_child].offset;
+ index = node.children_start + node.offsets[target_child];
depth++;
}
}
template<bool CHECK_X, bool CHECK_Y, bool CHECK_Z>
- inline ST lookup_neighbor(Vec3i& cell, Vec3i& offset, int res, int
depth, unsigned int* index_trace) const
+ inline ST lookup_neighbor(Vec3i& cell, Vec3i& offset, int
stop_depth, int depth, unsigned int* index_trace) const
{
Vec3i target_neighbor = cell + offset;
int child_bit;
@@ -318,47 +346,47 @@
//recurse up the tree until we find a parent that contains both
cell and target_neighbor
for(int up = depth; up >= 0; up--)
{
- child_bit = child_bit_depth[up+res];
+ child_bit = child_bit_depth[up];
if (CHECK_X && CHECK_Y && CHECK_Z)
{
if ( ((target_neighbor.data[0] & child_bit) ==
(cell.data[0] & child_bit)) &&
((target_neighbor.data[1] & child_bit) ==
(cell.data[1] & child_bit)) &&
((target_neighbor.data[2] & child_bit) ==
(cell.data[2] & child_bit)) )
- return lookup_node(target_neighbor, res, up,
index_trace[up]);
+ return lookup_node(target_neighbor, stop_depth, up,
index_trace[up]);
}
else if (CHECK_X && CHECK_Y)
{
if ( ((target_neighbor.data[0] & child_bit) ==
(cell.data[0] & child_bit)) &&
((target_neighbor.data[1] & child_bit) ==
(cell.data[1] & child_bit)) )
- return lookup_node(target_neighbor, res, up,
index_trace[up]);
+ return lookup_node(target_neighbor, stop_depth, up,
index_trace[up]);
}
else if (CHECK_X && CHECK_Z)
{
if ( ((target_neighbor.data[0] & child_bit) ==
(cell.data[0] & child_bit)) &&
((target_neighbor.data[2] & child_bit) ==
(cell.data[2] & child_bit)) )
- return lookup_node(target_neighbor, res, up,
index_trace[up]);
+ return lookup_node(target_neighbor, stop_depth, up,
index_trace[up]);
}
else if (CHECK_Y && CHECK_Z)
{
if ( ((target_neighbor.data[1] & child_bit) ==
(cell.data[1] & child_bit)) &&
((target_neighbor.data[2] & child_bit) ==
(cell.data[2] & child_bit)) )
- return lookup_node(target_neighbor, res, up,
index_trace[up]);
+ return lookup_node(target_neighbor, stop_depth, up,
index_trace[up]);
}
else if (CHECK_X)
{
if ((target_neighbor.data[0] & child_bit) ==
(cell.data[0] & child_bit))
- return lookup_node(target_neighbor, res, up,
index_trace[up]);
+ return lookup_node(target_neighbor, stop_depth, up,
index_trace[up]);
}
else if (CHECK_Y)
{
if ((target_neighbor.data[1] & child_bit) ==
(cell.data[1] & child_bit))
- return lookup_node(target_neighbor, res, up,
index_trace[up]);
+ return lookup_node(target_neighbor, stop_depth, up,
index_trace[up]);
}
else if (CHECK_Z)
{
if ((target_neighbor.data[2] & child_bit) ==
(cell.data[2] & child_bit))
- return lookup_node(target_neighbor, res, up,
index_trace[up]);
+ return lookup_node(target_neighbor, stop_depth, up,
index_trace[up]);
}
}
@@ -368,40 +396,40 @@
if ( ((target_neighbor.data[0] & child_bit) == (cell.data[0]
& child_bit)) &&
((target_neighbor.data[1] & child_bit) == (cell.data[1]
& child_bit)) &&
((target_neighbor.data[2] & child_bit) == (cell.data[2]
& child_bit)) )
- return lookup_node(target_neighbor, res, 0, 0);
+ return lookup_node(target_neighbor, stop_depth, 0, 0);
}
else if (CHECK_X && CHECK_Y)
{
if ( ((target_neighbor.data[0] & child_bit) == (cell.data[0]
& child_bit)) &&
((target_neighbor.data[1] & child_bit) == (cell.data[1]
& child_bit)) )
- return lookup_node(target_neighbor, res, 0, 0);
+ return lookup_node(target_neighbor, stop_depth, 0, 0);
}
else if (CHECK_X && CHECK_Z)
{
if ( ((target_neighbor.data[0] & child_bit) == (cell.data[0]
& child_bit)) &&
((target_neighbor.data[2] & child_bit) == (cell.data[2]
& child_bit)) )
- return lookup_node(target_neighbor, res, 0, 0);
+ return lookup_node(target_neighbor, stop_depth, 0, 0);
}
else if (CHECK_Y && CHECK_Z)
{
if ( ((target_neighbor.data[1] & child_bit) == (cell.data[1]
& child_bit)) &&
((target_neighbor.data[2] & child_bit) == (cell.data[2]
& child_bit)) )
- return lookup_node(target_neighbor, res, 0, 0);
+ return lookup_node(target_neighbor, stop_depth, 0, 0);
}
else if (CHECK_X)
{
if ((target_neighbor.data[0] & child_bit) == (cell.data[0] &
child_bit))
- return lookup_node(target_neighbor, res, 0, 0);
+ return lookup_node(target_neighbor, stop_depth, 0, 0);
}
else if (CHECK_Y)
{
if ((target_neighbor.data[1] & child_bit) == (cell.data[1] &
child_bit))
- return lookup_node(target_neighbor, res, 0, 0);
+ return lookup_node(target_neighbor, stop_depth, 0, 0);
}
else if (CHECK_Z)
{
if ((target_neighbor.data[2] & child_bit) == (cell.data[2] &
child_bit))
- return lookup_node(target_neighbor, res, 0, 0);
+ return lookup_node(target_neighbor, stop_depth, 0, 0);
}
return 0;
@@ -417,24 +445,12 @@
int ts,
double variance_threshold,
int kernelWidth,
- int mresLevels = 2,
ST isomin = 0,
ST isomax = 0);
bool read_file(char* filebase);
bool write_file(char* filebase);
-
- private:
-
- /*
- inline ST lookup_safe(SCIRun::Array3<ST>& originalData, int d1, int
d2, int d3) const
- {
- if (d1>=0 && d1<originalData.dim1() && d2>=0 &&
d2<originalData.dim2() && d3>=0 && d3<originalData.dim3())
- return originalData(d1, d2, d3);
- else
- return 0;
- }
- */
+
};
};
Modified: trunk/scenes/octisovol.cc
==============================================================================
--- trunk/scenes/octisovol.cc (original)
+++ trunk/scenes/octisovol.cc Sat May 27 15:09:52 2006
@@ -126,7 +126,7 @@
else if (strcmp(c_buildfrom_filename,"")!=0)
{
cerr << "Creating octree volume from original rectilinear grid
data." << endl;
- ov = new OctreeVolume(c_buildfrom_filename, tstart, tend, variance,
kernel_width, mres_levels, (double)isomin, (double)isomax);
+ ov = new OctreeVolume(c_buildfrom_filename, tstart, tend, variance,
kernel_width, (double)isomin, (double)isomax);
}
else
{
@@ -160,13 +160,13 @@
//scene->setCamera(camera);
LightSet *lights = new LightSet();
- lights->add( new PointLight(Vector(eye + 25.f*(up + Cross(eye-lookat,
up))), Color(RGBColor(.6,.1,.1))) );
- lights->setAmbientLight( new ConstantAmbient(
Color(RGBColor(0.91,0.9,0.9) ) ));
+ lights->add( new PointLight(Vector(eye + 25.f*(up + Cross(eye-lookat,
up))), Color(RGBColor(1,1,1))) );
+ lights->setAmbientLight( new ConstantAmbient(
Color(RGBColor(0.4,0.5,0.5) ) ));
scene->setLights(lights);
// Background.
- scene->setBackground( new ConstantBackground( Color(RGB(0.8, 0.8,
0.8)) ) );
+ scene->setBackground( new ConstantBackground( Color(RGB(0.0, 0.0, 0.0))
) );
return scene;
}
- [MANTA] r1088 - in trunk: Model/Primitives scenes, knolla, 05/27/2006
Archive powered by MHonArc 2.6.16.
