-
sparker
thiago
Manta Interactive Ray Tracer Development Mailing List
Text archives Help
- From: "Thiago Ize" < >
- To:
- Subject: [Manta] r2289 - in trunk: Interface Model/Groups Model/Groups/BSP
- Date: Tue, 24 Jun 2008 17:27:13 -0600 (MDT)
Author: thiago
Date: Tue Jun 24 17:27:12 2008
New Revision: 2289
Added:
trunk/Model/Groups/BSP/
trunk/Model/Groups/BSP/BSH.cc
trunk/Model/Groups/BSP/BSH.h
trunk/Model/Groups/BSP/BSP.cc
trunk/Model/Groups/BSP/BSP.h
trunk/Model/Groups/BSP/Geometry.cc
trunk/Model/Groups/BSP/Geometry.h
trunk/Model/Groups/BSP/Polytope.cc
trunk/Model/Groups/BSP/Polytope.h
trunk/Model/Groups/BSP/aip.h
trunk/Model/Groups/BSP/chainHull.h
Modified:
trunk/Interface/AccelerationStructure.h
trunk/Model/Groups/CMakeLists.txt
trunk/Model/Groups/KDTree.cc
trunk/Model/Groups/KDTree.h
Log:
Model/Groups/BSP:
-Adding in BSP acceleration from my RT08 paper.
-There's a bunch of helper code and structures in the BSP
directory. If someone finds them useful for use in other manta
code, then feel free to clean it up and move it elsewhere.
Model/Groups/KDTree:
-Improved the file load/writing code to use binary format. Of course
this isn't as portable, but it results in accurate trees and much
faster load times.
Model/Groups/CMakeLists.txt:
-Added BSP
Interface/AccelerationStructure.h:
-Added interface for loading and saving acceleration structures to
file. This is really only useful for slow to build structures such
as the BSP and KDTree (currently the only two that implement this).
Modified: trunk/Interface/AccelerationStructure.h
==============================================================================
--- trunk/Interface/AccelerationStructure.h (original)
+++ trunk/Interface/AccelerationStructure.h Tue Jun 24 17:27:12 2008
@@ -37,6 +37,13 @@
virtual void addToUpdateGraph(ObjectUpdateGraph* graph,
ObjectUpdateGraphNode* parent) = 0;
+
+ virtual bool buildFromFile(const std::string &fileName) {
+ return false; //not implemented
+ }
+ virtual bool saveToFile(const std::string &fileName) {
+ return false; //not implemented
+ }
};
}
#endif //Manta_Interface_AccelerationStructure_h
Added: trunk/Model/Groups/BSP/BSH.cc
==============================================================================
--- (empty file)
+++ trunk/Model/Groups/BSP/BSH.cc Tue Jun 24 17:27:12 2008
@@ -0,0 +1,906 @@
+#include <Model/Groups/BSP/BSH.h>
+#include <Model/Primitives/MeshTriangle.h>
+#include <assert.h>
+using namespace Manta;
+
+//initializing static data members
+vector<Vector> BSH::triCentroid;
+vector<int> BSH::tmpTriID1, BSH::tmpTriID2;
+vector<BSH::ClippedTriangle> BSH::tmpCTri;
+
+void BSH::build(Mesh *mesh, const BBox &bounds)
+{
+ this->mesh = mesh;
+ nodes.resize(mesh->size()*2);
+
+ tmpTriID1.resize(mesh->size());
+ tmpTriID2.resize(mesh->size());
+
+ triCentroid.resize(mesh->size());
+ for (size_t i=0; i < mesh->size(); ++i) {
+ tmpTriID1[i] = i;
+ triCentroid[i] = mesh->getBBox(i).center();
+ }
+
+ int nextFree = 1;
+ nodes[0].parent = -1;
+ //spatialMedianBuild(0, nextFree, 0, mesh->size(), bounds);
+ localMedianBuild(0, nextFree, 0, mesh->size());
+
+ triID.resize(mesh->size());
+ for (size_t i=0; i < mesh->size(); ++i) {
+ const MeshTriangle *triangle = mesh->get(tmpTriID1[i]);
+ triID[i] = ClippedTriangle(tmpTriID1[i],
+ triangle->getVertex(0),
+ triangle->getVertex(1),
+ triangle->getVertex(2));
+ }
+
+ refit(0);
+
+ tmpTriID1.clear();
+ tmpTriID2.clear();
+ triCentroid.clear();
+}
+
+void BSH::buildSubset(Mesh *mesh, const vector<int> &subsetObjects,
+ const BBox &bounds)
+{
+ this->mesh = mesh;
+ nodes.resize(subsetObjects.size()*2);
+ tmpTriID1.resize(subsetObjects.size());
+ tmpTriID2.resize(subsetObjects.size());
+
+ //TODO: try and build bvh over subsetObjects.size() things, not the
+ //entire mesh.
+ triCentroid.resize(mesh->size());
+ for (size_t i=0; i < subsetObjects.size(); ++i) {
+ tmpTriID1[i] = subsetObjects[i];
+ triCentroid[subsetObjects[i]] = mesh->getBBox(subsetObjects[i]).center();
+ }
+
+ int nextFree = 1;
+ nodes[0].parent = -1;
+ //spatialMedianBuild(0, nextFree, 0, subsetObjects.size(), bounds);
+ localMedianBuild(0, nextFree, 0, subsetObjects.size());
+
+ triID.resize(subsetObjects.size());
+ for (size_t i=0; i < subsetObjects.size(); ++i) {
+ const MeshTriangle *triangle = mesh->get(tmpTriID1[i]);
+ triID[i] = ClippedTriangle(tmpTriID1[i],
+ triangle->getVertex(0),
+ triangle->getVertex(1),
+ triangle->getVertex(2));
+ }
+
+ refit(0);
+
+ tmpTriID1.clear();
+ tmpTriID2.clear();
+ triCentroid.clear();
+
+}
+
+void BSH::localMedianBuild(int nodeID, int &nextFree,
+ const size_t begin, const size_t end)
+{
+ BBox box;
+ for (size_t i=begin; i < end; ++i)
+ box.extendByPoint(triCentroid[tmpTriID1[i]]);
+
+ if (end - begin < 5 || box[0] == box[1]) {
+ nodes[nodeID].makeLeaf(begin, end);
+ return;
+ }
+
+ const int dim = box.longestAxis();
+ const Real split = 0.5f*(box[1][dim] + box[0][dim]);
+
+ for (size_t i=begin; i < end; ++i)
+ tmpTriID2[i] = tmpTriID1[i];
+ int *l = &tmpTriID1[begin];
+ int *r = &tmpTriID1[end - 1];
+
+ for (size_t i=begin;i<end;i++) {
+ int t = tmpTriID2[i];
+ if (triCentroid[t][dim] >= split)
+ *r-- = t;
+ else
+ *l++ = t;
+ }
+ int bestSplit = l - &tmpTriID1[begin];
+
+ if (bestSplit == 0) {
+ nodes[nodeID].makeLeaf(begin, end);
+ return;
+ }
+ if (begin + bestSplit == end) {
+ nodes[nodeID].makeLeaf(begin, end);
+ return;
+ }
+
+ nodes[nodeID].makeInternal(nextFree);
+ nodes[nextFree ].parent = nodeID;
+ nodes[nextFree+1].parent = nodeID;
+
+ nextFree+=2;
+
+ localMedianBuild(nodes[nodeID].children, nextFree,
+ begin, begin+bestSplit);
+ localMedianBuild(nodes[nodeID].children+1, nextFree,
+ begin+bestSplit, end);
+}
+
+void BSH::spatialMedianBuild(int nodeID, int &nextFree,
+ const size_t begin, const size_t end,
+ const BBox &bounds)
+{
+ //TODO: try a build that might be slower but produces something more
+ //suited to our needs. local median might be better quality, but
+ //slower to build.
+
+ BBox box;
+ for (size_t i=begin; i < end; ++i)
+ box.extendByPoint(triCentroid[tmpTriID1[i]]);
+
+ if (end - begin < 5 || box[0] == box[1]) {
+ nodes[nodeID].makeLeaf(begin, end);
+ return;
+ }
+
+ const int dim = bounds.longestAxis();
+ const Real split = 0.5f*(bounds[1][dim] + bounds[0][dim]);
+
+ for (size_t i=begin; i < end; ++i)
+ tmpTriID2[i] = tmpTriID1[i];
+ int *l = &tmpTriID1[begin];
+ int *r = &tmpTriID1[end - 1];
+
+ for (size_t i=begin;i<end;i++) {
+ int t = tmpTriID2[i];
+ if (triCentroid[t][dim] >= split)
+ *r-- = t;
+ else
+ *l++ = t;
+ }
+
+ BBox lBounds = bounds;
+ BBox rBounds = bounds;
+ lBounds[1][dim] = rBounds[0][dim] = split;
+
+ int bestSplit = l - &tmpTriID1[begin];
+
+ if (bestSplit == 0) {
+ spatialMedianBuild(nodeID, nextFree, begin, end, rBounds);
+ return;
+ }
+ if (begin + bestSplit == end) {
+ spatialMedianBuild(nodeID, nextFree, begin, end, lBounds);
+ return;
+ }
+
+ nodes[nodeID].makeInternal(nextFree);
+ nodes[nextFree ].parent = nodeID;
+ nodes[nextFree+1].parent = nodeID;
+
+ nextFree+=2;
+
+ spatialMedianBuild(nodes[nodeID].children, nextFree,
+ begin, begin+bestSplit, lBounds);
+ spatialMedianBuild(nodes[nodeID].children+1, nextFree,
+ begin+bestSplit, end, rBounds);
+}
+
+void BSH::refit(size_t nodeID)
+{
+ BSHNode &node = nodes[nodeID];
+ if (node.isLeaf()) {
+ for (int i=0; i < node.size; ++i) {
+ const int tri = triID[node.objectsIndex+i].originalTriID;
+ node.bounds.extendByBox(mesh->getBBox(tri));
+ }
+ }
+ else {
+ refit(node.children);
+ refit(node.children+1);
+
+ node.size = nodes[node.children].size + nodes[node.children+1].size;
+
+ if (nodes[node.children].size > 0 && nodes[node.children+1].size == 0) {
+ int parent = node.parent;
+ node = nodes[node.children];
+ node.parent = parent;
+ }
+ else if (nodes[node.children].size == 0 && nodes[node.children+1].size >
0) {
+ int parent = node.parent;
+ node = nodes[node.children+1];
+ node.parent = parent;
+ }
+
+ node.bounds.extendByBox(nodes[node.children].bounds);
+ node.bounds.extendByBox(nodes[node.children+1].bounds);
+ }
+
+ //compute a loose bound.
+ //TODO: compute a tighter bound
+ node.center = node.bounds.center();
+ node.radius = 0.5f * node.bounds.diagonal().length();
+}
+
+
+void BSH::splitTriangle(const ClippedTriangle &tri,
+ const BuildSplitPlane &plane,
+ ClippedTriangle negSide[2],
+ ClippedTriangle posSide[2])
+{
+ //This is the Reentrant Polygon Clipping algorithm by Sutherland
+ //and Hodgman.
+
+ int negCount=0;
+ int posCount=0;
+
+ Point s, f;
+ for (size_t i=0; i < 3; ++i) {
+ const Point &p = tri.tri[i];
+
+ if (p != tri.tri[i]) {
+ cout << p<<endl<<tri.tri[i]<<endl;
+ exit(1);
+ }
+
+ double p_d = signedDistance(plane, p);
+ if (i==0) { //first point?
+ s = f = p;
+ }
+ else {
+ double s_d = signedDistance(plane, s);
+
+ if (s_d*p_d < 0) { //does line sp cross dividing plane?
+ //compute intersection I, of sp and plane.
+ const double length1 = fabs(p_d);
+ const double length2 = fabs(s_d);
+
+ const double a = length1 / (length1 + length2);
+
+ //We do the a>.5 test for numerical stability.
+ Point I;
+ if (a > .5)
+ I = p + a*(s-p);
+ else
+ I = s + (1-a)*(p-s);
+
+ if(plane.normal[0]==1)
+ I[0]=-plane.d;
+ else if(plane.normal[1]==1)
+ I[1]=-plane.d;
+ else if(plane.normal[2]==1)
+ I[2]=-plane.d;
+
+ assert(negCount < 3);
+ assert(posCount < 3);
+ negSide[0].tri[negCount++] = I;
+ posSide[0].tri[posCount++] = I;
+ }
+ }
+
+ //How does p relate to plane?
+ if (p_d == 0) { //on the plane
+ assert(negCount < 3);
+ assert(posCount < 3);
+ negSide[0].tri[negCount++] = p;
+ posSide[0].tri[posCount++] = p;
+ }
+ else if (p_d < 0)
+ if (negCount == 3) {
+ negSide[1].tri[0] = negSide[0].tri[0];
+ negSide[1].tri[1] = negSide[0].tri[2];
+ negSide[1].tri[2] = p;
+ ++negCount;
+ }
+ else
+ negSide[0].tri[negCount++] = p;
+ else if (p_d > 0)
+ if (posCount == 3) {
+ posSide[1].tri[0] = posSide[0].tri[0];
+ posSide[1].tri[1] = posSide[0].tri[2];
+ posSide[1].tri[2] = p;
+ ++posCount;
+ }
+ else
+ posSide[0].tri[posCount++] = p;
+ s=p;
+ }
+
+ //now we need to close the two polygons.
+ //TODO: once we know it works, do this closing as part of the for
+ //loop.
+ const Point &p = f;
+ double p_d = signedDistance(plane, p);
+ double s_d = signedDistance(plane, s);
+
+ if (s_d*p_d < 0) { //does line sp cross dividing plane?
+ //compute intersection I, of sp and plane.
+ const double length1 = fabs(p_d);
+ const double length2 = fabs(s_d);
+
+ const double a = length1 / (length1 + length2);
+ Point I;
+ if (a > .5)
+ I = p + a*(s-p);
+ else
+ I = s + (1-a)*(p-s);
+
+ if(plane.normal[0]==1)
+ I[0]=-plane.d;
+ else if(plane.normal[1]==1)
+ I[1]=-plane.d;
+ else if(plane.normal[2]==1)
+ I[2]=-plane.d;
+
+ if (negCount == 3) {
+ negSide[1].tri[0] = negSide[0].tri[0];
+ negSide[1].tri[1] = negSide[0].tri[2];
+ negSide[1].tri[2] = I;
+ ++negCount;
+ }
+ else
+ negSide[0].tri[negCount++] = I;
+
+ if (posCount == 3) {
+ posSide[1].tri[0] = posSide[0].tri[0];
+ posSide[1].tri[1] = posSide[0].tri[2];
+ posSide[1].tri[2] = I;
+ ++posCount;
+ }
+ else
+ posSide[0].tri[posCount++] = I;
+ }
+
+ assert(negCount > 0 && negCount < 5);
+ assert(posCount > 0 && posCount < 5);
+
+ if (negCount < 4)
+ negSide[1].originalTriID = -1;
+ if (posCount < 4)
+ posSide[1].originalTriID = -1;
+}
+
+
+void BSH::split(const BuildSplitPlane &plane, BSH &positiveBSH,
+ const Location coplanar_primitiveSide)
+{
+ size_t oldSize = size();
+ tmpCTri.clear();
+ positiveBSH.triID.clear();
+
+ split(plane, positiveBSH, coplanar_primitiveSide, 0);
+
+ triID = tmpCTri;
+
+ if (oldSize > 1) {
+ size_t i=0;
+ if (size()-numSplits() > 0 && size()-numSplits() < oldSize)
+ compact(i);
+ i=0;
+ if (positiveBSH.size()-positiveBSH.numSplits() > 0 &&
+ positiveBSH.size()-positiveBSH.numSplits() < oldSize)
+ positiveBSH.compact(i);
+ }
+}
+
+void BSH::nodeTriIDCopy(size_t nodeID, const vector<ClippedTriangle> &triID,
+ vector<ClippedTriangle> &tmpCTri)
+{
+ BSHNode &node = nodes[nodeID];
+ if (node.isLeaf()) {
+ tmpCTri.insert(tmpCTri.end(), triID.begin()+node.objectsIndex,
+ triID.begin()+node.objectsIndex+node.size);
+ node.objectsIndex = tmpCTri.size()-node.size;
+ }
+ else {
+ nodeTriIDCopy(node.children, triID, tmpCTri);
+ nodeTriIDCopy(node.children+1, triID, tmpCTri);
+ }
+}
+
+void BSH::split(const BuildSplitPlane &plane, BSH &positiveBSH,
+ const Location coplanar_primitiveSide, size_t nodeID)
+{
+ //assume positiveBSH is a copy of *this.
+ //negative is us.
+
+ BSHNode &posNode = positiveBSH.nodes[nodeID];
+ BSHNode &negNode = nodes[nodeID];
+
+ Location position = posNode.whichSide(plane);
+ if (position == negativeSide) {
+ posNode.makeLeaf(0,0);
+ posNode.bounds.reset();
+ nodeTriIDCopy(nodeID, triID, tmpCTri);
+ return;
+ }
+ else if (position == positiveSide) {
+ negNode.makeLeaf(0,0);
+ negNode.bounds.reset();
+ positiveBSH.nodeTriIDCopy(nodeID, triID, positiveBSH.triID);
+ return;
+ }
+ else if (posNode.isLeaf()) {
+
+ negNode.bounds.reset();
+ posNode.bounds.reset();
+
+ unsigned int size = posNode.size;
+ posNode.size=0;
+ negNode.size=0;
+
+ for (size_t i=0; i < size; ++i) {
+ ClippedTriangle &tri = triID[posNode.objectsIndex+i];
+ const Location pos = getLocation(tri, plane);
+ switch(pos) {
+ case negativeSide:
+ tmpCTri.push_back(tri);
+ negNode.size++;
+ break;
+ case positiveSide:
+ positiveBSH.triID.push_back(tri);
+ posNode.size++;
+ break;
+ case bothSides:
+ //need to clip the triangle
+ tmpCTri.push_back(tri);
+ negNode.size++;
+
+ positiveBSH.triID.push_back(tri);
+ posNode.size++;
+
+#if 1 //do triangle clipping
+ tmpCTri.push_back(tri);
+ positiveBSH.triID.push_back(tri);
+
+ splitTriangle(tri, plane, &tmpCTri[tmpCTri.size()-2],
+ &positiveBSH.triID[positiveBSH.triID.size()-2]);
+
+ if (tmpCTri.back().originalTriID == -1) {
+ tmpCTri.pop_back();
+ }
+ else {
+ negNode.size++;
+ }
+ if (positiveBSH.triID.back().originalTriID == -1) {
+ positiveBSH.triID.pop_back();
+ }
+ else {
+ posNode.size++;
+ }
+#endif
+ break;
+ case eitherSide:
+ if (!tri.fixedPlane and true){ //set to false to not use fixedPlane
+ tri.fixedPlane = true;
+ tri.plane = plane;
+ }
+ if (coplanar_primitiveSide==negativeSide) {
+ tmpCTri.push_back(tri);
+ negNode.size++;
+ }
+ else {
+ positiveBSH.triID.push_back(tri);
+ posNode.size++;
+ }
+ break;
+ }
+ }
+
+ negNode.objectsIndex = tmpCTri.size() - negNode.size;
+ posNode.objectsIndex = positiveBSH.triID.size() - posNode.size;
+
+ //need to find out how many split triangles there are now. We
+ //know that split triangles will always be next to each other.
+ negNode.numSplits=0;
+ for (int i=0; i < negNode.size; ++i) {
+ if (i!=0 &&
+ tmpCTri[negNode.objectsIndex+i-1].originalTriID ==
+ tmpCTri[negNode.objectsIndex+i].originalTriID)
+ negNode.numSplits++;
+
+ negNode.bounds.extendByPoint(tmpCTri[negNode.objectsIndex+i].tri[0]);
+ negNode.bounds.extendByPoint(tmpCTri[negNode.objectsIndex+i].tri[1]);
+ negNode.bounds.extendByPoint(tmpCTri[negNode.objectsIndex+i].tri[2]);
+ }
+ posNode.numSplits=0;
+ for (int i=0; i < posNode.size; ++i) {
+ if (i!=0 &&
+ positiveBSH.triID[posNode.objectsIndex+i-1].originalTriID ==
+ positiveBSH.triID[posNode.objectsIndex+i].originalTriID)
+ posNode.numSplits++;
+
+
posNode.bounds.extendByPoint(positiveBSH.triID[posNode.objectsIndex+i].tri[0]);
+
posNode.bounds.extendByPoint(positiveBSH.triID[posNode.objectsIndex+i].tri[1]);
+
posNode.bounds.extendByPoint(positiveBSH.triID[posNode.objectsIndex+i].tri[2]);
+ }
+ }
+ else {
+ split(plane, positiveBSH, coplanar_primitiveSide, negNode.children);
+ split(plane, positiveBSH, coplanar_primitiveSide, negNode.children+1);
+
+ posNode.size = positiveBSH.nodes[posNode.children].size +
+ positiveBSH.nodes[posNode.children+1].size;
+ negNode.size = nodes[negNode.children].size +
+ nodes[negNode.children+1].size;
+ posNode.numSplits = positiveBSH.nodes[posNode.children].numSplits +
+ positiveBSH.nodes[posNode.children+1].numSplits;
+ negNode.numSplits = nodes[negNode.children].numSplits +
+ nodes[negNode.children+1].numSplits;
+
+ //If only left child has items, make left child become child of
+ //its grand parent (i.e. remove current node).
+ if (nodes[negNode.children].size > 0 && nodes[negNode.children+1].size
== 0) {
+ int parent = negNode.parent;
+ negNode = nodes[negNode.children];
+ negNode.parent = parent;
+ }
+ else if (nodes[negNode.children].size == 0 &&
nodes[negNode.children+1].size > 0) {
+ int parent = negNode.parent;
+ negNode = nodes[negNode.children+1];
+ negNode.parent = parent;
+ }
+ else if (negNode.size == 0)
+ negNode.makeLeaf(0,0);
+
+ if (positiveBSH.nodes[posNode.children].size > 0 &&
positiveBSH.nodes[posNode.children+1].size == 0) {
+ int parent = posNode.parent;
+ posNode = positiveBSH.nodes[posNode.children];
+ posNode.parent = parent;
+ }
+ else if (positiveBSH.nodes[posNode.children].size == 0 &&
positiveBSH.nodes[posNode.children+1].size > 0) {
+ int parent = posNode.parent;
+ posNode = positiveBSH.nodes[posNode.children+1];
+ posNode.parent = parent;
+ }
+ else if (posNode.size == 0)
+ posNode.makeLeaf(0,0);
+
+ if (!posNode.isLeaf()) {
+ posNode.bounds.extendByBox(positiveBSH.nodes[posNode.children].bounds);
+
posNode.bounds.extendByBox(positiveBSH.nodes[posNode.children+1].bounds);
+ //TODO: compute a tighter bound
+ //compute a loose bound.
+ posNode.center = posNode.bounds.center();
+ posNode.radius = 0.5f * posNode.bounds.diagonal().length();
+ }
+
+ if (!negNode.isLeaf()) {
+ negNode.bounds.extendByBox(nodes[negNode.children+0].bounds);
+ negNode.bounds.extendByBox(nodes[negNode.children+1].bounds);
+
+ negNode.center = negNode.bounds.center();
+ negNode.radius = 0.5f * negNode.bounds.diagonal().length();
+ }
+ }
+}
+
+void BSH::compact(size_t &highestSeen, size_t node)
+{
+
+ BSHNode &currNode = nodes[node];
+
+ if (highestSeen < node) { //this is only true for left children
+ highestSeen = node+1;//+1 for right child.
+ }
+
+ if (currNode.isLeaf())
+ return;
+
+ if (currNode.children > highestSeen+2) {
+ //looks like we have some nodes that went unused.
+ nodes[++highestSeen] = nodes[currNode.children];
+ nodes[highestSeen].parent = node;
+ nodes[++highestSeen] = nodes[currNode.children+1];
+ nodes[highestSeen].parent = node;
+ currNode.children = highestSeen-1;
+ }
+ compact(highestSeen, currNode.children);
+ compact(highestSeen, currNode.children+1);
+
+ if (node == 0) {
+ nodes.resize(highestSeen+1);
+ }
+}
+
+BSH::Iterator BSH::begin() const
+{
+ Iterator iter(*this, 0);
+ while (!nodes[iter.currNode].isLeaf()) {
+ iter.currNode = nodes[iter.currNode].children;
+ }
+ return iter.next();
+}
+
+
+BSH::Iterator &BSH::Iterator::next()
+{
+ if (bsh.nodes[currNode].isLeaf()) {
+ if (++i < bsh.nodes[currNode].size)
+ return *this;
+ else if (currNode == 0) {
+ i=-1;
+ return *this;
+ }
+ }
+ bool retreating = true;
+ unsigned int prevNode = currNode;
+ currNode = bsh.nodes[currNode].parent;
+
+ while(true) {
+ if (!bsh.nodes[currNode].isLeaf()) {
+ if (retreating) {
+ if (bsh.nodes[currNode].children == prevNode) {
+ retreating = false;
+ currNode = bsh.nodes[currNode].children+1;
+ }
+ else if (currNode == 0) {
+ i=-1;
+ return *this;
+ }
+ else {
+ prevNode = currNode;
+ currNode = bsh.nodes[currNode].parent;
+ }
+ }
+ else
+ currNode = bsh.nodes[currNode].children; //go left.
+ }
+ else { //we're at a leaf
+ if (bsh.nodes[currNode].size > 0) {
+ i = 0;
+ return *this;
+ }
+ else {
+ if (currNode == 0) {
+ i=-1;
+ return *this;
+ }
+ retreating = true;
+ prevNode = currNode;
+ currNode = bsh.nodes[currNode].parent;
+ }
+ }
+ }
+}
+
+Location BSH::getLocation(const ClippedTriangle &ctri,
+ const BuildSplitPlane &plane) const
+{
+ if (ctri.plane.normal == plane.normal && ctri.plane.d == plane.d) {
+ return eitherSide;
+ }
+
+ double p1_dist = signedDistance(plane, ctri.tri[0]);
+ double p2_dist = signedDistance(plane, ctri.tri[1]);
+ double p3_dist = signedDistance(plane, ctri.tri[2]);
+
+ const double MIN_DIST = 1e-6;
+ if ( fabs(p1_dist) < MIN_DIST )
+ p1_dist = 0;
+ if ( fabs(p2_dist) < MIN_DIST )
+ p2_dist = 0;
+ if ( fabs(p3_dist) < MIN_DIST )
+ p3_dist = 0;
+
+ const double points1_2 = p1_dist * p2_dist;
+ const double points1_3 = p1_dist * p3_dist;
+ const double points2_3 = p2_dist * p3_dist;
+
+ //Check to see if the triangle intersects the plane. It's ok if
+ //the triangle lies on the plane, hence we allow distance to be 0,
+ //but need to check all combinations of points since a 0 would
+ //mask out the sign of the other point.
+ if (points1_2 >= 0 && points1_3 >= 0 && points2_3 >= 0) {
+ if (p1_dist > 0 || p2_dist > 0 || p3_dist > 0) {
+ //positive distance places tri on right child.
+ return positiveSide;
+ }
+ else if (p1_dist == 0 && p2_dist == 0 && p3_dist == 0) {
+ //all zero distance places tri on split, so it can go to either side//
on left child.
+ return eitherSide;
+ }
+ else {
+ //we have negative (and maybe zero) distances.
+ return negativeSide;
+ }
+ }
+ else {
+ //need to place triangle in both children.
+ return bothSides;
+ }
+}
+
+
+void BSH::pointsNearPlaneDist(vector<Point> &points,
+ const BuildSplitPlane &plane,
+ const int nodeID) const
+{
+ const BSHNode &node = nodes[nodeID];
+
+ const double MAX_DIST = 2.5e-6;
+
+ bool nearPlane = node.nearPlaneDist(plane, MAX_DIST);
+ if (!nearPlane)
+ return;
+
+ if (node.isLeaf()) {
+ for (int i=0; i < node.size; ++i) {
+
+ const Point &p0 = triID[node.objectsIndex+i].tri[0];
+ const Point &p1 = triID[node.objectsIndex+i].tri[1];
+ const Point &p2 = triID[node.objectsIndex+i].tri[2];
+
+ double dist[3] = {signedDistance(plane, p0),
+ signedDistance(plane, p1),
+ signedDistance(plane, p2)};
+
+ int pointOnPlaneID=-1;
+ for (int k=0; k < 3; ++k)
+ if (fabs(dist[k]) <= MAX_DIST) {
+ dist[k]=0;
+ pointOnPlaneID=k;
+ }
+
+ for (int k=0; k<3; ++k) {
+ if (pointOnPlaneID>=0) {
+ if (dist[k] == 0)
+ points.push_back(triID[node.objectsIndex+i].tri[k]);
+ else {
+ //if one of the points was close enough to be considered
+ //on the plane, but this point isn't, it's possible it
+ //should be on the plane, but was just slightly too far
+ //away. So instead we check to see if the edge formed by
+ //this point and the point known to be on the plane is
+ //orthogonal to the plane normal, and so must be on the
+ //plane.
+ const Point &p = triID[node.objectsIndex+i].tri[pointOnPlaneID];
+ Point edge = (p - triID[node.objectsIndex+i].tri[k]).normal();
+ double costheta = 1 - fabs(Dot(edge, plane.normal));
+ if (costheta > 0.999999) {
+#ifdef BUILD_DEBUG
+ cout <<"angle: " <<costheta<<endl;
+ cout <<"got one! it had a distance of: "
+ <<signedDistance(plane, triID[node.objectsIndex+i].tri[k])
+ <<" " <<dist[k]<<endl
+ <<acos(costheta)<<" with d of " <<plane.d<< endl
+ <<Dot(plane.normal,
triID[node.objectsIndex+i].tri[k])<<endl;
+#endif
+ points.push_back(triID[node.objectsIndex+i].tri[k]);
+ // exit(1);
+ }
+ }
+ }
+ }
+ }
+ }
+ else {
+ pointsNearPlaneDist(points, plane, node.children);
+ pointsNearPlaneDist(points, plane, node.children+1);
+ }
+}
+
+void BSH::pointNearPlaneAngle(const vector<Point> &points,
+ const BuildSplitPlane &plane,
+ Point &closePoint, double &minCostheta,
+ const int nodeID) const
+{
+ assert(points.size() == 2);
+
+ const BSHNode &node = nodes[nodeID];
+
+ bool nearPlane = node.nearPlaneAngle(plane, points[0], minCostheta);
+ if (!nearPlane)
+ return;
+
+ if (node.isLeaf()) {
+ const Point edge = points[1]-points[0];
+ for (int i=0; i < node.size; ++i) {
+ for (int k=0; k<3; ++k) {
+ const Point &p = triID[node.objectsIndex+i].tri[k];
+ if (p == points[0] || (points.size()>1 && p == points[1]))
+ continue;
+ //TODO: optimize this by passing in some of these vectors to the
function.
+ const Point normal2 = Cross(edge, p-points[0]);
+ const float normal2_length2 = static_cast<float>(normal2.length2());
+ if (normal2_length2 < 1e-12f)
+ continue;
+
+ const float inv_length = inverseSqrt(normal2_length2);
+ const double costheta = 1-fabs(Dot(normal2,
plane.normal)*inv_length);
+
+ if (costheta==minCostheta) {
+ if (closePoint == p)
+ continue;
+ if (normal2_length2 <= Cross(edge, closePoint-points[0]).length2())
+ continue;
+// cout <<"still equal! " << costheta<<endl;
+ }
+ if (costheta <= minCostheta) {
+ closePoint = p;
+ minCostheta = costheta;
+ }
+ }
+ }
+ }
+ else {
+ pointNearPlaneAngle(points, plane, closePoint,
+ minCostheta, node.children);
+ pointNearPlaneAngle(points, plane, closePoint,
+ minCostheta, node.children+1);
+ }
+}
+
+void BSH::countLocations(const BuildSplitPlane &plane,
+ int &nLeft, int &nRight, int &nEither,
+ const int nodeID) const
+{
+ const BSHNode &node = nodes[nodeID];
+
+ if (node.isLeaf()) {
+ if (node.numSplits > 0) {
+ Location location = getLocation(triID[node.objectsIndex],
+ plane);
+ switch (location) {
+ case eitherSide: ++nEither;
+ break;
+ case negativeSide: ++nLeft;
+ break;
+ case bothSides: ++nLeft;
+ case positiveSide: ++nRight;
+ break;
+ }
+
+ for (int i=1; i < node.size; ++i) {
+ if (triID[node.objectsIndex+i-1].originalTriID !=
+ triID[node.objectsIndex+i].originalTriID) {
+ Location location = getLocation(triID[node.objectsIndex+i],
+ plane);
+ switch (location) {
+ case eitherSide: ++nEither;
+ break;
+ case negativeSide: ++nLeft;
+ break;
+ case bothSides: ++nLeft;
+ case positiveSide: ++nRight;
+ break;
+ }
+ }
+ }
+ }
+ else {
+ for (int i=0; i < node.size; ++i) {
+ Location location = getLocation(triID[node.objectsIndex+i],
+ plane);
+ switch (location) {
+ case eitherSide: ++nEither;
+ break;
+ case negativeSide: ++nLeft;
+ break;
+ case bothSides: ++nLeft;
+ case positiveSide: ++nRight;
+ break;
+ }
+ }
+ }
+ }
+ else {
+ Location location = node.whichSide(plane);
+
+ switch (location) {
+ case negativeSide: nLeft+=node.size-node.numSplits;
+ break;
+ case positiveSide: nRight+=node.size-node.numSplits;
+ break;
+ default:
+ countLocations(plane, nLeft, nRight, nEither, node.children);
+ countLocations(plane, nLeft, nRight, nEither, node.children+1);
+ break;
+ }
+ }
+}
Added: trunk/Model/Groups/BSP/BSH.h
==============================================================================
--- (empty file)
+++ trunk/Model/Groups/BSP/BSH.h Tue Jun 24 17:27:12 2008
@@ -0,0 +1,200 @@
+#ifndef _MANTA_BSH_H
+#define _MANTA_BSH_H
+
+#include <Model/Groups/BSP/Geometry.h>
+#include <Core/Geometry/BBox.h>
+#include <Core/Geometry/Vector.h>
+#include <Model/Groups/Mesh.h>
+
+namespace Manta
+{
+ using namespace BSPs;
+ enum Location{negativeSide, positiveSide, bothSides, eitherSide};
+
+ class BSH
+ {
+ public:
+ struct BSHNode
+ {
+ BBox bounds;
+ Point center;
+ double radius;
+
+ union {
+ unsigned int children; //INTERNAL NODE:
+ //leftChild=children. rightChild=children+1.
+ unsigned int objectsIndex; //LEAF NODE:
+ };
+
+ int parent;
+
+ int size; //number of primitives contained by this node and its
+ //children
+ int numSplits; //number of triangles that were split by this node tree.
+
+ bool _isLeaf;
+
+ inline void makeLeaf(int begin, int end) {
+ objectsIndex = begin;
+ size = end - begin;
+ numSplits=0;
+ _isLeaf = true;
+ }
+
+ inline void makeInternal(int first_child) {
+ children = first_child;
+ _isLeaf = false;
+ numSplits = 0;
+ }
+
+ inline bool isLeaf() const {
+ return _isLeaf;
+ }
+
+ inline Location whichSide(const BuildSplitPlane& plane) const {
+ const double dist = signedDistance(plane, center);
+ if (dist < -radius)
+ return negativeSide;
+ else if (dist > radius)
+ return positiveSide;
+ else {
+ return bothSides;
+ }
+ }
+ inline bool nearPlaneDist(const BuildSplitPlane& plane,
+ double maxError) const {
+ const double dist = signedDistance(plane, center);
+ return (fabs(dist) <= radius+maxError);
+ }
+ inline bool nearPlaneAngle(const BuildSplitPlane& plane,
+ const Point &p,
+ double minAngle) const {
+ const double dist = signedDistance(plane, center);
+ if (fabs(dist) < radius)
+ return true;
+ const double r = dist < 0 ? -radius : radius;
+
+ const Point center_p = center-p;
+
+ const Point v1 = (center_p - plane.normal*dist);
+ const Point v2 = (center_p - plane.normal*r);
+
+ const double costheta = 1-fabs(Dot(v1, v2)) *
+ (inverseSqrt(static_cast<float>(v1.length2())) *
+ inverseSqrt(static_cast<float>(v2.length2())));
+ return costheta < minAngle;
+ }
+ };
+
+ vector<BSHNode> nodes;
+
+ struct ClippedTriangle {
+ int originalTriID;
+ Point tri[3]; //vertices of the clipped triangle
+ BuildSplitPlane plane;
+ bool fixedPlane;
+
+ ClippedTriangle() { } //Do not use. This is only for std::vector
+ ClippedTriangle(int originalTriID, const Point &v0,
+ const Point &v1, const Point &v2) :
+ originalTriID(originalTriID), fixedPlane(false) {
+ tri[0] = v0;
+ tri[1] = v1;
+ tri[2] = v2;
+
+ plane.normal = getTriangleNormal(tri[0], tri[1], tri[2]).normal();
+ plane.d = Dot(plane.normal, -(tri[0]+tri[1]+tri[2])/3.0);
+ }
+ };
+ vector<ClippedTriangle> triID;
+
+ Mesh* mesh;
+
+ BSH() { }
+ BSH(const BSH &bsh) {
+ nodes = bsh.nodes;
+ mesh = bsh.mesh;
+ }
+
+ //this becomes negative side of plane and positiveBSH is the
+ //positive side.
+ void split(const BuildSplitPlane &plane, BSH &positiveBSH,
+ const Location coplanar_primitiveSide);
+
+ void build(Mesh *mesh, const BBox &bounds);
+ void buildSubset(Mesh *mesh, const vector<int> &subsetObjects,
+ const BBox &bounds);
+
+ inline size_t size() const { return nodes[0].size; }
+ inline size_t numSplits() const { return nodes[0].numSplits; }
+
+ void pointsNearPlaneDist(vector<Point> &points,
+ const BuildSplitPlane &plane,
+ const int nodeID=0) const;
+ void pointNearPlaneAngle(const vector<Point> &points,
+ const BuildSplitPlane &plane,
+ Point &closePoint, double &minAngle,
+ const int nodeID=0) const;
+
+ void countLocations(const BuildSplitPlane &plane,
+ int &nLeft, int &nRight, int &nEither,
+ const int nodeID=0) const;
+
+ class Iterator {
+ public:
+ const BSH& bsh;
+ int currNode;
+ int i;
+
+ Iterator(const BSH& bsh, int currNode) :
+ bsh(bsh), currNode(currNode), i(-1) { }
+
+ const ClippedTriangle &get() const {
+ return bsh.triID[bsh.nodes[currNode].objectsIndex+i];
+ }
+ Iterator &next();
+ bool operator!=(const Iterator &iter) const {
+ return (currNode != iter.currNode) || (i != iter.i);
+ }
+ };
+ Iterator begin() const;
+ inline static const Iterator& end() {
+ const static BSH foo;
+ const static Iterator endIterator(foo, 0);
+ return endIterator;
+ }
+
+ private:
+ void localMedianBuild(int nodeID, int &nextFree,
+ const size_t begin, const size_t end);
+
+ void spatialMedianBuild(int nodeID, int &nextFree,
+ const size_t begin, const size_t end,
+ const BBox &bounds);
+ void refit(size_t nodeID);
+
+ void splitTriangle(const ClippedTriangle &tri,
+ const BuildSplitPlane &plane,
+ ClippedTriangle negSide[2],
+ ClippedTriangle posSide[2]);
+
+ void split(const BuildSplitPlane &plane, BSH &positiveBSH,
+ const Location coplanar_primitiveSide,
+ size_t nodeID);
+
+ void nodeTriIDCopy(size_t nodeID, const vector<ClippedTriangle> &triID,
+ vector<ClippedTriangle> &tmpCTri);
+
+ void compact(size_t &highestSeen, size_t node=0);
+
+ Location getLocation(const ClippedTriangle &ctri,
+ const BuildSplitPlane &plane) const;
+
+ //These are used only during the build (maybe I should just pass
+ //these as arguments?).
+ static vector<int> tmpTriID1, tmpTriID2;
+ static vector<ClippedTriangle> tmpCTri;
+ static vector<Vector> triCentroid;
+ };
+};
+#endif
Added: trunk/Model/Groups/BSP/BSP.cc
==============================================================================
--- (empty file)
+++ trunk/Model/Groups/BSP/BSP.cc Tue Jun 24 17:27:12 2008
@@ -0,0 +1,1801 @@
+#include <Model/Groups/BSP/BSP.h>
+#include <Model/Groups/KDTree.h> //for the mailbox
+#include <Model/Primitives/MeshTriangle.h>
+#include <Core/Geometry/Vector.h>
+#include <Core/Geometry/AffineTransformT.h>
+#include <Core/Math/MinMax.h>
+#include <Core/Util/Preprocessor.h>
+#include <Interface/Context.h>
+#include <Core/Thread/Time.h>
+#include <Interface/MantaInterface.h>
+
+#include <fstream>
+#include <iostream>
+#include <limits>
+#include <set>
+#include <iomanip>
+
+// if defined, we'll add _realtive_ epsilons to plane distance (i.e.,
compare with t_p +/- 0.0001*t_p
+#define SSE_REL_EPSILONS
+
+using namespace Manta;
+
+void BSP::preprocess(const PreprocessContext &context)
+{
+ if (mesh) {
+ static_cast<Group*>(mesh)->computeBounds(context, bounds);
+ mesh->preprocess(context);
+ }
+
+ context.manta_interface->registerSerialPreRenderCallback(
+ Callback::create(this, &BSP::update));
+}
+
+void BSP::setGroup(Group* new_group)
+{
+ mesh = dynamic_cast<Mesh*>(new_group);
+}
+
+void BSP::rebuild(int proc, int numProcs)
+{
+ if (proc != 0)
+ return;
+
+ double startTime = Time::currentSeconds();
+
+ bounds.reset();
+
+ //Something's wrong, let's bail. Since we reset the bounds right
+ //above, this BSP should never get intersected.
+ if (!mesh)
+ return;
+
+ PreprocessContext context;
+ static_cast<Group*>(mesh)->computeBounds(context, bounds);
+
+ objects.clear();
+ nodes.clear();
+ nodes.reserve(mesh->size()*4);
+
+ Point *vertexData = new Point[mesh->size()*3];
+ for (size_t i=0; i < mesh->size(); ++i) {
+ MeshTriangle *tri = mesh->get(i);
+ vertexData[i*3+0] = tri->getVertex(0);
+ vertexData[i*3+1] = tri->getVertex(1);
+ vertexData[i*3+2] = tri->getVertex(2);
+ vertices.push_back(&vertexData[i*3]);
+ }
+
+ int nextFreeNode = 1;
+
+ //During traversal we do a very fast bounding box test to quickly
+ //decide whether we need to actually intersect the BSP. So here we
+ //place "fake" initial splitting planes to simulate the planes that
+ //make up the bounding box test.
+ vector<BuildSplitPlane> splitPlanes;
+
+ cout <<"bounds are: "<<bounds[0] <<" , " <<bounds[1]<<endl;
+ cout <<"bounds SA: "<<bounds.computeArea()<<endl;
+
+ polytope.initialize(bounds);
+
+ float lArea = polytope.getSurfaceArea();
+ cout <<"root area is: "<<lArea<<endl;
+
+ BSH bsh;
+ bsh.build(mesh, bounds);
+ cout <<"BSH built\n";
+ build(0, nextFreeNode, splitPlanes, bsh, 1);
+
+ delete[] vertexData;
+ vertices.clear();
+
+ float buildTime = Time::currentSeconds() - startTime;
+ printf("BSP tree built in %f seconds\n", buildTime);
+ printStats();
+}
+
+void BSP::rebuildSubset(const vector<int> &subsetObjects, const BBox &bounds)
+{
+ isSubset = true;
+
+ nodes.clear();
+ nodes.reserve(subsetObjects.size()*4);
+
+ objects = subsetObjects;
+ this->bounds = bounds;
+
+ Point *vertexData = new Point[mesh->size()*3];
+ for (size_t i=0; i < mesh->size(); ++i) {
+ MeshTriangle *tri = mesh->get(i);
+ vertexData[i*3+0] = tri->getVertex(0);
+ vertexData[i*3+1] = tri->getVertex(1);
+ vertexData[i*3+2] = tri->getVertex(2);
+ vertices.push_back(&vertexData[i*3]);
+ }
+ int nextFreeNode = 1;
+
+ polytope.initialize(bounds);
+
+ BSH bsh;
+ bsh.buildSubset(mesh,subsetObjects, bounds);
+
+ vector<BuildSplitPlane> splitPlanes;
+ build(0, nextFreeNode, splitPlanes, bsh, 1);
+
+ delete[] vertexData;
+ vertices.clear();
+}
+
+#ifdef MANTA_SSE
+// packet BSP traversal
+// Object *split_primitive = NULL;
+
+typedef sse_t sse_arr[RayPacket::MaxSize/4];
+
+#define false4() cmp4_lt(zero4(),zero4())
+#endif
+
+void BSP::intersect(const RenderContext& context, RayPacket& rays) const {
+#ifndef MANTA_SSE
+ intersectSingleRays(context, rays);
+#else
+ //Can only do packet traversal if we have constant origins. Furthermore,
+ //it's only worth doing if there's more than one ray in the packet.
+ if (!rays.getFlag(RayPacket::ConstantOrigin) || rays.end()-rays.begin()==1)
+ intersectSingleRays(context, rays);
+ else
+ intersectPacket(context, rays);
+#endif
+}
+
+void BSP::intersectPacket(const RenderContext& context, RayPacket& rays)
const
+// void BSP::intersect_same_origin(const RenderContext& context, RayPacket&
rays) const
+{
+#ifdef MANTA_SSE
+
+// rays.normalizeDirections();
+// rays.computeSigns();
+ rays.computeInverseDirections();
+
+ const int ray_begin = rays.begin();
+ const int ray_end = rays.end();
+ const int sse_begin = ray_begin / 4;
+ const int sse_end = (ray_end+3) / 4;
+
+#ifdef COLLECT_STATS
+ nTotalRays += rays.end() - rays.begin();
+#endif
+
+#undef MAILBOX
+#ifdef MAILBOX
+ Mailbox mailbox;
+ mailbox.clear();
+#endif
+ assert(rays.getFlag(RayPacket::ConstantOrigin));
+
+ sse_t floatStackMem[RayPacket::MaxSize*80*3];
+ float *freeMem = (float *)floatStackMem;
+
+ sse_t *t_n = (sse_t*)freeMem; freeMem += RayPacket::MaxSize;
+ sse_t *t_f = (sse_t*)freeMem; freeMem += RayPacket::MaxSize;
+ sse_arr t_p;
+
+ // set all sse-blocks to full ray segment
+ for (int i=sse_begin;i<sse_end;i++) {
+ t_n[i] = set4(T_EPSILON);
+ t_f[i] = load44(&rays.getMinT(i*4));
+ }
+ // for incomplete border sse-blocks, invalidate ray segment
+ for (int i=4*sse_begin;i<ray_begin;i++) {
+ ((float*)t_n)[i] = MAXT;
+ ((float*)t_f)[i] = -MAXT;
+ }
+ for (int i=ray_end;i<4*sse_end;i++) {
+ ((float*)t_n)[i] = MAXT;
+ ((float*)t_f)[i] = -MAXT;
+ }
+
+ sse_t anyValid = false4();
+ for (int i=sse_begin;i<sse_end;i++) {
+MANTA_UNROLL(3)
+ for (int k=0;k<3;k++) {
+ const sse_t rcp_ki = oneOver(load44(&rays.getDirection(4*i,k)));
+ const sse_t dist_b0k = mul4(sub4(set4(bounds[0][k]),
+
load44(&rays.getOrigin(4*i,k))),rcp_ki);
+ const sse_t dist_b1k = mul4(sub4(set4(bounds[1][k]),
+
load44(&rays.getOrigin(4*i,k))),rcp_ki);
+ t_n[i] = max4(t_n[i],min4(dist_b0k,dist_b1k));
+ t_f[i] = min4(t_f[i],max4(dist_b0k,dist_b1k));
+ }
+ anyValid = or4(anyValid,cmp4_le(t_n[i],t_f[i]));
+ }
+
+ if (getmask4(anyValid) == 0)
+ // no active ray
+ return;
+
+ const Vector shared_origin = rays.getOrigin(rays.begin());
+
+ struct StackEntry
+ {
+ float *freeMem;
+ sse_t *t_n;
+ sse_t *t_f;
+ int nodeID;
+ };
+
+ StackEntry stackBase[80];
+ StackEntry *stackPtr = stackBase;
+
+ int nodeID = 0;
+ while (1) {
+ const BSPNode &node = nodes[nodeID];
+ if (node.isLeaf()) {
+ for (size_t i=0; i<node.numPrimitives(); ++i) {
+#ifdef MAILBOX
+ if (mailbox.testAndMark(objects[node.objectsIndex+i]))
+ continue;
+#endif
+ mesh->get(objects[node.objectsIndex+i])->intersect(context, rays);
+#ifdef COLLECT_STATS
+ ++nTriIntersects;
+#endif
+ // objects[node.objectsIndex+i]->intersect(context, rays);
+ }
+
+#if 1
+ {
+ // early ray termination
+ sse_t allDied = _mm_true;
+MANTA_UNROLL(4)
+ for (int i=sse_begin;i<sse_end;i++)
+ {
+ sse_t thisOneDied =// and4(cmp4_le(t_n[i],t_f[i]), // valid
+
cmp4_le(load44(&rays.getMinT(i*4)),t_f[i]);
+ allDied = and4(allDied,thisOneDied);
+ }
+ if (getmask4(allDied) == 0xf)
+ {
+ //(thiago) this isn't catching all the early outs.
+ //cout << "ERT" << endl;
+ return;
+ }
+ }
+#endif
+
+ pop_node:
+ if (stackPtr <= stackBase)
+ {
+ return;
+ }
+ --stackPtr;
+
+ freeMem = stackPtr->freeMem;
+ t_n = stackPtr->t_n;
+ t_f = stackPtr->t_f;
+ // val = stackPtr->val;
+
+MANTA_UNROLL(4)
+ for (int i=sse_begin;i<sse_end;i++)
+ {
+ t_f[i] = min4(t_f[i],load44(&rays.getMinT(i*4)));
+ // val[i] = cmp4_le(t_n[i],t_f[i]);
+ // val[i] = and4(val[i],cmp4_le(t_n[i],t_f[i]));
+ }
+ nodeID = stackPtr->nodeID;
+ } else {
+ // inner node
+// const SplitPlane &plane = node.split;
+
+ int originSideNodeID;
+ int otherSideNodeID;
+ sse_t originSideFlag = cast4_i2f(set4i(0));
+ sse_t otherSideFlag = cast4_i2f(set4i(0));
+
+ if (node.isKDTree()) {
+#ifdef COLLECT_STATS
+ nTraversals++;
+ nKDTraversals++;
+#endif
+ const int planeDim = node.kdtreePlaneDim();
+#define ASSUME_CONSTANT_ORIGIN
+ // NOTE: ASSUME CONSTANT ORIGIN RIGHT NOW!!!
+ const float origin_in_plane_eq = node.d - shared_origin[planeDim];
+ const int origin_halfspace_sign = (origin_in_plane_eq < 0.f);
+
+ // Dot(shared_origin,node.normal)+node.d > 0.f);
+ originSideNodeID = node.children + origin_halfspace_sign;
+ otherSideNodeID = node.children + 1-origin_halfspace_sign;
+
+#ifdef ASSUME_CONSTANT_ORIGIN
+ const sse_t org_in_plane = set4(origin_in_plane_eq);
+#endif
+ MANTA_UNROLL(8)
+ for (int i=sse_begin;i<sse_end;i++) {
+ const sse_t rcp_dir = load44(&rays.getInverseDirection(4*i,
planeDim));
+
+#ifndef ASSUME_CONSTANT_ORIGIN
+ const sse_t normal_dot_origin = dot4(org_x[i],org_y[i],org_z[i],
+ normal_0, normal_1, normal_2);
+ const sse_t negative_org_in_plane =
+ sub4(set4(-node.d),normal_dot_origin);
+#endif
+ t_p[i] = mul4(org_in_plane, rcp_dir);
+ }
+ }
+ else {
+#ifdef COLLECT_STATS
+ ++nTraversals;
+#endif
+
+#define ASSUME_CONSTANT_ORIGIN
+ // NOTE: ASSUME CONSTANT ORIGIN RIGHT NOW!!!
+ const float origin_in_plane_eq
+ = shared_origin[0]*node.normal[0]
+ + shared_origin[1]*node.normal[1]
+ + shared_origin[2]*node.normal[2]
+ + node.d;
+
+ const int origin_halfspace_sign = (origin_in_plane_eq > 0.f);
+
+ // Dot(shared_origin,node.normal)+node.d > 0.f);
+ originSideNodeID = node.children + origin_halfspace_sign;
+ otherSideNodeID = node.children + 1-origin_halfspace_sign;
+
+#ifdef ASSUME_CONSTANT_ORIGIN
+ const sse_t negative_org_in_plane = set4(-origin_in_plane_eq);
+#endif
+ const Vector normal = node.normal;
+ const sse_t normal_0 = set4(normal[0]);
+ const sse_t normal_1 = set4(normal[1]);
+ const sse_t normal_2 = set4(normal[2]);
+
+ MANTA_UNROLL(8)
+ for (int i=sse_begin;i<sse_end;i++) {
+ // const float t_p = -(Dot(node.normal, org) + node.d) /
normal_dot_direction;
+ const sse_t normal_dot_direction =
dot4(load44(&rays.getDirection(4*i,0)),
+
load44(&rays.getDirection(4*i,1)),
+
load44(&rays.getDirection(4*i,2)),
+ normal_0, normal_1,
normal_2);
+#ifndef ASSUME_CONSTANT_ORIGIN
+ const sse_t normal_dot_origin = dot4(org_x[i],org_y[i],org_z[i],
+ normal_0, normal_1, normal_2);
+ const sse_t negative_org_in_plane =
+ sub4(set4(-node.d),normal_dot_origin);
+#endif
+ t_p[i] = mul4(negative_org_in_plane,
oneOver4(normal_dot_direction));
+ }
+ }
+
+ const sse_t rel_err = set4(1.0001f);
+
+ MANTA_UNROLL(8)
+ for (int i=sse_begin;i<sse_end;i++) {
+#ifdef SSE_REL_EPSILONS
+ const sse_t thisIsValid = cmp4_le(t_n[i],mul4(t_f[i],rel_err));
+#else
+ const sse_t thisIsValid = cmp4_le(t_n[i],t_f[i]);
+#endif
+ const sse_t facing_away = cmp4_le(t_p[i],zero4());
+ const sse_t org_in_plane = cmp4_eq(t_p[i],zero4());
+#ifdef SSE_REL_EPSILONS
+ const sse_t thisOnOrigin
+ = or4(facing_away, // faces away (in which case it's always on org
side) ...
+ cmp4_le(t_n[i],mul4(t_p[i],rel_err)) // or it points to the
plane but doesn't
+ );
+#else
+ const sse_t thisOnOrigin
+ = or4(facing_away, // faces away (in which case it's always on org
side) ...
+ cmp4_le(t_n[i],t_p[i]) // or it points to the plane but
doesn't
+ );
+#endif
+ originSideFlag = or4(originSideFlag,
+ and4(thisIsValid,
+ thisOnOrigin));
+
+#ifdef SSE_REL_EPSILONS
+ const sse_t thisOnOther = or4(org_in_plane,
+ andnot4(facing_away,
+ // _not_ facing away, i.e.
facing towrad, _and_
+ cmp4_ge(mul4(rel_err,
+ t_f[i]),
+ t_p[i]) // ray seg
does not end before plane
+ ));
+#else
+ const sse_t thisOnOther = or4(org_in_plane,
+ andnot4(facing_away,
+ // _not_ facing away, i.e.
facing towrad, _and_
+ cmp4_ge(t_f[i],t_p[i]) // ray
seg does not end before plane
+ ));
+#endif
+// val[i] = cmp4_le(t_n[i],t_f[i]);
+// originSideFlag = or4(originSideFlag, and4(val[i],thisOnOrigin));
+// otherSideFlag = or4(otherSideFlag, and4(val[i],thisOnOther));
+ otherSideFlag = or4(otherSideFlag,
+ and4(thisIsValid,
+ thisOnOther));
+ }
+#if 0
+ // FORCE traversing both children -- essetially tests if tree
+ // contains all tris, and if tritests work ...
+ const int originSide = 0xf;
+ const int otherSide = 0xf;
+#else
+ const int originSide = getmask4(originSideFlag);
+ const int otherSide = getmask4(otherSideFlag);
+#endif
+ if (originSide) {
+ if (otherSide) {
+ // both sides
+
+ sse_t *const frnt_t_f = (sse_t*)freeMem; freeMem +=
RayPacket::MaxSize;
+// sse_t * frnt_val = (sse_t*)freeMem; freeMem +=
RayPacket::MaxSize;
+ sse_t *const back_t_n = (sse_t*)freeMem; freeMem +=
RayPacket::MaxSize;
+MANTA_UNROLL(4)
+ for (int i=sse_begin;i<sse_end;i++) {
+
+ /*
+ - ray facing away: ray interval (n,f) unchanged on
+ origin side, invalidated (in whatever way) on other
+ side
+ - ray facing towards plane: (n,min(f,d)) on front side,
+ and (max(n,d),f) on back side that is, assuming all
+ rays agree on what origin side and other side are ...
+
+ inverted:
+ - origin interval is: (n,f) for rays facing away, and
+ (n,min(f,d)) for those pointing towards
+ - other interval is (max(n,d),f) for those pointing
+ towards, and (inf,f) for those pointing away
+ */
+ const sse_t facing_away = cmp4_le(t_p[i],zero4());
+ const sse_t org_in_plane = cmp4_eq(t_p[i],zero4());
+ frnt_t_f[i] = mask4(facing_away,t_f[i],min4(t_p[i],t_f[i]));
+ back_t_n[i] = mask4(andnot4(org_in_plane,facing_away),
+ set4(1e20f),
+ max4(t_p[i],t_n[i]));
+ }
+
+// stackPtr->val = back_val;
+ stackPtr->freeMem = freeMem;
+ stackPtr->t_n = back_t_n;
+ stackPtr->t_f = t_f;
+ stackPtr->nodeID = otherSideNodeID;
+ stackPtr++;
+
+ nodeID = originSideNodeID;
+ t_f = frnt_t_f;
+// val = frnt_val;
+
+ } else {
+ // _only_ origin side
+ nodeID = originSideNodeID;
+ }
+ } else {
+ if (otherSide) {
+ // _only_ other side
+ nodeID = otherSideNodeID;
+ } else {
+ goto pop_node;
+ exit(1);
+ }
+ }
+ }
+ }
+#else
+ cerr<<"oops, shouldn't have gotten to intersectPacket()\n";
+#endif
+}
+
+void BSP::intersectSingleRays(const RenderContext& context, RayPacket& rays)
const
+{
+ //we normalize the rays so that we get consistent distances when we
+ //do the epsilon test on the t value. kdtrees don't need that test,
+ //so they also don't need to normalized rays.
+ rays.normalizeDirections();
+ TrvStack trvStack[80];
+
+ rays.computeInverseDirections();
+ for(int r=rays.begin(); r<rays.end(); r++) {
+
+ RayPacket singleRayPacket(rays, r, r+1);
+ Vector origin = rays.getOrigin( r );
+ Vector direction = rays.getDirection( r );
+ Vector inverse_direction = rays.getInverseDirection( r );
+
traverse(context,singleRayPacket,origin,direction,inverse_direction,trvStack);
+ }
+}
+
+void BSP::traverse(const RenderContext &context, RayPacket &ray,
+ const Vector &org, const Vector &dir, const Vector &rcp,
TrvStack *const stackBase)
+ const
+{
+#ifdef COLLECT_STATS
+ nTotalRays += ray.end() - ray.begin();
+#endif
+
+ float t_n = T_EPSILON;
+ float t_f = ray.getMinT( ray.begin() );
+
+ const int signs[3] = {dir[0] < 0, dir[1] < 0, dir[2] < 0};
+
+ for (int k=0;k<3;k++) {
+ if (signs[k]) {
+ const float k0 = (bounds[1][k] - org[k]) * rcp[k];
+ t_n = Max(t_n,k0);
+ const float k1 = (bounds[0][k] - org[k]) * rcp[k];
+ t_f = Min(t_f,k1);
+ } else {
+ const float k0 = (bounds[0][k] - org[k]) * rcp[k];
+ t_n = Max(t_n,k0);
+ const float k1 = (bounds[1][k] - org[k]) * rcp[k];
+ t_f = Min(t_f,k1);
+ }
+ if (t_n > t_f)
+ return;
+ };
+ TrvStack *stackPtr = stackBase;
+ int nodeID = 0;
+
+#ifdef BSP_INTERSECTION
+ int depth = 0;
+ int nearPlane = -1;
+ int farPlane = -1;
+
+#endif
+
+ while (1) {
+ const BSPNode &node = nodes[nodeID];
+ if (node.isLeaf()) {
+#ifdef COLLECT_STATS
+ nLeafs++;
+#endif
+ for (size_t i=0; i<node.numPrimitives(); ++i) {
+#ifdef BSP_INTERSECTION
+ //Is the closest splitting plane that the ray intersected the
+ //triangle face plane? If so, and we've also intersected the
+ //planes that go through the 3 triangle edges, then we must
+ //have intersected the triangle!
+
+ /*
+ //This is actually slower and results in some errors in the cylinder.
+ if (node.numPrimitives() == 1 && ((node.triangleBounds & 0x1) ==
0x1)) {
+
+#ifdef COLLECT_STATS
+ nBSPTriIntersects++;
+#endif
+ const MeshTriangle* const tri =
mesh->get(objects[node.objectsIndex]);
+ const Vector v0 = tri->getVertex(0);
+ const Vector v1 = tri->getVertex(1);
+ const Vector v2 = tri->getVertex(2);
+
+ const Vector normal = Cross(v1 - v0, v2 - v0);
+ const float normal_dot_direction = Dot(normal, dir);
+ const float inv_dot = 1.f/normal_dot_direction;
+ const float d = Dot(normal, -v0);
+ const float t_p = -(Dot(normal, org) + d) * inv_dot; ///
normal_dot_direction;
+
+ if (t_p+T_EPSILON >= t_n && t_p <= t_f+T_EPSILON) {
+ ray.hit(ray.begin(), t_p,
mesh->materials[mesh->face_material[0]],
+ tri, (TexCoordMapper*)tri);
+ return;
+ }
+ }
+ */
+
+ if (node.triangleBounds == 0x3) {
+ if (nearPlane == node.getTriPlaneDepth()) {
+#ifdef COLLECT_STATS
+ nBSPTriIntersects++;
+#endif
+ const MeshTriangle* const tri =
mesh->get(objects[node.objectsIndex]);
+ ray.hit(ray.begin(), t_n,
mesh->materials[mesh->face_material[objects[node.objectsIndex]]],
+ tri, (TexCoordMapper*)tri);
+ }
+ else if (farPlane == node.getTriPlaneDepth()) {
+#ifdef COLLECT_STATS
+ nBSPTriIntersects++;
+#endif
+ const MeshTriangle* const tri =
mesh->get(objects[node.objectsIndex]);
+ ray.hit(ray.begin(), t_f,
mesh->materials[mesh->face_material[objects[node.objectsIndex]]],
+ tri, (TexCoordMapper*)tri);
+ }
+ else if (1) {
+ //In an ideal world, if both the near and far planes don't
+ //correspond with the triangle, then we know the ray
+ //missed. But in practice there might be multiple planes
+ //with the same hit point (a plane gets used twice or
+ //numerical precision issues occur), so we still need to
+ //check the misses.
+#ifdef COLLECT_STATS
+ nTriIntersects++;
+#endif
+// const float oldT = ray.getMinT( ray.begin() );
+ mesh->get(objects[node.objectsIndex+i])->intersect(context, ray);
+// const float newT = ray.getMinT( ray.begin() );
+// if (newT != oldT) {
+// //we got a hit that we didn't expect!!!!
+// }
+ }
+ }
+
+ else {
+#ifdef COLLECT_STATS
+ nTriIntersects++;
+#endif
+ mesh->get(objects[node.objectsIndex+i])->intersect(context, ray);
+ }
+#else
+#ifdef COLLECT_STATS
+ nTriIntersects++;
+#endif
+
+ mesh->get(objects[node.objectsIndex+i])->intersect(context, ray);
+#endif //BSP_INTERSECTION
+ }
+ if (stackPtr <= stackBase)
+ return;
+ --stackPtr;
+ const float newT = ray.getMinT( ray.begin() );
+ if (newT <= t_f)
+ return;
+ nodeID = stackPtr->nodeID;
+
+ t_n = t_f;
+#ifdef BSP_INTERSECTION
+ nearPlane = farPlane;
+
+ if (newT < stackPtr->t) {
+ t_f = newT;
+ //t_f does not correspond to a plane (it hit a primitive not
+ //on a split)
+ farPlane = -1;
+ }
+ else {
+ t_f = stackPtr->t;
+ farPlane = stackPtr->farPlane;
+ }
+ depth = stackPtr->depth+1;
+#else
+ t_f = Min(newT,stackPtr->t);
+#endif //BSP_INTERSECTION
+ }
+ else if (node.isKDTree()) {
+#ifdef COLLECT_STATS
+ nTraversals++;
+ nKDTraversals++;
+#endif
+ const int planeDim = node.kdtreePlaneDim();
+ const float t_p = (node.d - org[planeDim]) * rcp[planeDim];
+
+ const int frontChild = node.children+signs[planeDim];
+ const int backChild = node.children+1-signs[planeDim];
+
+ //TODO: This T_EPSILON stuff is making the far/near plane
+ //triangle intersection not work well. Should see if there's a
+ //good way to fix this. The problem is that this causes us to
+ //have the incorrect far/near plane and we require that one of
+ //these two planes be the triangle plane in order for an
+ //intersection to occur.
+
+ //A solution would be to keep a list of nearPlanes and farPlanes
+ //that are close (within epsilon) to t_n and t_f.
+
+ if (t_p < t_n - 1e-5f) {
+ nodeID = backChild;
+ }
+ else if (t_p > t_f + 1e-5f) {
+ nodeID = frontChild;
+ }
+ else {
+ stackPtr->nodeID = backChild;
+ stackPtr->t = t_f;
+#ifdef BSP_INTERSECTION
+ stackPtr->depth = depth;
+ stackPtr->farPlane = farPlane;
+ farPlane = depth;
+#endif
+ ++stackPtr;
+ nodeID = frontChild;
+ t_f = t_p;
+ }
+#ifdef BSP_INTERSECTION
+ depth++;
+#endif
+ }
+ else {
+#ifdef COLLECT_STATS
+ nTraversals++;
+#endif
+ const float normal_dot_direction = Dot(node.normal, dir);
+ const float inv_dot = 1.f/normal_dot_direction;
+ const float t_p = -(Dot(node.normal, org) + node.d) * inv_dot; ///
normal_dot_direction;
+
+ const int side = ((unsigned int &)normal_dot_direction) >> 31;
+ const int frontChild = node.children+side; //normal_dot_direction >
0.f ? node.children : node.children+1;
+ const int backChild = node.children+1-side;//normal_dot_direction >
0.f ? node.children+1 : node.children;
+
+ if (t_p < t_n - 1e-5f) {
+ nodeID = backChild;
+ }
+ else if (t_p > t_f + 1e-5f) {
+ nodeID = frontChild;
+ }
+ else {
+ stackPtr->nodeID = backChild;
+ stackPtr->t = t_f;
+#ifdef BSP_INTERSECTION
+ stackPtr->depth = depth;
+ stackPtr->farPlane = farPlane;
+ farPlane = depth;
+#endif
+ ++stackPtr;
+ nodeID = frontChild;
+ t_f = t_p;
+ }
+#ifdef BSP_INTERSECTION
+ ++depth;
+#endif
+ }
+ }
+};
+
+void BSP::setTriangleBounds(vector<BuildSplitPlane> &splitPlanes, int
&nodeID,
+ int &nextFreeNode)
+{
+#ifdef BSP_INTERSECTION
+ assert(nodes[nodeID].numPrimitives() == 1);
+ const size_t triID = objects[nodes[nodeID].objectsIndex];
+
+ //The node needs to be contained in the triangle. In otherwords, it
+ //can be equal to the triangle, or it can be inside the triangle.
+ //To be more precise, if one of the faces of the node is completly
+ //contained within the triangle, then we can intersect that triangle
+ //purely through the traversal.
+
+ //This also means that it is ok to have very large nodes (high
+ //Surface Area) as long as one of the faces corresponds to a
+ //triangle, since any further subdividing to get empty nodes would
+ //not result in any fewer triangle or node tests (we aren't doing
+ //triangle tests after all).
+
+ //TODO: modify SAH criteria so that these kinds of nodes have
+ //priority. Perhaps make C_isec zero.
+
+ unsigned int triangleBounds = 0;
+
+#if 1
+ Polygon triangle;
+ triangle.push_back(vertices[triID][0]);
+ triangle.push_back(vertices[triID][1]);
+ triangle.push_back(vertices[triID][2]);
+ for (size_t f=0; f < polytope.faces.size(); ++f) {
+ float intersectedArea = 0;
+ if (overlaps(triangle, polytope.faces[f], intersectedArea, true)) {
+ if (almostEqualRelative(intersectedArea, getArea(polytope.faces[f]))) {
+ for (size_t i=0; i < splitPlanes.size(); ++i) {
+ bool pointsOnPlane[3];
+ for (int k=0; k < 3; ++k)
+ pointsOnPlane[k] = fabs(signedDistance(splitPlanes[i],
vertices[triID][k])) < 1e-6;
+ if (pointsOnPlane[0] && pointsOnPlane[1] && pointsOnPlane[2]) {
+ nodes[nodeID].setTriPlaneDepth(i);
+ nodes[nodeID].triangleBounds = 0x3;
+ triangleBounds = 0xf;
+// cerr << nodes[nodeID].triangleBounds<<endl;
+// return;
+ }
+ }
+// cerr << nodes[nodeID].triangleBounds<<endl;
+// return;
+ }
+// else
+// cout <<intersectedArea << " " << getArea(polytope.faces[f])<< "
"
+// <<intersectedArea / getArea(polytope.faces[f])<<endl<<endl;
+ }
+ }
+ //#else
+ //find which splits correspond with the triangle edges and plane.
+ for (size_t i=0; i < splitPlanes.size(); ++i) {
+ bool pointsOnPlane[3];
+ for (int k=0; k < 3; ++k)
+ pointsOnPlane[k] = fabs(signedDistance(splitPlanes[i],
vertices[triID][k])) < 1e-6;
+
+ if (pointsOnPlane[0] && pointsOnPlane[1] && pointsOnPlane[2]) {
+ //if all three edges are on the split plane, then the triangle
+ //is on the plane.
+
+ //This means during traversal we already calculated the t value.
+ //If all the other sides have been tested and this is the
+ //closest intersection point, then we know for sure we have an
+ //intersection.
+ if (true||(triangleBounds & 0x8) == 0) {
+ nodes[nodeID].setTriPlaneDepth(i);
+ triangleBounds |= 0x8;
+ }
+ }
+ else if (pointsOnPlane[0] && pointsOnPlane[1]) {
+ triangleBounds |= 0x4;
+ }
+ else if (pointsOnPlane[0] && pointsOnPlane[2]) {
+ triangleBounds |= 0x2;
+ }
+ else if (pointsOnPlane[1] && pointsOnPlane[2]) {
+ triangleBounds |= 0x1;
+ }
+ }
+
+#endif
+
+ unsigned int numSplits = splitPlanes.size();
+
+ if (triangleBounds == 0xe ||
+ triangleBounds == 0xd ||
+ triangleBounds == 0xb ||
+ triangleBounds == 0x7) {
+ numSplits++;
+ //We are missing just one split. It's definitely cheaper to just
+ //add this extra node traversal than doing a full triangle test
+
+ BuildSplitPlane plane;
+ const Point triNormal = getTriangleNormal(vertices[triID][0],
+ vertices[triID][1],
+ vertices[triID][2]).normal();
+ Point p1, p2, p3;
+ if (triangleBounds == 0xe) {
+ //missing 0x1 == vertices 1 and 2
+ p1 = vertices[triID][1];
+ p2 = vertices[triID][2];
+ p3 = vertices[triID][0];
+ }
+ else if (triangleBounds == 0xd) {
+ //missing 0x2 == vertices 0 and 2
+ p1 = vertices[triID][0];
+ p2 = vertices[triID][2];
+ p3 = vertices[triID][1];
+ }
+ else if (triangleBounds == 0xb) {
+ //missing 0x4 == vertices 0 and 1
+ p1 = vertices[triID][0];
+ p2 = vertices[triID][1];
+ p3 = vertices[triID][2];
+ }
+
+ if (triangleBounds == 0x7) {
+ //missing 0x8 == triangle autopartition
+ plane.normal = triNormal;
+ p3 = vertices[triID][0];
+ plane.d = Dot(plane.normal, -p3);
+ nodes[nodeID].setTriPlaneDepth(splitPlanes.size());
+ }
+ else {
+ plane.normal = Cross(triNormal, p2-p1).normal();
+ plane.d = Dot(plane.normal, -p1);
+ }
+
+ if (triangleBounds == 0xf)
+ nodes[nodeID].triangleBounds = 0x3;
+ else if (triangleBounds == 0x8)
+ nodes[nodeID].triangleBounds = 0x2;
+ else if (triangleBounds == 0x7)
+ nodes[nodeID].triangleBounds = 0x1;
+ else nodes[nodeID].triangleBounds = 0x0;
+
+ int leftChild = nextFreeNode++;
+ int rightChild = nextFreeNode++;
+ nodes.resize(nextFreeNode);
+ BSPNode leafNode = nodes[nodeID];
+ nodes[nodeID].makeInternal(leftChild, plane.normal, plane.d);
+
+ double p3_dist = signedDistance(plane, p3);
+ if (p3_dist > 0) {
+ swap(leftChild, rightChild);
+ }
+
+ nodes[leftChild] = leafNode;
+ nodes[leftChild].triangleBounds = 0x3;
+
+ nodes[rightChild].makeLeaf(objects.size(), 0);
+
+ nodeID = leftChild;
+ }
+
+ return;
+
+ if ((nodes[nodeID].triangleBounds == 0x3 ||
+ nodes[nodeID].triangleBounds == 0x1) &&
+ nodes[nodeID].getTriPlaneDepth() < (int) splitPlanes.size()-1) {
+// cout <<"hi\n";
+ //add triangle autopartition
+ BuildSplitPlane plane;
+ const Point triNormal = getTriangleNormal(vertices[triID][0],
+ vertices[triID][1],
+ vertices[triID][2]).normal();
+ plane.normal = triNormal;
+ plane.d = Dot(plane.normal, -vertices[triID][0]);
+ nodes[nodeID].setTriPlaneDepth(splitPlanes.size());
+
+ int leftChild = nextFreeNode++;
+ int rightChild = nextFreeNode++;
+ nodes.resize(nextFreeNode);
+ BSPNode leafNode = nodes[nodeID];
+ nodes[nodeID].makeInternal(leftChild, plane.normal, plane.d);
+
+ nodes[leftChild] = leafNode;
+ nodes[leftChild].triangleBounds = 0x3;
+ nodes[leftChild].setTriPlaneDepth(numSplits);
+
+ nodes[rightChild].makeLeaf(objects.size(), 0);
+
+ nodeID = leftChild;
+ }
+#endif
+// cerr << nodes[nodeID].triangleBounds<<endl;
+}
+
+void BSP::build(int nodeID, int &nextFreeNode,
+ vector<BuildSplitPlane> &splitPlanes, BSH &primitives, int
depth)
+{
+ nodes.resize(nextFreeNode);
+
+ bool forceCreateLeaf = false;
+ if (primitives.size()-primitives.numSplits() <=0){
+ createLeaf:
+ if (primitives.numSplits() > 0) {
+ nodes[nodeID].makeLeaf(objects.size(),
+ primitives.size()-primitives.numSplits());
+ BSH::Iterator iter = primitives.begin();
+ objects.push_back(iter.get().originalTriID);
+ for (iter.next(); iter != primitives.end(); iter.next()) {
+ if (iter.get().originalTriID != objects.back())
+ objects.push_back(iter.get().originalTriID);
+ }
+ }
+ else {
+ nodes[nodeID].makeLeaf(objects.size(), primitives.size());
+ for (BSH::Iterator iter = primitives.begin();
+ iter != primitives.end(); iter.next()) {
+ objects.push_back(iter.get().originalTriID);
+ }
+ }
+
+ //Doing this seems to improve performance for single ray
+ //traversal. Packet traversal on the other hand does better with
+ //shallow trees.
+#ifdef BSP_INTERSECTION
+ if (primitives.size()-primitives.numSplits() == 1) {
+ setTriangleBounds(splitPlanes, nodeID, nextFreeNode);
+
+ if (!forceCreateLeaf && nodes[nodeID].triangleBounds != 0x3) {
+ objects.pop_back();
+ goto createInteriorNode;
+ }
+ }
+#endif
+
+ }
+ else {
+ createInteriorNode:
+
+ Location coplanar_primitivesSide;
+ BuildSplitPlane plane = getBuildSplitPlane(primitives, splitPlanes,
coplanar_primitivesSide);
+ if (plane.normal == Point(0,0,0)) {
+ forceCreateLeaf = true;
+ goto createLeaf;
+ }
+
+ splitPlanes.push_back(plane);
+
+ BSH posPrimitives(primitives);
+ primitives.split(splitPlanes.back(), posPrimitives,
coplanar_primitivesSide);
+
+ int leftChild = nextFreeNode++;
+ int rightChild = nextFreeNode++;
+ nodes[nodeID].makeInternal(leftChild, plane.normal, plane.d);
+
+ Polytope side1, side2;
+ polytope.split(splitPlanes.back(), side1, side2);
+
+ polytope = side1;
+ polytope.updateVertices();
+
+ build(leftChild, nextFreeNode, splitPlanes, primitives, depth+1);
+
+ polytope = side2;
+ polytope.updateVertices();
+
+ splitPlanes.resize(depth);
+ splitPlanes.back().flip();
+
+ build(rightChild, nextFreeNode, splitPlanes, posPrimitives, depth+1);
+ }
+}
+
+BuildSplitPlane BSP::getBuildSplitPlane(BSH &primitives,
+ vector<BuildSplitPlane> &splitPlanes,
+ Location &coplanar_primitivesSide)
+{
+ //I am assuming we just have triangles!!!
+// return getRandomSplit(primitives, splitPlanes,
coplanar_primitivesSide);
+ return getSAHSplit(primitives, splitPlanes, coplanar_primitivesSide);
+}
+
+void BSP::printPrimitivesAsObj(const BSH &primitives) const
+{
+ for (BSH::Iterator iter = primitives.begin();
+ iter != primitives.end(); iter.next())
+ for (size_t k = 0; k < 3; ++k)
+ cout << "v "<<iter.get().tri[k]<<endl;
+ //cout << "v " << vertices[iter.get().originalTriID][k] <<"\n";
+ for (size_t i = 0; i < primitives.size(); ++i)
+ cout << "f " << i*3 +1<< " " << i*3+2 << " " <<i*3+3 <<endl;
+}
+void BSP::printOriginalPrimitivesAsObj(const BSH &primitives) const
+{
+ for (BSH::Iterator iter = primitives.begin();
+ iter != primitives.end(); iter.next())
+ for (size_t k = 0; k < 3; ++k)
+ //cout << "v "<<iter.get().tri[k]<<endl;
+ cout << "v " << vertices[iter.get().originalTriID][k] <<"\n";
+ for (size_t i = 0; i < primitives.size(); ++i)
+ cout << "f " << i*3 +1<< " " << i*3+2 << " " <<i*3+3 <<endl;
+}
+
+
+const float C_isec = 75;
+const float C_KDTree_trav = 10;
+const float C_BSP_trav = 20;
+
+inline bool isKDTreeSplit(const BuildSplitPlane &splitPlane) {
+ return (splitPlane.normal.minComponent() < -0.9999999 ||
+ splitPlane.normal.maxComponent() > 0.9999999);
+}
+
+BuildSplitPlane BSP::getSAHSplit(BSH &primitives,
+ vector<BuildSplitPlane> &splitPlanes,
+ Location &coplanar_primitivesSide)
+{
+ float parentArea;
+ parentArea = polytope.getSurfaceArea();
+// cout <<"parent Area: " << parentArea <<endl;
+
+ SplitData bestSplits[2]; //first is for bsp splits, second for kd.
+ bestSplits[0].cost = bestSplits[1].cost =
+ (primitives.size()-primitives.numSplits())*C_isec;
+
+ //if the node is really small, just stop. (is this really needed?)
+ if (parentArea < 1e-10) {
+ BuildSplitPlane plane;
+ plane.normal = Point(0,0,0);
+ plane.d = 0;
+ return plane;
+ }
+
+ splitPlanes.resize(splitPlanes.size()+1);
+
+ int previousTriID=-1;
+ const int lastPass = 16;
+
+ set<BuildSplitPlane> splitPlaneSet;
+
+ for (BSH::Iterator iter = primitives.begin();
+ iter!=primitives.end(); iter.next()) {
+
+ SplitData splitdata;
+ splitdata.splitPrimitive = iter.get().originalTriID;
+
+ for (int pass=0; pass < lastPass; ++pass) {
+ Polygon splitPoints;
+
+ splitdata.splitType = pass;
+
+ //set the split to use the triangle face.
+ splitdata.split = iter.get().plane;
+
+ if (pass == 0 && iter.get().originalTriID != previousTriID) {
+ //try triangle plane In theory, the originalTri normal's and
+ //the current tri normal are the same. In practice though, the
+ //current tri normal might have more error since the triangle
+ //is smaller, so we use the original normal.
+
+ splitPoints.push_back(vertices[iter.get().originalTriID][0]);
+ splitPoints.push_back(vertices[iter.get().originalTriID][1]);
+ splitPoints.push_back(vertices[iter.get().originalTriID][2]);
+ }
+ else if (pass > 0 && pass < 4 && iter.get().originalTriID !=
previousTriID
+ ) {
+ //use plane on triangle edge that is orthogonal to triangle
+ //normal Edges that are on the bounding polytope (and
+ //obviously those outside it), are useless.
+ Point p1, p2;
+ switch ((pass-1)%3) {
+ case 0:
+#if 1
+ p1 = vertices[iter.get().originalTriID][0];//iter.get().tri[0];
+ p2 = vertices[iter.get().originalTriID][1];//iter.get().tri[1];
+ break;
+ case 1:
+ p1 = vertices[iter.get().originalTriID][1];//iter.get().tri[1];
+ p2 = vertices[iter.get().originalTriID][2];//iter.get().tri[2];
+ break;
+ case 2:
+ p1 = vertices[iter.get().originalTriID][2];//iter.get().tri[2];
+ p2 = vertices[iter.get().originalTriID][0];//iter.get().tri[0];
+ break;
+#else
+ p1 = iter.get().tri[0];
+ p2 = iter.get().tri[1];
+ break;
+ case 1:
+ p1 = iter.get().tri[1];
+ p2 = iter.get().tri[2];
+ break;
+ case 2:
+ p1 = iter.get().tri[2];
+ p2 = iter.get().tri[0];
+ break;
+
+#endif
+ }
+
+ splitdata.split.normal = Cross(splitdata.split.normal, p2 -
p1).normal();
+ if (pass > 3) {
+ AffineTransformT<double> transform;
+ transform.initWithRotation(p2-p1, (.5-drand48())*2.5);
+ splitdata.split.normal =
transform.multiply_vector(splitdata.split.normal).normal();
+ }
+
+ splitdata.split.d = Dot(splitdata.split.normal, -p1);
+
+ splitPoints.push_back(p1);
+ splitPoints.push_back(p2);
+ }
+ else if (pass > 6 && pass < 16) {
+ const int k = (pass-7)%3;
+ const int v = (pass-7)/3;
+
+#if 1
+ //We are only interested in spliting planes that don't split
+ //the triangle (i.e. use those only on the triangle boundary).
+ const double split = vertices[iter.get().originalTriID][v][k];
+ if ((split - vertices[iter.get().originalTriID][(v+1)%3][k])*
+ (split - vertices[iter.get().originalTriID][(v+2)%3][k]) < 0)
+#else
+ const double split = iter.get().tri[v][k];
+ if ((split - iter.get().tri[(v+1)%3][k])*
+ (split - iter.get().tri[(v+2)%3][k]) < 0)
+#endif
+ continue;
+
+ splitdata.split.d = -split;
+ splitdata.split.normal = Point(0,0,0); //set it all to 0
+ splitdata.split.normal[k] = 1; //now set one of them to 1
+ }
+ else
+ continue;
+
+ //Check to see if we've already tried this split. If so, don't
+ //bother testing it again.
+ if (splitPlaneSet.insert(splitdata.split).second == false)
+ continue;
+
+ //set to false to not use getBuildSplitPlane
+ if (false and !(pass == 0 && iter.get().fixedPlane) ) //are we using
the fixed plane?
+ splitdata.split = getBestSplitPlane(primitives, splitdata.split,
+ splitPoints);
+
+ splitPlanes.back() = splitdata.split;
+
+ bool success = getSurfaceAreas(splitdata);
+ if (!success) {
+ //Either something went bad or the split was useless (cut off
+ //only a vertex or an edge of the polytope).
+ continue;
+ }
+ calculateSAH(parentArea, splitdata, primitives, polytope);
+
+// printf("%d - cost: %f areas: %f %f counts: %d %d
probabilities are %f %f \t%f\n",
+// pass, splitdata.cost, splitdata.lArea, splitdata.rArea,
splitdata.nLeft, splitdata.nRight, splitdata.lArea/parentArea,
splitdata.rArea/parentArea,
+// (splitdata.lArea+splitdata.rArea)/parentArea);
+
+
+ if (polytope.debugPrint) {
+ cout <<"\n\ndebug print:\n";
+ printPrimitivesAsObj(primitives);
+ cout <<endl;
+ Polytope side1, side2;
+ Polytope poly;
+ poly.initialize(bounds);
+ poly.printAsObj(primitives.size()*3);
+ for (size_t i=6; i < splitPlanes.size(); ++i) {
+ cout <<"\nsplit " << i<<endl;
+ side1.clear();
+ side2.clear();
+ poly.split(splitPlanes[i], side1, side2);
+ side1.printAsObj(primitives.size()*3);
+ side2.printAsObj(primitives.size()*3);
+ cout <<endl;
+ cout << "plane " << splitPlanes[i].normal<< " "
<<splitPlanes[i].d<<endl;
+ poly = side1;
+ }
+ cout <<"\n\ninput mesh:\n";
+ printPrimitivesAsObj(primitives);
+
+ exit(1);
+ }
+
+ if (splitdata.lArea*.98 > parentArea || splitdata.rArea*.98 >
parentArea) {
+ //This should never happen. If it does, the surface area
+ //calculation got screwed up. If we use this split plane, bad
+ //things might later happen.
+#ifdef BUILD_DEBUG
+ printf( "%f+%f = %f\n", splitdata.lArea/parentArea,
splitdata.rArea/parentArea, (splitdata.lArea+splitdata.rArea)/parentArea);
+ printf("%d Split has cost %f and #left/right are: %d %d\n",
+ iter.get().originalTriID, splitdata.cost, splitdata.nLeft,
splitdata.nRight);
+
+ cout.precision(20);
+ Polytope side1, side2;
+ polytope.split(splitdata.split, side1, side2);
+ polytope.printAsObj(primitives.size()*3);
+ cout << "plane " << splitdata.split.normal<< " "
<<splitdata.split.d<<endl;
+ side1.printAsObj(primitives.size()*3);
+ side2.printAsObj(primitives.size()*3);
+ cout << "parentArea: " <<parentArea<<endl;
+ cout << "side1 area: " <<lArea<<endl;
+ cout << "side2 area: " <<rArea<<endl;
+
+ Polytope poly;
+ poly.initialize(bounds);
+ poly.printAsObj(primitives.size()*3);
+ for (size_t i=6; i < splitPlanes.size(); ++i) {
+ cout <<"\nsplit " << i<<endl;
+ side1.clear();
+ side2.clear();
+ poly.split(splitPlanes[i], side1, side2);
+ cout<<"area: "<<side1.getSurfaceArea() << " and is
flat()="<<side1.isFlat()<<endl;
+ side1.printAsObj(primitives.size()*3);
+ cout<<"area: "<<side2.getSurfaceArea() << " and is
flat()="<<side2.isFlat()<<endl;
+ side2.printAsObj(primitives.size()*3);
+ cout <<endl;
+ cout << "plane " << splitPlanes[i].normal<< " "
<<splitPlanes[i].d<<endl;
+
+ poly = side1;
+ }
+
+ cout <<"\n\ninput mesh:\n";
+ printPrimitivesAsObj(primitives);
+ if (parentArea > 1e-7) //smaller than this and there's
+ //probably not much we can do. (maybe compute
convex hull??)
+ exit(1);
+#endif
+ continue;
+ }
+
+ if (splitdata.lArea + splitdata.rArea <= parentArea*0.99) {
+ //it is legitimate (but still useless) for one area to be 1
+ //and another 0 (split along polytope edge for instance).
+ if (splitdata.lArea < parentArea*.99 && splitdata.rArea <
parentArea*.99) {
+ //This should never happen. If it does, the surface area
+ //calculation got screwed up. If we use this split plane, bad
+ //things might later happen.
+ cout.precision(20);
+ printf( "%f+%f = %f (area too small)\n",
splitdata.lArea/parentArea, splitdata.rArea/parentArea,
(splitdata.lArea+splitdata.rArea)/parentArea);
+ printf("%d Split has cost %f and #left/right are: %d %d\n",
+ iter.get().originalTriID, splitdata.cost, splitdata.nLeft,
splitdata.nRight);
+
+#ifdef BUILD_DEBUG
+ Polytope side1, side2;
+ Polytope poly;
+ poly.initialize(bounds);
+
+ for (size_t i=6; i < splitPlanes.size(); ++i) {
+ cout <<"\nsplit " << i<<endl;
+ side1.clear();
+ side2.clear();
+ poly.split(splitPlanes[i], side1, side2);
+ side1.printAsObj(primitives.size()*3);
+ side2.printAsObj(primitives.size()*3);
+ cout <<endl;
+ cout << "plane " << splitPlanes[i].normal<< " "
<<splitPlanes[i].d<<endl;
+ poly = side1;
+ }
+
+ cout <<"\n\ndebug print:\n";
+ printPrimitivesAsObj(primitives);
+ cout <<endl;
+
+ if (parentArea > 1e-7)
+ exit(1);
+#endif
+ continue;
+ }
+ }
+
+
+ if (isKDTreeSplit(splitdata.split)) {
+ if ( splitdata.cost < bestSplits[1].cost ||
+ (splitdata.cost == bestSplits[1].cost &&
+ bestSplits[1].splitPrimitive < splitdata.splitPrimitive))
+ bestSplits[1] = splitdata;
+ }
+ else {
+ if (splitdata.cost < bestSplits[0].cost ||
+ (splitdata.cost == bestSplits[0].cost &&
+ bestSplits[0].splitPrimitive < splitdata.splitPrimitive))
+ bestSplits[0] = splitdata;
+ }
+ }
+ previousTriID = splitdata.splitPrimitive;
+ }
+
+
+// if (bestSplits[0].splitPrimitive >= 0)
+// printf("\nSplit has cost %f and #left/right are: %d %d and child
probabilities are %f %f. C_bsp=%f\n",
+// bestSplits[0].cost, bestSplits[0].nLeft, bestSplits[0].nRight,
bestSplits[0].lArea/parentArea, bestSplits[0].rArea/parentArea, C_BSP_trav);
+// if (bestSplits[1].splitPrimitive >= 0)
+// printf("Split has cost %f and #left/right are: %d %d and child
probabilities are %f %f\n",
+// bestSplits[1].cost, bestSplits[1].nLeft, bestSplits[1].nRight,
bestSplits[1].lArea/parentArea, bestSplits[1].rArea/parentArea);
+
+
+ //Now that we've tested all the split candidates, we need to figure
+ //out whether the best kdtree split is better than the best bsp
+ //split.
+
+ //First lets check whether we can use the varying BSP cost.
+ if (bestSplits[1].splitPrimitive >= 0) {
+ //We found a valid kd-tree split, so lets use the varying BSP cost.
+
+ const float C_trav = .1*C_isec*(primitives.size() - 2) + C_KDTree_trav;
+ const float lProb = bestSplits[0].lArea / parentArea;
+ const float rProb = bestSplits[0].rArea / parentArea;
+
+ bestSplits[0].cost = (bestSplits[0].nLeft*lProb +
+ bestSplits[0].nRight*rProb)*C_isec + C_trav;
+ }
+ else {
+ bestSplits[0].cost += C_BSP_trav;
+ if (bestSplits[0].cost >=
(primitives.size()-primitives.numSplits())*C_isec)
+ bestSplits[0].splitPrimitive = -1;
+ }
+
+ SplitData bestSplit;
+ if (bestSplits[0].cost <= bestSplits[1].cost)
+ bestSplit = bestSplits[0];
+ else
+ bestSplit = bestSplits[1];
+
+ splitPlanes.resize(splitPlanes.size()-1);
+
+ if (primitives.size()/(primitives.size()-primitives.numSplits()) >= 100) {
+ cout <<primitives.numSplits()<< " split triangles "
<<primitives.size()<<endl;
+ //I don't know why this is happening. But if it does, just get out.
+ if (primitives.size()/(primitives.size()-primitives.numSplits()) >= 600)
{
+ bestSplit.split.normal = Point(0,0,0);
+ return bestSplit.split;
+ }
+ }
+
+ if (bestSplit.splitPrimitive < 0) {
+// cout <<"found nothing!"<<endl;
+
+ if (primitives.size() - primitives.numSplits() >= 16) {
+ cout <<primitives.size() - primitives.numSplits()
+ <<"too many primitives in leaf!\n";
+// cout.precision(10);
+// printPrimitivesAsObj(primitives);
+
+// Polytope side1, side2;
+// polytope.split(splitdata.split, side1, side2);
+// polytope.printAsObj(primitives.size()*3);
+// side1.printAsObj(primitives.size()*3);
+// side2.printAsObj(primitives.size()*3);
+ //exit(1);
+ }
+
+ bestSplit.split.normal = Point(0,0,0);
+ return bestSplit.split;
+ }
+
+ coplanar_primitivesSide = bestSplit.coplanar_side;
+
+ //axis aligned splits need to have positive normals.
+ if (isKDTreeSplit(bestSplit.split) && //is it a kd-tree split?
+ bestSplit.split.normal.maxComponent() < 0.1) { //is it negative?
+ bestSplit.split.normal*=-1;
+ bestSplit.split.d*=-1;
+ if (coplanar_primitivesSide == negativeSide)
+ coplanar_primitivesSide = positiveSide;
+ else if (coplanar_primitivesSide == positiveSide)
+ coplanar_primitivesSide = negativeSide;
+ }
+
+// printf("\n%f,%f,%f Split has cost %f and #left/right are: %d %d and
child probabilities are %f %f\n",
+//
bestSplit.split.normal[0],bestSplit.split.normal[1],bestSplit.split.normal[2],
bestSplit.cost, bestSplit.nLeft, bestSplit.nRight,
bestSplit.lArea/parentArea, bestSplit.rArea/parentArea);
+
+ return bestSplit.split;
+ }
+
+void BSP::calculateSAH(float parentArea, SplitData &splitdata,
+ const BSH &primitives, const Polytope &polytope)
+{
+ int nLeft = 0;
+ int nRight = 0;
+ int nEither = 0;
+ primitives.countLocations(splitdata.split, nLeft, nRight, nEither);
+
+ //TODO: optimize this to remove all the useless computations.
+
+ const float lProb = splitdata.lArea / parentArea;
+ const float rProb = splitdata.rArea / parentArea;
+
+ const float C_trav = isKDTreeSplit(splitdata.split) ? C_KDTree_trav : 0;
+ //if it's a BSP split, we'll add in the traversal cost later.
+
+ float leftCost, rightCost;
+ leftCost = ((nLeft+nEither)*lProb + nRight*rProb)*C_isec + C_trav;
+ rightCost = (nLeft*lProb + (nRight+nEither)*rProb)*C_isec + C_trav;
+
+ if (leftCost < rightCost) {
+ nLeft+=nEither;
+ splitdata.cost = leftCost;
+ splitdata.coplanar_side = negativeSide;
+ }
+ else {
+ nRight+=nEither;
+ splitdata.cost = rightCost;
+ splitdata.coplanar_side = positiveSide;
+ }
+
+ splitdata.nLeft = nLeft;
+ splitdata.nRight = nRight;
+}
+
+bool BSP::getSurfaceAreas(SplitData &splitdata)
+{
+ Polytope side1, side2;
+ polytope.split(splitdata.split, side1, side2);
+ splitdata.lArea = side1.getSurfaceArea();
+ splitdata.rArea = side2.getSurfaceArea();
+
+ if (splitdata.lArea <= 0 || splitdata.rArea <= 0) {
+ return false;
+ }
+ return true;
+}
+
+
+ BuildSplitPlane BSP::getBestSplitPlane(const BSH &primitives,
+ BuildSplitPlane plane,
+ const vector<Point> &initialVertices)
+ {
+ //some splits, like axis aligned splits, don't do this.
+ if(initialVertices.empty())
+ return plane;
+
+ //static to reduce (de)allocations. Not thread-safe.
+ //static Polygon points;
+ Polygon points; //Can't do static Polygon due to getBestSplitPlane being
called recursively!
+ points = initialVertices;
+
+ if (initialVertices.size() < 3) {
+// primitives.pointsNearPlaneDist(points, plane);
+
+ //Also want the plane to exactly go through the vertices of the
+ //bounding polytope
+ for (size_t i=0; i < polytope.vertices.size(); ++i) {
+
+ Point normal2 = Cross(points[1]-points[0],
+ polytope.vertices[i]-points[0]);
+ if (normal2.length2() < 1e-12)
+ continue;
+ normal2.normalize();
+ double costheta = 1-fabs(Dot(normal2, plane.normal));
+
+ if( costheta < 0.001 ) {
+ points.push_back(polytope.vertices[i]);
+ break;
+ }
+ }
+
+ if (points.size() >= 3)
+ radialSort(points);
+ if (points.size() < 3) {
+ Point closePoint = initialVertices[0];
+ double minCostheta = .2;
+
+ primitives.pointNearPlaneAngle(points, plane, closePoint,
minCostheta);
+ if (closePoint == points[0])
+ return plane;
+ points.push_back(closePoint);
+ if (points.size() >= 3)
+ radialSort(points);
+ if (points.size() < 3)
+ return plane;
+ }
+
+ //Now that we have 3 points, lets try
+ plane.normal = getOrderedFaceNormal(points);
+ const Point center = getCenter(points);
+ plane.d = Dot(plane.normal, -center);
+ }
+
+ primitives.pointsNearPlaneDist(points, plane);
+ const double MAX_DIST = 2.5e-6; //TODO: ideally, make this a
+ //function of distance from vertex
+ //to triangle that defined plane.
+
+ //Also want the plane to exactly go through the vertices of the
+ //bounding polytope
+ for (size_t i=0; i < polytope.vertices.size(); ++i) {
+ double p_dist = signedDistance(plane, polytope.vertices[i]);
+// double errorScalar = distanceAlongPlane(plane, initialVertices,
+// polytope.vertices[i]);
+ if (fabs(p_dist) <= 1*MAX_DIST) {
+// if (p_dist != 0)
+// cout << "got some at: " <<p_dist<<endl;
+ points.push_back(polytope.vertices[i]);
+ }
+ }
+
+ if (points.size() >= 3)
+ radialSort(points);
+
+ if (points.size() < 3) {
+ return plane;
+ }
+
+ BuildSplitPlane optimizedPlane;
+ optimizedPlane.normal = getOrderedFaceNormal(points);
+
+ //we just want a small correction
+ double costheta = fabs(Dot(optimizedPlane.normal, plane.normal));
+ if( (costheta < 0.99 && initialVertices.size()>=3) ||
+ (costheta < 0.5 && initialVertices.size() < 3)) {
+#ifdef BUILD_DEBUG
+ if (costheta > .9) {
+ cout << "corrected too much!\n"
+ <<plane.normal <<endl
+ <<optimizedPlane.normal<<endl;
+ cout<<initialVertices.size()<< " and angle: " <<costheta<<endl;
+ // exit(0);
+ }
+#endif
+ return plane;
+ }
+
+ //Is it better to place the center at the initialVertices or using
+ //all the vertices that were used to calculate the optimized
+ //normal?
+ Point center = getCenter(points);//getCenter(initialVertices);
+ optimizedPlane.d = Dot(optimizedPlane.normal, -center);
+
+ //TODO: check if some points are still far away and if so, remove them
+ //and try again.
+
+ return optimizedPlane;
+ }
+
+
+bool BSP::buildFromFile(const string &file)
+{
+ double startTime = Time::currentSeconds();
+
+ bounds.reset();
+
+ //Something's wrong, let's bail. Since we reset the bounds right
+ //above, this KDTree should never get intersected.
+ if (!mesh)
+ return false;
+
+ PreprocessContext context;
+ static_cast<Group*>(mesh)->computeBounds(context, bounds);
+
+ bool isAscii = false;
+ bool isKDTree = false;
+ if (!strncmp(file.c_str()+file.length()-8, ".kdtreer", 8) ||
+ !strncmp(file.c_str()+file.length()-7, ".kdtree", 7)) {
+ isAscii = isKDTree = true;
+ }
+ if (!strncmp(file.c_str()+file.length()-4, ".bsp", 4)) {
+ isAscii = true;
+ }
+
+ if (isAscii) {
+ objects.clear();
+ nodes.clear();
+ nodes.reserve(mesh->size()*4);
+
+ 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;
+ objects.push_back(item);
+ }
+
+ while (in.good()) {
+ BSPNode node;
+ int i;
+
+ in >> i;
+ int isLeaf = i;
+
+ if (in.eof())
+ break;
+
+ if (isLeaf) {
+ int nObjects;
+ unsigned int children;
+ in >> nObjects;
+ in >> children;
+ node.makeLeaf(children, nObjects);
+ if (!isKDTree) {
+ int triangleBounds = 0;
+ int triPlaneDepth = 0;
+ in >> triangleBounds;
+ in >> triPlaneDepth;
+#ifdef BSP_INTERSECTION
+ node.triangleBounds = triangleBounds;
+ node.setTriPlaneDepth(triPlaneDepth);
+#endif
+ }
+ }
+ else {
+ int type;
+ in >> type;
+
+ float d;
+ in >> d;
+
+ Vector normal(0,0,0);
+ if (type == 0) { //kdtree format
+ d*=-1;
+
+ unsigned int planeDim;
+ in >> planeDim;
+ normal[planeDim] = 1;
+ }
+ else {
+ in >> normal;
+ }
+ unsigned int children;
+ in >> children;
+
+ node.makeInternal(children, normal, d);
+ }
+ nodes.push_back(node);
+ }
+ in.close();
+ }
+ else {
+ ifstream in(file.c_str(), ios::binary | ios::in);
+
+ unsigned int objects_size;
+ in.read((char*)&objects_size, sizeof(objects_size));
+ objects.resize(objects_size);
+
+ in.read((char*) &objects[0], sizeof(objects[0])*objects.size());
+
+ unsigned int nodes_size;
+ in.read((char*)&nodes_size, sizeof(nodes_size));
+ nodes.resize(nodes_size);
+
+ in.read((char*) &nodes[0], sizeof(BSPNode)*nodes.size());
+ in.close();
+ }
+
+ cout << "done building" << endl << flush;
+ float buildTime = Time::currentSeconds() - startTime;
+ printf("BSP tree built in %f seconds\n", buildTime);
+ printStats();
+
+ return true;
+}
+
+
+bool BSP::saveToFile(const string &file)
+{
+#if 0
+ //ascii write (human readable, but big files and inaccurate!)
+ ofstream out(file.c_str());
+ if (!out) return false;
+
+ out << objects.size() << endl;
+ for (unsigned int i=0; i < objects.size(); ++i) {
+ out << objects[i] <<" ";
+ }
+
+ out<<endl;
+
+ for (unsigned int i=0; i < nodes.size(); ++i) {
+ out << nodes[i].isLeaf() << " ";
+ if (nodes[i].isLeaf()) {
+ out << nodes[i].numPrimitives() << " " << nodes[i].objectsIndex << "
"
+ <<nodes[i].triangleBounds << " " <<nodes[i].triPlaneDepth<<endl;
+ }
+ else {
+ if (nodes[i].isKDTree())
+ out << 0 << " " << nodes[i].d << " " << nodes[i].kdtreePlaneDim()
<< " " << nodes[i].children<<endl;
+ else
+ out << 1 << " " << nodes[i].d << " " << nodes[i].normal << " " <<
nodes[i].children<<endl;
+ }
+ }
+#else
+ ofstream out(file.c_str(), ios::out | ios::binary);
+ if (!out) return false;
+
+ const unsigned int objects_size = objects.size();
+ const unsigned int nodes_size = nodes.size();
+ out.write((char*) &objects_size, sizeof(objects_size));
+ out.write((char*) &objects[0], sizeof(objects[0])*objects.size());
+ out.write((char*) &nodes_size, sizeof(nodes_size));
+ out.write((char*) &nodes[0], sizeof(BSPNode)*nodes.size());
+#endif
+ out.close();
+ return true;
+}
+
+
+void BSP::printStats()
+{
+ TreeStats stats = {0};
+ collectTreeStats(0, 0, stats);
+
+ printf("BSP tree made from %d primitives has: \n", (int)mesh->size());
+ printf(" - %d nodes\n", (int)nodes.size());
+ printf(" - %d primitive references\n", (int)objects.size());
+ printf(" - %d max depth\n", stats.maxDepth);
+ printf(" - %d max objects in a leaf\n", stats.maxObjectsInLeaf);
+
+ for (int i=0; i<stats.MAX; ++i)
+ if (stats.childrenLeafHist[i] > 0)
+ printf(" - %d leaves with %d children\n", stats.childrenLeafHist[i],
i);
+
+ float totalKD = 0;
+ float totalBSP = 0;
+ for (int i=0; i<stats.MAX; ++i)
+ if (stats.splitType[i][0] + stats.splitType[i][1] > 0) {
+ totalBSP += stats.splitType[i][0];
+ totalKD += stats.splitType[i][1];
+ float totalSplits = stats.splitType[i][0] + stats.splitType[i][1];
+ printf(" - %d, %d\t (%3.0f%% %3.0f%%) bsp/kdtree splitType\n",
+ stats.splitType[i][0], stats.splitType[i][1],
+ 100*(stats.splitType[i][0] / totalSplits),
+ 100*(stats.splitType[i][1] / totalSplits));
+ }
+
+ printf(" - %3.1f%% nodes kdtree\n", 100*totalKD / (totalKD+totalBSP));
+}
+
+void BSP::collectTreeStats(int nodeID, int depth, TreeStats &stats)
+{
+ if (depth >= stats.MAX) {
+ cout << "Tree is too deep!\n";
+ exit(1);
+ }
+
+ const BSPNode &node = nodes[nodeID];
+ if (node.isLeaf()) {
+ if ((int)node.numPrimitives() >= stats.MAX) {
+ cerr << "There are too many primitives in the leaf! " <<
node.numPrimitives()<<endl;
+ return;
+ }
+
+ stats.maxDepth = Max(stats.maxDepth, depth);
+ stats.maxObjectsInLeaf = Max((unsigned int)stats.maxObjectsInLeaf,
+ node.numPrimitives());
+ stats.childrenLeafHist[node.numPrimitives()]++;
+ }
+ else {
+ if (node.isKDTree())
+ stats.splitType[depth][1]++;
+ else
+ stats.splitType[depth][0]++;
+ collectTreeStats(node.children, depth+1, stats);
+ collectTreeStats(node.children+1, depth+1, stats);
+ }
+}
Added: trunk/Model/Groups/BSP/BSP.h
==============================================================================
--- (empty file)
+++ trunk/Model/Groups/BSP/BSP.h Tue Jun 24 17:27:12 2008
@@ -0,0 +1,292 @@
+#ifndef _MANTA_BSP_H_
+#define _MANTA_BSP_H_
+
+#include <Core/Geometry/BBox.h>
+#include <Core/Geometry/VectorT.h>
+#include <Interface/AccelerationStructure.h>
+#include <Interface/RayPacket.h>
+#include <Model/Groups/Group.h>
+#include <Model/Groups/Mesh.h>
+#include <Model/Groups/BSP/BSH.h>
+#include <Model/Groups/BSP/Polytope.h>
+#include <Model/Groups/BSP/Geometry.h>
+
+#include <assert.h>
+
+//BSP_INTERSECTION improves performance a little bit. It seems to
+//work, but considering how much more complicated it makes things for
+//just a small gain, I personally prefer to leave it off.
+//#define BSP_INTERSECTION
+
+//The BSP acceleration structure has a very slow build (best to save
+//the tree to a file for future use) and only handles triangle meshes.
+//It is however often faster than the other acceleration structures when
+//the scene is "complicated" or if there are non-coherent packets.
+
+
+namespace Manta
+{
+ class BSP : public AccelerationStructure
+ {
+ public:
+ struct BSPNode
+ {
+// //plane is: dot(normal, X) + d == 0
+// //Half-space is dot(normal, X) + d <= 0
+ Vector normal;
+ union {
+ Real d;
+#ifdef BSP_INTERSECTION
+ struct {
+ unsigned int objectsIndex : 30;
+ unsigned int triangleBounds : 2; //0x1 is all 3 edges bounded,
+ //0x2 is face bounded
+ };
+#else
+ unsigned int objectsIndex;
+#endif
+ };
+
+ unsigned int children; //leftChild=children. rightChild=children+1.
+
+ inline void makeLeaf(unsigned int index, int nObjects) {
+ setAsLeaf();
+ setNumPrimitives(nObjects);
+ objectsIndex = index;
+#ifdef BSP_INTERSECTION
+ triangleBounds = 0;
+#endif
+ }
+
+ inline void makeInternal(int first_child, const Vector &normal, Real
d) {
+ children = first_child;
+ this->normal = normal;
+ this->d = d;
+ checkAndMakeKDTree();
+ }
+
+ inline bool isLeaf() const {
+ return children == 0;
+ }
+ inline void setAsLeaf() {
+ children=0;
+ }
+
+ inline unsigned int numPrimitives() const {
+ assert(isLeaf());
+ union {
+ Real f;
+ unsigned int i;
+ };
+ f = normal[0];
+ return i;
+ }
+ inline void setNumPrimitives(unsigned int size) {
+ assert(isLeaf());
+ union {
+ Real f;
+ unsigned int i;
+ };
+ i = size;
+ normal[0] = f;
+ }
+
+ inline void checkAndMakeKDTree() {
+ assert(!isLeaf());
+
+ if (!(normal[0]*normal[1] == 0 &&
+ normal[0]*normal[2] == 0 &&
+ normal[1]*normal[2] == 0))
+ return;
+
+ if (normal.maxComponent() > 0.9999999) {
+ union {
+ Real f;
+ unsigned int i;
+ };
+ i = normal.indexOfMaxComponent();
+ normal[1] = f;
+ normal[0] = 10;
+ d*=-1;
+ }
+ }
+ inline bool isKDTree() const {
+ assert(!isLeaf());
+ return normal[0] == 10;
+ }
+
+ inline int kdtreePlaneDim() const {
+ assert(isKDTree());
+ union {
+ Real f;
+ unsigned int i;
+ };
+ f = normal[1];
+ return i;
+ }
+
+#ifdef BSP_INTERSECTION
+ inline void setTriPlaneDepth(int depth) {
+ union {
+ Real f;
+ int i;
+ };
+ i = depth;
+ normal[2] = f;
+ }
+ inline int getTriPlaneDepth() const {
+ union {
+ Real f;
+ int i;
+ };
+ f = normal[2];
+ return i;
+ }
+#endif
+ };
+
+ vector<BSPNode> nodes;
+ vector<int> objects; //vector of objects added to BSP tree. Might
+ //contain the same object multiple times
+ //(object split by plane)
+
+ //If this BSP is a subset of another acceleration structure, then
+ //if intersect() is called on this bsp we don't need to do a
+ //bounds intersection test since we implicitly know that already
+ //happened in order to get to this node.
+ bool isSubset;
+
+ //getting the triangle vertex data is expensive, especially
+ //converting from float to double. So here we make it easier to
+ //get. Each element of vertices contains 3 Points corresponding
+ //to a triangle.
+ vector<Point*>vertices;
+
+ BBox bounds;
+
+ Mesh *mesh;
+
+ Polytope polytope;
+
+ //for performance measurements. Can be removed.
+ mutable long nTriIntersects;
+ mutable long nBSPTriIntersects;
+ mutable long nTraversals;
+ mutable long nKDTraversals;
+ mutable long nLeafs;
+ mutable long nTotalRays;
+
+ void update(int proc, int numProcs)
+ {
+ //#define COLLECT_STATS
+ #ifdef COLLECT_STATS
+ printf("tri intersections per ray: %f\n",
(float)nTriIntersects/nTotalRays);
+ printf("BSP tri intersections per ray: %f\n",
(float)nBSPTriIntersects/nTotalRays);
+ printf("node traversals per ray: %f\n",
(float)nTraversals/nTotalRays);
+ printf("node bsp traversals per ray: %f\n",
(float)(nTraversals-nKDTraversals)/nTotalRays);
+ printf("node kdtree traversals per ray: %f\n",
(float)(nKDTraversals)/nTotalRays);
+ printf("leaf visits per ray: %f\n",
(float)(nLeafs)/nTotalRays);
+ printf("number of rays: %ld\n", nTotalRays);
+ nTriIntersects = 0;
+ nBSPTriIntersects = 0;
+ nTraversals = 0;
+ nKDTraversals = 0;
+ nLeafs = 0;
+ nTotalRays = 0;
+ #endif
+ }
+
+ BSP() : mesh(NULL), nTriIntersects(0), nBSPTriIntersects(0),
+ nTraversals(0), nKDTraversals(0),
+ nLeafs(0), nTotalRays(0)
+ {}
+
+ void preprocess(const PreprocessContext&);
+ void intersect(const RenderContext& context, RayPacket& rays) const;
+ inline void intersectSingleRays(const RenderContext& context, RayPacket&
rays) const;
+ inline void intersectPacket(const RenderContext& context, RayPacket&
rays) const;
+
+ struct TrvStack {
+ int nodeID;
+ float t;
+ int depth;
+ int farPlane;
+ };
+
+ void traverse(const RenderContext &context, RayPacket &ray,
+ const Vector &org, const Vector &dir, const Vector
&rcp,
+ TrvStack *const stackBase) const;
+
+ void computeBounds(const PreprocessContext&,
+ BBox& bbox) const
+ {
+ if (!nodes.empty())
+ bbox.extendByBox(bounds);
+ }
+
+ void setGroup(Group* new_group);
+ void groupDirty() {} // does nothing for now
+ Group* getGroup() const { return mesh; }
+
+ virtual void addToUpdateGraph(ObjectUpdateGraph* graph,
+ ObjectUpdateGraphNode* parent) { }
+
+ void rebuild(int proc=0, int numProcs=1);
+
+ void rebuildSubset(const vector<int> &subsetObjects, const BBox &bounds);
+
+ void setTriangleBounds(vector<BuildSplitPlane> &splitPlanes, int &nodeID,
+ int &nextFreeNode);
+
+ void build(int nodeID, int &nextFreeNode,
+ vector<BuildSplitPlane> &splitplanes, BSH &bsh, int depth =
0);
+
+ BuildSplitPlane getBuildSplitPlane(BSH &primitives,
+ vector<BuildSplitPlane> &splitPlanes,
+ Location &coplanar_side);
+ BuildSplitPlane getSAHSplit(BSH &primitives,
+ vector<BuildSplitPlane> &splitPlanes,
+ Location &coplanar_primitivesSide);
+
+ struct SplitData {
+ SplitData() : splitPrimitive(-1), splitType(-1) {
+ }
+ BuildSplitPlane split;
+ int splitPrimitive;
+ int splitType; //0 for face, 1, 2, 3 for the three edges, etc...
+ int nLeft, nRight;
+ float lArea, rArea;
+ float cost;
+ Location coplanar_side;
+ };
+
+ //assume last BuildSplitPlane in planes determines left and right SA.
+ bool getSurfaceAreas(SplitData &splitdata);
+
+ void calculateSAH(float parentArea, SplitData &splitdata,
+ const BSH &primitives, const Polytope &polytope);
+
+ BuildSplitPlane getBestSplitPlane(const BSH &primitives,
+ BuildSplitPlane plane,
+ const vector<Point> &initialVertices);
+
+ bool buildFromFile(const string &file);
+ bool saveToFile(const string &file);
+
+ void printPrimitivesAsObj(const BSH &primitives) const;
+ void printOriginalPrimitivesAsObj(const BSH &primitives) const;
+
+ void printStats();
+ struct TreeStats{
+ static const int MAX = 1024;
+ float total_SAH_cost;
+ int maxDepth;
+ int maxObjectsInLeaf;
+ int childrenLeafHist[MAX]; //Better not have more than MAX primitives
in a node!
+ int splitType[MAX][2]; //Tree better not be more than MAX deep!
+ };
+ void collectTreeStats(int nodeID, int depth, TreeStats &stats);
+ };
+};
+
+#endif
Added: trunk/Model/Groups/BSP/Geometry.cc
==============================================================================
--- (empty file)
+++ trunk/Model/Groups/BSP/Geometry.cc Tue Jun 24 17:27:12 2008
@@ -0,0 +1,393 @@
+#include <Model/Groups/BSP/Geometry.h>
+#include <Core/Geometry/AffineTransformT.h>
+#include <Model/Groups/BSP/aip.h>
+#include <set>
+using namespace std;
+
+namespace Manta {
+ namespace BSPs {
+ //http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
+ bool almostEqualRelative(double A, double B,
+ double maxRelativeError,
+ double maxAbsoluteError)
+
+ {
+ if (fabs(A - B) < maxAbsoluteError)
+ return true;
+ double relativeError;
+ if (fabs(B) > fabs(A))
+ relativeError = fabs((A - B) / B);
+ else
+ relativeError = fabs((A - B) / A);
+ if (relativeError <= maxRelativeError)
+ return true;
+ return false;
+ }
+
+ struct PointAngle
+ {
+ Point p;
+ double cos_theta;
+ PointAngle(Point p, double cos_theta):
+ p(p), cos_theta(cos_theta)
+ { }
+ bool operator<(const PointAngle &b) const
+ {
+ //Two points that are equal might have very different
+ //angles. i.e. 0 and 2PI. Because of this two points that are
+ //almost equal might also have very different angles (0.001 and
+ //2PI-.0001), so we might end up with duplicate points.
+ return cos_theta < b.cos_theta;
+ }
+ };
+
+ //if the points in cap form a convex planar polygon, but are out
+ //of order, a radial sort will make it a legit polygon again.
+ void radialSort(Polygon &cap)
+ {
+ //2 or less elements are already radially sorted (if you ignore
+ //the winding order).
+ if (cap.size() < 3)
+ return;
+
+ Point center = getCenter(cap);
+
+ //TODO: Optimize by replacing set with static vector + sort + unique.
+
+ //Rather than store angle, we optimize by using costheta. This
+ //lets us avoid calling expensive acos.
+ set<PointAngle> pointAngles;
+ pointAngles.insert(PointAngle(cap[0], 1));
+
+ const Point base_vec = (cap[0] - center).normal();
+
+ const Point base_direction = getFaceNormal(cap);
+
+ for (size_t i=1; i < cap.size(); ++i) {
+ if (!(cap[i] == cap[0])) {
+ const Point vec = (cap[i] - center).normal();
+ double costheta = Dot(base_vec, vec);
+ //costheta only defined [-1,1]. Precision issues can make it
+ //slightly outside this.
+ if (costheta > 1)
+ costheta = 1;
+ else if (costheta < -1)
+ costheta = -1;
+// double angle = acos(costheta);
+
+ const Point direction = Cross(base_vec, vec);
+ const double sameDirection = Dot(base_direction, direction);
+ if (sameDirection < 0) {
+ costheta = 2-costheta;
+// angle = M_PI*2 - angle;
+ }
+ pointAngles.insert(PointAngle(cap[i], costheta));
+ }
+ }
+
+ cap.clear();
+ int i=0;
+ for (set<PointAngle>::iterator iter = pointAngles.begin();
+ iter != pointAngles.end(); ++iter, ++i) {
+ bool okToAdd = true;
+
+ //ouch, this is expensive. Should check and see whether there's
+ //a faster way.
+ for (size_t j=0; j < cap.size(); ++j) {
+ if (almostEqualPoints(iter->p, cap[j])) {
+ okToAdd = false;
+ break;
+ }
+ }
+ if (okToAdd)
+ cap.push_back(iter->p);
+ }
+ }
+
+ bool verifyUnique(const Polygon &face)
+ {
+ set<Point, ltPoint> face_set(face.begin(), face.end());
+ if (face.size() != face_set.size()) {
+ cout << face.size() << " vs " <<face_set.size()<<endl;
+ exit(1);
+ }
+ return true;
+ }
+
+ void printPolygon(const Polygon &face)
+ {
+ for (size_t i=0; i < face.size(); ++i) {
+ cout <<"v " <<face[i]<<endl;
+ }
+ for (size_t i=0; i < face.size()-2; ++i) {
+ cout << "f " << 1<< " " << i+2 << " " <<i+3 <<endl;
+ }
+ }
+
+ bool isDegenerateFace(const Polygon &face)
+ {
+ //if the three points are each off by 1e-x, then the area will be
+ //of magnitude 1e-2x.
+ double area=0;
+ for (size_t i=2; i < face.size(); ++i) {
+ area += getTriangleArea(face[0], face[i-1], face[i]);
+ }
+ return area < 1e-12; //This would work better if I could compare
+ //with the areas of the other faces. (relative
+ //compare)
+ }
+
+ bool overlaps(const Polygon &face1, const Polygon &face2,
+ float &intersectedArea, bool
forceIntersectionAreaComputation)
+ {
+ Point normal1 = getOrderedFaceNormal(face1);
+ Point normal2 = getOrderedFaceNormal(face2);
+ double normalDot = Dot(normal1, normal2);
+ if (fabs(normalDot) < 0.9999) {
+#ifdef BUILD_DEBUG
+ if (fabs(normalDot) > .9999)
+ cout<<"normals are not parallel: " <<normal1 << " and "
<<normal2<<endl;
+#endif
+ return false;
+ }
+
+ if (normalDot < 0)
+ normal2*=-1;
+
+ Point center1 = getCenter(face1);
+ Point center2 = getCenter(face2);
+
+ if (!forceIntersectionAreaComputation &&
+ almostEqualPoints(center1, center2)) {
+// cout <<"almost same centers: " << center1 << " and "<<
center2<<endl;
+ return true;
+ }
+
+ Point averageNormal = 0.5*(normal1+normal2);
+ double d1 = Dot(averageNormal, -center1);
+ double d2 = Dot(averageNormal, -center2);
+ if (!almostEqualRelative(d1, d2, 5e-6)) {
+#ifdef BUILD_DEBUG
+ //I'm curious if it would have made it...
+ if (almostEqualRelative(d1, d2, 1e-4)) {
+ AffineTransformT<double> transform;
+ transform.initWithRotation(averageNormal, Point(0,0,1));
+
+ vector<Point2D> face1_t, face2_t;
+ for (size_t i=0; i < face1.size(); ++i) {
+ Point v_t = transform.multiply_point(face1[i]);
+ face1_t.push_back(Point2D(v_t[0], v_t[1]));
+ }
+ for (size_t i=0; i < face2.size(); ++i) {
+ Point v_t = transform.multiply_point(face2[i]);
+ face2_t.push_back(Point2D(v_t[0], v_t[1]));
+ }
+
+ if (normalDot < 0)
+ reverse(face2_t.begin(), face2_t.end());
+
+ intersectedArea = PolygonIntersect::intersectionArea(face1_t,
face2_t);
+
+ float area1 = getArea(face1);
+ float area2 = getArea(face2);
+
+ float areaRatio = intersectedArea/Min(area1, area2);
+
+ if (areaRatio > .8 ) {
+ cout <<"normalDot was: " <<normalDot<<endl;
+ cout <<"distances differ: " <<d1 <<" and "<<d2
+ <<" relative: " << ((d1-d2)/d2)<<endl;
+ }
+ }
+#endif
+ return false;
+ }
+
+ //they are pretty much on the same plane. Now lets find out if
+ //they overlap at all.
+ AffineTransformT<double> transform;
+ transform.initWithRotation(averageNormal, Point(0,0,1));
+
+ vector<Point2D> face1_t, face2_t;
+ for (size_t i=0; i < face1.size(); ++i) {
+ Point v_t = transform.multiply_point(face1[i]);
+ face1_t.push_back(Point2D(v_t[0], v_t[1]));
+ }
+ for (size_t i=0; i < face2.size(); ++i) {
+ Point v_t = transform.multiply_point(face2[i]);
+ face2_t.push_back(Point2D(v_t[0], v_t[1]));
+ }
+
+ if (normalDot < 0)
+ reverse(face2_t.begin(), face2_t.end());
+
+ intersectedArea = PolygonIntersect::intersectionArea(face1_t, face2_t);
+
+ float area1 = getArea(face1);
+ float area2 = getArea(face2);
+
+// cout <<"overlap area is: " << area<<endl;
+// cout <<"face 1 area is: " <<area1<<endl;
+// cout <<"face 2 area is: " <<area2<<endl;
+
+ float areaRatio = intersectedArea/Min(area1, area2);
+
+// cout <<"distance: " <<(d1-d2)/d1<<endl;
+#ifdef BUILD_DEBUG
+ if (areaRatio > 1e-2 && areaRatio < .8)
+ cout <<"area ratio: "<<areaRatio << endl;
+#endif
+// exit(1);
+ if (areaRatio > .2 ) return true;
+ else return false;
+
+ }
+
+ bool isSubset(const Polygon &smallFace, const Polygon &bigFace)
+ {
+ return overlaps(smallFace, bigFace);
+ }
+
+ Point getTriangleNormal(const Point &p0, const Point &p1,
+ const Point &p2)
+ {
+ const Point edge1 = p1 - p0;
+ const Point edge2 = p2 - p0;
+ return Cross(edge1, edge2);
+ }
+
+ //This tries to return a normalized "high quality" normal for the
+ //convex face.
+ Point getFaceNormal(const Polygon& face)
+ {
+ //Assume that the vertices in the face are not ordered. The normal
+ //would be defined by the cross product of the two main principal
+ //components (wouldn't that be the third componenent?). But doing
+ //PCA is a lot of work (both for the computer and for the
+ //programmer).
+
+ //find centroid
+ Point center = getCenter(face);
+
+ //If the vertices do not define a convex polygon, then the center
+ //and face[0] might end up lying on the same point (or very close
+ //to it). It is for this reason that we require face to be a
+ //convex poygon. If the vertices are not planar (i.e. we use all
+ //the vertices of the polytope), then base_vec might point in the
+ //direction of the normal, which would obviously prevent us from
+ //evee returning the correct normal.
+
+ //But what we can try is to find the vertex furthest from the
+ //center. This would give a fairly good base_vec.
+ Point base_vec;
+ double biggest=0;
+ for (size_t i=1; i < face.size(); ++i) {
+ const Point vec = (face[i] - center);
+ double length = vec.length2();
+ if (length > biggest) {
+ base_vec = vec;
+ biggest = length;
+ }
+ }
+
+ Point bestNormal;
+ biggest=0;
+ for (size_t i=1; i < face.size(); ++i) {
+ const Point vec = (face[i] - center);
+ Point normal = Cross(base_vec, vec);
+ double length = normal.length2();
+ if (length > biggest) {
+ bestNormal = normal;
+ biggest = length;
+ }
+ }
+
+ //If the face is degenerate, biggest will still be 0.
+ if (biggest == 0)
+ return Point(0,0,0);
+ else
+ return bestNormal / Sqrt(biggest);
+ }
+
+ Point getCenter(const Polygon& vertices)
+ {
+ Point center(0,0,0);
+ for (size_t i=0; i < vertices.size(); ++i)
+ center += vertices[i];
+ center /= vertices.size();
+ return center;
+ }
+
+ Point getOrderedFaceNormal(const Polygon& face)
+ {
+ //adapted from http://tog.acm.org/GraphicsGems/gemsiii/newell.c
+ //This assumes the vertices are in the correct order.
+ if (face.size() < 3)
+ return Point (0,0,0);
+
+ Point u, v;
+ Point normal(0,0,0);
+ for (size_t i=0; i < face.size(); ++i) {
+ u = face[i];
+ v = face[(i+1) % face.size()];
+ normal[0] += (u[1]-v[1]) * (u[2]+v[2]);
+ normal[1] += (u[2]-v[2]) * (u[0]+v[0]);
+ normal[2] += (u[0]-v[0]) * (u[1]+v[1]);
+ }
+ /* normalize the polygon normal to obtain the first
+ three coefficients of the plane equation
+ */
+ return normal.normal();
+ }
+
+ double getArea(const Polygon& face)
+ {
+ //TODO: look at http://jgt.akpeters.com/papers/Sunday02/FastArea.html
+ double area=0;
+ for (size_t i=2; i < face.size(); ++i) {
+ area += getTriangleArea(face[0], face[i-1], face[i]);
+ }
+ return area;
+ }
+
+ double getArea2D(vector<Point2D>& face)
+ {
+ //This comes from sunday's JGT paper on newell normal...
+ // guarantee the first two vertices are also at array end
+ face.push_back(face[0]);
+ face.push_back(face[1]);
+
+ double sum = 0.0;
+ // double *xptr = x+1, *ylow = y, *yhigh = y+2;
+ for (size_t i=1; i <= face.size()-2; i++) {
+ //sum += (*xptr++) * ( (*yhigh++) - (*ylow++) );
+ sum += face[i][0] * (face[i+1][1]-face[i-1][1]);
+ }
+ face.resize(face.size()-2);
+ return 0.5*sum;
+ }
+
+ double getTriangleArea(const Point &p0, const Point &p1,
+ const Point &p2)
+ {
+ return 0.5 * getTriangleNormal(p0, p1, p2).length();
+ }
+
+ void getOrthonormalBasis(const Point &w, Point &u, Point &v)
+ {
+ //From Hughes and Moller: Building an Orthonormal Basis from a
+ //Unit Vector (in JGT).
+ const double x = fabs(w[0]);
+ const double y = fabs(w[1]);
+ const double z = fabs(w[2]);
+ if (x <= y && x <= z)
+ u = Point(0, -w[2], w[1]);
+ else if (y <= x && y <= z)
+ u = Point(-w[2], 0, w[0]);
+ else //if (z <= x && z <= y)
+ u = Point(-w[1], w[0], 0);
+ u = u.normal();
+ v = Cross(w, u); //this is normalized if w is normalized.
+ }
+}
+}
Added: trunk/Model/Groups/BSP/Geometry.h
==============================================================================
--- (empty file)
+++ trunk/Model/Groups/BSP/Geometry.h Tue Jun 24 17:27:12 2008
@@ -0,0 +1,142 @@
+#ifndef Geometry_h
+#define Geometry_h
+
+#include <Core/Geometry/VectorT.h>
+#include <Core/Math/SSEDefs.h>
+
+#include <vector>
+
+#include <iostream>
+using namespace std;
+
+//#define BUILD_DEBUG
+
+typedef Manta::VectorT<double, 2> Point2D;
+
+namespace Manta {
+ namespace BSPs {
+
+ typedef VectorT<double, 3> Point;
+ typedef vector<Point> Polygon; //Assume this is a convex planar polygon.
+
+ class BuildSplitPlane
+ {
+ public:
+ //plane is: dot(normal, X) + d == 0
+ //Half-space is dot(normal, X) + d <= 0
+ Point normal;
+ double d;
+
+ //point to other half-space
+ void flip() {
+ normal = -normal;
+ d = -d;
+ }
+
+ bool operator<(const BuildSplitPlane &p) const
+ {
+ for (int i=0; i < 3; ++i) {
+ if (normal[i] < p.normal[i]) return true;
+ if (normal[i] > p.normal[i]) return false;
+ }
+ return d < p.d;
+ }
+ };
+
+
+ void printPolygon(const Polygon &face);
+
+ bool isDegenerateFace(const Polygon &face);
+
+ bool overlaps(const Polygon &face1, const Polygon &face2,
+ float &intersectedArea,
+ bool forceIntersectionAreaComputation=false);
+ inline bool overlaps(const Polygon &face1, const Polygon &face2) {
+ float intersectedArea=0;
+ return overlaps(face1, face2, intersectedArea);
+ }
+
+ bool isSubset(const Polygon &smallFace, const Polygon &bigFace);
+
+ Point getTriangleNormal(const Point &p0, const Point &p1,
+ const Point &p2);
+
+ //The points in face might not be ordered.
+ //Note that the returned normal could point in either direction.
+ Point getFaceNormal(const Polygon& face);
+
+ //The points in the face are ordered around the polygon.
+ Point getCenter(const Polygon& face);
+ //The points in the face are ordered around the polygon.
+ Point getOrderedFaceNormal(const Polygon& face);
+
+ void getOrthonormalBasis(const Point &v, Point &u, Point &v);
+
+ double getArea(const Polygon& face);
+
+ double getArea2D(vector<Point2D>& face);
+ double getTriangleArea(const Point &p1, const Point &p2,
+ const Point &p3);
+
+ void radialSort(Polygon &cap);
+
+ bool verifyUnique(const Polygon &face);
+
+ bool almostEqualRelative(double A, double B,
+ double maxRelativeError=1e-6,
+ double maxAbsoluteError=1e-7);
+
+ inline bool almostEqualPoints(const Point &l, const Point &r,
+ double maxAbsoluteError=1e-6)
+ {
+ return fabs(l[0]-r[0]) < maxAbsoluteError &&
+ fabs(l[1]-r[1]) < maxAbsoluteError &&
+ fabs(l[2]-r[2]) < maxAbsoluteError;
+ }
+
+ inline bool almostEqualPointsRelative(const Point &l, const Point &r,
+ double maxRelativeError=1e-6, double
maxAbsoluteError=1e-7)
+ {
+ return almostEqualRelative(l[0], r[0], maxRelativeError,
maxAbsoluteError) &&
+ almostEqualRelative(l[1], r[1], maxRelativeError, maxAbsoluteError) &&
+ almostEqualRelative(l[2], r[2], maxRelativeError, maxAbsoluteError);
+ }
+
+
+
+ inline double signedDistance(const BuildSplitPlane &plane, const Point &p)
+ {
+ return Dot(plane.normal, p) + plane.d;
+ }
+
+ double distanceAlongPlane(const Point basis[2], const Point ¢er,
+ const BuildSplitPlane &plane, const Point &p);
+
+
+ inline float inverseSqrt(float f) {
+#ifdef MANTA_SSE
+ sse_t x = _mm_load_ss(&f);
+ x = _mm_rsqrt_ss(x); //we don't need much precision here.
+ float rsqrt;
+ _mm_store_ss( &rsqrt, x );
+ //The compiler should keep the next line in xmm registers
+ return rsqrt*(1.5f - 0.5f*f*rsqrt*rsqrt);
+#else
+ return 1.0f / Sqrt(f);
+#endif
+ }
+
+ struct ltPoint {
+ bool operator()(const Point &a, const Point &b) const
+ {
+ for (int i=0; i < 3; ++i) {
+ if (a[i] < b[i]) return true;
+ if (a[i] > b[i]) return false;
+ }
+ return false; //they are equal
+ }
+ };
+ };
+};
+
+#endif //Geometry_h
Added: trunk/Model/Groups/BSP/Polytope.cc
==============================================================================
--- (empty file)
+++ trunk/Model/Groups/BSP/Polytope.cc Tue Jun 24 17:27:12 2008
@@ -0,0 +1,688 @@
+#include <Core/Geometry/VectorT.h>
+#include <Model/Groups/BSP/Polytope.h>
+#include <Model/Groups/BSP/BSP.h>
+#include <Model/Groups/BSP/aip.h>
+#include <Model/Groups/BSP/chainHull.h>
+#include <algorithm>
+#include <map>
+#include <iostream> //for debugging
+
+namespace Manta {
+
+ using namespace BSPs;
+
+ Polytope::Polytope() : debugPrint(false){
+ }
+
+ Polytope::~Polytope() {
+ }
+
+ void Polytope::printAsObj(int offset) const
+ {
+ cout <<endl;
+ for (size_t f=0; f < faces.size(); ++f) {
+ for (size_t v=0; v < faces[f].size(); ++v)
+ cout << "v " << faces[f][v] << "\n";
+ }
+
+ int v_count = 1;
+ for (size_t f=0; f < faces.size(); ++f) {
+ int v_start = v_count + offset;
+ cout <<endl;
+ for (size_t i=2; i < faces[f].size(); ++i) {
+ cout <<"f " << v_start << " " <<v_start+i-1 << " " <<v_start+i<<endl;
+ }
+ v_count += faces[f].size();
+ }
+
+ cout <<endl;
+ }
+
+ void Polytope::printAsPly()
+ {
+ int numVertex = 0;
+ for (size_t f=0; f < faces.size(); ++f)
+ numVertex += faces[f].size();
+
+ cout <<endl;
+ cout <<"element vertex " << numVertex <<endl;
+ cout <<"element face " << faces.size() <<endl;
+ cout <<"property list uchar int vertex_indices" <<endl;
+ cout <<"end_header\n";
+ for (size_t f=0; f < faces.size(); ++f) {
+ for (size_t v=0; v < faces[f].size(); ++v)
+// cout << "v "
+ cout<< faces[f][v]<<endl;
+ }
+
+ int v_count = 0;
+ for (size_t f=0; f < faces.size(); ++f) {
+ cout << "\n"<< faces[f].size()<< " ";
+ for (size_t v=0; v < faces[f].size(); ++v)
+ cout << v_count++ << " ";
+ }
+ cout <<endl;
+ }
+
+
+ void Polytope::clear()
+ {
+ faces.clear();
+ }
+
+ void Polytope::initialize(const BBox &bounds)
+ {
+ clear();
+ faces.resize(6);
+
+ faces[0].push_back(Point(bounds[0][0], bounds[0][1], bounds[0][2]));
+ faces[0].push_back(Point(bounds[0][0], bounds[0][1], bounds[1][2]));
+ faces[0].push_back(Point(bounds[0][0], bounds[1][1], bounds[1][2]));
+ faces[0].push_back(Point(bounds[0][0], bounds[1][1], bounds[0][2]));
+
+ faces[1].push_back(Point(bounds[0][0], bounds[0][1], bounds[0][2]));
+ faces[1].push_back(Point(bounds[0][0], bounds[1][1], bounds[0][2]));
+ faces[1].push_back(Point(bounds[1][0], bounds[1][1], bounds[0][2]));
+ faces[1].push_back(Point(bounds[1][0], bounds[0][1], bounds[0][2]));
+
+ faces[2].push_back(Point(bounds[0][0], bounds[0][1], bounds[0][2]));
+ faces[2].push_back(Point(bounds[0][0], bounds[0][1], bounds[1][2]));
+ faces[2].push_back(Point(bounds[1][0], bounds[0][1], bounds[1][2]));
+ faces[2].push_back(Point(bounds[1][0], bounds[0][1], bounds[0][2]));
+
+ faces[3].push_back(Point(bounds[1][0], bounds[0][1], bounds[0][2]));
+ faces[3].push_back(Point(bounds[1][0], bounds[0][1], bounds[1][2]));
+ faces[3].push_back(Point(bounds[1][0], bounds[1][1], bounds[1][2]));
+ faces[3].push_back(Point(bounds[1][0], bounds[1][1], bounds[0][2]));
+
+ faces[4].push_back(Point(bounds[0][0], bounds[0][1], bounds[1][2]));
+ faces[4].push_back(Point(bounds[0][0], bounds[1][1], bounds[1][2]));
+ faces[4].push_back(Point(bounds[1][0], bounds[1][1], bounds[1][2]));
+ faces[4].push_back(Point(bounds[1][0], bounds[0][1], bounds[1][2]));
+
+ faces[5].push_back(Point(bounds[0][0], bounds[1][1], bounds[0][2]));
+ faces[5].push_back(Point(bounds[0][0], bounds[1][1], bounds[1][2]));
+ faces[5].push_back(Point(bounds[1][0], bounds[1][1], bounds[1][2]));
+ faces[5].push_back(Point(bounds[1][0], bounds[1][1], bounds[0][2]));
+
+ //In case we have a flat BBox, we need to remove degenerate faces.
+ cleanup();
+
+ }
+
+ bool Polytope::verifyWindingOrder() const
+ {
+ for (size_t f=0; f < faces.size(); ++f) {
+ const Point v1 = getOrderedFaceNormal(faces[f]);
+ for (size_t i=2; i < faces[f].size(); ++i) {
+ const Point v2 = getTriangleNormal(faces[f][0], faces[f][i-1],
faces[f][i]).normal();
+ if (Dot(v1, v2) < .98) { // really it should be <1
+// if (getTriangleArea(faces[f][0], faces[f][i-1], faces[f][i]) <
1e-7)
+// continue;
+ cout << "winding order is jacked!"<<endl;
+ cout <<"face: "<<f<<endl;
+ cout << v1 << "\tvs\t"<<v2<<endl;
+ cout <<Dot(v1, v2)<<endl;
+
+ for (size_t i=2; i < faces[f].size(); ++i) {
+ cout <<" normal: "<< getTriangleNormal(faces[f][0],
faces[f][i-1], faces[f][i]) <<endl;
+ }
+
+ printAsObj();
+// exit(1);
+ return false;
+ }
+ }
+ }
+ return true;
+ }
+
+ bool Polytope::hasSimilarFace(const Polytope &poly, const Polygon &cap)
const
+ {
+ //We assume that if cap contains a face of the poly, then
+ //they are the same. This handles situations where the cap contains extra
+ //(redundant) vertices that came from different faces.
+ //It should not be possible to have a cap with less vertices and
+ //still be equal because those extra points would have been added
+ //to the cap.
+
+ for (size_t face=0; face < poly.faces.size(); ++face) {
+#ifdef BUILD_DEBUG
+ const size_t polySize = poly.faces[face].size();
+ if (polySize <= cap.size()) {
+ if (isSubset(poly.faces[face], cap)) {
+
+ if (polySize < cap.size()) {
+ cout <<"there were less.\n";
+// for (size_t i =0; i < poly.faces[face].size(); ++i)
+// cout <<poly.faces[face][i]<<endl;
+// cout <<endl;
+// for (size_t i =0; i < cap.size(); ++i)
+// cout <<cap[i]<<endl;
+// poly.printAsObj();
+// // exit(1);
+// debugPrint = true;
+ }
+
+// return true;
+ }
+ }
+ else
+ if (isSubset(cap, poly.faces[face])) {
+ if (polySize > cap.size()) {
+ cout <<"there were more\n";
+ // exit(1);
+ }
+// return true;
+ }
+#endif
+ if (overlaps(poly.faces[face], cap))
+ return true;
+ }
+ return false;
+ }
+
+ bool Polytope::findSplitEdge(PointDistMap &pointDists, const Polygon
&face) const
+ {
+ const double MIN_DIST = 2e-6;
+
+ int v1 = -1,
+ v2 = -1;
+
+ //TODO: This can probably be heavily optimized. Should do that.
+ for (size_t i1=0; i1 < face.size()-1; ++i1) {
+ PointDist &i1distOnPlane = pointDists.find(face[i1])->second;
+ if (fabs(i1distOnPlane.dist) >= MIN_DIST || i1distOnPlane.offPlane)
+ continue;
+ for (size_t i2=i1+1; i2 < face.size(); ++i2) {
+ PointDist &i2distOnPlane = pointDists.find(face[i2])->second;
+ if (fabs(i2distOnPlane.dist) >= MIN_DIST || i2distOnPlane.offPlane)
+ continue;
+
+ //We now have two vertices i1 and i2 that are supposed to be
+ //on the plane. Now lets find out if the edge formed by i1 and
+ //i2 cleanly divides the set of vertices into those that are
+ //above and those that are below the plane.
+ //TODO: How do we handle dist==0? Should we handle that
+ //specifically?
+ double insideSign=0;
+ size_t i;
+ for (i=i1+1; i < i2; ++i) {
+ const PointDist &distOnPlane = pointDists.find(face[i])->second;
+ if (insideSign == 0)
+ insideSign = distOnPlane.dist;
+ else if (insideSign*distOnPlane.dist < 0)
+ break; //signs disagree
+ }
+ if (i != i2) {//did they all agree?
+ continue; //no.
+ }
+
+ insideSign = 0;
+ for (i=(i2+1)%face.size(); i!=i1; i=(i+1)%face.size()) {
+ const PointDist &distOnPlane = pointDists.find(face[i])->second;
+ if (insideSign == 0)
+ insideSign = distOnPlane.dist;
+ else if (insideSign*distOnPlane.dist < 0)
+ break; //signs disagree
+ }
+ if (i != i1) { //did they all agree?
+ continue; //no.
+ }
+
+ if (v1<0 ||
+ (fabs(pointDists.find(face[i1])->second.dist) +
fabs(pointDists.find(face[i2])->second.dist))
+ <(fabs(pointDists.find(face[v1])->second.dist) +
fabs(pointDists.find(face[v2])->second.dist))) {
+ //We found our edges!
+ v1 = i1;
+ v2 = i2;
+ }
+ }
+ }
+ if (v1 >= 0) {
+ //cout <<"using "<<v1<<" and "<<v2<<endl;
+ //i1distOnPlane.dist = 0;
+ //i2distOnPlane.dist = 0;
+ //Set the other vertices as being off the plane.
+ for (int i=0; i < static_cast<int>(face.size()); ++i) {
+ if (i != v1 && i != v2) {
+ PointDist &distOnPlane = pointDists.find(face[i])->second;
+ distOnPlane.offPlane = true;
+ }
+ }
+ return true;
+ }
+
+ //Found no edge.
+ v1=-1;
+ for (size_t i=0; i < face.size(); ++i) {
+ PointDist &distOnPlane = pointDists.find(face[i])->second;
+ if (v1 < 0 && fabs(distOnPlane.dist) < MIN_DIST &&
!distOnPlane.offPlane) {
+ distOnPlane.dist=0;
+ v1 = i;
+ }
+ else
+// perhaps, it should be all the points are ON the plane...?
+ distOnPlane.offPlane = true;
+ }
+ return false;
+ }
+
+ void Polytope::findPointsOnPlane(const BuildSplitPlane& plane,
+ PointDistMap &pointDists) const
+ {
+ //TODO: change from using stl::map to something based on vector.
+
+ bool debug = false;
+
+ const double MIN_DIST = 2e-6;
+ for (size_t f=0; f < faces.size(); ++f) {
+ unsigned int pointsOnPlane = 0;
+ unsigned int pointsAlmostOnPlane = 0;
+ for (size_t i=0; i < faces[f].size(); ++i) {
+ const double dist = signedDistance(plane, faces[f][i]);
+
+ const PointDist &distOnPlane =
+ pointDists.insert(make_pair(faces[f][i],
+ PointDist(dist, false))).first->second;
+ if (fabs(distOnPlane.dist) < MIN_DIST && !distOnPlane.offPlane)
+ ++pointsOnPlane;
+ //TODO: VERY IMPORTANT: include plane normal so that big
+ //planes are treated fairly (otherwise, the 3*MIN_DIST thing
+ //is too restrictive).
+ else if (fabs(distOnPlane.dist) < 3*MIN_DIST &&
!distOnPlane.offPlane)
+ ++pointsAlmostOnPlane;
+ }
+ if (pointsOnPlane > 2 &&
+ pointsOnPlane + pointsAlmostOnPlane < faces[f].size()) {
+
+#ifdef BUILD_DEBUG
+ cout <<"looks like we need to minimize face "<< f <<" ! "
+ <<pointsOnPlane << " + "<<pointsAlmostOnPlane<<" < "
+ <<faces[f].size() <<"\n";
+
+ for (size_t i=0; i < faces[f].size(); ++i) {
+ cout << " i: " <<signedDistance(plane, faces[f][i]) <<endl;
+ }
+#endif
+ debug = true;
+
+ findSplitEdge(pointDists, faces[f]);
+ }
+ }
+
+ //Now we need to figure out which points can be placed exactly on the
plane.
+ for (size_t f=0; f < faces.size(); ++f) {
+ unsigned int pointsOnPlane = 0;
+ unsigned int pointsAlmostOnPlane = 0;
+ for (size_t i=0; i < faces[f].size(); ++i) {
+ const PointDist &distOnPlane = pointDists.find(faces[f][i])->second;
+ if (fabs(distOnPlane.dist) < MIN_DIST && !distOnPlane.offPlane)
+ ++pointsOnPlane;
+ else if (fabs(distOnPlane.dist) < 3*MIN_DIST &&
!distOnPlane.offPlane)
+ ++pointsAlmostOnPlane;
+ }
+
+#ifdef BUILD_DEBUG
+ if (debug) {
+ cout <<"face " << f << " has: " <<pointsOnPlane<< " + "
+ <<pointsAlmostOnPlane<<" < "<<faces[f].size()<<endl;
+ }
+#endif
+ //all on plane
+ if (pointsOnPlane >= 3 &&
+ pointsOnPlane+pointsAlmostOnPlane == faces[f].size())
+ for (size_t i=0; i < faces[f].size(); ++i) {
+// cout <<"1: " <<f<<" "
<<i<<"\t"<<pointDists.find(faces[f][i])->second.distUlps
+// <<"and "<<pointDists.find(faces[f][i])->second.dist<<endl;
+ pointDists.find(faces[f][i])->second.dist = 0;
+ }
+ //an edge or just a corner
+ else if (pointsOnPlane > 0 && pointsOnPlane < 3)
+ for (size_t i=0; i < faces[f].size(); ++i) {
+ PointDist &distOnPlane = pointDists.find(faces[f][i])->second;
+// cout <<"2: " <<f<<" "
<<i<<"\t"<<pointDists.find(faces[f][i])->second.distUlps
+// <<"and "<<pointDists.find(faces[f][i])->second.dist<<endl;
+ if (fabs(distOnPlane.dist) < MIN_DIST && !distOnPlane.offPlane)
+ distOnPlane.dist = 0;
+ }
+ //partial face on plane (bad!)
+ else if (pointsOnPlane > 2) {
+#ifdef BUILD_DEBUG
+ cout.precision(20);
+ cout <<f<<" oops, partial face on plane! " <<pointsOnPlane<< " + "
+ <<pointsAlmostOnPlane<<" < "<<faces[f].size()<<endl;
+
+ cout <<f<<"out of: " <<faces.size()<<endl;
+ for (size_t i=0; i < faces[f].size(); ++i) {
+// cout << i<<": " <<signedUlpsDistance(plane, faces[f][i]);
+// <<" vs "
<<pointDists.find(faces[f][i])->second.distUlps<<"\t";
+ cout << i<<": " <<signedDistance(plane, faces[f][i])
+ <<" vs " <<pointDists.find(faces[f][i])->second.dist
+ <<"offPlane:
"<<pointDists.find(faces[f][i])->second.offPlane<<endl;
+ }
+ printAsObj();
+ cout << "plane " << plane.normal<< " " <<plane.d<<endl;
+#endif
+
+ //exit(1);
+ //debugPrint = true;//exit(1);
+ }
+ else {
+// cout <<pointsOnPlane<<" and " <<pointsAlmostOnPlane<<endl;
+// for (size_t i=0; i < faces[f].size(); ++i)
+// cout <<"3: " <<f<<" "
<<i<<"\t"<<pointDists.find(faces[f][i])->second.distUlps
+// <<"and "<<pointDists.find(faces[f][i])->second.dist<<endl;
+ }
+ }
+ }
+
+ void Polytope::split(const BuildSplitPlane &plane, Polytope &side1,
Polytope &side2) const
+ {
+ PointDistMap pointDists;
+ findPointsOnPlane(plane, pointDists);
+
+ //We are making these Polygons (really stl::vectors) static so
+ //that the vector doesn't allocate/deallocate memory all the
+ //time. Of course this would need to be changed into a some sort
+ //of pool if we end up threading this code.
+ // static
+ Polygon cap;
+ cap.clear();
+
+ for (size_t f=0; f < faces.size(); ++f) {
+#ifdef BUILD_DEBUG
+ verifyUnique(faces[f]);
+#endif
+ static Polygon face1, face2;
+ face1.clear();
+ face2.clear();
+ splitPolygon(pointDists, plane, faces[f], face1, face2, cap);
+
+ if (face1.size() > 2)
+ side1.faces.push_back(face1);
+ if (face2.size() > 2) {
+ side2.faces.push_back(face2);
+ }
+ }
+
+ //radially sort points in cap.
+
+ //It's possible to get no edge. Imagine a split plane that only
+ //touches an edge or a vertex of the polytope.
+ if (true) {
+
+ radialSort(cap);
+ //A radial sort removes duplicate points using several
+ //comparison methods. Since we then need to compare the cap with
+ //these faces to see if they are similar, we perform a radial
+ //sort on these faces as well so that if they are similar, they
+ //will show it.
+
+
+ if (isDegenerateFace(cap))
+ cap.clear();
+
+// side1.verifyWindingOrder();
+ Polytope side1Copy = side1;
+ for (size_t f=0; f < side1.faces.size(); ++f) {
+ radialSort(side1.faces[f]);
+ }
+#ifdef BUILD_DEBUG
+ if (!side1.verifyWindingOrder()) {
+ printAsObj();
+
+ for (size_t f=0; f < side2.faces.size(); ++f) {
+ radialSort(side2.faces[f]);
+ }
+ side2.printAsObj();
+ exit(1);
+ }
+#endif
+ for (size_t f=0; f < side2.faces.size(); ++f) {
+ radialSort(side2.faces[f]);
+ }
+
+ if (cap.size() > 2) {
+
+ side1.cleanup();
+ side2.cleanup();
+// we want to remove the face which is larger or smaller?
+ if (!hasSimilarFace(side1, cap))
+ side1.faces.push_back(cap);
+ else {
+// cout <<"cap1:\n";
+// printPolygon(cap);
+// for (size_t i=0; i < cap.size(); ++i)
+// cout <<i<<" has dist of:" <<signedUlpsDistance(plane,
cap[i])<<" \t"<< signedDistance(plane, cap[i])<<endl;
+// for (size_t f=0; f < faces.size(); ++f)
+// for (size_t i=0; i < faces[f].size(); ++i)
+// cout <<f << " " <<i<<" has dist of:"
<<signedUlpsDistance(plane, faces[f][i])<<" \t" <<signedDistance(plane,
faces[f][i])<<endl;
+
+ }
+ if (!hasSimilarFace(side2, cap))
+ side2.faces.push_back(cap);
+ else {
+// cout <<"cap2:\n";
+// printPolygon(cap);
+ }
+ }
+ else {
+// size_t oldsize = side1.faces.size();
+ side1.cleanup();
+// if (oldsize != side1.faces.size()) {
+// cout <<oldsize << " " <<side1.faces.size()<<endl;
+// exit(1);
+// }
+// oldsize = side2.faces.size();
+ side2.cleanup();
+// if (oldsize != side2.faces.size()) {
+// cout <<oldsize << " " <<side2.faces.size()<<endl;
+// exit(1);
+// }
+ }
+ }
+
+#ifdef BUILD_DEBUG
+ side1.verifyWindingOrder();
+ side2.verifyWindingOrder();
+#endif
+ }
+
+
+ void Polytope::splitPolygon(const PointDistMap &pointDists,
+ const BuildSplitPlane &plane, const Polygon
&in,
+ Polygon &side1, Polygon &side2, Polygon &cap)
const
+ {
+ //This is the Reentrant Polygon Clipping algorithm by Sutherland
+ //and Hodgman.
+
+ Point s, f;
+ for (size_t i=0; i < in.size(); ++i) {
+ const Point &p = in[i];
+
+ const PointDist &distOnPlane = pointDists.find(p)->second;
+ double p_d = distOnPlane.dist;
+ if (i==0) { //first point?
+ s = f = p;
+ }
+ else {
+ const PointDist &distOnPlane = pointDists.find(s)->second;
+ double s_d = distOnPlane.dist;
+
+ //if s_d and p_d are supposed to be 0, but aren't, then we
+ //might end up seeing that as an intersection if the signs
+ //differ. This will result in the cap having two points which
+ //are very similar.
+
+ if (s_d*p_d < 0) { //does line sp cross dividing plane?
+ //compute intersection I, of sp and plane.
+ const double length1 = fabs(p_d);
+ const double length2 = fabs(s_d);
+
+ const double a = length1 / (length1 + length2);
+
+ Point I;
+ //This version will give more consistent results, but I
+ //don't know if it buys us anything really useful.
+ if (a > .5)
+ I = p + a*(s-p);
+ else
+ I = s + (1-a)*(p-s);
+
+ side1.push_back(I);
+ side2.push_back(I);
+ cap.push_back(I);
+ }
+ }
+
+ //How does p relate to plane?
+ if (p_d == 0) { //on the plane
+ side1.push_back(p);
+ side2.push_back(p);
+
+ //Note: if the polygon is on the plane then we will end up
+ //with side1, side2, and cap all being copies of the polygon.
+ //Since cap is then added to side1 and side2, both sides will
+ //now have two copies of the polygon and double the area. If
+ //this repeats, the number of polygons and of course the
+ //areas, will double each time, which can be bad!.
+ cap.push_back(p);
+ }
+ else if (p_d < 0)
+ side1.push_back(p);
+ else if (p_d > 0) {
+ side2.push_back(p);
+ }
+ s=p;
+ }
+
+ //now we need to close the two polygons.
+ //TODO: once we know it works, do this closing as part of the for
+ //loop.
+ const Point &p = f;
+
+ const PointDist &distPOnPlane = pointDists.find(p)->second;
+ double p_d = distPOnPlane.dist;
+ const PointDist &distSOnPlane = pointDists.find(s)->second;
+ double s_d = distSOnPlane.dist;
+
+ if (s_d*p_d < 0) { //does line sp cross dividing plane?
+ //compute intersection I, of sp and plane.
+ const double length1 = fabs(p_d);
+ const double length2 = fabs(s_d);
+
+ const double a = length1 / (length1 + length2);
+ Point I;
+
+ //This version will give more consistent results, but I
+ //don't know if it buys us anything really useful.
+ if (a > .5)
+ I = p + a*(s-p);
+ else
+ I = s + (1-a)*(p-s);
+
+ side1.push_back(I);
+ side2.push_back(I);
+ cap.push_back(I);
+ }
+ }
+
+ void Polytope::cleanup()
+ {
+ for (vector<Polygon>::iterator iter=faces.begin(); iter!=faces.end();
++iter) {
+// cout <<"testing face\n";
+ if (isDegenerateFace(*iter)) {
+ iter = faces.erase(iter);
+ if (faces.empty())
+ break;
+ if (iter != faces.begin())
+ --iter;
+ }
+
+ //TODO: remove colinear vertices.
+
+ //if dot product of two adjacent edges is >0.9999, then remove
+ //middle vertex.
+
+ }
+ }
+
+ bool Polytope::isFlat() const {
+ //Determine whether all the vertices lie on the plane determined
+ //by one of the faces.
+
+// return faces.size() < 4;
+
+ //This is a very simple safe test (assuming nothing bad happened!)
+ if (faces.size() < 1)
+ return true;
+
+ //We are making this Polygon (really stl::vectors) static so
+ //that the vector doesn't allocate/deallocate memory all the
+ //time. Of course this would need to be changed into a some sort
+ //of pool if we end up threading this code.
+ static Polygon vertices;
+ vertices.clear();
+ for (size_t f=0; f < faces.size(); ++f) {
+ vertices.insert(vertices.end(), faces[f].begin(),faces[f].end());
+ }
+
+ Point normal = getFaceNormal(vertices);
+
+ Point u, v;
+ getOrthonormalBasis(normal, u, v);
+
+ static vector<Point2D> vertices_t;
+ vertices_t.clear();
+ for (size_t i=0; i < vertices.size(); ++i) {
+ Point2D v_t(Dot(vertices[i], u), Dot(vertices[i], v));
+ vertices_t.push_back(v_t);
+ }
+
+ static vector<Point2D> hull;
+ convexHull(vertices_t, hull);
+
+ const float flattenedArea = getArea2D(hull);
+
+ double area = 0;
+ for (size_t f=0; f < faces.size(); ++f) {
+ for (size_t i=2; i < faces[f].size(); ++i) {
+ area += getTriangleArea(faces[f][0], faces[f][i-1], faces[f][i]);
+ }
+ }
+
+// cout <<flattenedArea <<" > " <<0.5*area<<endl;
+
+ return flattenedArea >= 0.4999*area;
+ }
+
+
+ float Polytope::getSurfaceArea()
+ {
+ //TODO: if this function becomes expensive, store the area in each
+ //face and then only recompute when the face has been broken up.
+ double area = 0;
+ for (size_t f=0; f < faces.size(); ++f) {
+ for (size_t i=2; i < faces[f].size(); ++i) {
+ area += getTriangleArea(faces[f][0], faces[f][i-1], faces[f][i]);
+ }
+ }
+
+ if (isFlat())
+ return area*2;
+ else
+ return area;
+ }
+
+ void Polytope::updateVertices()
+ {
+ vertices.clear();
+ set<Point, ltPoint> vertices_set;
+ for (size_t f=0; f < faces.size(); ++f)
+ for (size_t i=0; i < faces[f].size(); ++i)
+ vertices_set.insert(faces[f][i]);
+ vertices.insert(vertices.begin(), vertices_set.begin(),
vertices_set.end());
+ }
+
+};
Added: trunk/Model/Groups/BSP/Polytope.h
==============================================================================
--- (empty file)
+++ trunk/Model/Groups/BSP/Polytope.h Tue Jun 24 17:27:12 2008
@@ -0,0 +1,57 @@
+#ifndef Polytope_h
+#define Polytope_h
+
+#include <Model/Groups/BSP/Geometry.h>
+#include <Core/Geometry/BBox.h>
+#include <vector>
+#include <set>
+#include <iostream>
+using namespace std;
+
+namespace Manta {
+ using namespace BSPs;
+
+ class Polytope {
+ public:
+ Polytope();
+ ~Polytope();
+
+ mutable bool debugPrint;
+
+ float getSurfaceArea();
+
+ void printAsObj(int offset=0) const;
+ void printAsPly();
+ void clear();
+ void initialize(const BBox &bounds);
+ bool isFlat() const;
+
+ void split(const BuildSplitPlane &plane, Polytope &side1, Polytope
&side2) const;
+
+ void updateVertices();
+// private:
+ vector<Polygon> faces;
+ vector<Point> vertices;
+
+ struct PointDist {
+ double dist;
+ bool offPlane;
+ PointDist(double dist, bool offPlane) : dist(dist), offPlane(offPlane)
{ }
+ };
+ typedef map<Point, PointDist, ltPoint> PointDistMap;
+
+ bool findSplitEdge(PointDistMap &pointDists, const Polygon &face) const;
+ void findPointsOnPlane(const BuildSplitPlane& plane,
+ PointDistMap &pointDists) const;
+
+ void splitPolygon(const PointDistMap &pointDists,
+ const BuildSplitPlane &plane, const Polygon &in,
+ Polygon &side1, Polygon &side2, Polygon &cap) const;
+
+ bool verifyWindingOrder() const;
+ bool hasSimilarFace(const Polytope &poly, const Polygon &cap) const;
+ void cleanup();
+ };
+};
+
+#endif
Added: trunk/Model/Groups/BSP/aip.h
==============================================================================
--- (empty file)
+++ trunk/Model/Groups/BSP/aip.h Tue Jun 24 17:27:12 2008
@@ -0,0 +1,216 @@
+#ifndef aip_h
+#define aip_h
+
+#include <stdio.h>
+#define bigReal 1.E38 // FLT_MAX, DBL_MAX if double above
+#include <limits.h>
+#include <Core/Geometry/VectorT.h>
+
+using namespace Manta;
+
+/**
+ * Area of Intersection of Polygons
+ *
+ * Algorithm based on http://cap-lore.com/MathPhys/IP/
+ *
+ * Adapted 9-May-2006 by Lagado to java.
+ */
+
+typedef VectorT<double, 2> Point2D;
+typedef long long hp;
+
+class PolygonIntersect
+{
+public:
+ /**
+ * return the area of intersection of two polygons
+ *
+ * Note: the area result has little more accuracy than a float
+ * This is true even if the polygon is specified with doubles.
+ */
+ static double intersectionArea(const vector<Point2D> &a,
+ const vector<Point2D> &b)
+ {
+ PolygonIntersect polygonIntersect;
+ return polygonIntersect.inter(a, b);
+ }
+
+
//--------------------------------------------------------------------------
+
+ struct Point {
+ double x; double y;
+ Point() { };
+ Point(double x, double y) { this->x = x; this->y = y; }
+ };
+ struct Box {
+ Point min; Point max;
+ Box(Point min, Point max) { this->min = min; this->max = max; }
+ };
+ struct Rng {
+ int mn; int mx;
+ Rng() { }
+ Rng(int mn, int mx) { this->mn = mn; this->mx = mx; }
+ };
+ struct IPoint { int x; int y; };
+ struct Vertex { IPoint ip; Rng rx; Rng ry; int in; };
+
+// static const double gamut = 500000000.;
+// static const double mid = 500000000 / 2.;
+#define gamut 500000000
+#define mid 250000000
+
+
//--------------------------------------------------------------------------
+
+private:
+ static void range(const vector<Point2D> &points, int c, Box &bbox)
+ {
+ while (c-- > 0) {
+ bbox.min.x = Min(bbox.min.x, points[c][0]);
+ bbox.min.y = Min(bbox.min.y, points[c][1]);
+ bbox.max.x = Max(bbox.max.x, points[c][0]);
+ bbox.max.y = Max(bbox.max.y, points[c][1]);
+ }
+ }
+
+ static hp area(IPoint a, IPoint p, IPoint q) {
+ return (hp)p.x * q.y - (hp)p.y * q.x +
+ (hp)a.x * (p.y - q.y) + (hp)a.y * (q.x - p.x);
+ }
+
+ static bool ovl(Rng p, Rng q) {
+ return p.mn < q.mx && q.mn < p.mx;
+ }
+
+
//--------------------------------------------------------------------------
+
+ hp ssss;
+ double sclx;
+ double scly;
+
+ PolygonIntersect() :ssss(0), sclx(0), scly(0){}
+
+ void cntrib(int f_x, int f_y, int t_x, int t_y, int w) {
+ ssss += (hp)w * (t_x - f_x) * (t_y + f_y) / 2;
+ }
+
+ void fit(const vector<Point2D> &x, int cx, Vertex ix[], int fudge, Box B)
+ {
+ int c = cx;
+ while (c-- > 0) {
+// ix[c] = new Vertex();
+// ix[c].ip = new IPoint();
+ ix[c].ip.x = ((int)((x[c][0] - B.min.x) * sclx - mid) & ~7)
+ | fudge | (c & 1);
+ ix[c].ip.y = ((int)((x[c][1] - B.min.y) * scly - mid) & ~7)
+ | fudge;
+ }
+
+ ix[0].ip.y += cx & 1;
+ ix[cx] = ix[0];
+
+ c = cx;
+ while (c-- > 0) {
+ ix[c].rx = ix[c].ip.x < ix[c + 1].ip.x ?
+ Rng(ix[c].ip.x, ix[c + 1].ip.x) :
+ Rng(ix[c + 1].ip.x, ix[c].ip.x);
+ ix[c].ry = ix[c].ip.y < ix[c + 1].ip.y ?
+ Rng(ix[c].ip.y, ix[c + 1].ip.y) :
+ Rng(ix[c + 1].ip.y, ix[c].ip.y);
+ ix[c].in = 0;
+ }
+ }
+
+ void cross(Vertex &a, Vertex &b, Vertex &c, Vertex &d,
+ double a1, double a2, double a3, double a4)
+ {
+ double r1 = a1 / ((double) a1 + a2);
+ double r2 = a3 / ((double) a3 + a4);
+
+ cntrib((int)(a.ip.x + r1 * (b.ip.x - a.ip.x)),
+ (int)(a.ip.y + r1 * (b.ip.y - a.ip.y)),
+ b.ip.x, b.ip.y, 1);
+ cntrib(d.ip.x, d.ip.y,
+ (int)(c.ip.x + r2 * (d.ip.x - c.ip.x)),
+ (int)(c.ip.y + r2 * (d.ip.y - c.ip.y)),
+ 1);
+ ++a.in;
+ --c.in;
+ }
+
+ void inness(Vertex P[], int cP, Vertex Q[], int cQ)
+ {
+ int s = 0;
+ int c = cQ;
+ IPoint p = P[0].ip;
+
+ while (c-- > 0) {
+ if (Q[c].rx.mn < p.x && p.x < Q[c].rx.mx) {
+ bool sgn = 0 < area(p, Q[c].ip, Q[c + 1].ip);
+ s += (sgn != Q[c].ip.x < Q[c + 1].ip.x) ? 0 : (sgn ? -1 : 1);
+ }
+ }
+ for (int j = 0; j < cP; ++j) {
+ if (s != 0)
+ cntrib(P[j].ip.x, P[j].ip.y,
+ P[j + 1].ip.x, P[j + 1].ip.y, s);
+ s += P[j].in;
+ }
+ }
+
+
//-------------------------------------------------------------------------
+
+ double inter(const vector<Point2D> &a, const vector<Point2D> &b)
+ {
+ int na = a.size();
+ int nb = b.size();
+ Vertex ipa[na + 1];
+ Vertex ipb[nb + 1];
+ Box bbox = Box(Point(bigReal, bigReal),
+ Point(-bigReal, -bigReal));
+
+ if (na < 3 || nb < 3)
+ return 0;
+
+ range(a, na, bbox);
+ range(b, nb, bbox);
+
+ double rngx = bbox.max.x - bbox.min.x;
+ sclx = gamut / rngx;
+ double rngy = bbox.max.y - bbox.min.y;
+ scly = gamut / rngy;
+ double ascale = sclx * scly;
+
+ fit(a, na, ipa, 0, bbox);
+ fit(b, nb, ipb, 2, bbox);
+
+ for (int j = 0; j < na; ++j) {
+ for (int k = 0; k < nb; ++k) {
+ if (ovl(ipa[j].rx, ipb[k].rx) && ovl(ipa[j].ry, ipb[k].ry)) {
+ hp a1 = -area(ipa[j].ip, ipb[k].ip, ipb[k + 1].ip);
+ hp a2 = area(ipa[j + 1].ip, ipb[k].ip, ipb[k + 1].ip);
+ bool o = a1 < 0;
+ if (o == a2 < 0) {
+ hp a3 = area(ipb[k].ip, ipa[j].ip, ipa[j + 1].ip);
+ hp a4 = -area(ipb[k + 1].ip, ipa[j].ip,
+ ipa[j + 1].ip);
+ if (a3 < 0 == a4 < 0) {
+ if (o)
+ cross(ipa[j], ipa[j + 1], ipb[k], ipb[k + 1],
+ a1, a2, a3, a4);
+ else
+ cross(ipb[k], ipb[k + 1], ipa[j], ipa[j + 1],
+ a3, a4, a1, a2);
+ }
+ }
+ }
+ }
+ }
+
+ inness(ipa, na, ipb, nb);
+ inness(ipb, nb, ipa, na);
+
+ return ssss / ascale;
+ }
+};
+
+#endif //aip_h
Added: trunk/Model/Groups/BSP/chainHull.h
==============================================================================
--- (empty file)
+++ trunk/Model/Groups/BSP/chainHull.h Tue Jun 24 17:27:12 2008
@@ -0,0 +1,56 @@
+#ifndef chainhull_h
+#define chainhull_h
+#include <Core/Geometry/VectorT.h>
+#include <algorithm>
+
+namespace Manta {
+namespace BSPs
+{
+ typedef VectorT<double, 2> Point2D;
+
+ //lexicographic sort
+ struct lexicographic2DPoint {
+ bool operator()(const Point2D &a, const Point2D &b)
+ {
+ return a[0] < b[0] || (a[0] == b[0] && a[1] < b[1]);
+ }
+ };
+
+ //http://www.algorithmist.com/index.php/Monotone_Chain_Convex_Hull.cpp
+ // 2D cross product. Return a positive value, if OAB makes a
+ // counter-clockwise turn, negative for clockwise turn, and zero if
+ // the points are collinear.
+ inline double cross(const Point2D &O, const Point2D &A, const Point2D &B)
+ {
+ return (A[0] - O[0]) * (B[1] - O[1]) -
+ (A[1] - O[1]) * (B[0] - O[0]);
+ }
+
+ // Returns a list of points on the convex hull in counter-clockwise
+ // order. Note: the last point in the returned list is the same as
+ // the first one.
+ void convexHull(vector<Point2D> &P, vector<Point2D> &H)
+ {
+ int n = P.size(), k = 0;
+ H.resize(2*n);
+
+ // Sort points lexicographically
+ sort(P.begin(), P.end(), lexicographic2DPoint());
+
+ // Build lower hull
+ for (int i = 0; i < n; i++) {
+ while (k >= 2 && cross(H[k-2], H[k-1], P[i]) <= 0) k--;
+ H[k++] = P[i];
+ }
+
+ // Build upper hull
+ for (int i = n-2, t = k+1; i >= 0; i--) {
+ while (k >= t && cross(H[k-2], H[k-1], P[i]) <= 0) k--;
+ H[k++] = P[i];
+ }
+
+ H.resize(k);
+ }
+};
+};
+#endif //chainhull_h
Modified: trunk/Model/Groups/CMakeLists.txt
==============================================================================
--- trunk/Model/Groups/CMakeLists.txt (original)
+++ trunk/Model/Groups/CMakeLists.txt Tue Jun 24 17:27:12 2008
@@ -12,6 +12,12 @@
Groups/ObjGroup.h
Groups/ObjGroup.cc
Groups/RecursiveGrid
+ Groups/BSP/BSP
+ Groups/BSP/BSH
+ Groups/BSP/Geometry
+ Groups/BSP/Polytope
+ Groups/BSP/aip.h
+ Groups/BSP/chainHull.h
)
# Include private code if available
Modified: trunk/Model/Groups/KDTree.cc
==============================================================================
--- trunk/Model/Groups/KDTree.cc (original)
+++ trunk/Model/Groups/KDTree.cc Tue Jun 24 17:27:12 2008
@@ -145,7 +145,7 @@
bool bestCommonToTheLeft = false;
vector<Event> event[3];
- if (bounds.computeArea() > 1e-8)
+ if (true ||bounds.computeArea() > 1e-10)
for (unsigned int i=0;i<primitiveID.size();i++) {
// bounds.reset();
// if (currGroup) {
@@ -336,7 +336,7 @@
const int ray_end = rays.end();
if (!rays.getFlag(RayPacket::ConstantSigns)) {
- TrvStack trvStack[80];
+ TrvStack trvStack[512];
for (int first=ray_begin; first < ray_end; ++first) {
int last = first+1;
for (; last < ray_end; ++last) {
@@ -364,7 +364,7 @@
return;
}
else if (ray_begin+1 == ray_end) {
- TrvStack trvStack[80];
+ TrvStack trvStack[512];
const Vector origin = rays.getOrigin( rays.begin() );
const Vector direction = rays.getDirection( rays.begin() );
const Vector inverse_direction = rays.getInverseDirection(rays.begin());
@@ -985,6 +985,7 @@
mesh = dynamic_cast<Mesh*>(currGroup);
+#if 0
ifstream in(file.c_str());
if (!in) return false;
unsigned int itemListSize;
@@ -1032,7 +1033,22 @@
}
node.push_back(n);
}
-
+#else
+ ifstream in(file.c_str(), ios::binary | ios::in);
+
+ unsigned int itemList_size;
+ in.read((char*)&itemList_size, sizeof(itemList_size));
+ itemList.resize(itemList_size);
+
+ in.read((char*) &itemList[0], sizeof(itemList[0])*itemList.size());
+
+ unsigned int node_size;
+ in.read((char*)&node_size, sizeof(node_size));
+ node.resize(node_size);
+
+ in.read((char*) &node[0], sizeof(Node)*node.size());
+#endif
+
in.close();
cout << "done building" << endl << flush;
@@ -1044,10 +1060,11 @@
}
-void KDTree::saveToFile(const string &file)
+bool KDTree::saveToFile(const string &file)
{
+#if 0
ofstream out(file.c_str());
- if (!out) return;
+ if (!out) return false;
out << itemList.size() << endl;
for (unsigned int i=0; i < itemList.size(); ++i) {
@@ -1065,7 +1082,19 @@
out << 0 << " " << node[i].planePos << " " << node[i].planeDim << " "
<< node[i].childIdx<<endl;
}
}
+#else
+ ofstream out(file.c_str(), ios::out | ios::binary);
+ if (!out) return false;
+
+ const unsigned int itemList_size = itemList.size();
+ const unsigned int node_size = node.size();
+ out.write((char*) &itemList_size, sizeof(itemList_size));
+ out.write((char*) &itemList[0], sizeof(itemList[0])*itemList_size);
+ out.write((char*) &node_size, sizeof(node_size));
+ out.write((char*) &node[0], sizeof(Node)*node_size);
+#endif
out.close();
+ return true;
}
void KDTree::printStats()
Modified: trunk/Model/Groups/KDTree.h
==============================================================================
--- trunk/Model/Groups/KDTree.h (original)
+++ trunk/Model/Groups/KDTree.h Tue Jun 24 17:27:12 2008
@@ -170,7 +170,7 @@
const BBox &bounds);
bool buildFromFile(const string &file);
- void saveToFile(const string &file);
+ bool saveToFile(const string &file);
void printStats();
#ifndef SWIG
- [Manta] r2289 - in trunk: Interface Model/Groups Model/Groups/BSP, Thiago Ize, 06/24/2008
Archive powered by MHonArc 2.6.16.
