Biomechanical Joint Model  Author: Anderson Maciel

---

mesh.cpp

Go to the documentation of this file.

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

 FILE               : mesh.C
 CREATION DATE      : Thu Dec  2 12:22:33 MET 1999
 CREATION AUTHOR(S) : aubel 
 ------------------------------------------------------------------------------
 KEYWORD(S)         : marching cubes, voxelisation, volume
 DEFINED TYPE(S)    : Mesh class for marching cubes
 RELATED TYPE(S)    : 

 FUNCTIONALITY      : 
 ------------------------------------------------------------------------------
 LANGAGE            : Ansi C++
 SYSTEM             : UNIX V 6.x / SGI
 GRAPHIC LIBRARY    : SGI  
 ------------------------------------------------------------------------------
 PROJECT DIRECTOR   : D. THALMANN
 LAB                : EPFL- LIG
 --------------------------------------------------------------------------#@*/
#include <vector>
#include <list>
#include <algorithm>

#include "mesh.h"
#include <RayCasting/bbox.h>
#include <LinAlg/vec.h>

//namespace MarchingCubes {

Mesh::Mesh(const VECArray &vertlist, const VECArray &normlist, 
           const vector<int> &trilist, const vector<int> &normIndex) : 
   Primitive(),
   vertices(vertlist), normals(normlist), faces(trilist), normalIndex(normIndex),
   bbox_h(getBoundingBox()) 
{ 
   bbox_h.createHierarchy(vertices, faces, (int)sqrt(static_cast<float>(vertices.size())) );   
}

bool 
Mesh::intersectedBy(Ray &ray) const
{      

    int i;
    REAL u,v;
    return intersectedBy(ray, u, v, i);
}

bool
Mesh::intersectedBy(Ray &ray, REAL &u, REAL &v, int &tri_idx) const
{
   vector<int> triangles;
   bbox_h.getPotentialTriangles(ray, triangles);

   REAL t, u_, v_;
   ray.t = infinity; // used in the loop
   tri_idx = -1;     // used for final test
      
   for (vector<int>::const_iterator tri_it=triangles.begin(); tri_it!=triangles.end(); tri_it++)
   {
      int i = *tri_it;      
      if ( intersect_triangle(ray.O, ray.D, vertices[faces[i]],vertices[faces[i+1]],
                              vertices[faces[i+2]], t, u_, v_) == true )
      {
         // if intersection in segment & closer to origin than intersections previously found
         if ( t>=0.0 && t<ray.t)
         {
            ray.t = t;
            u = u_; v = v_;
            tri_idx = i;
         }
      }
   }

   if (tri_idx != -1)
   {
      ray.P = ray.O + ray.t * ray.D;      
      return (ray.segment==true)? (ray.t <= 1) : true;
   }
   return false;
}

REAL
Mesh::getDistance(const VEC &X) const
{
   Ray ray(X, VEC(3, "0 0 1"));
   vector<REAL> intersections;

   // list all intersections with ray starting from X and pointing (arbitrarily) in +Z
   vector<int> triangles;
   bbox_h.getPotentialTriangles(ray, triangles);

   for (vector<int>::const_iterator tri_it=triangles.begin(); tri_it!=triangles.end(); tri_it++)
   {
      int i = *tri_it;
      REAL u,v;

      if ( intersect_triangle(ray.O, ray.D, vertices[faces[i]],vertices[faces[i+1]],
                              vertices[faces[i+2]], ray.t, u, v) == true )
      {
         if ( fabs(ray.t) < epsilon ) return 0;
         else if (ray.t>0) intersections.push_back(ray.t);
      }
   }

   // if no intersection, point is outside of mesh
   // we don't need to compute how far it is from the mesh
   // as findIntersection does not rely on this value
   if (intersections.empty()) return -1;
   
   vector<REAL>::const_iterator it = min_element(intersections.begin(), intersections.end());
   REAL dist = *it;
   return (intersections.size()%2 == 0)? -dist : dist;
}

  
BBox
Mesh::getBoundingBox()
{
   bbox.makeEmpty();

   // loop thru all vertices
   for (int i=0; i<vertices.size(); i++) bbox.extendBy(vertices[i]);

   return bbox;
}

// the values valp1 and valp2 are useless here, as an analytic solution
// can be computed
void
Mesh::findIntersection(const VEC &P1, const VEC &P2, REAL valp1, REAL valp2,
                       VEC &P, VEC &N) const
{
   REAL u,v;
   int tri_idx;
 
   Ray ray(P1, P2-P1);
   ray.segment = true;

   if (!intersectedBy(ray, u, v, tri_idx))
   {
      cout << "Fatal error: no intersection in Mesh::findIntersection\n";
      return;
   }
   
   P = ray.P;
   N = (1-u-v) * normals[normalIndex[tri_idx]] 
       + u * normals[normalIndex[tri_idx+1]] 
       + v * normals[normalIndex[tri_idx+2]];
   N.normalize();
}

//}




 

---