Text archives Help

---

Author: roni
Date: Sun Feb 10 17:19:38 2008
New Revision: 2044

Modified:
   trunk/Model/Primitives/QuadFacedHexahedron.cc
   trunk/Model/Primitives/QuadFacedHexahedron.h
   trunk/scenes/primtest.cc
Log:
New, much faster, algorithm for hex intersection, based on slabs with
special cases (thanks to Andrew for the idea).

Can't be used with InstanceRT (and other instance classes that don't
preserve the scratchpad).

Still has some problems when eye point is inside the hex, and the code
needs to be cleaned up.

primtest extended to include "simplehex" and "hex" (the latter does
not work properly because of the problems with InstanceRT).


Modified: trunk/Model/Primitives/QuadFacedHexahedron.cc
==============================================================================
--- trunk/Model/Primitives/QuadFacedHexahedron.cc       (original)
+++ trunk/Model/Primitives/QuadFacedHexahedron.cc       Sun Feb 10 17:19:38 
2008
@@ -1,94 +1,304 @@
 #include <Model/Primitives/QuadFacedHexahedron.h>
+#include <Core/Exceptions/BadPrimitive.h>
+#include <Core/Geometry/AffineTransform.h>
 #include <Core/Geometry/BBox.h>
 #include <Core/Geometry/Vector.h>
 #include <Interface/Context.h>
 #include <Interface/RayPacket.h>
 
-using namespace Manta;
+#include <sstream>
 
