-
sparker
thiago
Manta Interactive Ray Tracer Development Mailing List
Text archives Help
- From: cgribble@sci.utah.edu
- To: manta@sci.utah.edu
- Subject: [MANTA] r1089 - in trunk: Model/Primitives Model/Readers scenes
- Date: Tue, 30 May 2006 08:41:02 -0600 (MDT)
Author: cgribble
Date: Tue May 30 08:41:00 2006
New Revision: 1089
Added:
trunk/Model/Primitives/GridSpheres.cc
trunk/Model/Primitives/GridSpheres.h
Modified:
trunk/Model/Primitives/CMakeLists.txt
trunk/Model/Readers/ParticleNRRD.cc
trunk/Model/Readers/ParticleNRRD.h
trunk/scenes/pnrrd.cc
Log:
Model/Primitives/CMakeLists.txt
Model/Primitives/GridSpheres.cc
Model/Primitives/GridSpheres.h
Added naive GridSpheres implementation (from RTRT)
Rays are traced one-at-a-time through the grid
Several optimizations could be added, not the least of which is an actual
packetized traversal
Exhibits small errors with shadowing; not sure why...
Model/Readers/ParticleNRRD.h
Model/Readers/ParticleNRRD.cc
Added nuke flag to constructor to indicate that data should be nuked by the
destructor; false by default
scenes/pnrrd.cc
Uses GridSpheres object for particles by default
Modified: trunk/Model/Primitives/CMakeLists.txt
==============================================================================
--- trunk/Model/Primitives/CMakeLists.txt (original)
+++ trunk/Model/Primitives/CMakeLists.txt Tue May 30 08:41:00 2006
@@ -10,6 +10,8 @@
Primitives/Cylinder.h
Primitives/Disk.cc
Primitives/Disk.h
+ Primitives/GridSpheres.cc
+ Primitives/GridSpheres.h
Primitives/HeavyTriangle.cc
Primitives/HeavyTriangle.h
Primitives/Heightfield.cc
Added: trunk/Model/Primitives/GridSpheres.cc
==============================================================================
--- (empty file)
+++ trunk/Model/Primitives/GridSpheres.cc Tue May 30 08:41:00 2006
@@ -0,0 +1,923 @@
+
+#include <Core/Math/MiscMath.h>
+#include <Core/Math/Trig.h>
+#include <Core/Math/Expon.h>
+#include <Interface/Light.h>
+#include <Interface/LightSet.h>
+#include <Interface/Primitive.h>
+#include <Interface/Packet.h>
+#include <Interface/RayPacket.h>
+#include <Interface/AmbientLight.h>
+#include <Interface/Context.h>
+#include <Interface/ShadowAlgorithm.h>
+#include <Model/Primitives/GridSpheres.h>
+#include <Model/Textures/Constant.h>
+#include <SCIRun/Core/Containers/Array1.h>
+#include <SCIRun/Core/Math/MinMax.h>
+#include <SCIRun/Core/Util/Timer.h>
+
+#include <iostream>
+using std::cerr;
+
+#include <float.h>
+
+using namespace Manta;
+using namespace SCIRun;
+
+void GridSpheres::init(void)
+{
+ if (radius <= 0) {
+ if (ridx <= 0)
+ cerr<<"Resetting default radius to 1\n";
+ radius=1;
+ }
+
+ // Compute inverses
+ inv_radius=1/static_cast<Real>(radius);
+ inv_ncells=1/static_cast<Real>(ncells);
+
+ cerr<<"Recomputing min/max for GridSpheres\n";
+
+ // Initialize min/max arrays
+ min=new float[nvars];
+ max=new float[nvars];
+
+ for (int j=0; j<nvars; ++j) {
+ min[j]=FLT_MAX;
+ max[j]=-FLT_MAX;
+ }
+
+ // Find min/max values
+ float* data=spheres;
+ for (int i=0; i<nspheres; ++i) {
+ for (int j=0; j<nvars; ++j) {
+ min[j]=SCIRun::Min(min[j], data[j]);
+ max[j]=SCIRun::Max(max[j], data[j]);
+ }
+
+ data += nvars;
+ }
+
+ counts=0;
+ cells=0;
+}
+
+GridSpheres::GridSpheres(const Color& color, float* spheres, int nspheres,
+ int nvars, int ncells, int depth, float radius,
+ int ridx) :
+ spheres(spheres), nspheres(nspheres), nvars(nvars), ncells(ncells),
+ depth(depth), radius(radius), inv_radius(1/radius), ridx(ridx)
+{
+ colortex=new Constant<Color>(color);
+ init();
+}
+
+GridSpheres::GridSpheres(const Texture<Color>* colorfn, float* spheres,
+ int nspheres, int nvars, int ncells, int depth,
+ float radius, int ridx) :
+ spheres(spheres), nspheres(nspheres), nvars(nvars), ncells(ncells),
+ depth(depth), radius(radius), inv_radius(1/radius), ridx(ridx),
+ colortex(colorfn)
+{
+ init();
+}
+
+GridSpheres::~GridSpheres()
+{
+ delete [] min;
+ delete [] max;
+
+ if (counts)
+ delete[] counts;
+
+ if (cells)
+ delete[] cells;
+
+ if (macrocells)
+ delete [] macrocells;
+}
+
+void GridSpheres::preprocess(const PreprocessContext& context)
+{
+ // Preprocess material
+ LitMaterial::preprocess(context);
+
+ // Build grid
+ cerr<<"Building GridSpheres\n";
+
+ WallClockTimer timer;
+ timer.start();
+
+ cerr<<" min: ("<<min[0]<<", "<<min[1]<<", "<<min[2]<<")\n";
+ cerr<<" max: ("<<max[0]<<", "<<max[1]<<", "<<max[2]<<")\n";
+
+ // Determine the maximum radius
+ float max_radius;
+ if (ridx>0) {
+ cerr<<" GridSpheres::preprocess - ridx="<<ridx<<'\n';
+ max_radius=max[ridx];
+
+ if (max_radius <= 0) {
+ cerr<<" max_radius ("<<max_radius<<") <= 0, setting to default radius
("
+ <<radius<<")\n";
+ max_radius=radius;
+ }
+ } else
+ max_radius=radius;
+
+ // Bound the spheres
+ computeBounds(context, bounds);
+ diagonal=bounds.diagonal();
+ inv_diagonal=Vector(1/diagonal.x(), 1/diagonal.y(), 1/diagonal.z());
+
+ // Determine grid size
+ totalcells=1;
+ for (int i=0; i<=depth; ++i)
+ totalcells *= ncells;
+
+ int totalsize=totalcells*totalcells*totalcells;
+ cerr<<" Computing "<<totalcells<<'x'<<totalcells<<'x'<<totalcells
+ <<" grid ("<<totalsize<<" cells total)\n";
+
+ // Clear grid data
+ if (counts)
+ delete[] counts;
+
+ if (cells)
+ delete[] cells;
+
+ // Allocate and initialize counts
+ counts=new int[2*totalsize];
+ bzero(counts, 2*totalsize*sizeof(int));
+
+ cerr<<" 0/6: Allocation took "<<timer.time()<<" seconds\n";
+
+ // Generate map
+ int* map=new int[totalsize];
+ int idx=0;
+ for (int x=0; x<totalcells; ++x) {
+ for (int y=0; y<totalcells; ++y) {
+ for (int z=0; z<totalcells; ++z) {
+ map[idx++]=mapIdx(x, y, z, depth);
+ }
+ }
+ }
+
+ Real stime=timer.time();
+ cerr<<" 1/6: Generating map took "<<stime<<" seconds\n";
+
+ // Compute cell counts
+ float* data=spheres;
+ int tc2=totalcells*totalcells;
+ for (int i=0; i<nspheres; ++i) {
+ Real ctime=timer.time();
+ if (ctime - stime>5.0) {
+ cerr<<i<<"/"<<nspheres<<'\n';
+ stime=ctime;
+ }
+
+ // Compute cell overlap
+ Vector current_radius;
+ if (ridx>0) {
+ if (data[ridx] <= 0)
+ continue;
+
+ current_radius=Vector(data[ridx], data[ridx], data[ridx]);
+ } else {
+ current_radius=Vector(radius, radius, radius);
+ }
+
+ Vector center(data[0], data[1], data[2]);
+ BBox box(center - current_radius, center + current_radius);
+ int sx, sy, sz, ex, ey, ez;
+ transformToLattice(box, sx, sy, sz, ex, ey, ez);
+
+ int idx_x=sx*tc2;
+ for (int x=sx; x<=ex; ++x) {
+ int idx_y=idx_x + sy*totalcells;
+ idx_x += tc2;
+ for (int y=sy; y<=ey; ++y) {
+ int idx=idx_y + sz;
+ idx_y += totalcells;
+ for (int z=sz; z<=ez; ++z) {
+ int aidx=map[idx++];
+ counts[2*aidx + 1]++;
+ }
+ }
+ }
+
+ data += nvars;
+ }
+
+ cerr<<" 2/6: Counting cells took "<<timer.time()<<" seconds\n";
+
+ int total=0;
+ for (int i=0; i<totalsize; ++i) {
+ int count=counts[2*i + 1];
+ counts[2*i]=total;
+ total += count;
+ }
+
+ // Allocate cells
+ cerr<<" Allocating "<<total<<" grid cells ("
+ <<static_cast<Real>(total)/nspheres<<" per object, "
+ <<static_cast<Real>(total)/totalsize<<" per cell)\n";
+
+ cells=new int[total];
+ for (int i=0; i<total; ++i)
+ cells[i]=-1;
+
+ stime=timer.time();
+ cerr<<" 3/6: Calculating offsets took "<<stime<<" seconds\n";
+
+ // Populate the grid
+ Array1<int> current(totalsize);
+ current.initialize(0);
+ data=spheres;
+ for (int i=0; i<nspheres; ++i) {
+ Real ctime=timer.time();
+ if (ctime - stime>5.0) {
+ cerr<<i<<"/"<<nspheres<<'\n';
+ stime=ctime;
+ }
+
+ // Compute cell overlap
+ Vector current_radius;
+ if (ridx>0) {
+ if (data[ridx] <= 0)
+ continue;
+
+ current_radius=Vector(data[ridx], data[ridx], data[ridx]);
+ } else {
+ current_radius=Vector(radius, radius, radius);
+ }
+
+ Vector center(data[0], data[1], data[2]);
+ BBox box(center - current_radius, center + current_radius);
+ int sx, sy, sz, ex, ey, ez;
+ transformToLattice(box, sx, sy, sz, ex, ey, ez);
+
+ for (int x=sx; x<=ex; ++x) {
+ for (int y=sy; y<=ey; ++y) {
+ int idx=totalcells*(totalcells*x + y) + sz;
+ for (int z=sz; z<=ez; ++z) {
+ int aidx=map[idx++];
+ int cur=current[aidx]++;
+ int pos=counts[2*aidx] + cur;
+ cells[pos]=nvars*i;
+ }
+ }
+ }
+
+ data += nvars;
+ }
+
+ delete [] map;
+
+ cerr<<" 4/6: Filling grid took "<<timer.time()<<" seconds\n";
+
+ // Verify the grid
+ for (int i=0; i<totalsize; ++i) {
+ if (current[i] != counts[2*i + 1]) {
+ cerr<<"Fatal error: current="<<current[i]<<", but counts="
+ <<counts[2*i + 1]<<'\n';
+ exit(1);
+ }
+ }
+
+ for (int i=0; i<total; ++i) {
+ if (cells[i]==-1) {
+ cerr<<"Fatal error: cells["<<i<<"] == -1\n";
+ exit(1);
+ }
+ }
+
+ cerr<<" 5/6: Verifying grid took "<<timer.time()<<" seconds\n";
+
+ // Build the macrocells
+ if (depth > 0) {
+ macrocells=new MCell*[depth + 1];
+ macrocells[0]=0;
+
+ int size=ncells*ncells*ncells;
+ for (int d=depth; d>=1; d--) {
+ MCell* mcell=new MCell[size];
+ macrocells[d]=mcell;
+
+ float* mm=new float[2*nvars*size];
+ for (int i=0; i<size; ++i) {
+ // Minimum
+ mcell->min=mm;
+ mm += nvars;
+
+ // Maximum
+ mcell->max=mm;
+ mm += nvars;
+
+ mcell++;
+ }
+
+ size *= ncells*ncells*ncells;
+ }
+
+ MCell top;
+ fillMCell(top, depth, 0);
+ if (top.nspheres != total) {
+ cerr<<"Fatal error: macrocell construction went wrong\n";
+ cerr<<" Top macrocell: "<<top.nspheres<<'\n';
+ cerr<<" Total nspheres: "<<total<<'\n';
+ exit(1);
+ }
+
+ cerr<<" 6/6: Calculating macrocells took "<<timer.time()<<"
seconds\n";
+ } else {
+ macrocells=0;
+ cerr<<" 6/6: Macrocell hierarchy not built (depth == 0)\n";
+ }
+
+ cerr<<"Done building GridSpheres\n";
+}
+
+void GridSpheres::computeBounds(const PreprocessContext& context,
+ BBox& bbox) const
+{
+ // Determine the maximum radius
+ float max_radius;
+ if (ridx>0) {
+ max_radius=max[ridx];
+
+ if (max_radius <= 0) {
+ cerr<<" max_radius ("<<max_radius<<") <= 0, setting to default radius
("
+ <<radius<<")\n";
+ max_radius=radius;
+ }
+ } else
+ max_radius=radius;
+
+ // Bound the spheres
+ Vector mr(max_radius, max_radius, max_radius);
+ bbox.reset(Vector(min[0], min[1], min[2]) - mr,
+ Vector(max[0], max[1], max[2]) + mr);
+
+ const Vector eps3(1.e-3, 1.e-3, 1.e-3);
+ bbox.extendByPoint(bbox.getMin() - eps3);
+ bbox.extendByPoint(bbox.getMax() + eps3);
+
+ Vector diag(bbox.diagonal());
+ bbox.extendByPoint(bbox.getMin() - diag*eps3);
+ bbox.extendByPoint(bbox.getMax() + diag*eps3);
+}
+
+void GridSpheres::intersect(const RenderContext& context, RayPacket& rays)
const
+{
+ Vector bmin=bounds.getMin();
+ Vector bmax=bounds.getMax();
+
+ rays.computeInverseDirections();
+ rays.computeSigns();
+
+ for (int i=rays.begin(); i<rays.end(); ++i) {
+ const Vector origin(rays.getOrigin(i));
+ const Vector direction(rays.getDirection(i));
+ const Vector inv_direction(rays.getInverseDirection(i));
+
+ // Intersect ray with bounding box
+#if 1
+ Real tnear, tfar;
+ int di_dx;
+ int ddx;
+ int didx_dx;
+ int stop_x;
+ if (direction.x()>0) {
+ tnear=(bmin.x() - origin.x())*inv_direction.x();
+ tfar=(bmax.x() - origin.x())*inv_direction.x();
+ di_dx=1;
+ ddx=1;
+ didx_dx=ncells*ncells;
+ stop_x=ncells;
+ } else {
+ tnear=(bmax.x() - origin.x())*inv_direction.x();
+ tfar=(bmin.x() - origin.x())*inv_direction.x();
+ di_dx=-1;
+ ddx=0;
+ didx_dx=-ncells*ncells;
+ stop_x=-1;
+ }
+
+ Real y0, y1;
+ int di_dy;
+ int ddy;
+ int didx_dy;
+ int stop_y;
+ if (direction.y()>0) {
+ y0=(bmin.y() - origin.y())*inv_direction.y();
+ y1=(bmax.y() - origin.y())*inv_direction.y();
+ di_dy=1;
+ ddy=1;
+ didx_dy=ncells;
+ stop_y=ncells;
+ } else {
+ y0=(bmax.y() - origin.y())*inv_direction.y();
+ y1=(bmin.y() - origin.y())*inv_direction.y();
+ di_dy=-1;
+ ddy=0;
+ didx_dy=-ncells;
+ stop_y=-1;
+ }
+
+ if (y0>tnear)
+ tnear=y0;
+ if (y1<tfar)
+ tfar=y1;
+ if (tfar<tnear)
+ continue;
+
+ Real z0, z1;
+ int di_dz;
+ int ddz;
+ int didx_dz;
+ int stop_z;
+ if (direction.z()>0) {
+ z0=(bmin.z() - origin.z())*inv_direction.z();
+ z1=(bmax.z() - origin.z())*inv_direction.z();
+ di_dz=1;
+ ddz=1;
+ didx_dz=1;
+ stop_z=ncells;
+ } else {
+ z0=(bmax.z() - origin.z())*inv_direction.z();
+ z1=(bmin.z() - origin.z())*inv_direction.z();
+ di_dz=-1;
+ ddz=0;
+ didx_dz=-1;
+ stop_z=-1;
+ }
+
+ if (z0>tnear)
+ tnear=z0;
+ if (z1<tfar)
+ tfar=z1;
+ if (tfar<tnear)
+ continue;
+ if (tfar<static_cast<Real>(1.e-6))
+ continue;
+ if (tnear < 0)
+ tnear=0;
+
+ // Compute lattice coordinates
+ Vector p=origin + tnear*direction;
+ Vector lattice=(p - bmin)*inv_diagonal;
+ int ix=Clamp(static_cast<int>(lattice.x()*ncells), 0, ncells - 1);
+ int iy=Clamp(static_cast<int>(lattice.y()*ncells), 0, ncells - 1);
+ int iz=Clamp(static_cast<int>(lattice.z()*ncells), 0, ncells - 1);
+
+ // Compute cell index
+ int idx=(ix*ncells + iy)*ncells + iz;
+
+ // Compute delta t in each direction
+ Real
dt_dx=static_cast<Real>(di_dx)*diagonal.x()*inv_ncells*inv_direction.x();
+ Real dt_dy=di_dy*diagonal.y()*inv_ncells*inv_direction.y();
+ Real dt_dz=di_dz*diagonal.z()*inv_ncells*inv_direction.z();
+
+ // Compute far edges of the cell
+ Vector next(ix + ddx, iy + ddy, iz + ddz);
+ Vector far=diagonal*next*inv_ncells + bmin;
+
+ // Compute t values at far edges of the cell
+ Vector tnext=(far - origin)*inv_direction;
+
+ // Compute cell corner and direction
+ Vector factor(ncells*inv_diagonal);
+ Vector cellcorner((origin - bmin)*factor);
+ Vector celldir(direction*factor);
+
+ // Traverse the hierarchy
+ traverse(i, rays, depth, tnear, ix, iy, iz, idx, dt_dx, dt_dy, dt_dz,
+ tnext.x(), tnext.y(), tnext.z(),
+ cellcorner, celldir,
+ di_dx, di_dy, di_dz, didx_dx, didx_dy, didx_dz,
+ stop_x, stop_y, stop_z);
+#else
+ Vector v1=(bmin - origin)*inv_direction;
+ Vector v2=(bmax - origin)*inv_direction;
+ Vector vmin=Min(v1, v2);
+ Vector vmax=Max(v1, v2);
+ Real tnear=vmin.maxComponent();
+ Real tfar=vmax.minComponent();
+ if (tnear >= tfar)
+ continue;
+ if (tfar < static_cast<Real>(1.e-6))
+ continue;
+ if (tnear < 0)
+ tnear=0;
+
+ // Compute lattice coordinates
+ Vector p=origin + tnear*direction;
+ Vector lattice=(p - bmin)*inv_diagonal;
+ Vector integer(Clamp(static_cast<int>(lattice.x()*ncells), 0, ncells -
1),
+ Clamp(static_cast<int>(lattice.y()*ncells), 0, ncells -
1),
+ Clamp(static_cast<int>(lattice.z()*ncells), 0, ncells -
1));
+
+ // Compute delta i and stopping criteria in each direction
+ Vector di=Vector(1, 1, 1) - 2*rays.getSigns(i);
+ Vector stop(di.x()>0?ncells:-1, di.y()>0?ncells:-1, di.z()>0?ncells:-1);
+
+ // Compute cell index
+ int idx=(integer.x()*ncells + integer.y())*ncells + integer.z();
+
+ // Compute delta index for each direction
+ Vector didx=di*Vector(ncells*ncells, ncells, 1);
+
+ // Compute delta t for each direction
+ Vector dt=di*diagonal*inv_ncells*inv_direction;
+
+ // Compute far edges of the cell
+ Vector next=integer + Vector(1, 1, 1) - rays.getSigns(i);
+ Vector far=diagonal*next*inv_ncells + bmin;
+
+ // Compute t values at far edges of the cell
+ Vector tnext=(far - origin)*inv_direction;
+
+ // Compute cell corner and direction
+ Vector factor(ncells*inv_diagonal);
+ Vector cellcorner((origin - bmin)*factor);
+ Vector celldir(direction*factor);
+
+ // Traverse the hierarchy
+ traverse(i, rays, depth, tnear, integer.x(), integer.y(), integer.z(),
+ idx, dt.x(), dt.y(), dt.z(), tnext.x(), tnext.y(), tnext.z(),
+ cellcorner, celldir, di.x(), di.y(), di.z(),
+ didx.x(), didx.y(), didx.z(), stop.x(), stop.y(), stop.z());
+#endif
+ }
+}
+
+void GridSpheres::computeNormal(const RenderContext& context,
+ RayPacket& rays) const
+{
+ rays.computeHitPositions();
+ for (unsigned int i=rays.begin(); i<rays.end(); ++i) {
+ float* data=spheres + rays.scratchpad<int>(i);
+ Vector n=rays.getHitPosition(i) - Vector(data[0], data[1], data[2]);
+
+ if (ridx>0) {
+ if (data[ridx] <= 0)
+ n *= inv_radius;
+ else
+ n /= data[ridx];
+ } else {
+ n *= inv_radius;
+ }
+
+ rays.setNormal(i, n);
+ }
+
+ rays.setFlag(RayPacket::HaveUnitNormals);
+}
+
+void GridSpheres::shade(const RenderContext& context, RayPacket& rays) const
+{
+ // Shade a bunch of rays that have intersected the same particle
+
+ // Compute normals
+ rays.computeNormals(context);
+
+ // Compute colors
+ Packet<Color> diffuse;
+ colortex->mapValues(diffuse, context, rays);
+
+ // Compute ambient contributions for all rays
+ ColorArray totalLight;
+ activeLights->getAmbientLight()->computeAmbient(context, rays, totalLight);
+
+ // Normalize directions for proper dot product computation
+ rays.normalizeDirections();
+
+ ShadowAlgorithm::StateBuffer stateBuffer;
+ bool firstTime=true;
+ bool done;
+ int count=0;
+ do {
+ RayPacketData shadowData;
+ RayPacket shadowRays(shadowData, RayPacket::UnknownShape, 0, 0,
+ rays.getDepth(), 0);
+
+ // Call the shadow algorithm (SA) to generate shadow rays. We may not be
+ // able to compute all of them, so we pass along a buffer for the SA
+ // object to store it's state. The firstTime flag tells the SA to fill
+ // in the state rather than using anything in the state buffer. Most
+ // SAs will only need to store an int or two in the statebuffer.
+ done=context.shadowAlgorithm->computeShadows(context, activeLights,
+ rays, shadowRays,
+ firstTime, stateBuffer);
+
+ // Normalize directions for proper dot product computation
+ shadowRays.normalizeDirections();
+
+ for (int i=shadowRays.begin(); i < shadowRays.end(); ++i) {
+ if (!shadowRays.wasHit(i)) {
+ // Not in shadow, so compute the direct and specular contributions
+ Vector normal=rays.getNormal(i);
+ Vector shadowdir=shadowRays.getDirection(i);
+ ColorComponent cos_theta=Dot(shadowdir, normal);
+ Color light=shadowRays.getColor(i);
+ for (int j=0; j < Color::NumComponents; ++j)
+ totalLight[j][i] += light[j]*cos_theta;
+ }
+ }
+
+ firstTime=false;
+ } while(!done);
+
+ // Sum up diffuse/specular contributions
+ for (int i=rays.begin(); i < rays.end(); ++i) {
+ Color result;
+ for (int j=0;j<Color::NumComponents; ++j)
+ result[j]=totalLight[j][i]*diffuse.colordata[j][i];
+ rays.setColor(i, result);
+ }
+}
+
+void GridSpheres::computeTexCoords2(const RenderContext& context,
+ RayPacket& rays) const
+{
+ rays.computeHitPositions();
+ rays.computeNormals(context);
+ for(int i=rays.begin();i<rays.end();i++){
+ Vector n=rays.getNormal(i);
+ Real angle=Clamp(n.z(), (Real)-1, (Real)1);
+ Real theta=Acos(angle);
+ Real phi=Atan2(n.y(), n.x());
+ Real x=phi*(Real)(0.5*M_1_PI);
+ if (x < 0)
+ x += 1;
+ Real y=theta*(Real)M_1_PI;
+ rays.setTexCoords(i, Vector(x, y, 0));
+ }
+
+ rays.setFlag(RayPacket::HaveTexture2|RayPacket::HaveTexture3);
+}
+
+void GridSpheres::computeTexCoords3(const RenderContext& context,
+ RayPacket& rays) const
+{
+ rays.computeHitPositions();
+ rays.computeNormals(context);
+ for(int i=rays.begin();i<rays.end();i++){
+ Vector n=rays.getNormal(i);
+ Real angle=Clamp(n.z(), (Real)-1, (Real)1);
+ Real theta=Acos(angle);
+ Real phi=Atan2(n.y(), n.x());
+ Real x=phi*(Real)(0.5*M_1_PI);
+ if (x < 0)
+ x += 1;
+ Real y=theta*(Real)M_1_PI;
+ rays.setTexCoords(i, Vector(x, y, 0));
+ }
+
+ rays.setFlag(RayPacket::HaveTexture2|RayPacket::HaveTexture3);
+}
+
+void GridSpheres::traverse(int ray_idx, RayPacket& rays, int depth, Real
tnear,
+ int ix, int iy, int iz, int idx,
+ Real dt_dx, Real dt_dy, Real dt_dz,
+ Real tnext_x, Real tnext_y, Real tnext_z,
+ const Vector& cellcorner, const Vector& celldir,
+ int di_dx, int di_dy, int di_dz,
+ int didx_dx, int didx_dy, int didx_dz,
+ int stop_x, int stop_y, int stop_z) const
+{
+ if (depth>0) {
+ // Traverse the macrocell layers
+ MCell* mcells=macrocells[depth];
+ while (tnear < rays.getMinT(ray_idx)) {
+ MCell& mcell=mcells[idx];
+ if (mcell.nspheres > 0) {
+ // XXX: Range checking would go here... Ignore for now
+
+ // Compute lattice coordinates
+ int new_ix=Clamp(static_cast<int>((cellcorner.x() +
tnear*celldir.x() - ix)*ncells), 0, ncells - 1);
+ int new_iy=Clamp(static_cast<int>((cellcorner.y() +
tnear*celldir.y() - iy)*ncells), 0, ncells - 1);
+ int new_iz=Clamp(static_cast<int>((cellcorner.z() +
tnear*celldir.z() - iz)*ncells), 0, ncells - 1);
+
+ // Compute new cell index
+ int new_idx=((idx*ncells + new_ix)*ncells + new_iy)*ncells + new_iz;
+
+ // Compute new delta t in each direction
+ Real new_dt_dx=dt_dx*inv_ncells;
+ Real new_dt_dy=dt_dy*inv_ncells;
+ Real new_dt_dz=dt_dz*inv_ncells;
+
+ // Compute new t values at far edges of the cell
+ Vector signs=rays.getSigns(ray_idx);
+ Real new_tnext_x=tnext_x + (1 - 2*signs.x())*new_ix*new_dt_dx +
+ (1 - signs.x())*(new_dt_dx - dt_dx);
+ Real new_tnext_y=tnext_y + (1 - 2*signs.y())*new_iy*new_dt_dy +
+ (1 - signs.y())*(new_dt_dy - dt_dy);
+ Real new_tnext_z=tnext_z + (1 - 2*signs.z())*new_iz*new_dt_dz +
+ (1 - signs.z())*(new_dt_dz - dt_dz);
+
+ // Compute new cell corner and direction
+ Vector new_cellcorner=(cellcorner - Vector(ix, iy, iz))*ncells;
+ Vector new_celldir=celldir*ncells;
+
+ // Descend to next depth in the hierarchy
+ traverse(ray_idx, rays, depth - 1, tnear,
+ new_ix, new_iy, new_iz,
+ new_idx,
+ new_dt_dx, new_dt_dy, new_dt_dz,
+ new_tnext_x, new_tnext_y, new_tnext_z,
+ new_cellcorner, new_celldir,
+ di_dx, di_dy, di_dz,
+ didx_dx, didx_dy, didx_dz,
+ stop_x, stop_y, stop_z);
+ }
+
+ // March to next macrocell at the current depth
+ if (tnext_x<tnext_y && tnext_x<tnext_z) {
+ ix += di_dx;
+ if (ix == stop_x)
+ break;
+ tnear=tnext_x;
+ tnext_x += dt_dx;
+ idx += didx_dx;
+ } else if (tnext_y<tnext_z) {
+ iy += di_dy;
+ if (iy == stop_y)
+ break;
+ tnear=tnext_y;
+ tnext_y += dt_dy;
+ idx += didx_dy;
+ } else {
+ iz += di_dz;
+ if (iz == stop_z)
+ break;
+ tnear=tnext_z;
+ tnext_z += dt_dz;
+ idx += didx_dz;
+ }
+ }
+ } else {
+ // Traverse cells, intersecting spheres along the way as necessary
+ while (tnear < rays.getMinT(ray_idx)) {
+ int start=counts[2*idx];
+ int nsph=counts[2*idx + 1];
+ for (int j=0; j<nsph; ++j) {
+ float* data=spheres + cells[start + j];
+
+ // XXX: Range checking would go here... Ignore for now
+
+ // Sphere is in range, determine it's radius (squared)
+ float radius2;
+ if (ridx>0) {
+ float current_radius=data[ridx];
+ if (current_radius <= 0)
+ continue;
+
+ radius2=current_radius*current_radius;
+ } else {
+ radius2=radius*radius;
+ }
+
+ intersectSphere(rays, ray_idx, start + j,
+ Vector(data[0], data[1], data[2]), radius2);
+ }
+
+ // March to the next cell
+ if (tnext_x < tnext_y && tnext_x < tnext_z) {
+ ix += di_dx;
+ if (ix == stop_x)
+ break;
+ tnear=tnext_x;
+ tnext_x += dt_dx;
+ idx += didx_dx;
+ } else if (tnext_y < tnext_z){
+ iy += di_dy;
+ if (iy == stop_y)
+ break;
+ tnear=tnext_y;
+ tnext_y += dt_dy;
+ idx += didx_dy;
+ } else {
+ iz += di_dz;
+ if (iz == stop_z)
+ break;
+ tnear=tnext_z;
+ tnext_z += dt_dz;
+ idx += didx_dz;
+ }
+ }
+ }
+}
+
+void GridSpheres::intersectSphere(RayPacket& rays, int ray_idx, int idx,
+ const Vector& center, float radius2) const
+{
+ Vector CO=center - rays.getOrigin(ray_idx);
+ Real tca=Dot(CO, rays.getDirection(ray_idx));
+ Real l2oc=CO.length2();
+ if (l2oc <= radius2) {
+ // Ray origin is inside the sphere
+ Real t2hc=radius2 - l2oc + tca*tca;
+ Real thc=Sqrt(t2hc);
+ Real t=tca + thc;
+ if (rays.hit(ray_idx, tca - thc, this, this, this))
+ rays.scratchpad<int>(ray_idx)=cells[idx];
+ } else {
+ if (tca>=0) {
+ Real t2hc=radius2 - l2oc + tca*tca;
+ if (t2hc>0) {
+ Real thc=sqrt(t2hc);
+ if (rays.hit(ray_idx, tca - thc, this, this, this))
+ rays.scratchpad<int>(ray_idx)=cells[idx];
+ } else {
+ // Ray misses, no intersection
+ }
+ } else {
+ // Sphere is behind the ray origin
+ }
+ }
+}
+
+void GridSpheres::transformToLattice(const BBox& box, int& sx, int& sy, int&
sz,
+ int& ex, int& ey, int& ez) const
+{
+ Vector s=(box.getMin() - bounds.getMin())*inv_diagonal;
+ sx=static_cast<int>(s.x()*totalcells);
+ sy=static_cast<int>(s.y()*totalcells);
+ sz=static_cast<int>(s.z()*totalcells);
+ Clamp(sx, 0, totalcells - 1);
+ Clamp(sy, 0, totalcells - 1);
+ Clamp(sz, 0, totalcells - 1);
+
+ Vector e=(box.getMax() - bounds.getMin())*inv_diagonal;
+ ex=static_cast<int>(e.x()*totalcells);
+ ey=static_cast<int>(e.y()*totalcells);
+ ez=static_cast<int>(e.z()*totalcells);
+ Clamp(ex, 0, totalcells - 1);
+ Clamp(ey, 0, totalcells - 1);
+ Clamp(ez, 0, totalcells - 1);
+}
+
+int GridSpheres::mapIdx(int ix, int iy, int iz, int depth)
+{
+ if (depth==0) {
+ return (ix*ncells + iy)*ncells + iz;
+ } else {
+ int idx=mapIdx(static_cast<int>(ix*inv_ncells),
+ static_cast<int>(iy*inv_ncells),
+ static_cast<int>(iz*inv_ncells),
+ depth - 1);
+ int nx=ix%ncells;
+ int ny=iy%ncells;
+ int nz=iz%ncells;
+
+ return ((idx*ncells + nx)*ncells + ny)*ncells + nz;
+ }
+}
+
+void GridSpheres::fillMCell(MCell& mcell, int depth, int startidx) const
+{
+ mcell.nspheres=0;
+ mcell.min=new float[2*nvars];
+ mcell.max=mcell.min + nvars;
+
+ // Initialize min/max
+ for (int i=0; i<nvars; ++i) {
+ mcell.min[i]=FLT_MAX;
+ mcell.max[i]=-FLT_MAX;
+ }
+
+ // Determine min/max
+ int ncells3=ncells*ncells*ncells;
+ if (depth>0) {
+ MCell* mcells=macrocells[depth];
+ for (int i=0; i<ncells3; ++i) {
+ int idx=startidx + i;
+ fillMCell(mcells[idx], depth - 1, idx*ncells*ncells*ncells);
+ mcell.nspheres += mcells[idx].nspheres;
+ for (int j=0; j<nvars; ++j) {
+ if (mcells[idx].min[j]<mcell.min[j])
+ mcell.min[j]=mcells[idx].min[j];
+ if (mcells[idx].max[j]>mcell.max[j])
+ mcell.max[j]=mcells[idx].max[j];
+ }
+ }
+ } else {
+ for (int i=0; i<ncells3; ++i) {
+ int idx=startidx + i;
+ int nsph=counts[2*idx + 1];
+ mcell.nspheres += nsph;
+ int s=counts[2*idx];
+ for (int j=0; j<nsph; ++j) {
+ float* data=spheres + cells[s + j];
+ for (int k=0; k<nvars; ++k) {
+ if (data[k]<mcell.min[k])
+ mcell.min[k]=data[k];
+ if (data[k]>mcell.max[k])
+ mcell.max[k]=data[k];
+ }
+ }
+ }
+ }
+}
Added: trunk/Model/Primitives/GridSpheres.h
==============================================================================
--- (empty file)
+++ trunk/Model/Primitives/GridSpheres.h Tue May 30 08:41:00 2006
@@ -0,0 +1,92 @@
+
+#ifndef Manta_Model_GridSpheres_h
+#define Manta_Model_GridSpheres_h
+
+#include <Core/Color/Color.h>
+#include <Core/Geometry/BBox.h>
+#include <Interface/TexCoordMapper.h>
+#include <Interface/Texture.h>
+#include <Model/Primitives/PrimitiveCommon.h>
+#include <Model/Materials/LitMaterial.h>
+
+namespace Manta {
+ class GridSpheres : public PrimitiveCommon, public LitMaterial,
+ public TexCoordMapper
+ {
+ public:
+ GridSpheres(const Color& color, float* spheres, int nspheres,
+ int nvars, int ncells,
+ int depth, float radius, int ridx);
+ GridSpheres(const Texture<Color>* colorfn, float* spheres, int nspheres,
+ int nvars, int ncells,
+ int depth, float radius, int ridx);
+ ~GridSpheres(void);
+
+ void preprocess(const PreprocessContext&);
+
+ void computeBounds(const PreprocessContext& context,
+ BBox& bbox) const;
+ void intersect(const RenderContext& context, RayPacket& rays) const;
+ void computeNormal(const RenderContext& context, RayPacket& rays) const;
+
+ void shade(const RenderContext& context, RayPacket& rays) const;
+
+ void computeTexCoords2(const RenderContext& context,
+ RayPacket& rays) const;
+ void computeTexCoords3(const RenderContext& context,
+ RayPacket& rays) const;
+
+ private:
+ struct MCell {
+ int nspheres;
+ float* max;
+ float* min;
+ };
+
+ void init(void);
+
+ void traverse(int i, RayPacket& rays, int depth,
+ Real tnear,
+ int ix, int iy, int iz,
+ int idx,
+ Real dt_dx, Real dt_dy, Real dt_dz,
+ Real tnext_x, Real tnext_y, Real tnext_z,
+ const Vector& corner, const Vector& celldir,
+ int di_dx, int di_dy, int di_dz,
+ int didx_dx, int didx_dy,
+ int didx_dz, int stop_x, int stop_y, int stop_z) const;
+ void intersectSphere(RayPacket& rays, int i, int idx,
+ const Vector& center, float radius2) const;
+ void transformToLattice(const BBox& box, int& sx, int& sy, int& sz,
+ int& ex, int& ey, int& ez) const;
+ int mapIdx(int ix, int iy, int iz, int depth);
+ void fillMCell(MCell& mcell, int depth, int startidx) const;
+
+ float* spheres;
+ int nspheres;
+ int nvars;
+ Real radius;
+ Real inv_radius;
+ int ridx;
+
+ Vector diagonal;
+ Vector inv_diagonal;
+ int totalcells;
+
+ int* counts;
+ int* cells;
+ float* min;
+ float* max;
+ BBox bounds;
+
+ int ncells;
+ int depth;
+ Real inv_ncells;
+
+ MCell** macrocells;
+
+ const Texture<Color>* colortex;
+ };
+}
+
+#endif
Modified: trunk/Model/Readers/ParticleNRRD.cc
==============================================================================
--- trunk/Model/Readers/ParticleNRRD.cc (original)
+++ trunk/Model/Readers/ParticleNRRD.cc Tue May 30 08:41:00 2006
@@ -9,13 +9,14 @@
using namespace Manta;
using namespace std;
-ParticleNRRD::ParticleNRRD(void) :
- pdata(0), nvars(0), nparticles(0)
+ParticleNRRD::ParticleNRRD(bool nuke) :
+ pdata(0), nvars(0), nparticles(0), nuke(nuke)
{
// Do nothing
}
-ParticleNRRD::ParticleNRRD(string const fname)
+ParticleNRRD::ParticleNRRD(string const fname, bool nuke) :
+ pdata(0), nvars(0), nparticles(0), nuke(nuke)
{
if (!readFile(fname)) {
cerr<<"ParticleNRRD::ParticleNRRD: fatal error\n";
@@ -25,7 +26,7 @@
ParticleNRRD::~ParticleNRRD(void)
{
- if (pdata)
+ if (nuke && pdata)
delete pdata;
}
Modified: trunk/Model/Readers/ParticleNRRD.h
==============================================================================
--- trunk/Model/Readers/ParticleNRRD.h (original)
+++ trunk/Model/Readers/ParticleNRRD.h Tue May 30 08:41:00 2006
@@ -8,8 +8,8 @@
namespace Manta {
class ParticleNRRD {
public:
- ParticleNRRD(void);
- ParticleNRRD(string const fname);
+ ParticleNRRD(bool nuke=false);
+ ParticleNRRD(string const fname, bool nuke=false);
~ParticleNRRD(void);
unsigned int getNVars(void) const { return nvars; }
@@ -23,6 +23,7 @@
unsigned int nvars;
unsigned int nparticles;
float* pdata;
+ bool nuke;
};
}
Modified: trunk/scenes/pnrrd.cc
==============================================================================
--- trunk/scenes/pnrrd.cc (original)
+++ trunk/scenes/pnrrd.cc Tue May 30 08:41:00 2006
@@ -10,6 +10,7 @@
#include <Model/Groups/Group.h>
#include <Model/Lights/PointLight.h>
#include <Model/Materials/Lambertian.h>
+#include <Model/Primitives/GridSpheres.h>
#include <Model/Primitives/Sphere.h>
#include <Model/Readers/ParticleNRRD.h>
@@ -25,19 +26,30 @@
{
int argc=static_cast<int>(args.size());
string fname="";
+ int depth=1;
+ int ncells=2;
double radius=1.;
int ridx=-1;
Group* world=0;
for(int i=0; i<argc; ++i) {
string arg=args[i];
+#if 0
if (arg=="-bv") {
string s;
if (!getStringArg(i, args, s))
throw IllegalArgument("scene pnrrd -bv", i, args);
world=context.manta_interface->makeGroup(s);
+ } else
+#endif
+ if (arg=="-depth") {
+ if (!getIntArg(i, args, depth))
+ throw IllegalArgument("scene pnrrd -depth", i, args);
} else if (arg=="-i") {
if (!getStringArg(i, args, fname))
throw IllegalArgument("scene pnrrd -i", i, args);
+ } else if (arg=="-ncells") {
+ if (!getIntArg(i, args, ncells))
+ throw IllegalArgument("scene pnrrd -ncells", i, args);
} else if (arg=="-radius") {
if (!getDoubleArg(i, args, radius))
throw IllegalArgument("scene pnrrd -radius", i, args);
@@ -46,7 +58,9 @@
throw IllegalArgument("scene pnrrd -ridx", i, args);
} else {
cerr<<"Valid options for scene pnrrd:\n";
- cerr<<" -bv <string> bounding volume {bvh|grid|group}\n";
+ // cerr<<" -bv <string> bounding volume {bvh|grid|group}\n";
+ cerr<<" -depth <int> grid depth\n";
+ cerr<<" -ncells <int> grid resolution\n";
cerr<<" -i <string> filename\n";
cerr<<" -radius <float> particle radius\n";
cerr<<" -ridx <int> radius index\n";
@@ -59,6 +73,13 @@
// Load particle data, add particles to group
ParticleNRRD pnrrd(fname);
+#if 1
+ GridSpheres* spheres=new GridSpheres(Color(RGB(1, 1, 1)),
+ pnrrd.getParticleData(),
+ pnrrd.getNParticles(),
pnrrd.getNVars(),
+ ncells, depth, radius, ridx);
+ world->add(spheres);
+#else
float* pdata=pnrrd.getParticleData();
for (unsigned int i=0; i<pnrrd.getNParticles(); ++i) {
unsigned int idx=i*pnrrd.getNVars();
@@ -72,6 +93,7 @@
Material* white=new Lambertian(Color(RGB(1, 1, 1)));
world->add(new Sphere(white, Vector(x, y, z), radius));
}
+#endif
// Create scene
Scene* scene=new Scene();
@@ -83,6 +105,14 @@
lights->setAmbientLight(new ConstantAmbient(Color(RGB(0.4, 0.4, 0.4))));
scene->setLights(lights);
+ scene->addBookmark("debug 0", Vector(0.06, 26.9721, 0.06),
+ Vector(0.06, 0.06, 0.06), Vector(0, 0, 1), 0.59);
+ scene->addBookmark("debug 1", Vector(-10.3318,-2.88089, 24.71),
+ Vector(0.0812459, 0.0297618, 0.0654057),
+ Vector(0.00678671, -0.993405, -0.114458), 0.0232642);
+ scene->addBookmark("shadow error", Vector(-4.5545, 25.4447, 7.71331),
+ Vector(0.0534778, 0.0592496, 0.0578184),
+ Vector(-0.2524, -0.321088, 0.912796), 0.0877432);
scene->addBookmark("view 0", Vector(0.02, 1.02, 0.20),
Vector(0.02, 0.02, 0.20), Vector(0, 0, 1),
0.59);
- [MANTA] r1089 - in trunk: Model/Primitives Model/Readers scenes, cgribble, 05/30/2006
Archive powered by MHonArc 2.6.16.
