-
sparker
thiago
Manta Interactive Ray Tracer Development Mailing List
Text archives Help
- From: thiago@sci.utah.edu
- To: manta@sci.utah.edu
- Subject: [Manta] r1838 - in trunk: Model/Groups scenes
- Date: Thu, 8 Nov 2007 17:01:29 -0700 (MST)
Author: thiago
Date: Thu Nov 8 17:01:29 2007
New Revision: 1838
Added:
trunk/Model/Groups/KDTree.cc
trunk/Model/Groups/KDTree.h
Modified:
trunk/Model/Groups/CMakeLists.txt
trunk/scenes/triangleSceneViewer.cc
Log:
Added what should be a fully working KDTree. Note that the triangle
clipping for the build is not enabled since the code for that was not
compatible with the manta license. Writing our own version would give
around a 20% speedup in render performance (and a much slower build).
Modified: trunk/Model/Groups/CMakeLists.txt
==============================================================================
--- trunk/Model/Groups/CMakeLists.txt (original)
+++ trunk/Model/Groups/CMakeLists.txt Thu Nov 8 17:01:29 2007
@@ -5,6 +5,8 @@
Groups/GriddedGroup.h
Groups/Group.cc
Groups/Group.h
+ Groups/KDTree.h
+ Groups/KDTree.cc
Groups/Mesh.cc
Groups/Mesh.h
Groups/MovingMesh.cc
Added: trunk/Model/Groups/KDTree.cc
==============================================================================
--- (empty file)
+++ trunk/Model/Groups/KDTree.cc Thu Nov 8 17:01:29 2007
@@ -0,0 +1,986 @@
+#include <Model/Groups/KDTree.h>
+#include <Core/Exceptions/InternalError.h>
+#include <Core/Geometry/Vector.h>
+#include <Core/Math/MinMax.h>
+#include <Interface/Context.h>
+#include <SCIRun/Core/Thread/Time.h>
+#include <Interface/MantaInterface.h>
+
+#include <fstream>
+#include <iostream>
+#include <limits>
+
+#include <emmintrin.h>
+
+#define DBG(a)
+
+using namespace Manta;
+
+static const int ISEC_COST = 20;
+static const int TRAV_COST = 15;
+
+#ifdef MANTA_SSE
+#define SSE
+#endif
+
+#ifdef ALLOW_RESTRICT
+#define RESTRICT restrict
+#else
+#define RESTRICT /* doesn't exist on this compiler */
+#endif
+
+#define false4() cmp4_lt(zero4(),zero4())
+
+namespace Manta {
+
+// bool ClipTriangle(BBox &clipped,
+// const Vector &a,const Vector &b,const Vector&c,
+// const BBox &bounds);
+};
+
+void KDTree::preprocess(const PreprocessContext &context)
+{
+ if (currGroup) {
+ currGroup->computeBounds(context, bounds);
+ currGroup->preprocess(context);
+ }
+
+ context.manta_interface->registerSerialPreRenderCallback(
+ Callback::create(this, &KDTree::update));
+}
+
+void KDTree::setGroup(Group* new_group)
+{
+ currGroup = new_group;
+}
+
+
+void KDTree::rebuild(int proc, int numProcs)
+{
+ if (proc != 0)
+ return;
+
+ double startTime = SCIRun::Time::currentSeconds();
+
+ bounds.reset();
+
+ //Something's wrong, let's bail. Since we reset the bounds right
+ //above, this KDTree should never get intersected.
+ if (!currGroup)
+ return;
+
+ PreprocessContext context;
+ currGroup->computeBounds(context, bounds);
+
+ node.clear();
+
+ Node emptyNode;
+ node.push_back(emptyNode);
+
+ mesh = dynamic_cast<Mesh*>(currGroup);
+
+ vector<int> primitives;
+ for (unsigned int i=0; i < currGroup->size(); ++i) {
+ primitives.push_back(i);
+ }
+
+ cout <<"bounds are: "<<bounds[0] <<" , " <<bounds[1]<<endl;
+ // cout <<"bounds SA: "<<bounds.computeArea()<<endl;
+
+ build(0, primitives, bounds);
+
+ cout << "done building" << endl << flush;
+ float buildTime = SCIRun::Time::currentSeconds() - startTime;
+ printf("KDTree tree built in %f seconds\n", buildTime);
+ printStats();
+}
+
+struct Event
+{
+ enum Type { tri_end, tri_planar, tri_begin };
+ float pos;
+ Type type;
+ int tri;
+
+ Event(Type t, float p, int i) : pos(p), type(t), tri(i) {};
+};
+
+bool operator<(const Event &a, const Event &b)
+{
+ return a.pos < b.pos || a.pos == b.pos && a.type < b.type;
+}
+
+
+BBox clip(MeshTriangle *tri, const BBox &bounds)
+{
+ //Need to implement ClipTriangle
+ return BBox();
+
+ /*
+ BBox bb;
+ bb.reset();
+
+ Vector a = tri->getVertex(0);
+ Vector b = tri->getVertex(1);
+ Vector c = tri->getVertex(2);
+ ClipTriangle(bb,a,b,c,bounds);
+ return bb;
+*/
+};
+
+void KDTree::build(int nodeID,
+ vector<int> &primitiveID,
+ const BBox &bounds)
+{
+ DBG(cout << "build " << nodeID << " #=" << primitiveID.size() << " bb="
<< bounds[0] << ":" << bounds[1] << endl);
+ int bestDim = -1;
+ float bestSplit = -1;
+ float bestCost = primitiveID.size() * ISEC_COST;
+ bool bestCommonToTheLeft = false;
+
+ vector<Event> event[3];
+ if (bounds.computeArea() > 1e-5)
+ for (unsigned int i=0;i<primitiveID.size();i++) {
+ // bounds.reset();
+ // if (currGroup) {
+ // PreprocessContext context;
+ // currGroup->computeBounds(context, bounds);
+
+ BBox box;
+
+ //until someone implements triangle clipping, can't call clip()!
+ if (false && mesh)
+ box = clip(mesh->get(primitiveID[i]),bounds);
+ else {
+ PreprocessContext context;
+ currGroup->get(primitiveID[i])->computeBounds(context, box);
+ }
+ // if (box.computeArea() <= 1e-4)
+ // continue;
+
+ if (box.getMin() == box.getMax()) {
+ cerr << "zero triangle after clipping -- should not happen ..." <<
endl;
+ cerr << " nodeID " << nodeID << endl;
+ if (mesh) {
+ cerr << " " << mesh->get(primitiveID[i])->getVertex(0) << endl;
+ cerr << " " << mesh->get(primitiveID[i])->getVertex(1) << endl;
+ cerr << " " << mesh->get(primitiveID[i])->getVertex(2) << endl;
+ }
+ cerr << endl;
+ cerr << bounds.getMin() << " " << bounds.getMax() << endl;
+ cerr << box.getMin() << " " << box.getMax() << endl;
+ exit(0);
+ }
+ else {
+ for (int k=0;k<3;k++) {
+ if (box.getMin()[k] == box.getMax()[k])
+ event[k].push_back(Event(Event::tri_planar,box.getMin()[k],i));
+ else {
+ event[k].push_back(Event(Event::tri_begin,box.getMin()[k],i));
+ event[k].push_back(Event(Event::tri_end, box.getMax()[k],i));
+ }
+ }
+ }
+ }
+ // cout << __PRETTY_FUNCTION__ << endl;
+ const float boundsArea = bounds.computeArea();
+ for (int k=0;k<3;k++) {
+ BBox lBounds = bounds;
+ BBox rBounds = bounds;
+
+ DBG(cout << "sorting " << k << " " << event[k].size() << endl);
+ sort(event[k].begin(),event[k].end());
+
+ int Nl = 0;
+ int Nr = primitiveID.size();
+ int Np = 0;
+
+ unsigned int seqStart = 0;
+ while (seqStart < event[k].size()) {
+ unsigned int idx = seqStart;
+
+ int numEnding = 0;
+ while (idx < event[k].size() &&
+ event[k][idx].pos == event[k][seqStart].pos &&
+ event[k][idx].type == Event::tri_end)
+ { ++idx; ++numEnding; }
+
+ int numPlanar = 0;
+ while (idx < event[k].size() &&
+ event[k][idx].pos == event[k][seqStart].pos &&
+ event[k][idx].type == Event::tri_planar)
+ { ++idx; ++numPlanar; }
+
+ int numStarting = 0;
+ while (idx < event[k].size() &&
+ event[k][idx].pos == event[k][seqStart].pos &&
+ event[k][idx].type == Event::tri_begin)
+ { ++idx; ++numStarting; }
+
+ DBG(cout << "event " << event[k][seqStart].pos << "@" << k << endl;
+ cout << " " << numStarting << " " << numPlanar << " " <<
numEnding << endl);
+ lBounds[1][k] = event[k][seqStart].pos;
+ rBounds[0][k] = event[k][seqStart].pos;
+
+ float lProb = lBounds.computeArea() / boundsArea;
+ float rProb = rBounds.computeArea() / boundsArea;
+
+ Nr -= numEnding;
+ Nr -= numPlanar;
+
+ Np = numPlanar;
+
+ DBG(cout << "prob " << lProb << " " << rProb << endl;
+ cout << " N " << Nl << " " << Np << " " << Nr << endl);
+
+ {
+ // try putting common on the left
+ float cost = TRAV_COST + ISEC_COST * (lProb * (Nl+Np) + rProb * Nr);
+
+ DBG(cout << cost << endl);
+ if (cost < bestCost) {
+ bestCost = cost;
+ bestCommonToTheLeft = true;
+ bestDim = k;
+ bestSplit = event[k][seqStart].pos;
+ DBG(cout << "new favorite " << bestDim << " " << bestSplit <<
endl);
+ }
+ }
+ {
+ // try putting common on the right
+ float cost = TRAV_COST + ISEC_COST * (lProb * Nl + rProb * (Nr+Np));
+ DBG(cout << cost << endl);
+ if (cost < bestCost) {
+ bestCost = cost;
+ bestCommonToTheLeft = false;
+ bestDim = k;
+ bestSplit = event[k][seqStart].pos;
+ DBG(cout << "new favorite " << bestDim << " " << bestSplit <<
endl);
+ }
+ }
+
+ Nl += numPlanar;
+ Nl += numStarting;
+
+ seqStart = idx;
+ }
+ }
+ DBG(cout << "node " << nodeID << " bestSplit " << bestSplit << " dim " <<
bestDim << endl);
+ if (bestSplit == -1) {
+ // no split found ...
+ makeLeaf(nodeID,primitiveID);
+ DBG(cout << "leaf : " << nodeID << " " << primitiveID.size() << endl);
+ }
+ else {
+ vector<int> lPrimitive;
+ vector<int> rPrimitive;
+
+ bool *is_l = new bool[primitiveID.size()];
+ bool *is_r = new bool[primitiveID.size()];
+
+ for (unsigned int i=0;i<primitiveID.size();i++)
+ is_l[i] = is_r[i] = true;
+
+ for (unsigned int i=0;i<event[bestDim].size();i++)
+ {
+ switch (event[bestDim][i].type) {
+ case Event::tri_planar:
+ if (event[bestDim][i].pos == bestSplit && bestCommonToTheLeft ||
+ event[bestDim][i].pos < bestSplit)
+ is_r[event[bestDim][i].tri] = false;
+ else
+ is_l[event[bestDim][i].tri] = false;
+ break;
+ case Event::tri_begin:
+ if (event[bestDim][i].pos >= bestSplit)
+ is_l[event[bestDim][i].tri] = false;
+ break;
+ case Event::tri_end:
+ if (event[bestDim][i].pos <= bestSplit)
+ is_r[event[bestDim][i].tri] = false;
+ break;
+ };
+ }
+ for (unsigned int i=0;i<primitiveID.size();i++) {
+ if (is_l[i]) lPrimitive.push_back(primitiveID[i]);
+ if (is_r[i]) rPrimitive.push_back(primitiveID[i]);
+ }
+ delete[] is_l;
+ delete[] is_r;
+
+ primitiveID.clear();
+ makeInner(nodeID,bestDim,bestSplit);
+ BBox lBounds = bounds;
+ BBox rBounds = bounds;
+ lBounds[1][bestDim] = bestSplit;
+ rBounds[0][bestDim] = bestSplit;
+ build(node[nodeID].childIdx+0,lPrimitive,lBounds);
+ build(node[nodeID].childIdx+1,rPrimitive,rBounds);
+ }
+};
+
+
+void KDTree::intersect(const RenderContext& context, RayPacket& rays) const
+{
+ rays.normalizeDirections();
+ rays.computeSigns();
+ rays.computeInverseDirections();
+
+ const int ray_begin = rays.begin();
+ const int ray_end = rays.end();
+
+ if (!rays.getFlag(RayPacket::ConstantSigns)) {
+ for (int first=ray_begin; first < ray_end; ++first) {
+ int last = first+1;
+ for (; last < ray_end; ++last) {
+ if (rays.getSign(first, 0) != rays.getSign(last, 0) ||
+ rays.getSign(first, 1) != rays.getSign(last, 1) ||
+ rays.getSign(first, 2) != rays.getSign(last, 2)) {
+ break;
+ }
+ }
+ RayPacket sameSignSubpacket(rays, first, last);
+ sameSignSubpacket.setFlag(RayPacket::ConstantSigns);
+ sameSignSubpacket.resetFlag(RayPacket::HaveCornerRays);
+ intersect(context, sameSignSubpacket);
+ first = last-1;
+ }
+ return;
+ }
+
+ const bool COMMON_ORIGIN = rays.getFlag(RayPacket::ConstantOrigin);
+
+#ifdef MAILBOX
+ Mailbox mailbox;
+ mailbox.clear();
+#endif
+
+#ifdef COLLECT_STATS
+ nTotalRays += rays.end() - rays.begin();
+#endif
+
+ MANTA_ALIGN(16) float t_in[RayPacketData::MaxSize+1]; //last element is
min of all t_in
+ MANTA_ALIGN(16) float t_out[RayPacketData::MaxSize+1];//last element is
max of all t_out
+ MANTA_ALIGN(16) int valid[RayPacketData::MaxSize];
+
+#ifdef SSE
+ const int sse_begin = (ray_begin + 3) & (~3);
+ const int sse_end = (ray_end) & (~3);
+
+ sse_t anyValid = false4();
+
+ const int sign_x = rays.getSign(ray_begin,0);
+ const int sign_y = rays.getSign(ray_begin,1);
+ const int sign_z = rays.getSign(ray_begin,2);
+ const sse_t near_x = set4(bounds[sign_x][0]);
+ const sse_t near_y = set4(bounds[sign_y][1]);
+ const sse_t near_z = set4(bounds[sign_z][2]);
+ const sse_t far_x = set4(bounds[1-sign_x][0]);
+ const sse_t far_y = set4(bounds[1-sign_y][1]);
+ const sse_t far_z = set4(bounds[1-sign_z][2]);
+
+ int numRays = 0;
+ if (ray_begin < sse_begin || sse_end < ray_end) {
+ Vector origin;
+
+ if (COMMON_ORIGIN)
+ origin = rays.getOrigin(ray_begin);
+
+ bool leading = true;
+ for (int ray = ray_begin; /*no test*/ ; ++ray) {
+ if (leading && ray >= sse_begin) {
+ leading = false;
+ ray = sse_end;
+ }
+ if (!leading && ray >= ray_end)
+ break;
+
+ //do bounding box test
+
+ if (!COMMON_ORIGIN)
+ origin = rays.getOrigin( ray );
+ const Vector inverse_direction = rays.getInverseDirection( ray );
+
+ t_in[ray] = T_EPSILON;
+ t_out[ray] = rays.getMinT(ray);
+
+ for (int k=0;k<3;k++) {
+ //TODO: Check to see if using rays.getOrigin(0,k) and
+ //getInverseDirection(ray,k) are faster.
+ float t0 = (bounds.getMin()[k] - origin[k]) * inverse_direction[k];
+ float t1 = (bounds.getMax()[k] - origin[k]) * inverse_direction[k];
+
+ if (t0 > t1) {
+ if (t1 > t_in[ray]) t_in[ray] = t1;
+ if (t0 < t_out[ray]) t_out[ray] = t0;
+ } else {
+ if (t0 > t_in[ray]) t_in[ray] = t0;
+ if (t1 < t_out[ray]) t_out[ray] = t1;
+ }
+ if (t_in[ray] > t_out[ray])
+ break;
+ }
+ valid[ray] = (t_in[ray] <= t_out[ray]);
+ if (valid[ray])
+ numRays++;
+ }
+ }
+
+ sse_t org_x = load44(&rays.getOrigin(sse_begin,0));
+ sse_t org_y = load44(&rays.getOrigin(sse_begin,1));
+ sse_t org_z = load44(&rays.getOrigin(sse_begin,2));
+
+#ifdef __INTEL_COMPILER
+#pragma unroll(4)
+#endif
+ for (int ray = sse_begin; ray < sse_end; ray+=4) {
+ if (!COMMON_ORIGIN) {
+ org_x = load44(&rays.getOrigin(ray,0));
+ org_y = load44(&rays.getOrigin(ray,1));
+ org_z = load44(&rays.getOrigin(ray,2));
+ }
+
+ sse_t t0 = set4(T_EPSILON);
+ sse_t t1 = load44(&rays.getMinT(ray));
+ const sse_t t0x =
mul4(sub4(near_x,org_x),load44(&rays.getInverseDirection(ray,0)));
+ t0 = max4(t0,t0x);
+ const sse_t t1x = mul4(sub4(far_x,
org_x),load44(&rays.getInverseDirection(ray,0)));
+ t1 = min4(t1,t1x);
+ const sse_t t0y =
mul4(sub4(near_y,org_y),load44(&rays.getInverseDirection(ray,1)));
+ t0 = max4(t0,t0y);
+ const sse_t t1y = mul4(sub4(far_y,
org_y),load44(&rays.getInverseDirection(ray,1)));
+ t1 = min4(t1,t1y);
+ const sse_t t0z =
mul4(sub4(near_z,org_z),load44(&rays.getInverseDirection(ray,2)));
+ t0 = max4(t0,t0z);
+ const sse_t t1z = mul4(sub4(far_z,
org_z),load44(&rays.getInverseDirection(ray,2)));
+ t1 = min4(t1,t1z);
+
+ store44(&t_in[ray], t0);
+ store44(&t_out[ray], t1);
+ store44(reinterpret_cast<float*>(&valid[ray]), cmp4_le(t0,t1));
+ anyValid = or4(anyValid,load44(reinterpret_cast<float*>(&valid[ray])));
+ }
+
+ // cout << __PRETTY_FUNCTION__ << endl;
+ if (getmask4(anyValid) || numRays > 0) {
+ if (COMMON_ORIGIN)
+ intersectNode<true>(0, context, rays,
+ t_in,t_out,valid
+#ifdef MAILBOX
+ , mailbox
+#endif
+ );
+ else
+ intersectNode<false>(0, context, rays,
+ t_in,t_out,valid
+#ifdef MAILBOX
+ , mailbox
+#endif
+ );
+
+ }
+#else
+ int numRays = 0;
+ // cout << "new rays -> " << (ray_end - rays.begin()) << endl;
+ for(int ray=ray_begin; ray<ray_end; ray++) {
+
+ //do bounding box test
+ Vector origin = rays.getOrigin( ray );
+ Vector inverse_direction = rays.getInverseDirection( ray );
+
+ t_in[ray] = T_EPSILON;
+ t_out[ray] = rays.getMinT(ray);
+
+ for (int k=0;k<3;k++) {
+ float t0 = (bounds.getMin()[k] - origin[k]) * inverse_direction[k];
+ float t1 = (bounds.getMax()[k] - origin[k]) * inverse_direction[k];
+
+ // cout << " k=" << k << " : " << t0 << " " << t1 << endl;
+ // hopefully everything all right with nan's etc....
+ if (t0 > t1) {
+ if (t1 > t_in[ray]) t_in[ray] = t1;
+ if (t0 < t_out[ray]) t_out[ray] = t0;
+ } else {
+ if (t0 > t_in[ray]) t_in[ray] = t0;
+ if (t1 < t_out[ray]) t_out[ray] = t1;
+ }
+ if (t_in[ray] > t_out[ray])
+ break;
+ }
+ valid[ray] = (t_in[ray] <= t_out[ray]);
+ if (valid[ray])
+ numRays++;
+ }
+ // cout << "num rays " << numRays << endl;
+ if (numRays > 0) {
+
+ if (COMMON_ORIGIN)
+ intersectNode<true>(0, context, rays,
+ t_in,t_out,valid
+#ifdef MAILBOX
+ , mailbox
+#endif
+ );
+ else
+ intersectNode<false>(0, context, rays,
+ t_in,t_out,valid
+#ifdef MAILBOX
+ , mailbox
+#endif
+ );
+ }
+#endif
+}
+
+template<bool COMMON_ORIGIN>
+void KDTree::intersectNode(int nodeID, const RenderContext& context,
+ RayPacket& rays,
+ const float *const t_in,
+ const float *const t_out,
+ const int *const valid
+#ifdef MAILBOX
+ , Mailbox &mailbox
+#endif
+ ) const
+{
+#ifdef COLLECT_STATS
+ ++nTraversals;
+#endif
+
+ const Node &node = this->node[nodeID];
+ if (node.isLeaf) {
+
+ // cout << "LEAF " << node.numPrimitives << endl;
+ int primOffset = node.childIdx;
+ for (int i=0; i<node.numPrimitives; ++i) {
+ // printf("testing object from node %d\n", nodeID);
+ const int triID = itemList[primOffset+i];
+#ifdef MAILBOX
+ if (mailbox.testAndMark(triID))
+ continue;
+#endif
+
+ currGroup->get(triID)->intersect(context, rays);
+#ifdef COLLECT_STATS
+ ++nIntersects;
+#endif
+ }
+ } else {
+
+ const int ray_begin = rays.begin();
+ const int ray_end = rays.end();
+
+#ifdef SSE
+ const int sse_begin = (ray_begin + 3) & (~3);
+ const int sse_end = (ray_end) & (~3);
+#endif
+
+ const int sign = rays.getSign(ray_begin,node.planeDim);
+
+#ifdef SSE
+
+ MANTA_ALIGN(16) float t_plane[RayPacketData::MaxSize];
+
+ sse_t frontMask = false4();
+ sse_t backMask = frontMask;
+
+ const sse_t planePos4 = set4(node.planePos);
+ int front=0;
+ int back=0;
+ if (ray_begin < sse_begin || sse_end < ray_end) {
+ bool leading = true;
+ for (int ray = ray_begin; /*no test*/ ; ++ray) {
+ if (leading && ray >= sse_begin) {
+ leading = false;
+ ray = sse_end;
+ }
+ if (!leading && ray >= ray_end)
+ break;
+
+ t_plane[ray]
+ = (node.planePos - rays.getOrigin(ray,node.planeDim)) //TODO:
optimize for common origin
+ * rays.getInverseDirection(ray,node.planeDim);
+ if (valid[ray]) {
+ if (t_out[ray] < t_plane[ray])
+ front++;
+ else if (t_in[ray] > t_plane[ray])
+ back++;
+ else {
+ front++;
+ back++;
+ }
+ }
+ }
+ }
+
+ sse_t org_k = load44(&rays.getOrigin(sse_begin,node.planeDim));
+ sse_t plane_org = sub4(planePos4, org_k);
+
+#ifdef __INTEL_COMPILER
+#pragma unroll(4)
+#endif
+ for (int ray = sse_begin; ray < sse_end; ray+=4) {
+ if (!COMMON_ORIGIN) {
+ org_k = load44(&rays.getOrigin(ray,node.planeDim));
+ plane_org = sub4(planePos4, org_k);
+ }
+ const sse_t rcp = load44(&rays.getInverseDirection(ray,node.planeDim));
+ const sse_t tp = mul4(plane_org,rcp);
+ store44(&t_plane[ray], tp);
+
+ const sse_t frontOnly = cmp4_lt(load44(&t_out[ray]),tp);
+ backMask = or4(backMask, andnot4(frontOnly,
+ load44(reinterpret_cast<const
float*>(&valid[ray]))));
+ const sse_t backOnly = cmp4_gt(load44(&t_in[ray]),tp);
+ frontMask = or4(frontMask, andnot4(backOnly,
+ load44(reinterpret_cast<const
float*>(&valid[ray]))));
+ }
+
+// float *t_plane = (float *)&t_plane4[0]; // to tie it in w/ code
below...
+ front += getmask4(frontMask);
+ back += getmask4(backMask);
+#else
+ float t_plane[RayPacketData::MaxSize];
+ int front=0;
+ int back=0;
+ for(int ray=ray_begin; ray<ray_end; ray++) {
+ t_plane[ray]
+ = (node.planePos - rays.getOrigin(ray,node.planeDim)) //TODO:
optimize for common origin
+ * rays.getInverseDirection(ray,node.planeDim);
+ if (valid[ray]) {
+ if (t_out[ray] < t_plane[ray])
+ front++;
+ else if (t_in[ray] > t_plane[ray])
+ back++;
+ else {
+ front++;
+ back++;
+ }
+ }
+ }
+#endif
+
+ const int backChildID = node.childIdx+1-sign;
+ const int frontChildID = node.childIdx+sign;
+
+ if (front == 0) {
+ intersectNode<COMMON_ORIGIN>(backChildID,context,rays,t_in,t_out,valid
+#ifdef MAILBOX
+ , mailbox
+#endif
+);
+ } else if (back == 0) {
+ intersectNode<COMMON_ORIGIN>(frontChildID,context,rays,t_in,t_out,valid
+#ifdef MAILBOX
+ , mailbox
+#endif
+);
+ } else {
+
+ MANTA_ALIGN(16) float new_t_out[RayPacketData::MaxSize+1];
+ MANTA_ALIGN(16) int new_valid[RayPacketData::MaxSize];
+
+#ifdef SSE
+ if (ray_begin < sse_begin || sse_end < ray_end) {
+ bool leading = true;
+ for (int ray = ray_begin; /*no test*/ ; ++ray) {
+ if (leading && ray >= sse_begin) {
+ leading = false;
+ ray = sse_end;
+ }
+ if (!leading && ray >= ray_end)
+ break;
+
+ if (!valid[ray]) {
+ new_valid[ray] = false;
+ } else {
+ if (t_in[ray] > t_plane[ray]) {
+ new_valid[ray] = false;
+ }
+ else {
+ new_valid[ray] = true;
+ new_t_out[ray] = SCIRun::Min(t_out[ray], t_plane[ray]);
+ }
+ }
+ }
+ }
+#ifdef __INTEL_COMPILER
+#pragma unroll(4)
+#endif
+ for (int ray = sse_begin; ray < sse_end; ray+=4) {
+ //no need to do t_in here
+ const sse_t nto = min4(load44(&t_out[ray]),load44(&t_plane[ray]));
+ store44(&new_t_out[ray], nto);
+ store44(reinterpret_cast<float*>(&new_valid[ray]),
+ and4(load44(reinterpret_cast<const float*>(&valid[ray])),
+ cmp4_le(load44(&t_in[ray]),nto)));
+ }
+
+ intersectNode<COMMON_ORIGIN>(frontChildID,context,rays,
+ t_in,new_t_out,new_valid
+#ifdef MAILBOX
+ , mailbox
+#endif
+ );
+#else
+
+ // both ... that's the tricky part...
+ for(int ray=ray_begin; ray<ray_end; ray++) {
+ if (!valid[ray]) {
+ new_valid[ray] = false;
+ } else {
+ if (t_in[ray] > t_plane[ray]) {
+ new_valid[ray] = false;
+ }
+ else {
+ new_valid[ray] = true;
+ new_t_out[ray] = SCIRun::Min(t_out[ray], t_plane[ray]);
+ }
+ }
+ }
+
+ intersectNode<COMMON_ORIGIN>(frontChildID,context,rays,
+ t_in,new_t_out,new_valid
+#ifdef MAILBOX
+ , mailbox
+#endif
+ );
+#endif //SSE
+
+ MANTA_ALIGN(16) float new_t_in[RayPacketData::MaxSize+1];
+ back = 0;
+
+#ifdef SSE
+ if (ray_begin < sse_begin || sse_end < ray_end) {
+ bool leading = true;
+ for (int ray = ray_begin; /*no test*/ ; ++ray) {
+ if (leading && ray >= sse_begin) {
+ leading = false;
+ ray = sse_end;
+ }
+ if (!leading && ray >= ray_end)
+ break;
+
+ if (!valid[ray]) {
+ new_valid[ray] = false;
+ } else {
+ if (t_in[ray] > t_plane[ray]) {
+ new_t_in[ray] = t_in[ray];
+ // might may have changed in prev isec step:
+ new_t_out[ray] = SCIRun::Min(t_out[ray],rays.getMinT(ray));
+ new_valid[ray] = (new_t_in[ray]<=new_t_out[ray]);
+ } else if (t_out[ray] < t_plane[ray]) {
+ new_valid[ray] = false;
+ } else {
+ new_t_in[ray] = t_plane[ray]; //t_in[ray];
+ new_t_out[ray] = SCIRun::Min(t_out[ray],rays.getMinT(ray));
+ new_valid[ray] = (new_t_in[ray]<=new_t_out[ray]);
+ }
+ back += new_valid[ray];
+ }
+ }
+ }
+
+
+ backMask = false4();
+#ifdef __INTEL_COMPILER
+#pragma unroll(4)
+#endif
+ for (int ray = sse_begin; ray < sse_end; ray+=4) {
+ const sse_t nti = max4(load44(&t_in[ray]),load44(&t_plane[ray]));
+ store44(&new_t_in[ray], nti);
+ const sse_t nto =
min4(load44(&t_out[ray]),load44(&rays.getMinT(ray)));
+ store44(&new_t_out[ray], nto);
+ const sse_t nvi = and4(load44(reinterpret_cast<const
float*>(&valid[ray])),
+ cmp4_le(nti,nto));
+ store44(reinterpret_cast<float*>(&new_valid[ray]), nvi);
+
+ backMask =
or4(backMask,load44(reinterpret_cast<float*>(&new_valid[ray])));
+ }
+ back += getmask4(backMask);
+
+ if (back) {
+#else
+
+ // now, clip to back plane ...
+ for(int ray=ray_begin; ray<ray_end; ray++) {
+ if (!valid[ray]) {
+ new_valid[ray] = false;
+ } else {
+ if (t_in[ray] > t_plane[ray]) {
+ new_t_in[ray] = t_in[ray];
+ // might may have changed in prev isec step:
+ new_t_out[ray] = std::min(t_out[ray],rays.getMinT(ray));
+ new_valid[ray] = (new_t_in[ray]<=new_t_out[ray]);
+ } else if (t_out[ray] < t_plane[ray]) {
+ new_valid[ray] = false;
+ } else {
+ new_t_in[ray] = t_plane[ray]; //t_in[ray];
+ new_t_out[ray] = std::min(t_out[ray],rays.getMinT(ray));
+ new_valid[ray] = (new_t_in[ray]<=new_t_out[ray]);
+ }
+ back += new_valid[ray];
+ }
+ }
+
+ if (back) {
+#endif //SSE
+
+ intersectNode<COMMON_ORIGIN>(backChildID,context,rays,
+ new_t_in,new_t_out,new_valid
+#ifdef MAILBOX
+ , mailbox
+#endif
+ );
+ }
+ }
+ }
+}
+
+
+
+/**
+ * Test whether any of the rays in a ray packet intersects this
+ * bounding box.
+ *
+ * @param rays[in] the set of rays to check.
+ * @return true if any rays touch, false if all of them miss.
+ */
+bool KDTree::intersectBounds(const RenderContext& context, RayPacket& rays)
const
+{
+ throw SCIRun::InternalError("hopefully done need this function", __FILE__,
__LINE__);
+ return false;
+}
+
+
+bool KDTree::buildFromFile(const string &file)
+{
+ double startTime = SCIRun::Time::currentSeconds();
+
+ bounds.reset();
+
+ //Something's wrong, let's bail. Since we reset the bounds right
+ //above, this KDTree should never get intersected.
+ if (!currGroup)
+ return false;
+
+ PreprocessContext context;
+ currGroup->computeBounds(context, bounds);
+
+ node.clear();
+
+ mesh = dynamic_cast<Mesh*>(currGroup);
+
+ ifstream in(file.c_str());
+ if (!in) return false;
+ unsigned int itemListSize;
+ in >> itemListSize;
+
+ for (unsigned int i=0; i < itemListSize; ++i) {
+ int item;
+ in >> item;
+ if (!in) return false;
+ itemList.push_back(item);
+ }
+
+ while (in.good()) {
+ Node n;
+ int i;
+ unsigned int u;
+
+ in >> i;
+ n.isLeaf = i;
+
+ if (in.eof())
+ break;
+
+ if (n.isLeaf) {
+ in >> n.numPrimitives;
+ in >> u;
+ n.childIdx = u;
+ }
+ else {
+ int type;
+ in >> type;
+
+ if (type != 0)
+ return false; //can only read KDTree format.
+
+ float f;
+ in >> f;
+ n.planePos = f;
+
+
+ in >> u;
+ n.planeDim = u;
+ in >> u;
+ n.childIdx = u;
+ }
+ node.push_back(n);
+ }
+
+ in.close();
+
+ cout << "done building" << endl << flush;
+ float buildTime = SCIRun::Time::currentSeconds() - startTime;
+ printf("KDTree tree built in %f seconds\n", buildTime);
+ printStats();
+
+ return true;
+}
+
+
+void KDTree::saveToFile(const string &file)
+{
+ ofstream out(file.c_str());
+ if (!out) return;
+
+ out << itemList.size() << endl;
+ for (unsigned int i=0; i < itemList.size(); ++i) {
+ out << itemList[i] <<" ";
+ }
+
+ out<<endl;
+
+ for (unsigned int i=0; i < node.size(); ++i) {
+ out << node[i].isLeaf << " ";
+ if (node[i].isLeaf) {
+ out << node[i].numPrimitives << " " << node[i].childIdx <<endl;
+ }
+ else {
+ out << 0 << " " << node[i].planePos << " " << node[i].planeDim << " "
<< node[i].childIdx<<endl;
+ }
+ }
+ out.close();
+}
+
+void KDTree::printStats()
+{
+ TreeStats stats = {0};
+ collectTreeStats(0, 0, stats);
+
+ printf("KDtree made from %d primitives has: \n", (int) currGroup->size());
+ printf(" - %d nodes\n", (int)node.size());
+ printf(" - %d primitive references\n", (int)itemList.size());
+ printf(" - %d max depth\n", stats.maxDepth);
+ printf(" - %d max objects in a leaf\n", stats.maxObjectsInLeaf);
+
+ for (int i=0; i<1024; ++i)
+ if (stats.childrenLeafHist[i] > 0)
+ printf(" - %d leaves with %d children\n", stats.childrenLeafHist[i],
i);
+}
+
+void KDTree::collectTreeStats(int nodeID, int depth, TreeStats &stats)
+{
+ const Node &node = this->node[nodeID];
+ if (node.isLeaf) {
+ stats.maxDepth = SCIRun::Max(stats.maxDepth, depth);
+ stats.maxObjectsInLeaf = SCIRun::Max(stats.maxObjectsInLeaf,
+ node.numPrimitives);
+ stats.childrenLeafHist[std::min(node.numPrimitives,1023)]++;
+ }
+ else {
+ collectTreeStats(node.childIdx+0, depth+1, stats);
+ collectTreeStats(node.childIdx+1, depth+1, stats);
+ }
+}
Added: trunk/Model/Groups/KDTree.h
==============================================================================
--- (empty file)
+++ trunk/Model/Groups/KDTree.h Thu Nov 8 17:01:29 2007
@@ -0,0 +1,166 @@
+#ifndef Manta_Model_Groups_KDTree_h
+#define Manta_Model_Groups_KDTree_h
+
+#include <Core/Geometry/BBox.h>
+#include <Core/Geometry/Vector.h>
+#include <Interface/AccelerationStructure.h>
+#include <Interface/MeshTriangle.h>
+#include <Interface/RayPacket.h>
+#include <Model/Groups/Group.h>
+#include <Model/Groups/Mesh.h>
+
+#include <assert.h>
+
+namespace Manta
+{
+
+ class KDTree : public AccelerationStructure
+ {
+ public:
+
+
+#ifdef MANTA_SSE
+# define MAILBOX
+
+// must be a power of two...
+# define MAILBOX_ENTRIES 64
+# define MAILBOX_MASK (MAILBOX_ENTRIES-1)
+
+ struct Mailbox
+ {
+ MANTA_ALIGN(16) int triID[MAILBOX_ENTRIES];
+
+ inline void clear()
+ {
+ sse_int_t m1 = set4i(-1);
+#ifdef __INTEL_COMPILER
+#pragma unroll(4)
+#endif
+ for (int i=0;i<MAILBOX_ENTRIES;i+=4)
+ store44i((sse_int_t*)(triID+i),m1);
+ }
+ inline bool testAndMark(int id)
+ {
+ const int hash = id & MAILBOX_MASK;
+ if (triID[hash] == id)
+ return true;
+ triID[hash] = id;
+ return false;
+ }
+ };
+#endif
+
+
+ struct Node
+ {
+ union {
+ float planePos;
+ int numPrimitives;
+ };
+ unsigned int isLeaf : 1;
+ unsigned int planeDim : 2;
+ unsigned int childIdx : 29;
+ };
+ vector <Node> node;
+ vector <int> itemList;
+
+ void makeLeaf(int nodeID, vector<int> &prim)
+ {
+ node[nodeID].isLeaf = 1;
+ node[nodeID].childIdx = itemList.size();
+ node[nodeID].numPrimitives = prim.size();
+ copy(prim.begin(),prim.end(),back_inserter(itemList));
+ }
+ void makeInner(int nodeID,
+ int dim, float pos)
+ {
+ node[nodeID].isLeaf = 0;
+ node[nodeID].childIdx = node.size();
+ node[nodeID].planePos = pos;
+ node[nodeID].planeDim = dim;
+ node.push_back(Node());
+ node.push_back(Node());
+ }
+
+ BBox bounds;
+
+ Group *currGroup;
+ Mesh *mesh;
+
+ //for performance measurements. Can be removed.
+ mutable long nIntersects;
+ mutable long nTraversals;
+ mutable long nTotalRays;
+ mutable long nFrontIABranch;
+ mutable long nBackIABranch;
+
+ void update(int proc, int numProcs)
+ {
+// #define COLLECT_STATS
+ #ifdef COLLECT_STATS
+ printf("primitive intersections per ray: %f\n",
(float)nIntersects/nTotalRays);
+ printf("node traversals per ray: %f\n",
(float)nTraversals/nTotalRays);
+ printf("IA front branches taken per ray: %f\n",
(float)nFrontIABranch/nTotalRays);
+ printf("IA back branches taken per ray: %f\n",
(float)nBackIABranch/nTotalRays);
+ printf("Traversals taken by IA: %f%%\n",
100*(float)(nFrontIABranch+nBackIABranch)/nTraversals);
+ printf("number of rays: %ld\n", nTotalRays);
+ nIntersects = 0;
+ nTraversals = 0;
+ nFrontIABranch = 0;
+ nBackIABranch = 0;
+ nTotalRays = 0;
+ #endif
+ }
+
+ KDTree() : currGroup(NULL), mesh(NULL), nIntersects(0), nTraversals(0),
nTotalRays(0)
+ {}
+
+ void preprocess(const PreprocessContext&);
+ void intersect(const RenderContext& context, RayPacket& rays) const;
+
+ template<bool COMMON_ORIGIN>
+ void intersectNode(int nodeID,
+ const RenderContext& context,
+ RayPacket& rays,
+ const float *const t_in,
+ const float *const t_out,
+ const int *const valid
+#ifdef MAILBOX
+ , Mailbox &mailbox
+#endif
+ ) const;
+ bool intersectBounds(const RenderContext& context, RayPacket& rays)
const;
+
+ void computeBounds(const PreprocessContext&,
+ BBox& bbox) const
+ {
+ if (!node.empty())
+ bbox.extendByBox(bounds);
+ }
+
+ void setGroup(Group* new_group);
+ Group* getGroup() const { return currGroup; }
+
+ virtual void addToUpdateGraph(ObjectUpdateGraph* graph,
+ ObjectUpdateGraphNode* parent) { }
+
+ void rebuild(int proc=0, int numProcs=1);
+
+ void build(int nodeID,
+ vector<int> &primitives,
+ const BBox &bounds);
+
+ bool buildFromFile(const string &file);
+ void saveToFile(const string &file);
+
+ void printStats();
+ struct TreeStats{
+ int maxDepth;
+ int maxObjectsInLeaf;
+ int childrenLeafHist[1024]; //hopefully 1024 is more than enough
+ };
+ void collectTreeStats(int nodeID, int depth, TreeStats &stats);
+ };
+};
+
+#endif
Modified: trunk/scenes/triangleSceneViewer.cc
==============================================================================
--- trunk/scenes/triangleSceneViewer.cc (original)
+++ trunk/scenes/triangleSceneViewer.cc Thu Nov 8 17:01:29 2007
@@ -9,6 +9,7 @@
#include <Model/Backgrounds/ConstantBackground.h>
#include <Model/MiscObjects/KeyFrameAnimation.h>
#include <Model/Groups/Group.h>
+#include <Model/Groups/KDTree.h>
#include <Model/Groups/ObjGroup.h>
#include <Model/Groups/DynBVH.h>
#include <Model/Lights/PointLight.h>
@@ -36,6 +37,7 @@
#ifdef USE_PRIVATE_CODE
cerr << " -CGT - use Coherent Grid Traversal acceleration
structure\n";
#endif
+ cerr << " -KDTree - use KDTree acceleration structure\n";
cerr << " -model - Required. The file to load (obj or ply file)\n";
cerr << " Can call this multiple times to load an
animation.\n";
cerr << " -animationLength - Number of seconds animation takes\n";
@@ -74,6 +76,9 @@
} else if (arg == "-DynBVH") {
delete as;
as = new DynBVH;
+ } else if (arg == "-KDTree") {
+ delete as;
+ as = new KDTree;
} else if (arg == "-CGT") {
#ifdef USE_PRIVATE_CODE
delete as;
- [Manta] r1838 - in trunk: Model/Groups scenes, thiago, 11/08/2007
Archive powered by MHonArc 2.6.16.