-const unsigned QuadFacedHexahedron::faceIndex[6][4] = { {0, 3, 2, 1},   
//bottom
-                                                        {0, 4, 7, 3},   
//left
-                                                        {2, 3, 7, 6},   
//back
-                                                        {1, 2, 6, 5},   
//right
-                                                        {0, 1, 5, 4},   
//front
-                                                        {4, 5, 6, 7} }; //top
+using namespace Manta;
 
 QuadFacedHexahedron::QuadFacedHexahedron(Material *matl,
                                          const Vector& v0, const Vector& v1, 
const Vector& v2, const Vector& v3,
-                                         const Vector& v4, const Vector& v5, 
const Vector& v6, const Vector& v7)
+                                         const Vector& v4, const Vector& v5, 
const Vector& v6, const Vector& v7,
+                                         bool coplanar_test)
   : PrimitiveCommon(matl)
 {
-  // Save the vertices.
-  v[0] = v0;
-  v[1] = v1;
-  v[2] = v2;
-  v[3] = v3;
-  v[4] = v4;
-  v[5] = v5;
-  v[6] = v6;
-  v[7] = v7;
+  // Save the lower left front corner, and the upper right back
+  // corner; these will serve as the anchor points for the six planes
+  // bounding the hex.
+  corner[0] = v0;
+  corner[1] = v6;
+
+  // Compute the plane normals.  
+  normal[Bottom] = Cross(v3 - v0, v1 - v0).normal();
+  normal[Left]   = Cross(v4 - v0, v3 - v0).normal();
+  normal[Front]  = Cross(v1 - v0, v4 - v0).normal();
+  normal[Right]  = Cross(v5 - v6, v2 - v6).normal();
+  normal[Back]   = Cross(v2 - v6, v7 - v6).normal();
+  normal[Top]    = Cross(v7 - v6, v5 - v6).normal();
+
+  // Test for coplanarity in each face.
+  if(coplanar_test){
+    bool ok = true;
+    std::stringstream s;
+    if(!coplanar(corner[0], normal[Bottom], v0, v1, v2, v3)){
+      s << "QuadFacedHexahedron: bottom face "
+        << "(v0 = " << v0
+        << ", v1 = " << v1
+        << ", v2 = " << v2
+        << ", v3 = " << v3 << ") is not co-planar";
+      ok = false;
+    }
+    if(!coplanar(corner[0], normal[Left], v0, v3, v7, v4)){
+      s << "QuadFacedHexahedron: left face "
+        << "(v0 = " << v0
+        << ", v3 = " << v3
+        << ", v7 = " << v7
+        << ", v4 = " << v4 << ") is not co-planar";
+      ok = false;
+    }
+    if(!coplanar(corner[0], normal[Front], v0, v1, v5, v4)){
+      s << "QuadFacedHexahedron: front face "
+        << "(v0 = " << v0
+        << ", v1 = " << v1
+        << ", v5 = " << v5
+        << ", v4 = " << v4 << ") is not co-planar";
+      ok = false;
+    }
+    if(!coplanar(corner[1], normal[Top], v4, v5, v6, v7)){
+      s << "QuadFacedHexahedron: top face "
+        << "(v4 = " << v4
+        << ", v5 = " << v5
+        << ", v6 = " << v6
+        << ", v7 = " << v7 << ") is not co-planar";
+      ok = false;
+    }
+    if(!coplanar(corner[1], normal[Right], v5, v1, v2, v6)){
+      s << "QuadFacedHexahedron: right face "
+        << "(v5 = " << v5
+        << ", v1 = " << v1
+        << ", v2 = " << v2
+        << ", v6 = " << v6 << ") is not co-planar";
+      ok = false;
+    }
+    if(!coplanar(corner[1], normal[Back], v3, v2, v6, v7)){
+      s << "QuadFacedHexahedron: back face "
+        << "(v3 = " << v3
+        << ", v2 = " << v2
+        << ", v6 = " << v6
+        << ", v7 = " << v7 << ") is not co-planar";
+      ok = false;
+    }
+
+    if(!ok)
+      throw BadPrimitive(s.str());
+  }
 }
 
 void QuadFacedHexahedron::computeBounds(const PreprocessContext& context,
                                         BBox& bbox) const {
+  Vector v[8];
+  v[0] = corner[0];
+
+  v[1] = intersectThreePlanes(corner[1], normal[Right],
+                              corner[0], normal[Front],
+                              corner[0], normal[Bottom]);
+
+  v[2] = intersectThreePlanes(corner[1], normal[Right],
+                              corner[1], normal[Back],
+                              corner[0], normal[Bottom]);
+
+  v[3] = intersectThreePlanes(corner[0], normal[Left],
+                              corner[1], normal[Back],
+                              corner[0], normal[Bottom]);
+
+  v[4] = intersectThreePlanes(corner[0], normal[Left],
+                              corner[0], normal[Front],
+                              corner[1], normal[Top]);
+
+  v[5] = intersectThreePlanes(corner[1], normal[Right],
+                              corner[0], normal[Front],
+                              corner[1], normal[Top]);
+
+  v[6] = corner[1];
+
+  v[7] = intersectThreePlanes(corner[0], normal[Left],
+                              corner[1], normal[Back],
+                              corner[1], normal[Top]);
+
   for(unsigned i=0; i<8; i++)
     bbox.extendByPoint(v[i]);
 }
 
 void QuadFacedHexahedron::intersect(const RenderContext& context,
                                     RayPacket& rays) const {
-  // Each face is a quad, so compute the intersection with each one
-  // and pick the smallest positive t value.
-  for(int i=rays.begin(); i<rays.end(); i++){
-    //Real smallest_t = std::numeric_limits<Real>::max();
-    Real smallest_t = rays.getMinT(i);
-    int hitface = -1;
-
-    for(unsigned j=0; j<6; j++){
-      const Vector *vv[] = {v + faceIndex[j][0],
-                            v + faceIndex[j][1],
-                            v + faceIndex[j][2],
-                            v + faceIndex[j][3]};
-
-      // Ray plane intersection, with the plane containing v0->v3,
-      const Vector normal = Cross(*(vv[1])-*(vv[0]), *(vv[3])-*(vv[0]));
+  const int debugFlag = rays.getAllFlags() & RayPacket::DebugPacket;
+  if(debugFlag)
+    std::cout << "QuadFacedHexadron::intersect()" << std::endl;
 
-      // Check if ray is parallel to plane.
-      Real dn = Dot(rays.getDirection(i), normal);
-      if(dn == 0)
-        continue;
+  Real t[6];
 
-      Real ao = Dot(*(vv[0])-rays.getOrigin(i), normal);
-      Real t = ao/dn;
+  for(int i=rays.begin(); i<rays.end(); i++){
+    // TODO(choudhury): Make the slab intersection computation
+    // cleaner.  Sort the results between slabs early and eliminate
+    // half of the test cases ("if(t[f] < t[f+1]) ... else if(t[f+1] <
+    // t[f]) ...") later on.  Also make the entire special case testing 
cleaner.
+
+    // TODO(choudhury): Rendering from inside the hex is broken for
+    // some reason.  Probably has to do with incomplete special case
+    // testing.
+
+    for(int f=0; f<6; f += 2){
+      Real dn0 = Dot(rays.getDirection(i), normal[f]);
+      Real dn1 = Dot(rays.getDirection(i), normal[f+1]);
+      Real ao0 = 0.0, ao1 = 0.0;
 
-      const Vector p = rays.getOrigin(i) + t*rays.getDirection(i);
+      if(dn0 == 0){
+        t[f] = std::numeric_limits<Real>::max();
+      }
+      else{
+        ao0 = Dot(corner[0] - rays.getOrigin(i), normal[f]);
+        t[f] = ao0/dn0;
+      }
 
-      // See if the hit position falls inside the quad.
-      bool inside = true;
-      for(int k=0; k<4; k++){
-        const Real d = Dot(Cross(*(vv[(k+1)%4]) - *(vv[k]), p - *(vv[k])), 
normal);
-        if(d < 0.0){
-          inside = false;
-          break;
-        }
+      if(dn1 == 0){
+        t[f+1] = std::numeric_limits<Real>::max();
+      }
+      else{
+        ao1 = Dot(corner[1] - rays.getOrigin(i), normal[f+1]);
+        t[f+1] = ao1/dn1;
       }
 
-      if(inside && (t < smallest_t) ){
-        hitface = j;
-        smallest_t = t;
+      // Check if we are inside the slab, on the same side as the hex
+      // ("cis").
+      if(-ao0 <= 0.0 && -ao1 <= 0.0){
+        // For both t's positive, entry occurs at negative infinity,
+        // and exit occurs at the least positive of the two forward
+        // intersections.
+        if(t[f] < t[f+1] && t[f] >= 0.0)
+          t[f+1] = -std::numeric_limits<Real>::max();
+        else if(t[f+1] < t[f] && t[f+1] >= 0.0)
+          t[f] = -std::numeric_limits<Real>::max();
+
+        // However if both are negative, then the entry is the max and
+        // the exit is infinity.
+        if(t[f] > t[f+1] && t[f] <= 0.0)
+          t[f+1] = std::numeric_limits<Real>::max();
+        else if(t[f+1] > t[f] && t[f+1] <= 0.0)
+          t[f] = std::numeric_limits<Real>::max();
+      }
+      else if(-ao0 >= 0.0 && -ao1 >= 0.0){
+        // If we are inside the slab but on the other side from the
+        // hex ("trans"), the situation is converse of the above case.
+        if(t[f] > t[f+1] && t[f+1] >= 0.0)
+          t[f+1] = std::numeric_limits<Real>::max();
+        else if(t[f+1] > t[f] && t[f] >= 0.0)
+          t[f] = std::numeric_limits<Real>::max();
+
+        // Ditto for t's both negative.
+        if(t[f] < t[f+1] && t[f+1] <= 0.0)
+          t[f+1] = -std::numeric_limits<Real>::max();
+        else if(t[f+1] < t[f] && t[f] <= 0.0)
+          t[f] = -std::numeric_limits<Real>::max();
+      }
+      else{
+        // If we're outside there are special cases too.  If one t is
+        // negative and the other positive, really we are entering the
+        // slab at the positive t, and then never leaving, so convert
+        // the negative t to positive infinity.
+        if(t[f] < 0.0 && t[f+1] > 0.0)
+          t[f] = std::numeric_limits<Real>::max();
+        else if(t[f+1] < 0.0 && t[f] > 0.0)
+          t[f+1] = std::numeric_limits<Real>::max();
       }
     }
 
-    if(hitface >= 0){
-      rays.hit(i, smallest_t, getMaterial(), this, getTexCoordMapper());
-      rays.getScratchpad<int>(0)[i] = hitface;
+    if(debugFlag){
+      std::cout << "Bottom: " << t[0] << std::endl
+                << "Left  : " << t[1] << std::endl
+                << "Front : " << t[2] << std::endl
+                << "Right : " << t[3] << std::endl
+                << "Back  : " << t[4] << std::endl
+                << "Top   : " << t[5] << std::endl;
+    }
+
+    // Compute the max t over the set of minimum t, one from each
+    // opposing pair of planes, and the "min over max".
+    Real exit_t = Min(Min(Max(t[Front], t[Back]),
+                          Max(t[Left], t[Right])),
+                      Max(t[Bottom], t[Top]));
+
+    int nearest[] = {Back, Right, Top};
+
+    Real vals[3] = {t[Back], t[Right], t[Top]};
+    if(t[Front] < t[Back]){
+      nearest[0] = Front;
+      vals[0] = t[Front];
+    }
+    if(t[Left] < t[Right]){
+      nearest[1] = Left;
+      vals[1] = t[Left];
+    }
+    if(t[Bottom] < t[Top]){
+      nearest[2] = Bottom;
+      vals[2] = t[Bottom];
+    }
+
+    int faceIdx=-1;
+    if(vals[0] > vals[1] && vals[0] > vals[2])
+      faceIdx = 0;
+    else if(vals[1] > vals[2])
+      faceIdx = 1;
+    else
+      faceIdx = 2;
+
+    Real entry_t = vals[faceIdx];
+    int hitface = nearest[faceIdx];
+
+    if(debugFlag){
+      std::cout << "entry_t = " << entry_t << std::endl
+                << "exit_t  = " << exit_t << std::endl;
+    }
+
+    // If the exit_t is less than the entry_t, then there was no 
intersection.
+    if(exit_t < entry_t)
+      continue;
+
+    // If the entry_t is negative, use the exit_t.
+    if(entry_t < 0.0){
+      if(rays.hit(i, exit_t, getMaterial(), this, getTexCoordMapper()))
+        rays.getScratchpad<int>(0)[i] = hitface;
     }
+    else if(rays.hit(i, entry_t, getMaterial(), this, getTexCoordMapper()))
+      rays.getScratchpad<int>(0)[i] = hitface;
   }
 }
 
 void QuadFacedHexahedron::computeNormal(const RenderContext& context,
                                         RayPacket& rays) const {
-  for(int i=rays.begin(); i<rays.end(); i++)
-    rays.setNormal(i, Cross(v[faceIndex[rays.getScratchpad<int>(0)[i]][1]] - 
v[faceIndex[rays.getScratchpad<int>(0)[i]][0]],
-                            v[faceIndex[rays.getScratchpad<int>(0)[i]][3]] - 
v[faceIndex[rays.getScratchpad<int>(0)[i]][0]]));
+  // NOTE(choudhury): using this class with an RTInstance (and
+  // possibly other types of instances) corrupts the scratchpad and
+  // this section of code becomes seemingly buggy.  Don't use this
+  // class with RTInstance and other offending instance classes.
+  for(int i=rays.begin(); i<rays.end(); i++){
+    const int which = rays.getScratchpad<int>(0)[i];
+    rays.setNormal(i, normal[which]);
+  }
+}
+
+Vector QuadFacedHexahedron::intersectThreePlanes(const Vector& p1, const 
Vector& N1,
+                                                 const Vector& p2, const 
Vector& N2,
+                                                 const Vector& p3, const 
Vector& N3){
+  AffineTransform A_inv;
+  A_inv.initWithBasis(Vector(N1[0], N2[0], N3[0]),
+                      Vector(N1[1], N2[1], N3[1]),
+                      Vector(N1[2], N2[2], N3[2]),
+                      Vector(0.0,0.0,0.0));
+  A_inv.invert();
+
+  Vector b(Dot(p1, N1), Dot(p2, N2), Dot(p3, N3));
+
+  return A_inv.multiply_vector(b);
+}
+
+bool QuadFacedHexahedron::coplanar(const Vector& p, const Vector& N,
+                                   const Vector& v0, const Vector& v1, const 
Vector& v2, const Vector& v3){
+  static const Real eps = 1.0e-3;
+  return (Dot(v0-p, N) < eps &&
+          Dot(v1-p, N) < eps &&
+          Dot(v2-p, N) < eps &&
+          Dot(v3-p, N) < eps);
 }

Modified: trunk/Model/Primitives/QuadFacedHexahedron.h
==============================================================================
--- trunk/Model/Primitives/QuadFacedHexahedron.h        (original)
+++ trunk/Model/Primitives/QuadFacedHexahedron.h        Sun Feb 10 17:19:38 
2008
@@ -10,7 +10,8 @@
   public:
     QuadFacedHexahedron(Material *matl,
                         const Vector& v0, const Vector& v1, const Vector& 
v2, const Vector& v3,
-                        const Vector& v4, const Vector& v5, const Vector& 
v6, const Vector& v7);
+                        const Vector& v4, const Vector& v5, const Vector& 
v6, const Vector& v7,
+                        bool coplanar_test = false);
 
     void computeBounds(const PreprocessContext& context,
                       BBox& bbox_) const;
@@ -22,10 +23,23 @@
                       RayPacket& rays) const;
 
   private:
-    static const unsigned faceIndex[6][4];
+    static Vector intersectThreePlanes(const Vector& p1, const Vector& N1,
+                                       const Vector& p2, const Vector& N2,
+                                       const Vector& p3, const Vector& N3);
+
+    static bool coplanar(const Vector& p, const Vector& N,
+                         const Vector& v0, const Vector& v1, const Vector& 
v2, const Vector& v3);
+
+    enum Face {
+      Left = 0,
+      Right,
+      Bottom,
+      Top,
+      Front,
+      Back
+    };
 
-    Vector v[8];
+    Vector corner[2], normal[6];
   };
 }
-
 #endif

Modified: trunk/scenes/primtest.cc
==============================================================================
--- trunk/scenes/primtest.cc    (original)
+++ trunk/scenes/primtest.cc    Sun Feb 10 17:19:38 2008
@@ -38,6 +38,7 @@
 #include <Model/Primitives/Heightfield.h>
 #include <Model/Primitives/Hemisphere.h>
 #include <Model/Primitives/Parallelogram.h>
+#include <Model/Primitives/QuadFacedHexahedron.h>
 #include <Model/Primitives/Sphere.h>
 #include <Model/Primitives/SuperEllipsoid.h>
 #include <Model/Primitives/Torus.h>
@@ -306,6 +307,30 @@
         group->add( prim );
       }
     }
