Biomechanical Joint Model  Author: Anderson Maciel

---

bbox.cpp

Go to the documentation of this file.

/*@#---------------------------------------------------------------------------
 PROJECT            : PHD THESIS
 MODULE             : PHYSICALLY-BASED DEFORMATION

 FILE               : bbox.C
 CREATION DATE      : Fri Dec 10 11:58:19 MET 1999
 CREATION AUTHOR(S) : aubel 
 ------------------------------------------------------------------------------
 KEYWORD(S)         : bounding box
 DEFINED TYPE(S)    : 
 RELATED TYPE(S)    : 

 FUNCTIONALITY      : Bounding Box class to avoid relying on Inventor's bboxes
 ------------------------------------------------------------------------------
 LANGAGE            : Ansi C++
 SYSTEM             : UNIX V 6.x / SGI
 GRAPHIC LIBRARY    : SGI  
 ------------------------------------------------------------------------------
 PROJECT DIRECTOR   : D. THALMANN
 LAB                : EPFL- LIG
 --------------------------------------------------------------------------#@*/

#include "bbox.h"


void 
BBox::extendBy(const VEC &X)
{
   assert(X.dim()==n_);
   
   if (empty==true)
   {
      UR = LL = X;
      empty = false;
   }
   else
   {
      for (int k=0; k<n_; k++) 
      {
         if (X[k]<LL[k]) LL[k] = X[k];
         else if (X[k]>UR[k]) UR[k] = X[k];
      }
   }
}

void
BBox::extendBy(const BBox &bb)
{
   if (bb.isEmpty()) return;

   VEC V1(bb.getMin()), V2(bb.getMax());
   VEC Y(n_);
 
   // for each vertex of BBox
   for (int p=0; p<(0x2<<n_); p++)
   {     
      for (int k=0; k<n_; k++) 
         Y[k] = ( p &(0x1<<k) )? V1[k] : V2[k];      
      extendBy(Y);
   }
}

void
BBox::scale(REAL factor)
{
   // get center
   VEC C = 0.5*(LL + UR);
   // get dim with scaling
   REAL mid_scale = 0.5 * factor;
   VEC D = mid_scale*(UR - LL);
   LL = C-D;
   UR = C+D;
}

bool
BBox::intersectedBy(const BBox &bb) const
{
   VEC D(bb.getDim());
   for (int p=0; p<(2<<n_); p++)
   {
      VEC P(bb.getMin());           
      for (int k=0; k<n_; k++) 
      {
         if (p & (0x1<<k)) P[k] += D[k];
      }
      if ( contains(P) ) return true;
   }
   return false;
}

// modified from Graphics Gems I pp. 736
bool
BBox::intersectedBy(Ray &ray) const
{
   static const int RIGHT = 0;
   static const int LEFT = 1;
   static const int MIDDLE = 2;
    
   bool inside = true;
   REAL candidatePlane[dimspace]; // n_ = 2 or 3, so is always <= dimspace
   int  quadrant[dimspace];
   int  whichPlane, i;
   REAL maxT[dimspace];

   REAL O[3], D[3],P[3];
   for (i=0;i<3;i++) { O[i] = ray.O[i]; D[i] = ray.D[i]; };

   for (i=0; i<n_; i++)
   {
      if (O[i] < LL[i])
      {
         quadrant[i] = LEFT;
         candidatePlane[i] = LL[i];
         inside = false;
      }
      else if (O[i] > UR[i])
      {
         quadrant[i] = RIGHT;
         candidatePlane[i] = UR[i];
         inside = false;
      }
      else quadrant[i] = MIDDLE;
   }

   if (inside) 
   {
      if ( ray.segment == true)
      {
         VEC newO = ray.O + ray.D;
         if ( contains(newO) ) return false;

         // test with reverse ray        
         Ray rray(newO, -ray.D);
         rray.segment = true;
         if (intersectedBy(rray))
         {
            ray.t = 1 - rray.t;
            ray.P = rray.P;
            return true;
         }
         else cout << "Fatal error: could not find intersection in BBox::intersectedBy\n";
         return false;
      }
      else
      {
         // test with reverse ray
         REAL Dnorm = ray.D.norm();
         REAL coeff = 2*getDim().norm() / Dnorm;
         Ray rray(ray.O + coeff*ray.D, -ray.D); 
         if (!intersectedBy(rray))
         {
            // last remaining possibility, ray.O is on one of the bbox faces
            bool found  = false;
            for (i=0; !found && i<n_; i++)
            {
               if (fabs(O[i]-LL[i]) < epsilon)
                  found = true;
               if (fabs(O[i]-UR[i]) < epsilon)
                  found = true;
            }
            if (!found)
            {
               cout << "Fatal error: couldn't find intersection in BBox::intersectedBy\n";
               return false;
            }
            rray.P = ray.O;
         }
         ray.P = rray.P;
         ray.t = VEC(ray.P - ray.O).norm() / Dnorm;
         return ray.t >= 0;
      }
   }
   
   // Calculate T distances to candidate planes
   for (i = 0; i < n_; i++)
      if (quadrant[i] != MIDDLE && D[i] !=0.)
         maxT[i] = (candidatePlane[i]-O[i]) / D[i];
      else
         maxT[i] = -1.;

   // Get largest of the maxT's for final choice of intersection
   whichPlane = 0;
   for (i = 1; i < n_; i++)
      if (maxT[whichPlane] < maxT[i])
         whichPlane = i;


   // Check final candidate actually inside box
   if (maxT[whichPlane] < 0.) return (false);
   if (ray.segment == true && maxT[whichPlane] > 1.) return false;
   ray.t = maxT[whichPlane];

   for (i = 0; i < n_; i++)
      if (whichPlane != i) {
         P[i] = O[i] + maxT[whichPlane] *D[i];
         if ( P[i] > UR[i] || P[i] < LL[i] )
            return (false);      // outside box
      }else {
         P[i] = candidatePlane[i];
      }
   ray.P = VEC(3,P);
   return (true);                          // ray hits box
}       


