-
sparker
thiago
Manta Interactive Ray Tracer Development Mailing List
Text archives Help
[MANTA] r1079 - in trunk: . Model/Intersections Model/MiscObjects Model/Primitives scenes
Chronological Thread
- From: knolla@sci.utah.edu
- To: manta@sci.utah.edu
- Subject: [MANTA] r1079 - in trunk: . Model/Intersections Model/MiscObjects Model/Primitives scenes
- Date: Tue, 23 May 2006 11:00:01 -0600 (MDT)
Author: knolla
Date: Tue May 23 10:59:53 2006
New Revision: 1079
Added:
trunk/Model/MiscObjects/BrickArray3.h
trunk/Model/Primitives/IsosurfaceGridVolume.cc
trunk/Model/Primitives/IsosurfaceGridVolume.h
trunk/scenes/gridisovol.cc
Modified:
trunk/CMakeLists.txt
trunk/Model/Intersections/IsosurfaceImplicit.cc
trunk/Model/Intersections/IsosurfaceImplicit.h
trunk/Model/MiscObjects/CMakeLists.txt
trunk/Model/Primitives/CMakeLists.txt
trunk/Model/Primitives/IsosurfaceOctreeVolume.cc
trunk/Model/Primitives/OctreeVolume.cc
trunk/Model/Primitives/OctreeVolume.h
trunk/scenes/CMakeLists.txt
trunk/scenes/octisovol.cc
Log:
Implemented grid isosurface volumes, a la Parker et al 98
Modified: trunk/CMakeLists.txt
==============================================================================
--- trunk/CMakeLists.txt (original)
+++ trunk/CMakeLists.txt Tue May 23 10:59:53 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" 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" STRING "Optimized Flags")
+ FIRST_TIME_SET(CMAKE_CXX_FLAGS_RELEASE "-DSCI_ASSERTION_LEVEL=0 -fgcse-sm
-funroll-loops -fstrict-aliasing -fsched-interblock -ffast-math
-freorder-blocks -fpermissive -march=prescott -mtune=prescott -msse -msse2
-msse3 -mfpmath=sse" STRING "Optimized Flags")
ENDIF (APPLE)
##################################################################
Modified: trunk/Model/Intersections/IsosurfaceImplicit.cc
==============================================================================
--- trunk/Model/Intersections/IsosurfaceImplicit.cc (original)
+++ trunk/Model/Intersections/IsosurfaceImplicit.cc Tue May 23 10:59:53
2006
@@ -5,7 +5,7 @@
using namespace Manta;
//From Steven Parker's 1997 RTRT isosurface intersection
-bool IsosurfaceImplicit::single_intersect(RayPacket& rays, int which_one,
+bool IsosurfaceImplicit::single_intersect(const Vector& orig, const Vector&
dir,
const Vector& pmin, const Vector& pmax, float
rho[2][2][2],
float iso, float tmin, float tmax, float& t)
{
@@ -18,9 +18,6 @@
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];
Modified: trunk/Model/Intersections/IsosurfaceImplicit.h
==============================================================================
--- trunk/Model/Intersections/IsosurfaceImplicit.h (original)
+++ trunk/Model/Intersections/IsosurfaceImplicit.h Tue May 23 10:59:53
2006
@@ -9,7 +9,7 @@
{
struct IsosurfaceImplicit
{
- static bool single_intersect(RayPacket& rays, int which_one,
+ static bool single_intersect(const Vector& orig, const Vector& dir,
const Vector& pmin, const Vector& pmax, float
rho[2][2][2],
float iso, float tmin, float tmax, float& t);
Added: trunk/Model/MiscObjects/BrickArray3.h
==============================================================================
--- (empty file)
+++ trunk/Model/MiscObjects/BrickArray3.h Tue May 23 10:59:53 2006
@@ -0,0 +1,263 @@
+
+/*
+ * BrickArray3.h: Interface to dynamic 3D array class
+ *
+ * Written by:
+ * Steven G. Parker
+ * Department of Computer Science
+ * University of Utah
+ * March 1994
+ *
+ * Copyright (C) 1994 SCI Group
+ *
+ * Ported into Manta by Aaron Knoll.
+ * XXX - this should ideally be in Core or SCIRun/Core/Containers. Move if
necessary.
+ */
+
+#ifndef _MANTA_BRICKARRAY3_H_
+#define _MANTA_BRICKARRAY3_H_
+
+#include <math.h>
+
+#define CACHESIZE 64
+#define PAGESIZE 4096
+
+template<class T>
+class BrickArray3 {
+ //this is the data
+ T* objs;
+
+ char* data;
+
+ //use this to make sure the data is deleted only after all the referers
are done with it
+ int* refcnt;
+
+public:
+ //these are index arrays, the data is reorganized so that the memory for
elem[x,y,z] is close to elem[x+-q,y+-q,z+-q]
+ int* idx1;
+ int* idx2;
+ int* idx3;
+
+ int dm1;
+ int dm2;
+ int dm3;
+ int totaldm1;
+ int totaldm2;
+ int totaldm3;
+
+ int L1, L2;
+ void allocate();
+ BrickArray3<T>& operator=(const BrickArray3&);
+public:
+ BrickArray3();
+ BrickArray3(int, int, int);
+ ~BrickArray3();
+
+ //when you refer to an element, use the index arrays to find it's
optimized memory location first
+ inline T& operator()(int d1, int d2, int d3) const
+ {
+ return objs[idx1[d1]+idx2[d2]+idx3[d3]];
+ }
+
+ //accessors
+ inline int dim1() const {return dm1;}
+ inline int dim2() const {return dm2;}
+ inline int dim3() const {return dm3;}
+
+ void resize(int, int, int);
+
+ void initialize(const T&);
+
+ //allow the user to manipulate the data directly
+ inline T* get_dataptr() {return objs;}
+
+ //how big is it?
+ inline unsigned long get_datasize()
+ {
+ return totaldm1*totaldm2*totaldm3*sizeof(T);
+ }
+
+ //allow more than one BrickArray to share the same data
+ void share(const BrickArray3<T>& copy);
+};
+
+template<class T>
+BrickArray3<T>::BrickArray3()
+{
+ objs=0;
+ data=0;
+ dm1=dm2=dm3=0;
+ totaldm1=totaldm2=totaldm3=0;
+ idx1=idx2=idx3=0;
+ L1=L2=0;
+}
+
+template<class T>
+void BrickArray3<T>::allocate()
+{
+ if(dm1==0 || dm2==0 || dm3==0){
+ objs=0;
+ data=0;
+ refcnt=0;
+ dm1=dm2=dm3=0;
+ totaldm1=totaldm2=totaldm3=0;
+ idx1=idx2=idx3=0;
+ L1=L2=0;
+ return;
+ }
+
+ //CACHESIZE byte cache line
+ //CACHESIZE = sizeof(T) * L1^3
+ //L1 = 3rdrt(CACHESIZE/sizeof(T))
+ L1=(int)(pow((CACHESIZE*1.0)/(double)(sizeof(T)), 1./3.)+.1);
+
+ // 16K page size
+ //bricksize = (sizeof(T)*L1)^3
+ //16K = bricksize * L2^3
+ //L2 = 3rdrt(16K/bricksize)
+ L2=(int)(pow((PAGESIZE*1.0)/(sizeof(T)*L1*L1*L1), 1./3.)+.1);
+
+ //pad to make an even number of bricks in each dimension
+ int totalx=(dm1+L2*L1-1)/(L2*L1);
+ int totaly=(dm2+L2*L1-1)/(L2*L1);
+ int totalz=(dm3+L2*L1-1)/(L2*L1);
+
+ //create the x address index table
+ idx1=new int[dm1];
+ for(int x=0;x<dm1;x++){
+ int m1x=x%L1;
+ int xx=x/L1;
+ int m2x=xx%L2;
+ int m3x=xx/L2;
+ idx1[x]=
+ m3x*totaly*totalz*L2*L2*L2*L1*L1*L1+
+ m2x*L2*L2*L1*L1*L1+
+ m1x*L1*L1;
+ }
+ //create the y address index table
+ idx2=new int[dm2];
+ for(int y=0;y<dm2;y++){
+ int m1y=y%L1;
+ int yy=y/L1;
+ int m2y=yy%L2;
+ int m3y=yy/L2;
+ idx2[y]=
+ m3y*totalz*L2*L2*L2*L1*L1*L1+
+ m2y*L2*L1*L1*L1+
+ m1y*L1;
+ }
+ //create the z address index table
+ idx3=new int[dm3];
+ for(int z=0;z<dm3;z++){
+ int m1z=z%L1;
+ int zz=z/L1;
+ int m2z=zz%L2;
+ int m3z=zz/L2;
+ idx3[z]=
+ m3z*L2*L2*L2*L1*L1*L1+
+ m2z*L1*L1*L1+
+ m1z;
+ }
+
+ //these are the padded dimensions
+ totaldm1=totalx*L2*L1;
+ totaldm2=totaly*L2*L1;
+ totaldm3=totalz*L2*L1;
+ int totalsize=totaldm1*totaldm2*totaldm3;
+
+ //create the space to hold the data
+ //objs=new T[totalsize];
+ //the padding by CACHESIZE+PAGESIZE followed by
+ //the %CACHESIZE aligns start on a double word boundary.
+ data=new char[totalsize*sizeof(T)+CACHESIZE+4096]; //original used
4096, not 16384, why?
+ unsigned long off=(unsigned long)data%CACHESIZE;
+
+ if (off)
+ objs=(T*)(data+CACHESIZE-off);
+ else
+ objs=(T*)data;
+
+ //there is at least one referer to the data
+ refcnt = new int;
+ *refcnt=1;
+}
+
+template<class T>
+ void BrickArray3<T>::resize(int d1, int d2, int d3)
+{
+ if(objs && dm1==d2 && dm2==d2 && dm3==d3)return;
+ dm1=d1;
+ dm2=d2;
+ dm3=d3;
+ if (objs) {
+ (*refcnt--);
+ if (*refcnt == 0) {
+ delete[] data;
+ delete[] idx1;
+ delete[] idx2;
+ delete[] idx3;
+ delete refcnt;
+ }
+ }
+ allocate();
+}
+
+template<class T>
+ BrickArray3<T>::BrickArray3(int dm1, int dm2, int dm3)
+ : dm1(dm1), dm2(dm2),dm3(dm3)
+{
+ allocate();
+}
+
+template<class T>
+ BrickArray3<T>::~BrickArray3()
+{
+ if(objs){
+ if(*refcnt == 0){
+ delete[] data;
+ delete[] idx1;
+ delete[] idx2;
+ delete[] idx3;
+ delete refcnt;
+ }
+ }
+}
+
+template<class T>
+ void BrickArray3<T>::initialize(const T& t)
+{
+ int n=dm1*dm2*dm3;
+ for(int i=0;i<n;i++)
+ objs[i]=t;
+}
+
+template<class T>
+ void BrickArray3<T>::share(const BrickArray3<T>& copy)
+{
+ if(objs){
+ (*refcnt--);
+ if(*refcnt == 0){
+ delete[] data;
+ delete[] idx1;
+ delete[] idx2;
+ delete[] idx3;
+ }
+ }
+ objs=copy.objs;
+ data=copy.data;
+ refcnt=copy.refcnt;
+ dm1=copy.dm1;
+ dm2=copy.dm2;
+ dm3=copy.dm3;
+ totaldm1=copy.totaldm1;
+ totaldm2=copy.totaldm2;
+ totaldm3=copy.totaldm3;
+ idx1=copy.idx1;
+ idx2=copy.idx2;
+ idx3=copy.idx3;
+ L1=copy.L1;
+ L2=copy.L2;
+ (*refcnt)++;
+}
+
+#endif
Modified: trunk/Model/MiscObjects/CMakeLists.txt
==============================================================================
--- trunk/Model/MiscObjects/CMakeLists.txt (original)
+++ trunk/Model/MiscObjects/CMakeLists.txt Tue May 23 10:59:53 2006
@@ -1,4 +1,5 @@
SET (Manta_MiscObjects_SRCS
+ MiscObjects/BrickArray3.h
MiscObjects/CuttingPlane.h
MiscObjects/CuttingPlane.cc
MiscObjects/Difference.h
Modified: trunk/Model/Primitives/CMakeLists.txt
==============================================================================
--- trunk/Model/Primitives/CMakeLists.txt (original)
+++ trunk/Model/Primitives/CMakeLists.txt Tue May 23 10:59:53 2006
@@ -18,6 +18,8 @@
Primitives/Hemisphere.h
Primitives/OctreeVolume.cc
Primitives/OctreeVolume.h
+ Primitives/IsosurfaceGridVolume.cc
+ Primitives/IsosurfaceGridVolume.h
Primitives/IsosurfaceOctreeVolume.cc
Primitives/IsosurfaceOctreeVolume.h
Primitives/Parallelogram.cc
Added: trunk/Model/Primitives/IsosurfaceGridVolume.cc
==============================================================================
--- (empty file)
+++ trunk/Model/Primitives/IsosurfaceGridVolume.cc Tue May 23 10:59:53
2006
@@ -0,0 +1,672 @@
+
+#include <Model/Primitives/IsosurfaceGridVolume.h>
+#include <Model/Primitives/OctreeVolume.h>
+#include <Interface/RayPacket.h>
+#include <Model/Intersections/IsosurfaceImplicit.h>
+
+using namespace std;
+using namespace Manta;
+using namespace SCIRun;
+
+IsosurfaceGridVolume::IsosurfaceGridVolume(char* filebase, int _mc_depth,
Material* _matl)
+: PrimitiveCommon(_matl)
+{
+ isoval = 22;
+ mc_depth = _mc_depth;
+
+ this->filebase=strdup(filebase);
+ if (mc_depth<=0)
+ mc_depth=1;
+ char filename[2048];
+
+ sprintf(filename, "%s.hdr", filebase);
+ ifstream in(filename);
+
+ if(!in){
+ cerr << "Error opening header: " << filename << endl;
+ exit(1);
+ }
+ in >> nx >> ny >> nz;
+ double x,y,z;
+ in >> x >> y >> z;
+ bounds[0]=Vector(x,y,z);
+ in >> x >> y >> z;
+ bounds[1]=Vector(x,y,z);
+ double val; //do it this way, if ST is char, this will just read the
first char's instead of the whole number
+ in >> val;
+ datamin = (ST)val;
+ in >> val;
+ datamax = (ST)val;
+
+ if(!in)
+ {
+ cerr << "Error reading header: " << filename << '\n';
+ exit(1);
+ }
+ datadiag=bounds[1] - bounds[0];
+ sdiag=datadiag/Vector(nx-1,ny-1,nz-1);
+
+ blockdata.resize(nx, ny, nz);
+
+ sprintf(filename, "%s.raw", filebase);
+ FILE* din = fopen(filename, "rb");
+ if (!din)
+ {
+ cerr << "Error opening original data file: " << filename << '\n';
+ exit(1);
+ }
+
+ Array3<ST> indata;
+ indata.resize(nx, ny, nz);
+ cerr << "Reading " << filename << "...";
+ fread((char*)(&indata.get_dataptr()[0][0][0]), 1, indata.get_datasize(),
din);
+ cerr << "done.\n";
+ fclose(din);
+
+ cerr << "Copying original data to BrickArray3...";
+ for(int i=0; i<nx; i++)
+ for(int j=0; j<ny; j++)
+ for(int k=0; k<nz; k++)
+ blockdata(i,j,k) = indata(i,j,k);
+ cerr << "done!";
+
+ xsize=new int[mc_depth];
+ ysize=new int[mc_depth];
+ zsize=new int[mc_depth];
+ int tx=nx-1;
+ int ty=ny-1;
+ int tz=nz-1;
+ for(int i=mc_depth-1;i>=0;i--)
+ {
+ int nx=(int)(pow(tx, 1./(i+1))+.9);
+ tx=(tx+nx-1)/nx;
+ xsize[mc_depth-i-1]=nx;
+ int ny=(int)(pow(ty, 1./(i+1))+.9);
+ ty=(ty+ny-1)/ny;
+ ysize[mc_depth-i-1]=ny;
+ int nz=(int)(pow(tz, 1./(i+1))+.9);
+ tz=(tz+nz-1)/nz;
+ zsize[mc_depth-i-1]=nz;
+ }
+ ixsize=new double[mc_depth];
+ iysize=new double[mc_depth];
+ izsize=new double[mc_depth];
+ cerr << "Calculating mc_depths...\n";
+ for(int i=0;i<mc_depth;i++){
+ cerr << "xsize=" << xsize[i] << ", ysize=" << ysize[i] << ", zsize="
<< zsize[i] << '\n';
+ ixsize[i]=1./xsize[i];
+ iysize[i]=1./ysize[i];
+ izsize[i]=1./zsize[i];
+ }
+ cerr << "X: ";
+ tx=1;
+ for(int i=mc_depth-1;i>=0;i--){
+ cerr << xsize[i] << ' ';
+ tx*=xsize[i];
+ }
+ cerr << "(" << tx << ")\n";
+ if(tx<nx-1)
+ {
+ cerr << "TX TOO SMALL!\n";
+ exit(1);
+ }
+ cerr << "Y: ";
+ ty=1;
+ for(int i=mc_depth-1;i>=0;i--)
+ {
+ cerr << ysize[i] << ' ';
+ ty*=ysize[i];
+ }
+ cerr << "(" << ty << ")\n";
+ if(ty<ny-1)
+ {
+ cerr << "TY TOO SMALL!\n";
+ exit(1);
+ }
+ cerr << "Z: ";
+ tz=1;
+ for(int i=mc_depth-1;i>=0;i--){
+ cerr << zsize[i] << ' ';
+ tz*=zsize[i];
+ }
+ cerr << "(" << tz << ")\n";
+ if(tz<nz-1){
+ cerr << "TZ TOO SMALL!\n";
+ exit(1);
+ }
+ hierdiag=datadiag*Vector(tx,ty,tz)/Vector(nx-1,ny-1,nz-1);
+ ihierdiag=Vector(1.,1.,1.)/hierdiag;
+
+ if(mc_depth==1){
+ macrocells=0;
+ }
+ else
+ {
+ macrocells=new BrickArray3<VMCell>[mc_depth+1];
+ int xs=1;
+ int ys=1;
+ int zs=1;
+ for(int d=mc_depth-1;d>=1;d--){
+ xs*=xsize[d];
+ ys*=ysize[d];
+ zs*=zsize[d];
+ macrocells[d].resize(xs, ys, zs);
+ cerr << "Depth " << d << ": " << xs << "x" << ys << "x" << zs <<
'\n';
+ }
+ cerr << "Creating macrocells...";
+
+ VMCell top;
+ calc_mcell(mc_depth-1, 0, 0, 0, top);
+ cerr << "Min: " << (int)top.min << ", Max: " << (int)top.max << endl;
+ cerr << "\ndone!\n";
+ }
+
+#ifdef TEST_OCTDATA
+ //build the octree data from the single-res data.
+ double threshold = 0.0;
+ int max_depth = -1;
+ cerr << "Building multi-resolution octree data. (var thresh = "<<
threshold << ", max depth = " << max_depth <<")...";
+ octdata.make_from_single_res_data(indata, threshold, 2, 0);
+ octdata.write_file(filebase);
+#endif
+
+}
+
+IsosurfaceGridVolume::~IsosurfaceGridVolume()
+{
+}
+
+void IsosurfaceGridVolume::preprocess( PreprocessContext const &context )
+{
+ PrimitiveCommon::preprocess(context);
+}
+
+void IsosurfaceGridVolume::computeBounds( PreprocessContext const &context,
BBox &box ) const
+{
+ box = bounds;
+}
+
+void IsosurfaceGridVolume::computeNormal(const RenderContext& context,
RayPacket& rays) const
+{
+}
+
+BBox IsosurfaceGridVolume::getBounds() const
+{
+ return bounds;
+}
+
+void IsosurfaceGridVolume::intersect(RenderContext const &context, RayPacket
&packet) const
+{
+ for ( int i = packet.rayBegin; i < packet.rayEnd; i++ )
+ single_intersect(packet, i);
+}
+
+
+void IsosurfaceGridVolume::single_intersect(RayPacket& rays, int which_one)
const
+{
+ int start_depth = this->mc_depth-1;
+ //find tenter, texit, penter using ray-AABB intersection
+ Vector dir = rays.getDirection(which_one);
+ Vector inv_dir = dir.inverse();
+ Vector orig = rays.getOrigin(which_one);
+ Vector min = bounds[0];
+ Vector max = bounds[1];
+
+ int dix_dx;
+ int ddx;
+ int diy_dy;
+ int ddy;
+ int diz_dz;
+ int ddz;
+
+ float tenter, texit;
+ if(dir.x() > 0)
+ {
+ tenter = inv_dir.x() * (min.x()-orig.x());
+ texit = inv_dir.x() * (max.x()-orig.x());
+ dix_dx=1;
+ ddx=1;
+ }
+ else
+ {
+ tenter = inv_dir.x() * (max.x()-orig.x());
+ texit = inv_dir.x() * (min.x()-orig.x());
+ dix_dx=-1;
+ ddx=0;
+ }
+
+ float y0, y1;
+ if(dir.y() > 0)
+ {
+ y0 = inv_dir.y() * (min.y()-orig.y());
+ y1 = inv_dir.y() * (max.y()-orig.y());
+ diy_dy=1;
+ ddy=1;
+ }
+ else
+ {
+ y0 = inv_dir.y() * (max.y()-orig.y());
+ y1 = inv_dir.y() * (min.y()-orig.y());
+ diy_dy=-1;
+ ddy=0;
+ }
+ tenter = MAX(tenter, y0);
+ texit = MIN(texit, y1);
+ if (tenter > texit)
+ return;
+
+ float z0, z1;
+ if(dir.z() > 0)
+ {
+ z0 = inv_dir.z() * (min.z()-orig.z());
+ z1 = inv_dir.z() * (max.z()-orig.z());
+ diz_dz=1;
+ ddz=1;
+ }
+ else
+ {
+ z0 = inv_dir.z() * (max.z()-orig.z());
+ z1 = inv_dir.z() * (min.z()-orig.z());
+ diz_dz=-1;
+ ddz=0;
+ }
+ tenter = MAX(tenter, z0);
+ texit = MIN(texit, z1);
+ if (tenter > texit)
+ return;
+
+ if (texit < 0.)
+ return;
+
+ tenter = MAX(0., tenter);
+ Vector p = orig + dir*tenter;
+ Vector s((p-min)*ihierdiag);
+ int cx=xsize[start_depth];
+ int cy=ysize[start_depth];
+ int cz=zsize[start_depth];
+ int ix=(int)(s.x()*cx);
+ int iy=(int)(s.y()*cy);
+ int iz=(int)(s.z()*cz);
+ if(ix>=cx)
+ ix--;
+ if(iy>=cy)
+ iy--;
+ if(iz>=cz)
+ iz--;
+ if(ix<0)
+ ix++;
+ if(iy<0)
+ iy++;
+ if(iz<0)
+ iz++;
+
+
+ float next_x, next_y, next_z;
+ float dtdx, dtdy, dtdz;
+ float icx=ixsize[start_depth];
+ float x=min.x()+hierdiag.x()*float(ix+ddx)*icx;
+ next_x=(x-orig.x())*inv_dir.x();
+ dtdx=dix_dx*hierdiag.x()*icx*inv_dir.x();
+ float icy=iysize[start_depth];
+ float y=min.y()+hierdiag.y()*float(iy+ddy)*icy;
+ next_y=(y-orig.y())*inv_dir.y();
+ dtdy=diy_dy*hierdiag.y()*icy*inv_dir.y();
+ float icz=izsize[start_depth];
+ float z=min.z()+hierdiag.z()*float(iz+ddz)*icz;
+ next_z=(z-orig.z())*inv_dir.z();
+ dtdz=diz_dz*hierdiag.z()*icz*inv_dir.z();
+
+ Vector cellsize(cx,cy,cz);
+ Vector cellcorner((orig-min)*ihierdiag*cellsize);
+ Vector celldir(dir*ihierdiag*cellsize);
+
+ isect(rays, which_one, orig, dir, inv_dir,
+ start_depth, tenter, dtdx, dtdy, dtdz, next_x, next_y, next_z,
+ ix, iy, iz, dix_dx, diy_dy, diz_dz,
+ 0, 0, 0, cellcorner, celldir);
+}
+
+void IsosurfaceGridVolume::isect(RayPacket& rays, int which_one,
+ const Vector& orig, const Vector& dir, const Vector& inv_dir,
+ int depth, float t,
+ float dtdx, float dtdy, float dtdz,
+ float next_x, float next_y, float next_z,
+ int ix, int iy, int iz,
+ int dix_dx, int diy_dy, int diz_dz,
+ int startx, int starty, int startz,
+ const Vector& cellcorner, const Vector& celldir) const
+{
+ int cx=xsize[depth];
+ int cy=ysize[depth];
+ int cz=zsize[depth];
+ if(depth==0)
+ {
+ for(;;)
+ {
+ int gx=startx+ix;
+ int gy=starty+iy;
+ int gz=startz+iz;
+ if(gx<nx-1 && gy<ny-1 && gz<nz-1)
+ {
+ //alog( << "Doing cell: " << gx << ", " << gy << ", " << gz
+ //<< " (" << startx << "+" << ix << ", " << starty << "+" <<
iy << ", " << startz << "+" << iz << ")\n");
+ ST rhos[8];
+#ifdef USE_OCTREE_DATA
+#ifdef USE_OTD_NEIGHBOR_FIND
+ //neighbor find
+ int index_trace[octdata.max_depth];
+ ST min_rho, max_rho, this_rho;
+ Vec3i child_cell(gx, gy, gz);
+ rhos[0]=octdata.lookup_node_trace(child_cell, 0, 0, 0,
index_trace);
+ int target_child = ((gx & 1) << 2) | ((gy & 1) << 1) | (gz &
1);
+ OctCap& cap =
octdata.caps[0][index_trace[octdata.max_depth-1]];
+ //0,0,1
+ Vec3i offset(0,0,1);
+ if (target_child & 1)
+ rhos[1] = octdata.lookup_neighbor<0,0,1>(child_cell,
offset, 0, depth, index_trace);
+ else
+ rhos[1] = cap.values[target_child | 1];
+ //0,1,1
+ offset.y = 1;
+ if (target_child & 3)
+ rhos[3] = octdata.lookup_neighbor<0,1,1>(child_cell,
offset, 0, depth, index_trace);
+ else
+ rhos[3] = cap.values[target_child | 3];
+ //1,1,1
+ offset.x = 1;
+ if (target_child & 7)
+ rhos[7] = octdata.lookup_neighbor<1,1,1>(child_cell,
offset, 0, depth, index_trace);
+ else
+ rhos[7] = cap.values[target_child | 7];
+ //1,1,0
+ offset.z = 0;
+ if (target_child & 6)
+ rhos[6] = octdata.lookup_neighbor<1,1,0>(child_cell,
offset, 0, depth, index_trace);
+ else
+ rhos[6] = cap.values[target_child | 6];
+ //1,0,0
+ offset.y = 0;
+ if (target_child & 4)
+ rhos[4] = octdata.lookup_neighbor<1,0,0>(child_cell,
offset, 0, depth, index_trace);
+ else
+ rhos[4] = cap.values[target_child | 4];
+ //1,0,1
+ offset.z = 1;
+ if (target_child & 5)
+ rhos[5] = octdata.lookup_neighbor<1,0,1>(child_cell,
offset, 0, depth, index_trace);
+ else
+ rhos[5] = cap.values[target_child | 5];
+ //0,1,0
+ offset.x = 0;
+ offset.y = 1;
+ offset.z = 0;
+ if (target_child & 2)
+ rhos[2] = octdata.lookup_neighbor<0,1,0>(child_cell,
offset, 0, depth, index_trace);
+ else
+ rhos[2] = cap.values[target_child | 2];
+#else
+ //pure point location
+ rhos[0]=octdata(gx, gy, gz);
+ rhos[1]=octdata(gx, gy, gz+1);
+ rhos[2]=octdata(gx, gy+1, gz);
+ rhos[3]=octdata(gx, gy+1, gz+1);
+ rhos[4]=octdata(gx+1, gy, gz);
+ rhos[5]=octdata(gx+1, gy, gz+1);
+ rhos[6]=octdata(gx+1, gy+1, gz);
+ rhos[7]=octdata(gx+1, gy+1, gz+1);
+#endif
+#else //lookup using the grid
+ rhos[0]=blockdata(gx, gy, gz);
+ rhos[1]=blockdata(gx, gy, gz+1);
+ rhos[2]=blockdata(gx, gy+1, gz);
+ rhos[3]=blockdata(gx, gy+1, gz+1);
+ rhos[4]=blockdata(gx+1, gy, gz);
+ rhos[5]=blockdata(gx+1, gy, gz+1);
+ rhos[6]=blockdata(gx+1, gy+1, gz);
+ rhos[7]=blockdata(gx+1, gy+1, gz+1);
+#endif
+ ST min=rhos[0];
+ ST max=rhos[0];
+ for(int i=1;i<8;i++)
+ {
+ if(rhos[i]<min)
+ min=rhos[i];
+ if(rhos[i]>max)
+ max=rhos[i];
+ }
+
+ if(min < isoval && max>isoval)
+ {
+ //alog(<<"might have a hit" << endl);
+ float hit_t;
+ Vector p0(bounds[0]+sdiag*Vector(gx,gy,gz));
+ Vector p1(p0+sdiag);
+ float tmax=next_x;
+ if(next_y<tmax)
+ tmax=next_y;
+ if(next_z<tmax)
+ tmax=next_z;
+ float rho[2][2][2];
+ rho[0][0][0]=rhos[0];
+ rho[0][0][1]=rhos[1];
+ rho[0][1][0]=rhos[2];
+ rho[0][1][1]=rhos[3];
+ rho[1][0][0]=rhos[4];
+ rho[1][0][1]=rhos[5];
+ rho[1][1][0]=rhos[6];
+ rho[1][1][1]=rhos[7];
+
+ if (IsosurfaceImplicit::single_intersect(orig, dir, p0,
p1, rho,
+ isoval, t, tmax, 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, p0,
p1, phit, rho);
+ normal.normalize();
+ rays.setNormal(which_one, normal);
+ }
+ }
+ }
+ }
+
+ if(next_x < next_y && next_x < next_z)
+ {
+ // Step in x...
+ t=next_x;
+ next_x+=dtdx;
+ ix+=dix_dx;
+ if(ix<0 || ix>=cx)
+ break;
+ }
+ else if(next_y < next_z){
+ t=next_y;
+ next_y+=dtdy;
+ iy+=diy_dy;
+ if(iy<0 || iy>=cy)
+ break;
+ }
+ else {
+ t=next_z;
+ next_z+=dtdz;
+ iz+=diz_dz;
+ if(iz<0 || iz>=cz)
+ break;
+ }
+ }
+ }
+ else
+ {
+ BrickArray3<VMCell>& mcells=macrocells[depth];
+ for(;;)
+ {
+ int gx=startx+ix;
+ int gy=starty+iy;
+ int gz=startz+iz;
+ VMCell& mcell=mcells(gx,gy,gz);
+ //alog( << "doing macrocell: " << gx << ", " << gy << ", " << gz
<< ": " << mcell.min << ", " << mcell.max << endl);
+ if(mcell.max>isoval && mcell.min<isoval)
+ {
+ // Do this cell...
+ int new_cx=xsize[depth-1];
+ int new_cy=ysize[depth-1];
+ int new_cz=zsize[depth-1];
+ int
new_ix=(int)((cellcorner.data[0]+t*celldir.data[0]-ix)*new_cx);
+ int
new_iy=(int)((cellcorner.data[1]+t*celldir.data[1]-iy)*new_cy);
+ int
new_iz=(int)((cellcorner.data[2]+t*celldir.data[2]-iz)*new_cz);
+ //alog( << "new: " <<
(cellcorner.x()+t*celldir.x()-ix)*new_cx
+ //<< " " << (cellcorner.y()+t*celldir.y()-iy)*new_cy
+ //<< " " << (cellcorner.z()+t*celldir.z()-iz)*new_cz
+ //<< '\n');
+ if(new_ix<0)
+ new_ix=0;
+ else if(new_ix>=new_cx)
+ new_ix=new_cx-1;
+ if(new_iy<0)
+ new_iy=0;
+ else if(new_iy>=new_cy)
+ new_iy=new_cy-1;
+ if(new_iz<0)
+ new_iz=0;
+ else if(new_iz>=new_cz)
+ new_iz=new_cz-1;
+
+ float new_dtdx=dtdx*ixsize[depth-1];
+ float new_dtdy=dtdy*iysize[depth-1];
+ float new_dtdz=dtdz*izsize[depth-1];
+ const Vector dir = rays.getDirection(which_one);
+ float new_next_x;
+ if(dir.data[0] > 0){
+ new_next_x=next_x-dtdx+new_dtdx*(new_ix+1);
+ } else {
+ new_next_x=next_x-new_ix*new_dtdx;
+ }
+ float new_next_y;
+ if(dir.data[1] > 0){
+ new_next_y=next_y-dtdy+new_dtdy*(new_iy+1);
+ } else {
+ new_next_y=next_y-new_iy*new_dtdy;
+ }
+ float new_next_z;
+ if(dir.data[2] > 0){
+ new_next_z=next_z-dtdz+new_dtdz*(new_iz+1);
+ } else {
+ new_next_z=next_z-new_iz*new_dtdz;
+ }
+ int new_startx=gx*new_cx;
+ int new_starty=gy*new_cy;
+ int new_startz=gz*new_cz;
+ //alog( << "startz=" << startz << '\n');
+ //alog( << "iz=" << iz << '\n');
+ //alog( <
+ Vector cellsize(new_cx, new_cy, new_cz);
+ isect(rays, which_one,
+ orig, dir, inv_dir,
+ depth-1, t,
+ new_dtdx, new_dtdy, new_dtdz,
+ new_next_x, new_next_y, new_next_z,
+ new_ix, new_iy, new_iz,
+ dix_dx, diy_dy, diz_dz,
+ new_startx, new_starty, new_startz,
+ (cellcorner-Vector(ix, iy, iz))*cellsize,
celldir*cellsize);
+ }
+
+ if(next_x < next_y && next_x < next_z)
+ {
+ // Step in x...
+ t=next_x;
+ next_x+=dtdx;
+ ix+=dix_dx;
+ if(ix<0 || ix>=cx)
+ break;
+ }
+ else if(next_y < next_z)
+ {
+ t=next_y;
+ next_y+=dtdy;
+ iy+=diy_dy;
+ if(iy<0 || iy>=cy)
+ break;
+ }
+ else
+ {
+ t=next_z;
+ next_z+=dtdz;
+ iz+=diz_dz;
+ if(iz<0 || iz>=cz)
+ break;
+ }
+ if(rays.getMinT(which_one) < t)
+ break;
+ }
+ }
+ }
+
+void IsosurfaceGridVolume::calc_mcell(int depth, int startx, int starty, int
startz, VMCell& mcell)
+{
+ mcell.min=std::numeric_limits<ST>::max();
+ mcell.max=std::numeric_limits<ST>::min();
+ int endx=startx+xsize[depth];
+ int endy=starty+ysize[depth];
+ int endz=startz+zsize[depth];
+ if(endx>nx-1)
+ endx=nx-1;
+ if(endy>ny-1)
+ endy=ny-1;
+ if(endz>nz-1)
+ endz=nz-1;
+ if(startx>=endx || starty>=endy || startz>=endz){
+ /* This cell won't get used... */
+ mcell.min=datamax;
+ mcell.max=datamin;
+ return;
+ }
+ if(depth==0){
+ for(int ix=startx;ix<endx;ix++){
+ for(int iy=starty;iy<endy;iy++){
+ for(int iz=startz;iz<endz;iz++){
+ ST rhos[8];
+ rhos[0]=blockdata(ix, iy, iz);
+ rhos[1]=blockdata(ix, iy, iz+1);
+ rhos[2]=blockdata(ix, iy+1, iz);
+ rhos[3]=blockdata(ix, iy+1, iz+1);
+ rhos[4]=blockdata(ix+1, iy, iz);
+ rhos[5]=blockdata(ix+1, iy, iz+1);
+ rhos[6]=blockdata(ix+1, iy+1, iz);
+ rhos[7]=blockdata(ix+1, iy+1, iz+1);
+ ST min=rhos[0];
+ ST max=rhos[0];
+ for(int i=1;i<8;i++){
+ if(rhos[i]<min)
+ min=rhos[i];
+ if(rhos[i]>max)
+ max=rhos[i];
+ }
+ if(min<mcell.min)
+ mcell.min=min;
+ if(max>mcell.max)
+ mcell.max=max;
+ }
+ }
+ }
+ } else {
+ int nx=xsize[depth-1];
+ int ny=ysize[depth-1];
+ int nz=zsize[depth-1];
+ BrickArray3<VMCell>& mcells=macrocells[depth];
+ for(int x=startx;x<endx;x++){
+ for(int y=starty;y<endy;y++){
+ for(int z=startz;z<endz;z++){
+ VMCell tmp;
+ calc_mcell(depth-1, x*nx, y*ny, z*nz, tmp);
+ if(tmp.min < mcell.min)
+ mcell.min=tmp.min;
+ if(tmp.max > mcell.max)
+ mcell.max=tmp.max;
+ mcells(x,y,z)=tmp;
+ }
+ }
+ }
+ }
+}
Added: trunk/Model/Primitives/IsosurfaceGridVolume.h
==============================================================================
--- (empty file)
+++ trunk/Model/Primitives/IsosurfaceGridVolume.h Tue May 23 10:59:53
2006
@@ -0,0 +1,111 @@
+
+#ifndef _MANTA_ISOSURFACEGRIDVOLUME_H_
+#define _MANTA_ISOSURFACEGRIDVOLUME_H_
+
+#include <Model/Primitives/PrimitiveCommon.h>
+#include <SCIRun/Core/Containers/Array3.h>
+#include <Model/MiscObjects/BrickArray3.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>
+
+//generate the octree data
+//#define TEST_OCTDATA
+
+//test the octree data
+//#define USE_OCTREE_DATA
+//#define USE_OTD_NEIGHBOR_FIND
+
+typedef unsigned char ST;
+
+namespace Manta
+{
+
+ class IsosurfaceGridVolume : public PrimitiveCommon
+ {
+ public:
+ struct VMCell
+ {
+ ST max;
+ ST min;
+ };
+
+ BBox bounds;
+ Vector datadiag;
+ Vector hierdiag;
+ Vector ihierdiag;
+ Vector sdiag;
+ int nx,ny,nz;
+
+ #ifdef TEST_OCTDATA
+ OctreeData<ST> octdata;
+ #endif
+
+ BrickArray3<ST> blockdata;
+ ST datamin, datamax;
+ int mc_depth;
+ int* xsize;
+ int* ysize;
+ int* zsize;
+ double* ixsize;
+ double* iysize;
+ double* izsize;
+ ST isoval;
+
+ BrickArray3<VMCell>* macrocells;
+ char* filebase;
+
+ IsosurfaceGridVolume(char* filebase, int _mc_depth, Material* _matl);
+ ~IsosurfaceGridVolume();
+
+ 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;
+
+ inline double get_isovalue() const
+ {
+ return (double)isoval;
+ }
+
+ inline void set_isovalue(double val)
+ {
+ isoval = (ST)val;
+ }
+
+ inline void increase_isovalue()
+ {
+ if (isoval < std::numeric_limits<ST>::max())
+ isoval++;
+ }
+
+ inline void decrease_isovalue()
+ {
+ if (isoval > std::numeric_limits<ST>::min())
+ isoval--;
+ }
+
+ private:
+
+ inline void single_intersect(RayPacket& rays, int which_one) const;
+
+ void isect(RayPacket& rays, int which_one,
+ const Vector& orig, const Vector& dir, const Vector& inv_dir,
+ int depth, float t,
+ float dtdx, float dtdy, float dtdz,
+ float next_x, float next_y, float next_z,
+ int ix, int iy, int iz,
+ int dix_dx, int diy_dy, int diz_dz,
+ int startx, int starty, int startz,
+ const Vector& cellcorner, const Vector& celldir) const;
+
+ void calc_mcell(int depth, int startx, int starty, int startz,
VMCell& mcell);
+ };
+
+};
+
+#endif
Modified: trunk/Model/Primitives/IsosurfaceOctreeVolume.cc
==============================================================================
--- trunk/Model/Primitives/IsosurfaceOctreeVolume.cc (original)
+++ trunk/Model/Primitives/IsosurfaceOctreeVolume.cc Tue May 23 10:59:53
2006
@@ -11,7 +11,7 @@
#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};
@@ -27,6 +27,7 @@
void IsosurfaceOctreeVolume::preprocess( PreprocessContext const &context )
{
+ PrimitiveCommon::preprocess(context);
}
void IsosurfaceOctreeVolume::computeBounds( PreprocessContext const
&context, BBox &box ) const
@@ -59,7 +60,7 @@
Vector dir = rays.getDirection(which_one);
//Vector inv_dir = rays.getInverseDirection(which_one);
Vector inv_dir = dir.inverse();
-
+
#pragma unroll(3)
for(int axis=0; axis<3; axis++)
{
@@ -68,22 +69,22 @@
t1p.data[axis] = (octdata->padded_dims[axis] - orig.data[axis]) *
inv_dir.data[axis];
}
float tenter, texit, tenter_padded, texit_padded;
- if (rays.getSign(which_one,0))
- {
- tenter = t1.data[0];
- tenter_padded = t1p.data[0];
- texit = texit_padded = t0.data[0];
- }
- else
+ if (dir.data[0] > 0.f)
{
tenter = tenter_padded = t0.data[0];
texit = t1.data[0];
texit_padded = t1p.data[0];
}
+ else
+ {
+ tenter = t1.data[0];
+ tenter_padded = t1p.data[0];
+ texit = texit_padded = t0.data[0];
+ }
#pragma unroll(2)
for(int axis=1; axis<3; axis++)
{
- if (rays.getSign(which_one,axis))
+ if (dir.data[axis] > 0.f)
{
tenter = MAX(tenter, t0.data[axis]);
tenter_padded = MAX(tenter_padded, t0.data[axis]);
@@ -104,14 +105,26 @@
if (texit < 0.f)
return;
- //tenter = MAX(0.f, tenter);
+ //tenter_padded = MAX(0.f, tenter_padded);
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_padded);
+ #if 0
+ if (rays.hit(which_one, tenter_padded, PrimitiveCommon::getMaterial(),
this, 0))
+ {
+ rays.setNormal(which_one, Vector(-1,0,0));
+ return;
+ }
+ else
+ {
+ bool breakhere = true;
+ }
+ #else
+ single_traverse_node(rays, which_one, orig, dir, inv_dir, 0,
+ 0, 0, index_trace, cell, tenter_padded,
texit_padded);
+ #endif
}
bool IsosurfaceOctreeVolume::single_traverse_node(RayPacket& rays, int
which_one,
@@ -129,7 +142,7 @@
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;
+ Vector penter = orig + (dir*tenter);
Vec3i child_cell = cell;
Vec3i tc( penter.x() >= center.x(), penter.y() >= center.y(), penter.z()
>= center.z() );
@@ -140,36 +153,36 @@
Vec3i axis_isects;
if (tcenter.data[0] < tcenter.data[1] && tcenter.data[0] <
tcenter.data[2]){
- axis_isects[0] = 0;
+ axis_isects.data[0] = 0;
if (tcenter.data[1] < tcenter.data[2]){
- axis_isects[1] = 1;
- axis_isects[2] = 2;
+ axis_isects.data[1] = 1;
+ axis_isects.data[2] = 2;
}
else{
- axis_isects[1] = 2;
- axis_isects[2] = 1;
+ axis_isects.data[1] = 2;
+ axis_isects.data[2] = 1;
}
}
else if (tcenter.data[1] < tcenter.data[2]){
- axis_isects[0] = 1;
+ axis_isects.data[0] = 1;
if (tcenter.data[0] < tcenter.data[2]){
- axis_isects[1] = 0;
- axis_isects[2] = 2;
+ axis_isects.data[1] = 0;
+ axis_isects.data[2] = 2;
}
else{
- axis_isects[1] = 2;
- axis_isects[2] = 0;
+ axis_isects.data[1] = 2;
+ axis_isects.data[2] = 0;
}
}
else{
- axis_isects[0] = 2;
+ axis_isects.data[0] = 2;
if (tcenter.data[0] < tcenter.data[1]){
- axis_isects[1] = 0;
- axis_isects[2] = 1;
+ axis_isects.data[1] = 0;
+ axis_isects.data[2] = 1;
}
else{
- axis_isects[1] = 1;
- axis_isects[2] = 0;
+ axis_isects.data[1] = 1;
+ axis_isects.data[2] = 0;
}
}
@@ -187,7 +200,7 @@
child_texit = texit;
else
{
- axis = axis_isects[i];
+ axis = axis_isects.data[i];
child_texit = tcenter.data[axis];
if (child_texit < tenter)
continue;
@@ -211,7 +224,7 @@
{
unsigned int child_idx = node.children_start +
node.children[target_child].offset;
if (next_depth+1 == octdata->get_max_depth()) //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;
@@ -234,12 +247,12 @@
if (target_child & trueaxisbit) //going from true to
false
{
target_child &= ~trueaxisbit;
- child_cell[axis] &= ~child_bit;
+ child_cell.data[axis] &= ~child_bit;
}
else
//going from false to true
{
target_child |= trueaxisbit;
- child_cell[axis] |= child_bit;
+ child_cell.data[axis] |= child_bit;
}
}
return false;
@@ -257,7 +270,7 @@
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;
+ Vector penter = orig + (dir*tenter);
Vec3i child_cell = cell;
Vec3i tc( penter.x() >= center.x(), penter.y() >= center.y(), penter.z()
>= center.z() );
@@ -268,36 +281,36 @@
Vec3i axis_isects;
if (tcenter.data[0] < tcenter.data[1] && tcenter.data[0] <
tcenter.data[2]){
- axis_isects[0] = 0;
+ axis_isects.data[0] = 0;
if (tcenter.data[1] < tcenter.data[2]){
- axis_isects[1] = 1;
- axis_isects[2] = 2;
+ axis_isects.data[1] = 1;
+ axis_isects.data[2] = 2;
}
else{
- axis_isects[1] = 2;
- axis_isects[2] = 1;
+ axis_isects.data[1] = 2;
+ axis_isects.data[2] = 1;
}
}
else if (tcenter.data[1] < tcenter.data[2]){
- axis_isects[0] = 1;
+ axis_isects.data[0] = 1;
if (tcenter.data[0] < tcenter.data[2]){
- axis_isects[1] = 0;
- axis_isects[2] = 2;
+ axis_isects.data[1] = 0;
+ axis_isects.data[2] = 2;
}
else{
- axis_isects[1] = 2;
- axis_isects[2] = 0;
+ axis_isects.data[1] = 2;
+ axis_isects.data[2] = 0;
}
}
else{
- axis_isects[0] = 2;
+ axis_isects.data[0] = 2;
if (tcenter.data[0] < tcenter.data[1]){
- axis_isects[1] = 0;
- axis_isects[2] = 1;
+ axis_isects.data[1] = 0;
+ axis_isects.data[2] = 1;
}
else{
- axis_isects[1] = 1;
- axis_isects[2] = 0;
+ axis_isects.data[1] = 1;
+ axis_isects.data[2] = 0;
}
}
@@ -315,7 +328,7 @@
child_texit = texit;
else
{
- axis = axis_isects[i];
+ axis = axis_isects.data[i];
child_texit = tcenter.data[axis];
if (child_texit < tenter)
continue;
@@ -351,7 +364,7 @@
}
else
this_rho = scalar;
- rho[offset.data[0]][offset.data[1]][offset.data[2]] =
static_cast<float>(this_rho);
+ rho[0][0][1] = static_cast<float>(this_rho);
//0,1,1
offset.data[1] = 1;
@@ -363,7 +376,7 @@
}
else
this_rho = scalar;
- rho[offset.data[0]][offset.data[1]][offset.data[2]] =
static_cast<float>(this_rho);
+ rho[0][1][1] = static_cast<float>(this_rho);
//1,1,1
offset.data[0] = 1;
@@ -375,7 +388,7 @@
}
else
this_rho = scalar;
- rho[offset.data[0]][offset.data[1]][offset.data[2]] =
static_cast<float>(this_rho);
+ rho[1][1][1] = static_cast<float>(this_rho);
//1,1,0
offset.data[2] = 0;
@@ -387,7 +400,7 @@
}
else
this_rho = scalar;
- rho[offset.data[0]][offset.data[1]][offset.data[2]] =
static_cast<float>(this_rho);
+ rho[1][1][0] = static_cast<float>(this_rho);
//1,0,0
offset.data[1] = 0;
@@ -399,7 +412,7 @@
}
else
this_rho = scalar;
- rho[offset.data[0]][offset.data[1]][offset.data[2]] =
static_cast<float>(this_rho);
+ rho[1][0][0] = static_cast<float>(this_rho);
//1,0,1
offset.data[2] = 1;
@@ -411,7 +424,7 @@
}
else
this_rho = scalar;
- rho[offset.data[0]][offset.data[1]][offset.data[2]] =
static_cast<float>(this_rho);
+ rho[1][0][1] = static_cast<float>(this_rho);
//0,1,0
offset.data[0] = 0;
@@ -425,19 +438,19 @@
}
else
this_rho = scalar;
- rho[offset.data[0]][offset.data[1]][offset.data[2]] =
static_cast<float>(this_rho);
+ rho[0][1][0] = 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);
+#define MYDATA octdata->indata //toggle this to octdata if you want to test
pure point location (no neighbor finding)
+ min_rho = max_rho = lookup_safe(MYDATA, child_cell.data[0],
child_cell.data[1], child_cell.data[2]);
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);
+ ST this_rho = lookup_safe(MYDATA, neighboridx.data[0],
neighboridx.data[1], neighboridx.data[2]);
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);
@@ -445,20 +458,28 @@
#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(rays,
which_one, cmin, cmax, rho,
+ 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 = rays.getOrigin(which_one) +
rays.getDirection(which_one)*hit_t;
+ Vector phit = orig + dir*hit_t;
IsosurfaceImplicit::single_normal(normal, cmin,
cmax, phit, rho);
normal.normalize();
rays.setNormal(which_one, normal);
return true;
}
}
+ #endif
}
}
else //not at cap-level depth
@@ -479,12 +500,12 @@
if (target_child & trueaxisbit) //going from true to
false
{
target_child &= ~trueaxisbit;
- child_cell[axis] &= ~child_bit;
+ child_cell.data[axis] &= ~child_bit;
}
else
//going from false to true
{
target_child |= trueaxisbit;
- child_cell[axis] |= child_bit;
+ child_cell.data[axis] |= child_bit;
}
}
return false;
@@ -498,7 +519,7 @@
OctCap& cap = octdata->get_cap(cap_index);
index_trace[depth] = cap_index;
- Vector penter = orig + dir*tenter;
+ 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);
@@ -511,36 +532,36 @@
Vec3i axis_isects;
if (tcenter.data[0] < tcenter.data[1] && tcenter.data[0] <
tcenter.data[2]){
- axis_isects[0] = 0;
+ axis_isects.data[0] = 0;
if (tcenter.data[1] < tcenter.data[2]){
- axis_isects[1] = 1;
- axis_isects[2] = 2;
+ axis_isects.data[1] = 1;
+ axis_isects.data[2] = 2;
}
else{
- axis_isects[1] = 2;
- axis_isects[2] = 1;
+ axis_isects.data[1] = 2;
+ axis_isects.data[2] = 1;
}
}
else if (tcenter.data[1] < tcenter.data[2]){
- axis_isects[0] = 1;
+ axis_isects.data[0] = 1;
if (tcenter.data[0] < tcenter.data[2]){
- axis_isects[1] = 0;
- axis_isects[2] = 2;
+ axis_isects.data[1] = 0;
+ axis_isects.data[2] = 2;
}
else{
- axis_isects[1] = 2;
- axis_isects[2] = 0;
+ axis_isects.data[1] = 2;
+ axis_isects.data[2] = 0;
}
}
else{
- axis_isects[0] = 2;
+ axis_isects.data[0] = 2;
if (tcenter.data[0] < tcenter.data[1]){
- axis_isects[1] = 0;
- axis_isects[2] = 1;
+ axis_isects.data[1] = 0;
+ axis_isects.data[2] = 1;
}
else{
- axis_isects[1] = 1;
- axis_isects[2] = 0;
+ axis_isects.data[1] = 1;
+ axis_isects.data[2] = 0;
}
}
@@ -556,7 +577,7 @@
child_texit = texit;
else
{
- axis = axis_isects[i];
+ axis = axis_isects.data[i];
child_texit = tcenter.data[axis];
if (child_texit < tenter)
continue;
@@ -586,7 +607,7 @@
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);
+ rho[0][0][1] = static_cast<float>(this_rho);
//0,1,1
offset.data[1] = 1;
@@ -596,7 +617,7 @@
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);
+ rho[0][1][1] = static_cast<float>(this_rho);
//1,1,1
offset.data[0] = 1;
@@ -606,7 +627,7 @@
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);
+ rho[1][1][1] = static_cast<float>(this_rho);
//1,1,0
offset.data[2] = 0;
@@ -616,7 +637,7 @@
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);
+ rho[1][1][0] = static_cast<float>(this_rho);
//1,0,0
offset.data[1] = 0;
@@ -626,7 +647,7 @@
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);
+ rho[1][0][0] = static_cast<float>(this_rho);
//1,0,1
offset.data[2] = 1;
@@ -636,7 +657,7 @@
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);
+ rho[1][0][1] = static_cast<float>(this_rho);
//0,1,0
offset.data[0] = 0;
@@ -648,19 +669,19 @@
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);
+ rho[0][1][0] = 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);
+#define MYDATA octdata->indata
+ min_rho = max_rho = lookup_safe(MYDATA, child_cell.data[0],
child_cell.data[1], child_cell.data[2]);
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);
+ ST this_rho = lookup_safe(MYDATA, neighboridx.data[0],
neighboridx.data[1], neighboridx.data[2]);
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);
@@ -668,20 +689,28 @@
#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(rays, which_one, cmin,
cmax, rho,
+ 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 = rays.getOrigin(which_one) +
rays.getDirection(which_one)*hit_t;
+ Vector phit = orig + dir*hit_t;
IsosurfaceImplicit::single_normal(normal, cmin, cmax,
phit, rho);
normal.normalize();
rays.setNormal(which_one, normal);
return true;
}
}
+ #endif
}
if (axis==-1)
@@ -693,12 +722,12 @@
if (target_child & trueaxisbit) //going from true to
false
{
target_child &= ~trueaxisbit;
- child_cell[axis] &= ~1;
+ child_cell.data[axis] &= ~1;
}
else
//going from false to true
{
target_child |= trueaxisbit;
- child_cell[axis] |= 1;
+ child_cell.data[axis] |= 1;
}
}
return false;
Modified: trunk/Model/Primitives/OctreeVolume.cc
==============================================================================
--- trunk/Model/Primitives/OctreeVolume.cc (original)
+++ trunk/Model/Primitives/OctreeVolume.cc Tue May 23 10:59:53 2006
@@ -158,7 +158,10 @@
exit(1);
}
- Array3<ST> indata(nx, ny, nz);
+#ifndef TEST_INDATA
+ Array3<ST> indata;
+#endif
+ indata.resize(nx, ny, nz);
cerr << "indata(150,150,150)=" << indata(150,150,150) << endl;
cerr << "Reading " << filebase_numbered << "...";
cerr << "(" << indata.get_datasize() << " read/" << nx * ny * nz
* sizeof(ST) << " expected)...";
@@ -512,6 +515,7 @@
{
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));
@@ -670,6 +674,22 @@
}
}
fclose(in2);
+
+#ifdef TEST_INDATA
+ sprintf(filename, "%s.raw", filebase);
+ FILE* din = fopen(filename, "rb");
+ if (!din)
+ {
+ cerr << "Error opening original data file: " << filename << '\n';
+ exit(1);
+ }
+ 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);
+ cerr << "done.\n";
+ fclose(din);
+#endif
cerr << "done!" << endl;
Modified: trunk/Model/Primitives/OctreeVolume.h
==============================================================================
--- trunk/Model/Primitives/OctreeVolume.h (original)
+++ trunk/Model/Primitives/OctreeVolume.h Tue May 23 10:59:53 2006
@@ -15,6 +15,11 @@
#ifndef _MANTA_OCTREEVOLUME_H_
#define _MANTA_OCTREEVOLUME_H_
+#define lookup_safe(data, x, y, z) \
+(x > 0 && x < data.dim1() && y > 0 \
+&& y < data.dim2() && z > 0 && z < data.dim3()) ? data(x,y,z) : 0 \
+\
+
typedef unsigned char ST;
//AARONBAD - this is already defined in Abe's evil.h.
@@ -27,6 +32,9 @@
#define MAX(a,b) (((a)>(b))?(a):(b))
#endif
+//enable this only when testing octree data and original data side by side
+//#define TEST_INDATA
+
#include <Core/Geometry/BBox.h>
#include <SCIRun/Core/Containers/Array3.h>
#include <Core/Geometry/vecdefs.h>
@@ -87,6 +95,9 @@
{
public:
+#ifdef TEST_INDATA
+ SCIRun::Array3<ST> indata;
+#endif
//the maximum depth of the octree, finest resolution
int max_depth;
@@ -414,6 +425,8 @@
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())
@@ -421,6 +434,7 @@
else
return 0;
}
+ */
};
};
Modified: trunk/scenes/CMakeLists.txt
==============================================================================
--- trunk/scenes/CMakeLists.txt (original)
+++ trunk/scenes/CMakeLists.txt Tue May 23 10:59:53 2006
@@ -73,6 +73,13 @@
TARGET_LINK_LIBRARIES(scene_octisovol ${MANTA_SCENE_LINK})
ENDIF(SCENE_OCTISOVOL)
+# grid isosurface
+SET(SCENE_GRIDISOVOL 0 CACHE BOOL "grid isosurface")
+IF(SCENE_GRIDISOVOL)
+ ADD_LIBRARY(scene_gridisovol gridisovol.cc)
+ TARGET_LINK_LIBRARIES(scene_gridisovol ${MANTA_SCENE_LINK})
+ENDIF(SCENE_GRIDISOVOL)
+
SET(SCENE_OBJVIEWER 0 CACHE BOOL "Wavefront Obj file viewer.")
IF(SCENE_OBJVIEWER)
ADD_LIBRARY(scene_objviewer objviewer.cc)
Added: trunk/scenes/gridisovol.cc
==============================================================================
--- (empty file)
+++ trunk/scenes/gridisovol.cc Tue May 23 10:59:53 2006
@@ -0,0 +1,109 @@
+
+// Manta scene for grid-based 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/IsosurfaceGridVolume.h>
+#include <Model/Lights/PointLight.h>
+#include <Model/Lights/HeadLight.h>
+#include <Model/Materials/Phong.h>
+#include <Model/Materials/Lambertian.h>
+#include <Model/Cameras/PinholeCamera.h>
+#include <Core/Geometry/AffineTransform.h>
+#include <Model/AmbientLights/ConstantAmbient.h>
+#include <Model/MiscObjects/CuttingPlane.h>
+#include <SCIRun/Core/Containers/Array3.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 = "";
+ int macrocells = 3;
+
+ 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("gridisovol -file
<filename>", i, args);
+ }
+ else if (args[i] == "-macrocells") {
+ if (!getIntArg(i, args, macrocells))
+ throw IllegalArgument("octisovol -mres_levels
<mres_levels>", i, args);
+ }
+ else {
+ cerr << "Read built grid volume:" << endl;
+ cerr << "-file <filename>" << endl;
+ cerr << "-macrocells <#macrocells>" << endl;
+ throw IllegalArgument( "octisovol", i, args );
+ }
+ }
+
+ char c_filename[2048];
+ strcpy(c_filename, filename.c_str());
+
+ // Create the scene.
+ Scene *scene = new Scene();
+ Group *group = new Group();
+
+ Manta::BBox bounds;
+
+ //Material* mat1 = new Phong(Color(RGBColor(0.05f, 0.3f, 0.6f)),
Color(RGBColor(1.f, 1.f, 1.f)), 50);
+ Material* mat1 = new Lambertian(Color(RGBColor(0.05f, 0.3f, 0.6f)));
+ IsosurfaceGridVolume* igv = new IsosurfaceGridVolume(c_filename,
macrocells, mat1);
+ bounds = igv->getBounds();
+ group->add(igv);
+
+ // Add the tree to the scene.
+ scene->setObject( group );
+
+ float fov = 45.f;
+ float view_dist = 0.5f * bounds[1].data[0] * (1.0f + 1.0f / tanf(0.5f *
fov));
+ Vector lookat = (bounds[1] - bounds[0]) * 0.5f;
+ Vector eye = lookat;
+ Vector up = Vector(0,0,1);
+ eye.data[0] -= view_dist;
+
+ // Set other important scene properties.
+ PinholeCamera *camera = new PinholeCamera(eye, lookat, up, fov);
+ //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) ) ));
+ scene->setLights(lights);
+
+
+ // Background.
+ scene->setBackground( new ConstantBackground( Color(RGB(0.8, 0.8,
0.8)) ) );
+
+ return scene;
+}
+
Modified: trunk/scenes/octisovol.cc
==============================================================================
--- trunk/scenes/octisovol.cc (original)
+++ trunk/scenes/octisovol.cc Tue May 23 10:59:53 2006
@@ -47,7 +47,7 @@
int tend = 0;
double variance = 0.1;
int kernel_width = 2;
- int mres_levels = 2;
+ int mres_levels = 0;
double isomin = 0;
double isomax = 255;
@@ -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, isomin, isomax);
+ ov = new OctreeVolume(c_buildfrom_filename, tstart, tend, variance,
kernel_width, mres_levels, (double)isomin, (double)isomax);
}
else
{
- [MANTA] r1079 - in trunk: . Model/Intersections Model/MiscObjects Model/Primitives scenes, knolla, 05/23/2006
Archive powered by MHonArc 2.6.16.