+  } else if(primtype == "simplehex"){
+    Real x = scale/max;
+
+    for(int i=0;i<numx;i++){
+      for(int j=0;j<numy;j++){
+        Vector p((numx>1 ? i/static_cast<Real>(numx-1)-(Real)0.5 : 
0)*scale*2,
+                 (numy>1 ? j/static_cast<Real>(numy-1)-(Real)0.5 : 
0)*scale*2,
+                 0);
+        Primitive* prim = new QuadFacedHexahedron(matl,
+                                                  p + Vector(-1.5*x, -1.5*x, 
-x),
+                                                  p + Vector( 1.5*x, -1.5*x, 
-x),
+                                                  p + Vector( 1.5*x,  1.5*x, 
-x),
+                                                  p + Vector(-1.5*x,  1.5*x, 
-x),
+
+                                                  p + Vector(-0.5*x, -0.5*x, 
 x),
+                                                  p + Vector( 0.5*x, -0.5*x, 
 x),
+                                                  p + Vector( 0.5*x,  0.5*x, 
 x),
+                                                  p + Vector(-0.5*x,  0.5*x, 
 x),
+                                                  true); // "true" means to 
test each face for co-planarity.
+        if ( mapr )
+          prim->setTexCoordMapper( mapr );
+        group->add( prim );
+      }
+    }
   } else if(primtype == "sphere"){
     Primitive* prim = new Sphere(matl, Vector(0,0,0), scale/max);
     if ( mapr )
@@ -314,6 +339,20 @@
   } else if(primtype == "box"){
     Vector p2(scale/max/1.732, scale/max/1.732, scale/max/1.732);
     Primitive* prim = new Cube(matl, -p2, p2);
+    if ( mapr )
+      prim->setTexCoordMapper( mapr );
+    spinprim = prim;
+  } else if(primtype == "hex"){
+    Real x = scale/max/1.732;
+    Primitive* prim = new QuadFacedHexahedron(matl,
+                                              Vector(-2*x, -x, -1.5*x),
+                                              Vector( x,  x, -x),
+                                              Vector( x,  x, -x),
+                                              Vector(-x,  x, -x),
+                                              Vector(-x, -x,  x),
+                                              Vector( x, -x,  x),
+                                              Vector( x,  x,  x),
+                                              Vector(-x,  x,  x));
     if ( mapr )
       prim->setTexCoordMapper( mapr );
     spinprim = prim;

---

• [Manta] r2044 - in trunk: Model/Primitives scenes, roni, 02/10/2008

---

Archive powered by MHonArc 2.6.16.