// intersection with a triangle
// algorithm simplified from Graphics Gems III pp.236
bool 
BBox::intersectedBy(const VEC &V0, const VEC &V1, const VEC &V2) const
{
   // start with a simple acceptance test
   if (contains(V0) || contains(V1) || contains(V2) ) return true;

   // test if all 3 vertices of triangle not on the same "side" of bbox
   // simple rejection test
   for (int i=0; i<n_; i++)
   {
      if ( (V0[i] - LL[i] < 0) && (V1[i] - LL[i] < 0) && (V2[i] - LL[i] < 0) )
         return false;

      if ( (V0[i] - UR[i] > 0) && (V1[i] - UR[i] > 0) && (V2[i] - UR[i] > 0) )
         return false;
   }
      
   // ok, it's not a simple case
   Ray ray(V0, V1-V0);
   ray.segment = true;
   if ( intersectedBy(ray) ) return true;

   ray.O = V1; ray.D = V2-V1;
   if ( intersectedBy(ray) ) return true;

   ray.O = V2; ray.D = V0-V2;
   if ( intersectedBy(ray) ) return true;

   // it's really not a simple case!
   // we still have to test whether the diagonals of the cube intersect the triangle
   
   REAL t,u,v;
   VEC P0 (LL), P1 (UR);
   if ( intersect_triangle(P0, P1-P0, V0, V1, V2, t, u, v) && t>=0 && t<=1) return true;

   swap(P0[1], P1[1]);
   if ( intersect_triangle(P0, P1-P0, V0, V1, V2, t, u, v) && t>=0 && t<=1) return true;

   swap(P0[0], P1[0]);
   if ( intersect_triangle(P0, P1-P0, V0, V1, V2, t, u, v) && t>=0 && t<=1) return true;

   swap(P0[1], P1[1]);
   return ( intersect_triangle(P0, P1-P0, V0, V1, V2, t, u, v) && t>=0 && t<=1);
}

BBoxHierarchy::~BBoxHierarchy()
{
   for (vector<BBoxHierarchy*>::iterator it=bboxlist.begin(); it!=bboxlist.end(); ++it)
      delete (*it);
}

// handle triangles only as input
void
BBoxHierarchy::createHierarchy(const VECArray &coords, const vector<int> &tris, 
                               int num_tris_max, int level)
{   
   // check which triangles are inside
   for (int i=0; i<tris.size(); i+=3)
   {
      if (bb.intersectedBy(coords[tris[i+0]], coords[tris[i+1]], coords[tris[i+2]]))
      {  
         triangles.push_back(i);
      }
   }

//    cout << " level = "<< level << " with bbox " << bb.getMin() << " " << bb.getDim()
//      << " has " << triangles.size() << " tris. Max allowed = " 
//      << num_tris_max << endl;

   // if too many triangles, subdivise
   // do not permit more than 6 levels of recursion as a security
   if ( level<5 && triangles.size() > num_tris_max )
   {
      VEC O = bb.getMin();
      VEC D = bb.getDim()/2;

      // for each octree
      for (int p=0;p<8;p++)
      {
         VEC newO = O;
         for (int k=0; k<3; k++) 
         {
            if (p & (0x1<<k)) newO[k] += D[k];
         }

         // recursive process
         BBoxHierarchy *smallerbb = new BBoxHierarchy(newO, newO+D, level+1);            
         smallerbb->createHierarchy(coords, tris, num_tris_max, level+1);

         if (smallerbb->getNumTriangles() > 0)
            bboxlist.push_back(smallerbb);
         else 
            delete smallerbb;    
      }
   }     
}


void
BBoxHierarchy::getPotentialTriangles(Ray &ray, vector<int> &trilist) const
{   
   bool force = bb.contains(ray.O) && bb.contains(ray.O+ray.D);

   if ( force || bb.intersectedBy(ray) )
   {
      if (bboxlist.empty() == false)
      {
         for (vector<BBoxHierarchy*>::const_iterator bb_it = bboxlist.begin();
              bb_it != bboxlist.end(); bb_it++)
         {
            (*bb_it)->getPotentialTriangles(ray, trilist);
         }
      }
      else 
      {
         if (trilist.empty() == true)
         {
            trilist.insert(trilist.end(), triangles.begin(), triangles.end());
            return;
         }
         
         int prev_size = trilist.size();

         // for each tri. check it hasn't been inserted before
         for (vector<int>::const_iterator tri_it = triangles.begin();
              tri_it != triangles.end(); tri_it++)
         {
            int i=0;
            for (; i<prev_size && *tri_it != trilist[i]; i++);
            if (i == prev_size) trilist.push_back(*tri_it);
         }
      }
   }
}

const BBoxHierarchy*
BBoxHierarchy::getSmallestBBox(const VEC &X, int m) const
{
   if ( !bb.contains(X) ) return NULL;
   
   for (vector<BBoxHierarchy*>::const_iterator bb_it=bboxlist.begin();
        bb_it != bboxlist.end(); ++bb_it)
   {
      if ( (*bb_it)->bb.contains(X) && (*bb_it)->getNumTriangles()>=m ) 
         return (*bb_it)->getSmallestBBox(X, m);
   }

   return this;      
}

---