-
sparker
thiago
Manta Interactive Ray Tracer Development Mailing List
Text archives Help
[MANTA] r1063 - in trunk: . Core Core/Geometry Core/Math Engine/PixelSamplers Model/Groups Model/Intersections Model/Primitives Model/Readers/BART scenes
Chronological Thread
- From: knolla@sci.utah.edu
- To: manta@sci.utah.edu
- Subject: [MANTA] r1063 - in trunk: . Core Core/Geometry Core/Math Engine/PixelSamplers Model/Groups Model/Intersections Model/Primitives Model/Readers/BART scenes
- Date: Tue, 16 May 2006 16:59:59 -0600 (MDT)
Author: knolla
Date: Tue May 16 16:59:57 2006
New Revision: 1063
Added:
trunk/Core/Geometry/varray.h (contents, props changed)
trunk/Core/Geometry/vecdefs.h
trunk/Core/Math/CubicSolver.h
trunk/Model/Intersections/IsosurfaceImplicit.cc
trunk/Model/Intersections/IsosurfaceImplicit.h
trunk/Model/Primitives/IsosurfaceOctreeVolume.cc
trunk/Model/Primitives/IsosurfaceOctreeVolume.h
trunk/Model/Primitives/OctreeVolume.cc
trunk/Model/Primitives/OctreeVolume.h
trunk/scenes/octisovol.cc
Removed:
trunk/Model/Groups/varray.h
Modified:
trunk/CMakeLists.txt
trunk/Core/CMakeLists.txt
trunk/Core/Geometry/VectorT.h
trunk/Engine/PixelSamplers/JitterSampler.h
trunk/Model/Groups/CMakeLists.txt
trunk/Model/Groups/GriddedGroup.cc
trunk/Model/Groups/GriddedGroup.h
trunk/Model/Groups/KDTree.cc
trunk/Model/Groups/KDTree.h
trunk/Model/Groups/KDTreeLoader.cc
trunk/Model/Groups/SSEKDTree.cc
trunk/Model/Groups/SSEKDTree.h
trunk/Model/Groups/TransparentKDTree.cc
trunk/Model/Groups/TransparentKDTree.h
trunk/Model/Groups/VerticalKDTree.cc
trunk/Model/Groups/VerticalKDTree.h
trunk/Model/Groups/VolumeGrid.h
trunk/Model/Intersections/CMakeLists.txt
trunk/Model/Primitives/CMakeLists.txt
trunk/Model/Readers/BART/parse.cc
trunk/scenes/CMakeLists.txt
Log:
added cubic solver; typedef support for common vectors (e.g. Vec3i); and
octree volume isosurfacing options
Modified: trunk/CMakeLists.txt
==============================================================================
--- trunk/CMakeLists.txt (original)
+++ trunk/CMakeLists.txt Tue May 16 16:59:57 2006
@@ -163,7 +163,7 @@
# Check for Mac OS
IF (APPLE)
#FIRST_TIME_SET(CMAKE_CXX_FLAGS_RELEASE "-DSCI_ASSERTION_LEVEL=0 -O3 -g
-fgcse-sm -funroll-loops -fstrict-aliasing -fsched-interblock
-falign-loops=16 -falign-jumps=16 -falign-functions=16
-falign-jumps-max-skip=15 -falign-loops-max-skip=15 -ffast-math
-freorder-blocks -mpowerpc-gpopt -force_cpusubtype_ALL -mtune=G5 -mcpu=G5
-mpowerpc64 -faltivec -mabi=altivec -mpowerpc-gfxopt -malign-natural" STRING
"Optimized Flags")
- FIRST_TIME_SET(CMAKE_CXX_FLAGS_RELEASE "-DSCI_ASSERTION_LEVEL=0 -O3 -g
-fgcse-sm -funroll-loops -fstrict-aliasing -fsched-interblock -ffast-math
-freorder-blocks -march=prescott -mtune=prescott -msse -msse2 -msse3
-mfpmath=sse -malign-double" STRING "Optimized Flags")
+ FIRST_TIME_SET(CMAKE_CXX_FLAGS_RELEASE "-DSCI_ASSERTION_LEVEL=0 -O3 -g
-fgcse-sm -funroll-loops -fstrict-aliasing -fsched-interblock -ffast-math
-freorder-blocks -fpermissive -march=prescott -mtune=prescott -msse -msse2
-msse3 -mfpmath=sse -malign-double" STRING "Optimized Flags")
ENDIF (APPLE)
##################################################################
Modified: trunk/Core/CMakeLists.txt
==============================================================================
--- trunk/Core/CMakeLists.txt (original)
+++ trunk/Core/CMakeLists.txt Tue May 16 16:59:57 2006
@@ -32,6 +32,8 @@
Geometry/Vector.h
Geometry/VectorT.cc
Geometry/VectorT.h
+ Geometry/varray.h
+ Geometry/vecdefs.h
)
SET (CORE_SOURCES ${CORE_SOURCES}
Math/HaltonSequence.cc
Modified: trunk/Core/Geometry/VectorT.h
==============================================================================
--- trunk/Core/Geometry/VectorT.h (original)
+++ trunk/Core/Geometry/VectorT.h Tue May 16 16:59:57 2006
@@ -65,6 +65,10 @@
typedef char unnamed[ Dim == 3 ? 1 : 0 ];
data[0] = x; data[1] = y; data[2] = z;
}
+ VectorT(T x, T y, T z, T w) {
+ typedef char unnamed[ Dim == 4 ? 1 : 0 ];
+ data[0] = x; data[1] = y; data[2] = z; data[3] = w;
+ }
VectorT(const Vector& copy) {
typedef char unnamed[ Dim == 3 ? 1 : 0 ];
@@ -304,7 +308,10 @@
return result;
}
+#if 0
+ //un-privatized for speed.
private:
+#endif
T data[Dim];
// Do not use this function!!!! Use getDataPtr() instead.
@@ -384,7 +391,7 @@
result[i] = Interpolate(p1[i], p2[i], weight);
return result;
}
-
+
template<typename T, int Dim>
std::ostream& operator<< (std::ostream& os, const VectorT<T,Dim>& v);
Added: trunk/Core/Geometry/varray.h
==============================================================================
--- (empty file)
+++ trunk/Core/Geometry/varray.h Tue May 16 16:59:57 2006
@@ -0,0 +1,231 @@
+/********************************************************************\
+File: varray.h
+Author: Hansong Zhang
+Department of Computer Science
+UNC-Chapel Hill
+
+$Id: varray.h,v 1.2 1997/09/10 16:00:59 zhangh Exp $
+\********************************************************************/
+
+#ifndef __VARRAY_H__
+#define __VARRAY_H__
+
+
+// #include "alloc.h"
+#include <stdlib.h>
+
+#include <assert.h>
+#include <stdio.h>
+
+//#define DEBUG_ARRAY
+
+//#define INT long long
+#define INT unsigned int
+
+template <class T>
+class VArray {
+protected:
+ INT curLen;
+ INT curMaxLen;
+ T* array;
+public:
+ VArray() {
+ curLen = curMaxLen = 0;
+ array = NULL;
+ }
+
+ VArray(INT initMaxLen) {
+ init(initMaxLen);
+ }
+ VArray(INT initMaxLen, INT initVal) {
+ init(initMaxLen, initVal);
+ }
+
+ void init(INT initMaxLen) {
+ curLen = 0;
+ curMaxLen = initMaxLen;
+ array = (T*)malloc(curMaxLen * sizeof(T));
+ assert (array != NULL);
+ }
+
+ void init(INT initMaxLen, INT initVal) {
+ curLen = 0;
+ curMaxLen = initMaxLen;
+ array = (T*)malloc(curMaxLen * sizeof(T));
+ assert (array != NULL);
+ memset(array, initVal, curMaxLen*sizeof(T));
+ }
+
+ inline INT getLen() { return curLen; }
+ inline INT getNum() { return curLen; }
+ inline void setLen(INT len) {
+ curLen = len;
+ if (curMaxLen < len) {
+ curMaxLen = len;
+ if (array)
+ array = (T*)realloc(array,
curMaxLen*sizeof(T));
+ else
+ array = (T*)malloc(curMaxLen*sizeof(T));
+ }
+ assert (array != NULL);
+ }
+ inline void setLen(INT len,const T &val) {
+ curLen = len;
+ if (curMaxLen < len) {
+ curMaxLen = len;
+ if (array)
+ array = (T*)realloc(array,
curMaxLen*sizeof(T));
+ else
+ array = (T*)malloc(curMaxLen*sizeof(T));
+ }
+
+ for (INT i=0;i<len;++i) {
+ array[i] = val;
+ }
+
+ assert (array != NULL);
+ }
+ inline void resize(INT len) {
+ curMaxLen = len;
+ array = (T*)realloc(array, curMaxLen*sizeof(T));
+ assert (array != NULL);
+ }
+ inline void reset() { curLen=0; }
+ inline T& append() {
+ if (curLen < curMaxLen)
+ curLen ++;
+ else {
+ curMaxLen = (INT)(curMaxLen*1.3f);
+ if (array)
+ array = (T*)realloc(array,
curMaxLen*sizeof(T));
+ else
+ array = (T*)malloc(curMaxLen*sizeof(T));
+
+ assert(array != NULL);
+ curLen ++;
+ }
+ return array[curLen-1];
+ }
+
+ inline void append(const T& data) {
+ if (curMaxLen == 0) {
+ curLen = 0;
+ curMaxLen = 8;
+ array = (T*)malloc(curMaxLen * sizeof(T));
+ }
+
+ if (curLen < curMaxLen) {
+ array[curLen ++] = data;
+ } else {
+ // fprintf(stderr, "cur:%d, max:%d\n", curLen,
curMaxLen);
+ curMaxLen = (INT)(curMaxLen*1.3f);
+ if (array)
+ array = (T*)realloc(array,
curMaxLen*sizeof(T));
+ else
+ array = (T*)malloc(curMaxLen*sizeof(T));
+ assert(array != NULL);
+ array[curLen ++] = data;
+ }
+ }
+
+ void fit() {
+ if (curLen != curMaxLen) {
+ array = (T*)realloc(array, curLen*sizeof(T));
+ assert(array != NULL);
+ curMaxLen = curLen;
+ }
+ }
+
+ inline T &operator[](INT i) {
+ //if (array && i < curLen)
+ return array[i];
+ //else {
+ // fprintf(stderr, "VArray(%p)::[] reference out of
bound (idx:%d size:%d)\n", this, i, curLen);
+ // assert(0);
+ // return array[0];
+ //}
+ }
+
+ inline T & get(INT i) {
+#ifdef DEBUG_ARRAY
+ if (i >= curLen) {
+ fprintf(stderr, "VArray(%p)::get() reference out of
bound (%d).\n", this, i);
+ assert(0);
+ return array[0];
+ }
+#endif
+ return (*this)[i];
+ }
+
+ inline T & _get(INT i) {
+#ifdef DEBUG_ARRAY
+ if (i >= curLen) {
+ fprintf(stderr, "VArray(%p)::_get() reference out of
bound (%d).\n", this, i);
+ assert(0);
+ return array[0];
+ }
+#endif
+ return array[i];
+ }
+
+ inline void set(INT i, const T& elem) {
+#ifdef DEBUG_ARRAY
+ if (array && i >= curLen) {
+ fprintf(stderr, "VArray(%p)::set() out of range
(%d).\n", this, i);
+ }
+#endif
+ array[i] = elem;
+ }
+
+ inline void _set(INT i, const T& elem) {
+ array[i] = elem;
+ }
+
+ inline T & last() {
+ return array[curLen-1];
+ }
+
+ inline T *getArray() {
+ return array;
+ }
+ inline T *dupArray() {
+ T *ret = (T*)malloc(curLen * sizeof(T));
+ if (! ret) {
+ fprintf(stderr, "Cannot alloc mem in dupArray().\n");
+ return NULL;
+ }
+ memcpy(ret, array, curLen * sizeof(T));
+ return ret;
+ }
+
+ ~VArray() { if (array) free(array); }
+ void freemem() {
+ if (array) {
+ free(array);
+ curLen = curMaxLen = 0;
+ }
+ }
+
+ void clear() {
+ if (array)
+ memset(array, 0, sizeof(T)*curMaxLen);
+ }
+
+};
+
+/*
+template <class T>
+class VPool : public VArray<T> {
+ public:
+ VPool(INT size):VArray<T>(size) {}
+ T& alloc(INT * offset=NULL) {
+ append();
+ if (offset) *offset = curLen-1;
+ return array[curLen-1];
+ }
+};
+*/
+
+
+
+#endif
Added: trunk/Core/Geometry/vecdefs.h
==============================================================================
--- (empty file)
+++ trunk/Core/Geometry/vecdefs.h Tue May 16 16:59:57 2006
@@ -0,0 +1,49 @@
+//vecdefs: common typedefs for the vector class
+
+
+#ifndef Manta_Core_vecdefs_h
+#define Manta_Core_vecdefs_h
+
+#include <MantaTypes.h>
+#include <Core/Geometry/VectorT.h>
+
+namespace Manta {
+ typedef VectorT<Real, 2> Vector2D;
+
+ typedef VectorT<float, 2> Vec2f;
+ typedef VectorT<float, 3> Vec3f;
+ typedef VectorT<float, 4> Vec4f;
+
+ typedef VectorT<int, 2> Vec2i;
+ typedef VectorT<int, 3> Vec3i;
+ typedef VectorT<int, 4> Vec4i;
+
+} // end namespace Manta
+
+#endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Added: trunk/Core/Math/CubicSolver.h
==============================================================================
--- (empty file)
+++ trunk/Core/Math/CubicSolver.h Tue May 16 16:59:57 2006
@@ -0,0 +1,147 @@
+/*
+ * Utility functions to find cubic roots,
+ * coefficients are passed like this:
+ *
+ * c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 = 0
+ *
+ * The functions return the number of non-complex roots and
+ * put the values into the s array.
+ *
+ * Author: Jochen Schwarze (schwarze@isa.de)
+ *
+ * Jan 26, 1990 Version for Graphics Gems
+ * Oct 11, 1990 Fixed sign problem for negative q's in SolveQuartic
+ * (reported by Mark Podlipec),
+ * Old-style function definitions,
+ * IsZero() as a macro
+ * Nov 23, 1990 Some systems do not declare acos() and cbrt() in
+ * <math.h>, though the functions exist in the library.
+ * If large coefficients are used, EQN_EPS should be
+ * reduced considerably (e.g. to 1E-30), results will be
+ * correct but multiple roots might be reported more
+ * than once.
+ *
+ * Aaron Knoll note: this is code really requires double precision.
+ * This enormously affects speed. In general, this should ONLY
+ * be used as reference code.
+ */
+
+#ifndef __MANTA_CORE_CUBICSOLVER_H__
+#define __MANTA_CORE_CUBICSOLVER_H__
+
+#include <math.h>
+#include <stdio.h>
+
+namespace Manta
+{
+
+#define EQN_EPS 1e-30
+#define IsZero(x) ((x) > -EQN_EPS && (x) < EQN_EPS)
+
+#ifdef NOCBRT
+#define cbrt(x) ((x) > 0.0 ? pow((x), 1.0/3.0) : \
+ ((x) < 0.0 ? -pow(-(x), 1.0/3.0) : 0.0))
+#endif
+
+ inline int SolveLinear(double *c, double *s) {
+ if (c[ 1 ] < 0.0000001 && c[ 1 ] > -0.0000001)
+ return 0;
+ s[0] = -c[0]/c[1];
+ return 1;
+ }
+
+ inline int SolveQuadratic(double *c, double *s) {
+ double A, B, C;
+ if (c[ 2 ] < 0.0000001 && c[ 2 ] > -0.0000001)
+ return SolveLinear(c, s);
+
+ /* Ax^2 + Bx + C = 0 */
+ A = c[2];
+ B = c[1];
+ C = c[0];
+
+ double disc = B*B - 4*A*C;
+ if (disc < 0) return 0;
+
+ disc = sqrt(disc);
+ s[0] = (-B - disc) / (2*A);
+ s[1] = (-B + disc) / (2*A);
+ return 2;
+ }
+
+ inline int SolveCubic(double c[4], double s[3])
+ {
+ int i, num;
+ double sub;
+ double A, B, C;
+ double sq_A, p, q;
+ double cb_p, D;
+
+ if (c[ 3 ] < 0.0000001 && c[ 3 ] > -0.0000001)
+ return SolveQuadratic(c, s);
+
+ /* normal form: x^3 + Ax^2 + Bx + C = 0 */
+
+ A = c[ 2 ] / c[ 3 ];
+ B = c[ 1 ] / c[ 3 ];
+ C = c[ 0 ] / c[ 3 ];
+
+ /* substitute x = y - A/3 to eliminate quadric term:
+ x^3 +px + q = 0 */
+
+ sq_A = A * A;
+ p = 1.0/3 * (- 1.0/3 * sq_A + B);
+ q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);
+
+ /* use Cardano's formula */
+
+ cb_p = p * p * p;
+ D = q * q + cb_p;
+
+ if (IsZero(D))
+ {
+ if (IsZero(q)) /* one triple solution */
+ {
+ s[ 0 ] = 0;
+ num = 1;
+ }
+ else /* one single and one double solution */
+ {
+ double u = cbrt(-q);
+ s[ 0 ] = 2 * u;
+ s[ 1 ] = - u;
+ num = 2;
+ }
+ }
+ else if (D < 0) /* Casus irreducibilis: three real solutions */
+ {
+ double phi = 1.0/3 * acos(-q / sqrt(-cb_p));
+ double t = 2 * sqrt(-p);
+
+ s[ 0 ] = t * cosf(phi);
+ s[ 1 ] = - t * cosf(phi + M_PI / 3);
+ s[ 2 ] = - t * cosf(phi - M_PI / 3);
+ num = 3;
+ }
+ else /* one real solution */
+ {
+ double sqrt_D = sqrt(D);
+ double u = cbrt(sqrt_D - q);
+ double v = - cbrt(sqrt_D + q);
+
+ s[ 0 ] = u + v;
+ num = 1;
+ }
+
+ /* resubstitute */
+
+ sub = 1.0/3 * A;
+
+ for (i = 0; i < num; ++i)
+ s[ i ] -= sub;
+
+ return num;
+ }
+};
+
+#endif
Modified: trunk/Engine/PixelSamplers/JitterSampler.h
==============================================================================
--- trunk/Engine/PixelSamplers/JitterSampler.h (original)
+++ trunk/Engine/PixelSamplers/JitterSampler.h Tue May 16 16:59:57 2006
@@ -4,7 +4,7 @@
#include <Interface/PixelSampler.h>
#include <Core/Math/MT_RNG.h>
-#include <Core/Math/vector2d.h>
+#include <Core/Geometry/vecdefs.h>
#include <sgi_stl_warnings_off.h>
#include <string>
Modified: trunk/Model/Groups/CMakeLists.txt
==============================================================================
--- trunk/Model/Groups/CMakeLists.txt (original)
+++ trunk/Model/Groups/CMakeLists.txt Tue May 16 16:59:57 2006
@@ -31,7 +31,6 @@
Groups/TransparentKDTree.cc
Groups/TransparentKDTree.h
Groups/VolumeGrid.h
- Groups/varray.h
)
# Determine if VerticalKDTree can be included
Modified: trunk/Model/Groups/GriddedGroup.cc
==============================================================================
--- trunk/Model/Groups/GriddedGroup.cc (original)
+++ trunk/Model/Groups/GriddedGroup.cc Tue May 16 16:59:57 2006
@@ -122,7 +122,7 @@
cerr << "4. Allocating " << sum << " list entries\n";
lists = new const Object*[sum];
- Array3<int> offset(nx, ny, nz);
+ GridArray3<int> offset(nx, ny, nz);
offset.initialize(0);
cerr << "5. Filling in lists\n";
for(int i=0;i<static_cast<int>(objs.size());i++){
Modified: trunk/Model/Groups/GriddedGroup.h
==============================================================================
--- trunk/Model/Groups/GriddedGroup.h (original)
+++ trunk/Model/Groups/GriddedGroup.h Tue May 16 16:59:57 2006
@@ -10,19 +10,19 @@
namespace Manta {
template<typename T>
- class Array3 {
+ class GridArray3 {
public:
- Array3() {
+ GridArray3() {
nx = ny = nz = extra = 0;
xidx = yidx = zidx = 0;
data = 0;
}
- Array3(int nx, int ny, int nz, int extra = 0) {
+ GridArray3(int nx, int ny, int nz, int extra = 0) {
xidx = yidx = zidx = 0;
data = 0;
resize(nx, ny, nz, extra);
}
- ~Array3() {
+ ~GridArray3() {
resize(0, 0, 0);
}
@@ -104,8 +104,8 @@
private:
- Array3(const Array3&);
- Array3& operator=(const Array3&);
+ GridArray3(const GridArray3&);
+ GridArray3& operator=(const GridArray3&);
int nx, ny, nz, extra;
int* xidx;
@@ -139,7 +139,7 @@
Vector diag;
Vector cellsize;
Vector inv_cellsize;
- Array3<int> cells;
+ GridArray3<int> cells;
const Object** lists;
void transformToLattice(const BBox& box,
Modified: trunk/Model/Groups/KDTree.cc
==============================================================================
--- trunk/Model/Groups/KDTree.cc (original)
+++ trunk/Model/Groups/KDTree.cc Tue May 16 16:59:57 2006
@@ -29,7 +29,7 @@
*/
#include <Model/Groups/KDTree.h>
-#include <Model/Groups/varray.h>
+#include <Core/Geometry/varray.h>
#include <Model/Intersections/AxisAlignedBox.h>
#include <Model/Intersections/TriangleEdge.h>
Modified: trunk/Model/Groups/KDTree.h
==============================================================================
--- trunk/Model/Groups/KDTree.h (original)
+++ trunk/Model/Groups/KDTree.h Tue May 16 16:59:57 2006
@@ -43,7 +43,7 @@
#include <Model/Groups/KDTreeLoader.h>
#include <Model/Groups/KDTreeLoaderIW.h>
-#include <Model/Groups/varray.h>
+#include <Core/Geometry/varray.h>
#include <Model/Readers/V3C1.h>
Modified: trunk/Model/Groups/KDTreeLoader.cc
==============================================================================
--- trunk/Model/Groups/KDTreeLoader.cc (original)
+++ trunk/Model/Groups/KDTreeLoader.cc Tue May 16 16:59:57 2006
@@ -38,7 +38,7 @@
#include <stdio.h>
-#include "varray.h"
+#include <Core/Geometry/varray.h>
#include "KDTree.h"
#include <Model/Intersections/AxisAlignedBox.h>
Modified: trunk/Model/Groups/SSEKDTree.cc
==============================================================================
--- trunk/Model/Groups/SSEKDTree.cc (original)
+++ trunk/Model/Groups/SSEKDTree.cc Tue May 16 16:59:57 2006
@@ -29,7 +29,7 @@
*/
#include <Model/Groups/SSEKDTree.h>
-#include <Model/Groups/varray.h>
+#include <Core/Geometry/varray.h>
#include <Model/Intersections/AxisAlignedBox.h>
#include <Interface/Context.h>
Modified: trunk/Model/Groups/SSEKDTree.h
==============================================================================
--- trunk/Model/Groups/SSEKDTree.h (original)
+++ trunk/Model/Groups/SSEKDTree.h Tue May 16 16:59:57 2006
@@ -44,7 +44,7 @@
#include <Interface/Texture.h>
#include <Model/Groups/KDTree.h>
-#include <Model/Groups/varray.h>
+#include <Core/Geometry/varray.h>
#include <KdtreeParameters.h>
#include <Model/Groups/VerticalKDTree.h>
Modified: trunk/Model/Groups/TransparentKDTree.cc
==============================================================================
--- trunk/Model/Groups/TransparentKDTree.cc (original)
+++ trunk/Model/Groups/TransparentKDTree.cc Tue May 16 16:59:57 2006
@@ -29,7 +29,7 @@
*/
#include <Model/Groups/TransparentKDTree.h>
-#include <Model/Groups/varray.h>
+#include <Core/Geometry/varray.h>
#include <Model/Intersections/AxisAlignedBox.h>
Modified: trunk/Model/Groups/TransparentKDTree.h
==============================================================================
--- trunk/Model/Groups/TransparentKDTree.h (original)
+++ trunk/Model/Groups/TransparentKDTree.h Tue May 16 16:59:57 2006
@@ -38,7 +38,7 @@
#include <Model/Groups/Group.h>
#include <Model/Primitives/PrimitiveCommon.h>
-#include <Model/Groups/varray.h>
+#include <Core/Geometry/varray.h>
#include <Model/Groups/KDTree.h>
namespace Manta {
Modified: trunk/Model/Groups/VerticalKDTree.cc
==============================================================================
--- trunk/Model/Groups/VerticalKDTree.cc (original)
+++ trunk/Model/Groups/VerticalKDTree.cc Tue May 16 16:59:57 2006
@@ -29,7 +29,7 @@
*/
#include <Model/Groups/VerticalKDTree.h>
-#include <Model/Groups/varray.h>
+#include <Core/Geometry/varray.h>
#include <Model/Intersections/AxisAlignedBox.h>
#include <Interface/Context.h>
Modified: trunk/Model/Groups/VerticalKDTree.h
==============================================================================
--- trunk/Model/Groups/VerticalKDTree.h (original)
+++ trunk/Model/Groups/VerticalKDTree.h Tue May 16 16:59:57 2006
@@ -44,7 +44,7 @@
#include <Interface/Texture.h>
#include <Model/Groups/KDTree.h>
-#include <Model/Groups/varray.h>
+#include <Core/Geometry/varray.h>
#include <KdtreeParameters.h>
#include <iostream>
Modified: trunk/Model/Groups/VolumeGrid.h
==============================================================================
--- trunk/Model/Groups/VolumeGrid.h (original)
+++ trunk/Model/Groups/VolumeGrid.h Tue May 16 16:59:57 2006
@@ -15,14 +15,13 @@
// The grid is a 3-D array of uniformly sized rectangular cells.
//
-#include "varray.h"
+#include <Core/Geometry/varray.h>
#include "VolumeDesc.h"
#include "PsiGammaTable.h"
#define GAMMA_TABLE_SIZE 512
namespace Manta {
-#define Vec3f VectorT<float,3>
-#define Vec4f VectorT<float,4>
+
struct AABox3f {
Vec3f min, max;
Vec3f &operator[] (int i) const {
@@ -474,7 +473,7 @@
//VoxelType _minVal, _maxVal;
Vec3f _voxelScale;
VolumeRenderingParams *_vrp;
- CopyColorMaterial *_ccm;
+ *_ccm;
public:
VolumeGrid() : __Grid< MultiResVolumeCube<VoxelType>
> () {
_ccm = new CopyColorMaterial;
Modified: trunk/Model/Intersections/CMakeLists.txt
==============================================================================
--- trunk/Model/Intersections/CMakeLists.txt (original)
+++ trunk/Model/Intersections/CMakeLists.txt Tue May 16 16:59:57 2006
@@ -1,5 +1,7 @@
SET(Manta_Intersections_SRCS
Intersections/AxisAlignedBox.h
+ Intersections/IsosurfaceImplicit.cc
+ Intersections/IsosurfaceImplicit.h
Intersections/TriangleEdge.h
Intersections/Plane.h
)
Added: trunk/Model/Intersections/IsosurfaceImplicit.cc
==============================================================================
--- (empty file)
+++ trunk/Model/Intersections/IsosurfaceImplicit.cc Tue May 16 16:59:57
2006
@@ -0,0 +1,118 @@
+
+#include <Model/Intersections/IsosurfaceImplicit.h>
+#include <Core/Math/CubicSolver.h>
+
+using namespace Manta;
+
+//From Steven Parker's 1997 RTRT isosurface intersection
+bool IsosurfaceImplicit::single_intersect(RayPacket& rays, int which_one,
+ const Vector& pmin, const Vector& pmax, float
rho[2][2][2],
+ float iso, float tmin, float tmax, float& t)
+{
+ double c[4];
+ double s[3];
+ double ua[2];
+ double ub[2];
+ double va[2];
+ double vb[2];
+ double wa[2];
+ double wb[2];
+ int i,j,k;
+
+ Vector orig = rays.getOrigin(which_one);
+ Vector dir = rays.getDirection(which_one);
+
+ ua[1] = (orig.x() - pmin.x()) / (pmax.x() - pmin.x());
+ ua[0] = 1 - ua[1];
+ ub[1] = dir.x() / (pmax.x() - pmin.x());
+ ub[0] = - ub[1];
+
+ va[1] = (orig.y() - pmin.y()) / (pmax.y() - pmin.y());
+ va[0] = 1 - va[1];
+ vb[1] = dir.y() / (pmax.y() - pmin.y());
+ vb[0] = - vb[1];
+
+ wa[1] = (orig.z() - pmin.z()) / (pmax.z() - pmin.z());
+ wa[0] = 1 - wa[1];
+ wb[1] = dir.z() / (pmax.z() - pmin.z());
+ wb[0] = - wb[1];
+
+
+ c[3] = c[2] = c[1] = c[0] = 0;
+ for (i=0; i < 2; i++)
+ for (j=0; j < 2; j++)
+ for (k=0; k < 2; k++) {
+
+ // cubic term
+ c[3] += ub[i]*vb[j]*wb[k] * rho[i][j][k];
+
+ // square term
+ c[2] += (ua[i]*vb[j]*wb[k] + ub[i]*va[j]*wb[k] +
ub[i]*vb[j]*wa[k]) * rho[i][j][k];
+
+ // linear term
+ c[1] += (ub[i]*va[j]*wa[k] + ua[i]*vb[j]*wa[k] +
ua[i]*va[j]*wb[k]) * rho[i][j][k];
+
+ // constant term
+ c[0] += ua[i]*va[j]*wa[k] * rho[i][j][k];
+ }
+ c[0] -= iso;
+
+ //gives us the roots of the cubic. s are the roots, n are the number of
roots.
+ int n = SolveCubic( c, s);
+
+
+ t = tmax;
+ for (i = 0; i < n; i++)
+ {
+ if (s[i] >= tmin && s[i] < t)
+ t = s[i];
+ }
+
+ return (t < tmax);
+
+}
+
+//From Steven Parker's 1997 RTRT isosurface intersection
+void IsosurfaceImplicit::single_normal(Vector& outNormal, const Vector&
pmin, const Vector& pmax, const Vector& p, float rho[2][2][2])
+{
+ float x = p.x();
+ float y = p.y();
+ float z = p.z();
+
+ float x_0 = pmin.x();
+ float y_0 = pmin.y();
+ float z_0 = pmin.z();
+
+ float x_1 = pmax.x();
+ float y_1 = pmax.y();
+ float z_1 = pmax.z();
+
+ outNormal.data[0] = - (y_1-y)*(z_1-z)*rho[0][0][0]
+ + (y_1-y)*(z_1-z)*rho[1][0][0]
+ - (y-y_0)*(z_1-z)*rho[0][1][0]
+ - (y_1-y)*(z-z_0)*rho[0][0][1]
+ + (y-y_0)*(z_1-z)*rho[1][1][0]
+ + (y_1-y)*(z-z_0)*rho[1][0][1]
+ - (y-y_0)*(z-z_0)*rho[0][1][1]
+ + (y-y_0)*(z-z_0)*rho[1][1][1];
+
+ outNormal.data[1] = - (x_1-x)*(z_1-z)*rho[0][0][0]
+ - (x-x_0)*(z_1-z)*rho[1][0][0]
+ + (x_1-x)*(z_1-z)*rho[0][1][0]
+ - (x_1-x)*(z-z_0)*rho[0][0][1]
+ + (x-x_0)*(z_1-z)*rho[1][1][0]
+ - (x-x_0)*(z-z_0)*rho[1][0][1]
+ + (x_1-x)*(z-z_0)*rho[0][1][1]
+ + (x-x_0)*(z-z_0)*rho[1][1][1];
+
+ outNormal.data[2] = - (x_1-x)*(y_1-y)*rho[0][0][0]
+ - (x-x_0)*(y_1-y)*rho[1][0][0]
+ - (x_1-x)*(y-y_0)*rho[0][1][0]
+ + (x_1-x)*(y_1-y)*rho[0][0][1]
+ - (x-x_0)*(y-y_0)*rho[1][1][0]
+ + (x-x_0)*(y_1-y)*rho[1][0][1]
+ + (x_1-x)*(y-y_0)*rho[0][1][1]
+ + (x-x_0)*(y-y_0)*rho[1][1][1];
+}
+
+
Added: trunk/Model/Intersections/IsosurfaceImplicit.h
==============================================================================
--- (empty file)
+++ trunk/Model/Intersections/IsosurfaceImplicit.h Tue May 16 16:59:57
2006
@@ -0,0 +1,22 @@
+#ifndef _MANTA_ISOSURFACEINTERSECTOR_H_
+#define _MANTA_ISOSURFACEINTERSECTOR_H_
+
+#include <Interface/RayPacket.h>
+#include <Core/Geometry/Vector.h>
+#include <Interface/Material.h>
+
+namespace Manta
+{
+ struct IsosurfaceImplicit
+ {
+ static bool single_intersect(RayPacket& rays, int which_one,
+ const Vector& pmin, const Vector& pmax, float
rho[2][2][2],
+ float iso, float tmin, float tmax, float& t);
+
+ static void single_normal(Vector& outNormal, const Vector& pmin,
+ const Vector& pmax, const Vector& p, float
rho[2][2][2]);
+ };
+};
+
+#endif
+
Modified: trunk/Model/Primitives/CMakeLists.txt
==============================================================================
--- trunk/Model/Primitives/CMakeLists.txt (original)
+++ trunk/Model/Primitives/CMakeLists.txt Tue May 16 16:59:57 2006
@@ -16,6 +16,10 @@
Primitives/Heightfield.h
Primitives/Hemisphere.cc
Primitives/Hemisphere.h
+ Primitives/OctreeVolume.cc
+ Primitives/OctreeVolume.h
+ Primitives/IsosurfaceOctreeVolume.cc
+ Primitives/IsosurfaceOctreeVolume.h
Primitives/Parallelogram.cc
Primitives/Parallelogram.h
Primitives/ParticleBVH.cc
Added: trunk/Model/Primitives/IsosurfaceOctreeVolume.cc
==============================================================================
--- (empty file)
+++ trunk/Model/Primitives/IsosurfaceOctreeVolume.cc Tue May 16 16:59:57
2006
@@ -0,0 +1,755 @@
+/*
+ Aaron Knoll
+ Single-ray, non-SSE octree isosurface intersector
+*/
+
+#include <Model/Primitives/IsosurfaceOctreeVolume.h>
+#include <Interface/RayPacket.h>
+#include <Model/Intersections/IsosurfaceImplicit.h>
+
+#define USE_OCTREE_DATA
+
+#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;
+
+static const int axis_table[] = {4, 2, 1};
+
+IsosurfaceOctreeVolume::IsosurfaceOctreeVolume(OctreeVolume* _octdata,
Material* _matl)
+: PrimitiveCommon(_matl), octdata(_octdata)
+{
+}
+
+IsosurfaceOctreeVolume::~IsosurfaceOctreeVolume()
+{
+}
+
+void IsosurfaceOctreeVolume::preprocess( PreprocessContext const &context )
+{
+}
+
+void IsosurfaceOctreeVolume::computeBounds( PreprocessContext const
&context, BBox &box ) const
+{
+ box = octdata->get_bounds();
+}
+
+void IsosurfaceOctreeVolume::computeNormal(const RenderContext& context,
RayPacket& rays) const
+{
+}
+
+BBox IsosurfaceOctreeVolume::getBounds() const
+{
+ return octdata->get_bounds();
+}
+
+bool IsosurfaceOctreeVolume::bbox_intersects(const BBox& box, const
RayPacket& rays) const
+{
+ for (int ray = rays.begin(); ray < rays.end(); ray++ )
+ {
+ float maximum_minimum = 1e-5;
+ float minimum_maximum = rays.getMinT(ray);
+
+ float x_minimum = (box[rays.getSign(ray,0)][0] -
rays.getOrigin(ray,0)) * rays.getInverseDirection(ray,0);
+ float x_maximum = (box[1-rays.getSign(ray,0)][0] -
rays.getOrigin(ray,0)) * rays.getInverseDirection(ray,0);
+
+ float y_minimum = (box[rays.getSign(ray,1)][1] -
rays.getOrigin(ray,1)) * rays.getInverseDirection(ray,1);
+ float y_maximum = (box[1-rays.getSign(ray,1)][1] -
rays.getOrigin(ray,1)) * rays.getInverseDirection(ray,1);
+
+ float z_minimum = (box[rays.getSign(ray,2)][2] -
rays.getOrigin(ray,2)) * rays.getInverseDirection(ray,2);
+ float z_maximum = (box[1-rays.getSign(ray,2)][2] -
rays.getOrigin(ray,2)) * rays.getInverseDirection(ray,2);
+
+ if ( minimum_maximum < x_minimum ||
+ maximum_minimum > x_maximum )
+ continue;
+ if ( minimum_maximum > x_maximum )
+ minimum_maximum = x_maximum;
+ if ( maximum_minimum < x_minimum )
+ maximum_minimum = x_minimum;
+ if ( minimum_maximum < y_minimum ||
+ maximum_minimum > y_maximum )
+ continue;
+ if ( minimum_maximum > y_maximum )
+ minimum_maximum = y_maximum;
+ if ( maximum_minimum < y_minimum )
+ maximum_minimum = y_minimum;
+
+ if ( minimum_maximum >= z_minimum &&
+ maximum_minimum <= z_maximum )
+ return true; // found a hit
+ }
+ return false;
+}
+
+void IsosurfaceOctreeVolume::intersect(RenderContext const &context,
RayPacket &packet) const
+{
+ if(bbox_intersects(octdata->get_bounds(), packet))
+ {
+ for ( int i = packet.rayBegin; i < packet.rayEnd; i++ )
+ single_intersect(packet, i);
+ }
+}
+
+void IsosurfaceOctreeVolume::single_intersect(RayPacket& rays, int
which_one) const
+{
+ //find tenter, texit, penter using ray-AABB intersection
+ Vector ws_dir = rays.getDirection(which_one);
+ Vector ws_inv_dir = rays.getInverseDirection(which_one);
+ Vector ws_orig = rays.getOrigin(which_one);
+ Vector min = octdata->get_bounds().getMin();
+ Vector max = octdata->get_bounds().getMax();
+
+ //first, transform ray from world space to octree space on [0, 1 <<
octdata->get_max_depth()]^3
+ Vector orig =
octdata->get_inv_node_width(octdata->get_max_depth()-1-octdata->get_multires_level())
* (ws_orig - min);
+ Vector dir =
octdata->get_inv_node_width(octdata->get_max_depth()-1-octdata->get_multires_level())
* ws_dir;
+ Vector inv_dir =
octdata->get_node_width(octdata->get_max_depth()-1-octdata->get_multires_level())
* ws_inv_dir;
+
+ //XXX - implement improved ray/AABB for this? (precompute signs?)
+ float tenter, texit;
+ float max_width = (float)(octdata->get_child_bit_depth(0)<<1);
+ if(dir.x() > 0)
+ {
+ tenter = inv_dir.x() * (-orig.x());
+ texit = inv_dir.x() * (max_width-orig.x());
+ }
+ else
+ {
+ tenter = inv_dir.x() * (max_width-orig.x());
+ texit = inv_dir.x() * (-orig.x());
+ }
+
+ float y0, y1;
+ if(dir.y() > 0)
+ {
+ y0 = inv_dir.y() * (-orig.y());
+ y1 = inv_dir.y() * (max_width-orig.y());
+ }
+ else
+ {
+ y0 = inv_dir.y() * (max_width-orig.y());
+ y1 = inv_dir.y() * (-orig.y());
+ }
+ tenter = MAX(tenter, y0);
+ texit = MIN(texit, y1);
+ if (tenter > texit)
+ return;
+
+ float z0, z1;
+ if(dir.z() > 0)
+ {
+ z0 = inv_dir.z() * (-orig.z());
+ z1 = inv_dir.z() * (max_width-orig.z());
+ }
+ else
+ {
+ z0 = inv_dir.z() * (max_width-orig.z());
+ z1 = inv_dir.z() * (-orig.z());
+ }
+ tenter = MAX(tenter, z0);
+ texit = MIN(texit, z1);
+ if (tenter > texit)
+ return;
+
+ if (texit < 0.)
+ return;
+
+ tenter = MAX(0., tenter);
+
+ unsigned int index_trace[octdata->get_max_depth()+1];
+
+ Vec3i cell(0,0,0);
+
+ single_traverse_node(rays, which_one, orig, dir, inv_dir,
octdata->get_multires_level(), 0, 0, index_trace, cell, tenter, texit);
+}
+
+bool IsosurfaceOctreeVolume::single_traverse_node(RayPacket& rays, int
which_one,
+ const Vector& orig, const Vector& dir, const
Vector& inv_dir, int res,
+ int depth, int node_index, unsigned int
index_trace[], Vec3i& cell, const float tenter,
+ const float texit) const
+{
+ OctNode& node = octdata->get_node(depth, node_index);
+
+ index_trace[depth] = node_index;
+
+ int child_bit = octdata->get_child_bit_depth(depth+res);
+ 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);
+
+ Vector penter = orig + dir*tenter;
+
+ Vec3i child_cell = cell;
+ Vec3i tc( penter.x() >= center.x(), penter.y() >= center.y(), penter.z()
>= center.z() );
+ int target_child = (tc.data[0] << 2) | (tc.data[1] << 1) | tc.data[2];
+ child_cell.data[0] |= tc.data[0] ? child_bit : 0;
+ child_cell.data[1] |= tc.data[1] ? child_bit : 0;
+ child_cell.data[2] |= tc.data[2] ? child_bit : 0;
+
+ Vec3i axis_isects;
+ if (tcenter.data[0] < tcenter.data[1] && tcenter.data[0] <
tcenter.data[2]){
+ axis_isects[0] = 0;
+ if (tcenter.data[1] < tcenter.data[2]){
+ axis_isects[1] = 1;
+ axis_isects[2] = 2;
+ }
+ else{
+ axis_isects[1] = 2;
+ axis_isects[2] = 1;
+ }
+ }
+ else if (tcenter.data[1] < tcenter.data[2]){
+ axis_isects[0] = 1;
+ if (tcenter.data[0] < tcenter.data[2]){
+ axis_isects[1] = 0;
+ axis_isects[2] = 2;
+ }
+ else{
+ axis_isects[1] = 2;
+ axis_isects[2] = 0;
+ }
+ }
+ else{
+ axis_isects[0] = 2;
+ if (tcenter.data[0] < tcenter.data[1]){
+ axis_isects[1] = 0;
+ axis_isects[2] = 1;
+ }
+ else{
+ axis_isects[1] = 1;
+ axis_isects[2] = 0;
+ }
+ }
+
+ const int axis_table[] = {4, 2, 1};
+
+ float child_tenter = tenter;
+ float child_texit;
+
+ int next_depth = depth+1;
+
+ for(int i=0;; i++)
+ {
+ int axis = -1;
+ if (i>=3)
+ child_texit = texit;
+ else
+ {
+ axis = axis_isects[i];
+ child_texit = tcenter.data[axis];
+ if (child_texit < tenter)
+ continue;
+ if (child_texit >= texit){
+ child_texit = texit;
+ axis = -1;
+ }
+ }
+
+ if (octdata->get_isovalue() >= node.child_mins[target_child] &&
octdata->get_isovalue() <= node.child_maxs[target_child])
+ {
+ if (node.scalar_leaf & (1 << target_child))
+ {
+ if (single_traverse_leaf(rays, which_one, orig, dir,
inv_dir, res,
+ next_depth, depth, node.children[target_child].scalar,
+ child_cell, index_trace, child_cell, child_tenter,
child_texit))
+ return true;
+ }
+ else
+ {
+ int child_idx = node.children_start +
node.children[target_child].offset;
+ if (next_depth+1 == octdata->get_max_depth()-res) //cap
+ {
+ if (single_traverse_cap(rays, which_one, orig, dir,
inv_dir,
+ res, 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,
+ index_trace, child_cell, child_tenter, child_texit))
+ return true;
+ }
+ }
+ }
+
+ if (axis==-1)
+ return false;
+
+ //move to the next target_child, update tenter and texit
+ child_tenter = child_texit;
+ int trueaxisbit = axis_table[axis];
+ if (target_child & trueaxisbit) //going from true to
false
+ {
+ target_child &= ~trueaxisbit;
+ child_cell[axis] &= ~child_bit;
+ }
+ else
//going from false to true
+ {
+ target_child |= trueaxisbit;
+ child_cell[axis] |= child_bit;
+ }
+ }
+ return false;
+}
+
+bool IsosurfaceOctreeVolume::single_traverse_leaf(RayPacket& rays, int
which_one,
+ const Vector& orig, const Vector& dir, const Vector&
inv_dir, int res,
+ 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);
+
+ 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);
+ Vector penter = orig + dir*tenter;
+
+ Vec3i child_cell = cell;
+ Vec3i tc( penter.x() >= center.x(), penter.y() >= center.y(), penter.z()
>= center.z() );
+ int target_child = (tc.data[0] << 2) | (tc.data[1] << 1) | tc.data[2];
+ child_cell.data[0] |= tc.data[0] ? child_bit : 0;
+ child_cell.data[1] |= tc.data[1] ? child_bit : 0;
+ child_cell.data[2] |= tc.data[2] ? child_bit : 0;
+
+ Vec3i axis_isects;
+ if (tcenter.data[0] < tcenter.data[1] && tcenter.data[0] <
tcenter.data[2]){
+ axis_isects[0] = 0;
+ if (tcenter.data[1] < tcenter.data[2]){
+ axis_isects[1] = 1;
+ axis_isects[2] = 2;
+ }
+ else{
+ axis_isects[1] = 2;
+ axis_isects[2] = 1;
+ }
+ }
+ else if (tcenter.data[1] < tcenter.data[2]){
+ axis_isects[0] = 1;
+ if (tcenter.data[0] < tcenter.data[2]){
+ axis_isects[1] = 0;
+ axis_isects[2] = 2;
+ }
+ else{
+ axis_isects[1] = 2;
+ axis_isects[2] = 0;
+ }
+ }
+ else{
+ axis_isects[0] = 2;
+ if (tcenter.data[0] < tcenter.data[1]){
+ axis_isects[1] = 0;
+ axis_isects[2] = 1;
+ }
+ else{
+ axis_isects[1] = 1;
+ axis_isects[2] = 0;
+ }
+ }
+
+ const int axis_table[] = {4, 2, 1};
+
+ float child_tenter = tenter;
+ float child_texit;
+
+ int next_depth = depth+1;
+
+ for(int i=0;; i++)
+ {
+ int axis = -1;
+ if (i>=3)
+ child_texit = texit;
+ else
+ {
+ axis = axis_isects[i];
+ child_texit = tcenter.data[axis];
+ if (child_texit < tenter)
+ continue;
+ if (child_texit >= texit){
+ child_texit = texit;
+ axis = -1;
+ }
+ }
+
+ 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)
+ {
+ //try isosurface intersection
+ Vector child_cell_float(child_cell.data[0],
child_cell.data[1], child_cell.data[2]);
+ Vector cmin = (child_cell_float *
octdata->get_node_width(depth)) + octdata->get_bounds().getMin();
+ Vector cmax = cmin + octdata->get_node_width(depth);
+
+#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);
+
+ //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);
+ min_rho = MIN(min_rho, this_rho);
+ max_rho = MAX(max_rho, this_rho);
+ }
+ else
+ this_rho = scalar;
+ rho[offset.data[0]][offset.data[1]][offset.data[2]] =
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, leaf_depth, index_trace);
+ min_rho = MIN(min_rho, this_rho);
+ max_rho = MAX(max_rho, this_rho);
+ }
+ else
+ this_rho = scalar;
+ rho[offset.data[0]][offset.data[1]][offset.data[2]] =
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, leaf_depth, index_trace);
+ min_rho = MIN(min_rho, this_rho);
+ max_rho = MAX(max_rho, this_rho);
+ }
+ else
+ this_rho = scalar;
+ rho[offset.data[0]][offset.data[1]][offset.data[2]] =
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, leaf_depth, index_trace);
+ min_rho = MIN(min_rho, this_rho);
+ max_rho = MAX(max_rho, this_rho);
+ }
+ else
+ this_rho = scalar;
+ rho[offset.data[0]][offset.data[1]][offset.data[2]] =
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, leaf_depth, index_trace);
+ min_rho = MIN(min_rho, this_rho);
+ max_rho = MAX(max_rho, this_rho);
+ }
+ else
+ this_rho = scalar;
+ rho[offset.data[0]][offset.data[1]][offset.data[2]] =
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, leaf_depth, index_trace);
+ min_rho = MIN(min_rho, this_rho);
+ max_rho = MAX(max_rho, this_rho);
+ }
+ else
+ this_rho = scalar;
+ rho[offset.data[0]][offset.data[1]][offset.data[2]] =
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, leaf_depth, index_trace);
+ min_rho = MIN(min_rho, this_rho);
+ max_rho = MAX(max_rho, this_rho);
+ }
+ else
+ this_rho = scalar;
+ rho[offset.data[0]][offset.data[1]][offset.data[2]] =
static_cast<float>(this_rho);
+#else
+ //use original grid data
+ float rho[2][2][2];
+ ST min_rho, max_rho;
+#define MYDATA indata //toggle this to octdata if you want to test pure
point location (no neighbor finding)
+ min_rho = max_rho = MYDATA(child_cell);
+ rho[0][0][0] = static_cast<float>(min_rho);
+ for(int c=1; c<8; c++)
+ {
+ Vec3i offset((c&4)!=0, (c&2)!=0, c&1);
+ Vec3i neighboridx = child_cell + offset;
+ ST this_rho = MYDATA(neighboridx);
+ rho[offset.data[0]][offset.data[1]][offset.data[2]] =
static_cast<float>(this_rho);
+ min_rho = MIN(this_rho, min_rho);
+ max_rho = MAX(this_rho, max_rho);
+ }
+#endif
+ if (octdata->get_isovalue() >= min_rho &&
octdata->get_isovalue() <= max_rho)
+ {
+ float hit_t;
+ if (IsosurfaceImplicit::single_intersect(rays,
which_one, cmin, cmax, rho,
+ octdata->get_isovalue(), child_tenter, child_texit,
hit_t))
+ {
+ if (rays.hit(which_one, hit_t,
PrimitiveCommon::getMaterial(), 0, 0))
+ {
+ Vector normal;
+ Vector phit = rays.getOrigin(which_one) +
rays.getDirection(which_one)*hit_t;
+ IsosurfaceImplicit::single_normal(normal, cmin,
cmax, phit, rho);
+ normal.normalize();
+ rays.setNormal(which_one, normal);
+ return true;
+ }
+ }
+ }
+ }
+ else //not at cap-level depth
+ {
+ if (single_traverse_leaf(rays, which_one, orig, dir,
inv_dir, res,
+ next_depth, leaf_depth, scalar, leaf_base_cell,
index_trace, child_cell, child_tenter, child_texit))
+ return true;
+ }
+ }
+
+ if (axis==-1)
+ return false;
+
+ //move to the next target_child, update tenter and texit
+ child_tenter = child_texit;
+ int trueaxisbit = axis_table[axis];
+ if (target_child & trueaxisbit) //going from true to
false
+ {
+ target_child &= ~trueaxisbit;
+ child_cell[axis] &= ~child_bit;
+ }
+ else
//going from false to true
+ {
+ target_child |= trueaxisbit;
+ child_cell[axis] |= child_bit;
+ }
+ }
+ return false;
+}
+
+bool IsosurfaceOctreeVolume::single_traverse_cap(RayPacket& rays, int
which_one,
+ const Vector& orig, const Vector& dir, const
Vector& inv_dir, int res,
+ int depth, int cap_index, unsigned int
index_trace[], Vec3i& cell,
+ const float tenter, const float texit) const
+{
+ OctCap& cap = octdata->get_cap(cap_index);
+ index_trace[depth] = cap_index;
+
+ Vector penter = orig + dir*tenter;
+ Vector center(static_cast<float>(cell.data[0] | 1),
static_cast<float>(cell.data[1] | 1), static_cast<float>(cell.data[2] | 1));
+ Vector tcenter = inv_dir * (center - orig);
+
+ Vec3i child_cell = cell;
+ Vec3i tc( penter.x() >= center.x(), penter.y() >= center.y(), penter.z()
>= center.z() );
+ int target_child = (tc.data[0] << 2) | (tc.data[1] << 1) | tc.data[2];
+ child_cell.data[0] |= tc.data[0];
+ child_cell.data[1] |= tc.data[1];
+ child_cell.data[2] |= tc.data[2];
+
+ Vec3i axis_isects;
+ if (tcenter.data[0] < tcenter.data[1] && tcenter.data[0] <
tcenter.data[2]){
+ axis_isects[0] = 0;
+ if (tcenter.data[1] < tcenter.data[2]){
+ axis_isects[1] = 1;
+ axis_isects[2] = 2;
+ }
+ else{
+ axis_isects[1] = 2;
+ axis_isects[2] = 1;
+ }
+ }
+ else if (tcenter.data[1] < tcenter.data[2]){
+ axis_isects[0] = 1;
+ if (tcenter.data[0] < tcenter.data[2]){
+ axis_isects[1] = 0;
+ axis_isects[2] = 2;
+ }
+ else{
+ axis_isects[1] = 2;
+ axis_isects[2] = 0;
+ }
+ }
+ else{
+ axis_isects[0] = 2;
+ if (tcenter.data[0] < tcenter.data[1]){
+ axis_isects[1] = 0;
+ axis_isects[2] = 1;
+ }
+ else{
+ axis_isects[1] = 1;
+ axis_isects[2] = 0;
+ }
+ }
+
+ const int axis_table[] = {4, 2, 1};
+
+ float child_tenter = tenter;
+ float child_texit;
+
+ for(int i=0;; i++)
+ {
+ int axis = -1;
+ if (i>=3)
+ child_texit = texit;
+ else
+ {
+ axis = axis_isects[i];
+ child_texit = tcenter.data[axis];
+ if (child_texit < tenter)
+ continue;
+ if (child_texit >= texit){
+ child_texit = texit;
+ axis = -1;
+ }
+ }
+
+ //try isosurface intersection in this node
+ Vector child_cell_float(child_cell.data[0], child_cell.data[1],
child_cell.data[2]);
+ Vector cmin = (child_cell_float * octdata->get_node_width(depth)) +
octdata->get_bounds().getMin();
+ Vector cmax = cmin + octdata->get_node_width(depth);
+
+#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[offset.data[0]][offset.data[1]][offset.data[2]] =
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[offset.data[0]][offset.data[1]][offset.data[2]] =
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[offset.data[0]][offset.data[1]][offset.data[2]] =
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[offset.data[0]][offset.data[1]][offset.data[2]] =
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[offset.data[0]][offset.data[1]][offset.data[2]] =
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[offset.data[0]][offset.data[1]][offset.data[2]] =
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[offset.data[0]][offset.data[1]][offset.data[2]] =
static_cast<float>(this_rho);
+#else
+ //use original grid data
+ float rho[2][2][2];
+ ST min_rho, max_rho;
+#define MYDATA indata
+ min_rho = max_rho = MYDATA(child_cell);
+ rho[0][0][0] = static_cast<float>(min_rho);
+ for(int c=1; c<8; c++)
+ {
+ Vec3i offset((c&4)!=0, (c&2)!=0, c&1);
+ Vec3i neighboridx = child_cell + offset;
+ ST this_rho = MYDATA(neighboridx);
+ rho[offset.data[0]][offset.data[1]][offset.data[2]] =
static_cast<float>(this_rho);
+ min_rho = MIN(this_rho, min_rho);
+ max_rho = MAX(this_rho, max_rho);
+ }
+#endif
+ if (octdata->get_isovalue() >= min_rho && octdata->get_isovalue() <=
max_rho)
+ {
+ float hit_t;
+ if (IsosurfaceImplicit::single_intersect(rays, which_one, cmin,
cmax, rho,
+ octdata->get_isovalue(), child_tenter, child_texit, hit_t))
+ {
+ if (rays.hit(which_one, hit_t,
PrimitiveCommon::getMaterial(), 0, 0))
+ {
+ Vector normal;
+ Vector phit = rays.getOrigin(which_one) +
rays.getDirection(which_one)*hit_t;
+ IsosurfaceImplicit::single_normal(normal, cmin, cmax,
phit, rho);
+ normal.normalize();
+ rays.setNormal(which_one, normal);
+ return true;
+ }
+ }
+ }
+
+ if (axis==-1)
+ return false;
+
+ //move to the next target_child, update tenter and texit
+ child_tenter = child_texit;
+ int trueaxisbit = axis_table[axis];
+ if (target_child & trueaxisbit) //going from true to
false
+ {
+ target_child &= ~trueaxisbit;
+ child_cell[axis] &= ~1;
+ }
+ else
//going from false to true
+ {
+ target_child |= trueaxisbit;
+ child_cell[axis] |= 1;
+ }
+ }
+ return false;
+}
+
+
Added: trunk/Model/Primitives/IsosurfaceOctreeVolume.h
==============================================================================
--- (empty file)
+++ trunk/Model/Primitives/IsosurfaceOctreeVolume.h Tue May 16 16:59:57
2006
@@ -0,0 +1,53 @@
+
+#ifndef _MANTA_ISOSURFACE_OCTREEVOLUME_HXX_
+#define _MANTA_ISOSURFACE_OCTREEVOLUME_HXX_
+
+#include <Model/Primitives/PrimitiveCommon.h>
+#include <Core/Geometry/Vector.h>
+#include <Core/Geometry/BBox.h>
+#include <Core/Color/Color.h>
+#include <Interface/Texture.h>
+#include <Model/Primitives/OctreeVolume.h>
+#include <Core/Math/simd.h>
+
+namespace Manta
+{
+ class IsosurfaceOctreeVolume : public PrimitiveCommon
+ {
+ const OctreeVolume* octdata;
+
+ public:
+ IsosurfaceOctreeVolume(OctreeVolume* _octdata, Material* _matl);
+ ~IsosurfaceOctreeVolume();
+
+ void preprocess( PreprocessContext const &context );
+ void computeBounds( PreprocessContext const &context, BBox &box
) const;
+ void intersect( RenderContext const &context, RayPacket &rays )
const;
+ void computeNormal(const RenderContext& context, RayPacket&
rays) const;
+
+ BBox getBounds() const;
+
+ private:
+
+ bool bbox_intersects(const BBox& box, const RayPacket& rays)
const;
+
+ void single_intersect(RayPacket& rays, int which_one) const;
+
+ bool single_traverse_node(RayPacket& rays, int which_one,
+ const Vector& orig, const Vector& dir, const
Vector& inv_dir, int res,
+ int depth, int node_index, unsigned int
index_trace[], Vec3i& cell,
+ const float tenter, const float texit) const;
+
+ bool single_traverse_leaf(RayPacket& rays, int which_one,
+ const Vector& orig, const Vector& dir, const Vector&
inv_dir, int res,
+ int depth, int leaf_depth, ST scalar, Vec3i&
leaf_base_cell,
+ unsigned int index_trace[], Vec3i& cell, const float
tenter, const float texit) const;
+
+ bool single_traverse_cap(RayPacket& rays, int which_one,
+ const Vector& orig, const Vector& dir, const
Vector& inv_dir, int res,
+ int depth, int cap_index, unsigned int
index_trace[], Vec3i& cell, const float tenter,
+ const float texit) const;
+ };
+};
+
+#endif
Added: trunk/Model/Primitives/OctreeVolume.cc
==============================================================================
--- (empty file)
+++ trunk/Model/Primitives/OctreeVolume.cc Tue May 16 16:59:57 2006
@@ -0,0 +1,715 @@
+/*
+* OctreeVolume.h : A multiresolution octree data
+* wrapper for scalar volumes, such as BrickArray3.
+*
+* Written by:
+* Aaron M Knoll
+* Department of Computer Science
+* University of Utah
+* October 2005
+*
+* Copyright (C) 2005 SCI Group
+*/
+
+#include <Model/Primitives/OctreeVolume.h>
+
+using namespace Manta;
+using namespace SCIRun;
+
+#define narrow_range(val, rmin, rmax) \
+ val = MAX(rmin, val);\
+ val = MIN(rmax, val);\
+\
+
+OctreeVolume::OctreeVolume(std::string filebase,
+ int startTimestep,
+ int endTimestep,
+ double variance_threshold,
+ int kernelWidth,
+ int mresLevels,
+ ST isomin,
+ ST isomax)
+{
+ if (isomin >= isomax)
+ {
+ isomin = std::numeric_limits<ST>::min();
+ isomax = std::numeric_limits<ST>::max();
+ }
+
+ current_timestep = 0;
+
+ isoval = 22;
+
+#pragma unroll(3)
+ for(int i=0; i<3; i++)
+ {
+ pmin.data[i] = bounds[0].data[i];
+ pmax.data[i] = bounds[1].data[i];
+ diag.data[i] = bounds[1].data[i] - bounds[0].data[i];
+ }
+
+ cout << "diag = " << diag.data[0] << "," << diag.data[1] << "," <<
diag.data[2] << endl;
+
+ //read in the octree data from .oth, .otd files
+ {
+ if (startTimestep >= endTimestep)
+ endTimestep = startTimestep;
+
+ num_timesteps = endTimestep - startTimestep + 1;
+ steps = new Timestep[num_timesteps];
+
+ char buf[4000];
+ sprintf(buf, "%s.hdr", filebase);
+ ifstream in(buf);
+
+ cerr<< "Reading original data: " << buf << endl;
+ int nx, ny, nz;
+
+ if(!in)
+ {
+ cerr << "Error opening header: " << buf << endl;
+ exit(1);
+ }
+
+ in >> actual_dims.data[0] >> actual_dims.data[1] >>
actual_dims.data[2];
+ double x,y,z;
+ in >> x >> y >> z;
+ in >> x >> y >> z;
+ double val;
+ in >> val;
+ //datamin = (ST)val;
+ in >> val;
+ //datamax = (ST)val;
+
+ datamin = isomin;
+ datamax = isomax;
+ use_adaptive_res = false;
+
+ //determine the max_depth
+ kernel_width = kernelWidth;
+ cerr << "(Original data has dimensions " << actual_dims.data[0] <<
", " << actual_dims.data[1] << ", " << actual_dims.data[2] << ")\n";
+ double ddepth_1 = log((double)nx) / log(2.0);
+ double ddepth_2 = log((double)ny) / log(2.0);
+ double ddepth_3 = log((double)nz) / log(2.0);
+ double ddepth = MAX(MAX(ddepth_1, ddepth_2), ddepth_3);
+ max_depth = (int)(ceil(ddepth));
+
+ //make the child_bit_depth arrays
+ child_bit_depth = new int[max_depth+1];
+ inv_child_bit_depth = new float[max_depth+1];
+ for(int d=0; d<=max_depth; d++)
+ {
+ child_bit_depth[d] = 1 << (max_depth - d - 1);
+ inv_child_bit_depth[d] = 1.0f /
static_cast<float>(child_bit_depth[d]);
+ }
+ max_depth_bit = 1 << max_depth;
+ inv_max_depth_bit = 1.0f / static_cast<float>(max_depth_bit);
+
+ if(!in)
+ {
+ cerr << "Error reading header: " << buf << '\n';
+ exit(1);
+ }
+
+ diag = pmax - pmin;
+
+ //read in the single-resolution data
+ char filebase_numbered[4000];
+ for(int fts = startTimestep, ts = 0; fts <= endTimestep; fts++, ts++)
+ {
+ sprintf(filebase_numbered,"%s_%04i.raw",filebase,fts);
+
+ int din = open(filebase_numbered, O_RDONLY);
+ if (din == -1)
+ {
+ cerr << "Error opening original data file: " <<
filebase_numbered << '\n';
+ exit(1);
+ }
+
+ Array3<ST> indata;
+ indata.resize(nx, ny, nz);
+ cerr << "Reading " << filebase_numbered << "...";
+ cerr << "(" << indata.get_datasize() << " read/" << nx * ny * nz
* sizeof(ST) << "expected)...";
+ read(din, indata.get_dataptr(), indata.get_datasize());
+ cerr << "done.\n";
+ close(din);
+
+ //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 << "\nDone making from single res data\n";
+ }
+
+ write_file(filebase);
+
+ multires_level = 0;
+ current_timestep = 0;
+
+ pad_scale_factor.data[0] = (float)max_depth_bit /
(float)actual_dims.data[0];
+ pad_scale_factor.data[1] = (float)max_depth_bit /
(float)actual_dims.data[1];
+ pad_scale_factor.data[2] = (float)max_depth_bit /
(float)actual_dims.data[2];
+
+ padded_diag = diag * pad_scale_factor;
+ padded_pmax = pmin + padded_diag;
+
+ node_width = new Vector[max_depth];
+ inv_node_width = new Vector[max_depth];
+ for(int d=0; d<max_depth; d++)
+ {
+ node_width[d] = padded_diag * (1.0f / (float)(1 << (d+1)));
+ inv_node_width[d] = node_width[d].inverse();
+ }
+ }
+
+ cout << "\nDone creating OctISOVolume!" << endl;
+}
+
+OctreeVolume::~OctreeVolume()
+{
+ for(int ts=0; ts<num_timesteps; ts++)
+ {
+ for(int r=0; r<mres_levels+1; r++)
+ {
+ delete[] steps[ts].caps[r];
+
+ 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;
+ }
+ delete[] steps;
+
+ delete[] child_bit_depth;
+ delete[] inv_child_bit_depth;
+
+ delete[] node_width;
+ delete[] inv_node_width;
+}
+
+
+void OctreeVolume::make_from_single_res_data(Array3<ST>& originalData,
+ int ts, double
variance_threshold,
+ int kernelWidth, int
mresLevels,
+ ST isomin, ST isomax)
+{
+ if (isomin >= isomax)
+ {
+ isomin = std::numeric_limits<ST>::min();
+ 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;
+
+ 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];
+
+ //1. start at the bottom. Create 2^(depth-1) nodes in each dim
+ //NOTE: d=0 is root node (coarsest depth)
+ // for depth 0, 1, 2, 3, ... max_depth-1
+ // we have dims 1, 2, 4, 8, ... 2^(max_depth-1)
+ for(int d=0; d<max_depth; d++)
+ {
+ int dim = 1 << d;
+ buildnodes[d].resize(dim, dim, dim);
+ }
+
+ 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 =
originalData.lookup_safe(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(2*i, 2*j, 2*k);
+ Vec3i kmax(kmin);
+ kmin -= (kernel_width-1)/2;
+ kmax += (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 =
originalData.lookup_safe(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;
+ 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";
+ }
+
+ delete[] buildnodes;
+}
+
+bool OctreeVolume::read_file(std::string filebase)
+{
+ cout << "\nReading octree volume data from \"" << filebase << "\".oth,
.otd...";
+
+ string filename = filebase + ".oth";
+ ifstream in(filename.c_str());
+ if (!in)
+ return false;
+
+ in >> max_depth;
+ in >> actual_dims.data[0] >> actual_dims.data[1] >> actual_dims.data[2];
+ in >> kernel_width;
+ in >> datamin;
+ in >> datamax;
+ in >> num_timesteps;
+ in >> mres_levels;
+
+ child_bit_depth = new int[max_depth+1];
+ inv_child_bit_depth = new float[max_depth+1];
+ for(int d=0; d<=max_depth; d++)
+ {
+ child_bit_depth[d] = 1 << (max_depth - d - 1);
+ inv_child_bit_depth[d] = 1.0f /
static_cast<float>(child_bit_depth[d]);
+ }
+ max_depth_bit = 1 << max_depth;
+ inv_max_depth_bit = 1.0f / static_cast<float>(max_depth_bit);
+
+ pad_scale_factor.data[0] = (float)max_depth_bit /
(float)actual_dims.data[0];
+ pad_scale_factor.data[1] = (float)max_depth_bit /
(float)actual_dims.data[1];
+ pad_scale_factor.data[2] = (float)max_depth_bit /
(float)actual_dims.data[2];
+
+ padded_diag = diag * pad_scale_factor;
+ padded_pmax = pmin + padded_diag;
+
+ node_width = new Vector[max_depth];
+ inv_node_width = new Vector[max_depth];
+ for(int d=0; d<max_depth; d++)
+ {
+ node_width[d] = padded_diag * (1.0f / (float)(1 << (d+1)));
+ inv_node_width[d] = node_width[d].inverse();
+ }
+
+ current_timestep = 0;
+ multires_level = 0;
+
+ 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];
+ }
+ }
+
+ in.close();
+
+ filename = filebase + ".otd";
+
+ FILE *in2 = fopen(filename.c_str(),"rb");
+ if (!in2)
+ return false;
+
+ for(int ts = 0; ts < num_timesteps; ts++)
+ {
+ //read in the nodes
+ for(int r=0; r<mres_levels+1; r++)
+ {
+ 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);
+ 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]];
+ // read(in2,
reinterpret_cast<char*>(steps[ts].caps[r]), steps[ts].num_caps[r] *
sizeof(OctCap));
+ int num_read =
fread((char*)steps[ts].caps[r],sizeof(OctCap),steps[ts].num_caps[r],in2);
+ if (num_read != steps[ts].num_caps[r])
+ {
+ cerr << "OctreeVolume: error reading octcaps" << endl;
+ exit(0);
+ }
+ }
+ }
+ fclose(in2);
+ cout << "done!\n";
+
+ return true;
+}
+
+bool OctreeVolume::write_file(std::string filebase)
+{
+ cout << "\nWriting octree volume data to \"" << filebase << "\".oth,
.otd...";
+
+ string filename = filebase + ".oth";
+ ofstream out(filename.c_str());
+ if (!out)
+ return false;
+
+ out << max_depth << endl;
+ out << actual_dims.data[0] << " " << actual_dims.data[1] << " " <<
actual_dims.data[2] << endl;
+ 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;
+ }
+ }
+ out.close();
+
+ filename = filebase + ".otd";
+ FILE* out2 = fopen(filename.c_str(), "wb");
+ if (!out2)
+ return false;
+
+ 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-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);
+ }
+ }
+ }
+ fclose(out2);
+
+ cout << "done!\n";
+
+ return true;
+}
+
Added: trunk/Model/Primitives/OctreeVolume.h
==============================================================================
--- (empty file)
+++ trunk/Model/Primitives/OctreeVolume.h Tue May 16 16:59:57 2006
@@ -0,0 +1,443 @@
+
+/*
+* OctreeVolume.h : A multiresolution octree data
+* wrapper for scalar volumes, such as Array3.
+*
+* Written by:
+* Aaron M Knoll
+* Department of Computer Science
+* University of Utah
+* October 2005
+*
+* Copyright (C) 2005 SCI Group
+*/
+
+#ifndef _MANTA_OCTREEVOLUME_H_
+#define _MANTA_OCTREEVOLUME_H_
+
+typedef unsigned char ST;
+
+//AARONBAD - this is already defined in Abe's evil.h.
+// But it's not that evil, really...
+// maybe should be in Core/Math?
+#ifndef MIN
+#define MIN(a,b) (((a)<(b))?(a):(b))
+#endif
+#ifndef MAX
+#define MAX(a,b) (((a)>(b))?(a):(b))
+#endif
+
+#include <Core/Geometry/BBox.h>
+#include <SCIRun/Core/Containers/Array3.h>
+#include <Core/Geometry/vecdefs.h>
+#include <Core/Geometry/Vector.h>
+#include <stdlib.h>
+#include <fstream>
+#include <fcntl.h>
+#include <math.h>
+
+using namespace std;
+
+namespace Manta{
+
+ typedef unsigned char ST;
+
+ union SO
+ {
+ ST scalar;
+ unsigned char offset;
+ };
+
+ //XXX: OctCap isn't really a leaf.
+ // rather, it's a "bottom node" which by nature CANNOT have children; it
MUST have scalar values
+ struct OctCap
+ {
+ 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
+ };
+
+ //node used in build
+ struct BuildNode
+ {
+ unsigned int myStart;
+ unsigned char myIndex;
+ unsigned int children_start;
+ SO children[8];
+ 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)
+
+ class OctreeVolume
+ {
+
+ public:
+ //the maximum depth of the octree, finest resolution
+ int max_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
+ // Ex: for isosurface rendering with forward differences, use width
= 2
+ // For full central differences, use width=3
+ int kernel_width;
+
+ //number of multiresolution levels.
+ int mres_levels;
+
+ //total number of timesteps in the data
+ int num_timesteps;
+
+ int* child_bit_depth;
+ float* inv_child_bit_depth;
+ int max_depth_bit;
+ float inv_max_depth_bit;
+
+ ST datamin, datamax;
+
+ 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;
+ };
+
+ Timestep* steps;
+
+ int* refcnt;
+
+ /* THESE VARS SHOULD BE IN OCTISOVOLUME */
+ BBox bounds;
+ Vector pmin;
+ Vector pmax;
+ Vector padded_pmax; //when actual dims are not exact
powers of 2
+
+ Vector diag;
+ Vector padded_diag;
+ Vector pad_scale_factor;
+ Vec3i actual_dims; //the actual dims, padded to the
lowest power of 2 that contains the largest one.
+
+ //tells us the width of each node (or voxel) in world space
+ Vector* node_width;
+ Vector* inv_node_width;
+
+ int multires_level;
+ int current_timestep;
+ ST isoval;
+ bool use_adaptive_res;
+ /* END VARS SHOULD BE IN OCTISOVOLUME */
+
+ inline float getMaxTime() const
+ {
+ return num_timesteps - 1;
+ }
+
+ inline bool setTime(const int t)
+ {
+ if (t >= num_timesteps || t < 0)
+ return false;
+
+ current_timestep = t;
+ return true;
+ }
+
+ inline void increment_timestep()
+ {
+ if (current_timestep < num_timesteps - 1)
+ current_timestep++;
+ }
+
+ inline int get_timestep() const
+ {
+ return current_timestep;
+ }
+
+ inline void decrement_isovalue()
+ {
+ if (isoval > datamin+1)
+ isoval--;
+ }
+ inline void increment_isovalue()
+ {
+ if (isoval < datamax-1)
+ isoval++;
+ }
+
+ inline int get_isovalue() const
+ {
+ return isoval;
+ }
+
+ inline int get_multires_level() const
+ {
+ return multires_level;
+ }
+
+ inline void decrease_resolution()
+ {
+ if (multires_level < mres_levels)
+ multires_level++;
+ }
+ inline void increase_resolution()
+ {
+ if (multires_level > 0)
+ multires_level--;
+ }
+
+ inline const BBox& get_bounds() const
+ {
+ return bounds;
+ }
+
+ inline Vector& get_node_width(int d) const
+ {
+ return node_width[d];
+ }
+
+ inline Vector& get_inv_node_width(int d) const
+ {
+ return inv_node_width[d];
+ }
+
+ inline int get_child_bit_depth(int d) const
+ {
+ return child_bit_depth[d];
+ }
+
+ inline int get_max_depth() const
+ {
+ return max_depth;
+ }
+
+
+ OctreeVolume(){}
+
+ OctreeVolume(std::string filebase)
+ {
+ read_file(filebase);
+ }
+
+ /*!
+ Multiple timestep constructor
+ */
+ OctreeVolume(std::string filebase,
+ int startTimestep,
+ int endTimestep,
+ double variance_threshold,
+ int kernelWidth,
+ int mresLevels,
+ ST isomin,
+ ST isomax);
+
+ ~OctreeVolume();
+
+ inline OctNode& get_node(int depth, unsigned int index) const
+ {
+ return
steps[current_timestep].nodes[multires_level][depth][index];
+ }
+
+ inline OctCap& get_cap(unsigned int index) const
+ {
+ return steps[current_timestep].caps[multires_level][index];
+ }
+
+ inline ST operator()(Vec3i& v) const
+ {
+ return lookup_node(v, 0, 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);
+ }
+
+ 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);
+ }
+
+ inline ST lookup_node_trace(Vec3i& p, int res, int depth, unsigned
int index, unsigned int* index_trace) const
+ {
+ for(;;)
+ {
+ OctNode& node =
steps[current_timestep].nodes[res][depth][index];
+ index_trace[depth] = index;
+
+ int child_bit = child_bit_depth[depth+res];
+ 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)
+ {
+ index = node.children_start +
node.children[target_child].offset;
+ 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];
+ }
+
+ index = node.children_start +
node.children[target_child].offset;
+ depth++;
+ }
+ }
+
+ inline ST lookup_node(Vec3i& p, int res, int depth, unsigned int
index) const
+ {
+ for(;;)
+ {
+ OctNode& node =
steps[current_timestep].nodes[res][depth][index];
+
+ int child_bit = child_bit_depth[depth+res];
+ 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)
+ {
+ index = node.children_start +
node.children[target_child].offset;
+ 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];
+ }
+
+ index = node.children_start +
node.children[target_child].offset;
+ 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
+ {
+ Vec3i target_neighbor = cell + offset;
+ int child_bit;
+
+ //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];
+
+ 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]);
+ }
+ 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]);
+ }
+ 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]);
+ }
+ 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]);
+ }
+ 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]);
+ }
+ 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]);
+ }
+ 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]);
+ }
+ }
+
+ child_bit = max_depth_bit;
+ 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, 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);
+ }
+ 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);
+ }
+ 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);
+ }
+ else if (CHECK_X)
+ {
+ if ((target_neighbor.data[0] & child_bit) == (cell.data[0] &
child_bit))
+ return lookup_node(target_neighbor, res, 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);
+ }
+ 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 0;
+ }
+
+ /*!
+ The octree data building routine. AARONBAD - this is highly
unoptimized. While the resulting data is fine,
+ the build process could be made faster (and maybe parallel?) by
just unrolling loops.
+ */
+
+ void make_from_single_res_data(SCIRun::Array3<ST>& originalData,
+ int ts,
+ double variance_threshold,
+ int kernelWidth,
+ int mresLevels = 2,
+ ST isomin = 0,
+ ST isomax = 0);
+
+ bool read_file(std::string filebase);
+ bool write_file(std::string filebase);
+ };
+
+};
+
+#endif
Modified: trunk/Model/Readers/BART/parse.cc
==============================================================================
--- trunk/Model/Readers/BART/parse.cc (original)
+++ trunk/Model/Readers/BART/parse.cc Tue May 16 16:59:57 2006
@@ -11,6 +11,8 @@
#include <Model/Lights/PointLight.h>
#include <Model/AmbientLights/ConstantAmbient.h>
#include <Model/Backgrounds/ConstantBackground.h>
+#include <Core/Geometry/vecdefs.h>
+#include <Core/Geometry/varray.h>
#include <Model/Groups/Group.h>
#include <Interface/Material.h>
#include <Model/Materials/Flat.h>
@@ -33,10 +35,6 @@
#ifndef M_PI
#define M_PI 3.1415926535897932384626433
#endif
-
-typedef float Vec2f[2];
-typedef float Vec3f[3];
-typedef float Vec4f[4];
#define X 0
#define Y 1
Modified: trunk/scenes/CMakeLists.txt
==============================================================================
--- trunk/scenes/CMakeLists.txt (original)
+++ trunk/scenes/CMakeLists.txt Tue May 16 16:59:57 2006
@@ -66,6 +66,13 @@
TARGET_LINK_LIBRARIES(scene_volume ${MANTA_SCENE_LINK})
ENDIF(SCENE_VOLUME)
+# octree isosurface
+SET(SCENE_OCTISOVOL 0 CACHE BOOL "octree isosurface")
+IF(SCENE_OCTISOVOL)
+ ADD_LIBRARY(scene_octisovol octisovol.cc)
+ TARGET_LINK_LIBRARIES(scene_octisovol ${MANTA_SCENE_LINK})
+ENDIF(SCENE_OCTISOVOL)
+
SET(SCENE_OBJVIEWER 0 CACHE BOOL "Wavefront Obj file viewer.")
IF(SCENE_OBJVIEWER)
ADD_LIBRARY(scene_objviewer objviewer.cc)
Added: trunk/scenes/octisovol.cc
==============================================================================
--- (empty file)
+++ trunk/scenes/octisovol.cc Tue May 16 16:59:57 2006
@@ -0,0 +1,156 @@
+
+// Manta scene for octree isosurface volume rendering
+// Aaron Knoll
+// knolla@sci.utah.edu
+// May 2006
+
+
+#include <Core/Geometry/Vector.h>
+#include <Core/Geometry/BBox.h>
+#include <Core/Exceptions/IllegalArgument.h>
+#include <Core/Exceptions/InputError.h>
+#include <Core/Util/Args.h>
+#include <Interface/Context.h>
+#include <Interface/LightSet.h>
+#include <Interface/MantaInterface.h>
+#include <Interface/Scene.h>
+#include <Model/AmbientLights/ArcAmbient.h>
+#include <Model/Backgrounds/ConstantBackground.h>
+#include <Model/Groups/Group.h>
+#include <Model/Primitives/OctreeVolume.h>
+#include <Model/Primitives/IsosurfaceOctreeVolume.h>
+#include <Model/Lights/PointLight.h>
+#include <Model/Lights/HeadLight.h>
+#include <Model/Materials/Phong.h>
+#include <Model/Cameras/PinholeCamera.h>
+#include <Core/Geometry/AffineTransform.h>
+#include <Model/AmbientLights/ConstantAmbient.h>
+#include <Model/MiscObjects/CuttingPlane.h>
+
+#include <SCIRun/Core/Thread/Time.h>
+
+using namespace Manta;
+using namespace SCIRun;
+
+enum CuttingPlaneType { CUTTING_NONE, CUTTING_DEFAULT, CUTTING_SPECIFIED };
+
+///////////////////////////////////////////////////////////////////////////
+extern "C"
+Scene* make_scene(const ReadContext& context, const vector<string>& args) {
+
+ string filename = "";
+ string buildfrom_filename;
+ int tstart = 0;
+ int tend = 0;
+ double variance = 0.1;
+ int kernel_width = 2;
+ int mres_levels = 2;
+ double isomin = 0;
+ double isomax = 255;
+
+ Vector plane_point;
+ Vector plane_normal;
+
+ CuttingPlaneType cutting_type = CUTTING_NONE;
+
+ // Parse args.i
+ for (int i=0;i<args.size();++i) {
+ if (args[i] == "-file") {
+ if (!getStringArg(i, args, filename))
+ throw IllegalArgument("octisovol -file
<filename>", i, args);
+ }
+ else if (args[i] == "-buildfrom") {
+ if (!getStringArg(i, args, buildfrom_filename))
+ throw IllegalArgument("octisovol -buildfrom
<filename>", i, args);
+ }
+ else if (args[i] == "-tstart") {
+ if (!getIntArg(i, args, tstart))
+ throw IllegalArgument("octisovol -tstart
<tstart>", i, args);
+ }
+ else if (args[i] == "-tend") {
+ if (!getIntArg(i, args, tend))
+ throw IllegalArgument("octisovol -tend
<tend>", i, args);
+ }
+ else if (args[i] == "-variance") {
+ if (!getDoubleArg(i, args, variance))
+ throw IllegalArgument("octisovol -variance
<variance>", i, args);
+ }
+ else if (args[i] == "-kernel_width") {
+ if (!getIntArg(i, args, kernel_width))
+ throw IllegalArgument("octisovol
-kernel_width <kernel_width>", i, args);
+ }
+ else if (args[i] == "-mres_levels") {
+ if (!getIntArg(i, args, mres_levels))
+ throw IllegalArgument("octisovol -mres_levels
<mres_levels>", i, args);
+ }
+ else if (args[i] == "-isomin") {
+ if (!getDoubleArg(i, args, isomin))
+ throw IllegalArgument("octisovol -isomin
<isomin>", i, args);
+ }
+ else if (args[i] == "-isomax") {
+ if (!getDoubleArg(i, args, isomin))
+ throw IllegalArgument("octisovol -isomax
<isomax>", i, args);
+ }
+ else {
+ cerr << "Read built octree volume:" << endl;
+ cerr << "-file <filename>" << endl;
+ cerr << "Octree build:" << endl;
+ cerr << " -buildfrom <filename> " << endl;
+ cerr << "Octree build options:" << endl;
+ cerr << " -tstart <tstart>" << endl;
+ cerr << " -tend <tend>" << endl;
+ cerr << " -variance <variance>" << endl;
+ cerr << " -kernel_width <kernel_width>" << endl;
+ cerr << " -mres_levels <mres_levels>" << endl;
+ cerr << " -isomin <isomin>" << endl;
+ cerr << " -isomax <isomax>" << endl;
+ throw IllegalArgument( "octisovol", i, args );
+ }
+ }
+
+ // Create the scene.
+ Scene *scene = new Scene();
+ Group *group = new Group();
+
+ Manta::BBox bounds;
+ OctreeVolume* ov;
+
+ if (filename != "")
+ {
+ ov = new OctreeVolume(filename);
+ }
+ else if (buildfrom_filename != "")
+ {
+ ov = new OctreeVolume(buildfrom_filename, tstart, tend, variance,
kernel_width, mres_levels, isomin, isomax);
+ }
+ else
+ {
+ throw InputError( "octisovol: No volume source specified. Use either
-file or -buildfrom.");
+ }
+
+ Material* mat1 = new Phong(Color(RGBColor(0.05f, 0.3f, 0.6f)),
Color(RGBColor(1.f, 1.f, 1.f)), 50);
+ IsosurfaceOctreeVolume* iov = new IsosurfaceOctreeVolume(ov, mat1);
+ bounds = iov->getBounds();
+ group->add(iov);
+
+ LightSet *lights = new LightSet();
+
+ lights->add( new HeadLight( 2.0, Color(RGB(1.0,1.0,1.0)) ));
+ lights->setAmbientLight( new ConstantAmbient(
Color(RGBColor(0.2,0.2,0.2) ) ));
+ scene->setLights(lights);
+
+ // Add the tree to the scene.
+ scene->setObject( group );
+
+ // Set other important scene properties.
+ PinholeCamera *camera = new PinholeCamera(
+ Vector(bounds[0].data[0]/2, bounds[0].data[1]/2,
bounds[0].data[2]*2),
+ Vector(bounds[1].data[0]/2, bounds[1].data[1]/2, 0),
+ Vector(0,1,0), 50.0 );
+
+ // Background.
+ scene->setBackground( new ConstantBackground( Color(RGB(0.8, 0.8,
0.8)) ) );
+
+ return scene;
+}
+
- [MANTA] r1063 - in trunk: . Core Core/Geometry Core/Math Engine/PixelSamplers Model/Groups Model/Intersections Model/Primitives Model/Readers/BART scenes, knolla, 05/16/2006
Archive powered by MHonArc 2.6.16.
