Biomechanical Joint Model  Author: Anderson Maciel

---

intersect_tri.cpp

Go to the documentation of this file.

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

 FILE               : intersect_tri.C
 CREATION DATE      : Thu Sep  2 09:52:30 MDT 1999
 CREATION AUTHOR(S) : aubel
 ------------------------------------------------------------------------------
 KEYWORD(S)         : Ray-triangle intersection
 DEFINED TYPE(S)    : 
 RELATED TYPE(S)    : 

 FUNCTIONALITY      : Compute the intersection of a ray and a triangle
 ------------------------------------------------------------------------------
 LANGAGE            : Ansi C++
 SYSTEM             : UNIX V 6.x / SGI
 GRAPHIC LIBRARY    : SGI  
 ------------------------------------------------------------------------------
 PROJECT DIRECTOR   : D. THALMANN
 LAB                : EPFL- LIG
 --------------------------------------------------------------------------#@*/

#include <iostream>
//#include <vector>
#include <algorithm>

#include <LinAlg/vecarray.h>
#include "ray.h"
#include <list>


// Calculates the intersection of the line defined by an origin orig and a
// direction dir with triangle [v0,v1,v2].
// If the line intercepts the triangle, true is returned.
// When true is returned,  t,u,v will satisfy:
//  (1-u-v).v0 + u.v1 + v.v2 = (1-t).orig + t.(orig+dir) = intersection point,
// so that u,v indicates a parametric distance from the vertices and t a
// parametric distance that can be used to determine if the segment
// [orig,orig+dist] intersects in fact the triangle.
// Note that t can be negative.
// Original code : http://www.acm.org/jgt/papers/MollerTrumbore97/
bool 
intersect_triangle ( const VEC &orig, const VEC &dir,
                     const VEC &v0, const VEC &v1, const VEC &v2,
                     REAL &t, REAL &u, REAL &v )
{
   VEC edge1 = v1 - v0;        // find vectors for two edges sharing v0 
   VEC edge2 = v2 - v0;

   // begin calculating determinant - also used to calculate U parameter 
   VEC pvec = cross_prod(dir, edge2);

   REAL det = dot_prod(edge1, pvec);  // if determinant is near zero, ray lies in plane of 

   if (fabs(det) < epsilon) 
   { 
       return false; 
   }
   REAL inv_det = 1.0 / det;

   VEC tvec = orig - v0;                     // calculate distance from v0 to ray origin 
   u = dot_prod(tvec, pvec) * inv_det;       // calculate U parameter and test bounds 
   if ( u<0.0 || u>1.0 ) 
   { 
       return false; 
   }

   VEC qvec = cross_prod(tvec, edge1);            // prepare to test V parameter 
   v = dot_prod(dir, qvec) * inv_det;        // calculate V parameter and test bounds 
   if ( v<0.0 || u+v>1.0 ) 
   { 
       return false; 
   }
   t = dot_prod(edge2,qvec) * inv_det;       // calculate t, ray intersects triangle 

   return true;
}

//
// For a point X, look for m closest points from list Y within a specified radius.  
// The accepted template types are containers of VECs : VECArray, vector<VEC>, list<VEC> etc.
//
vector<int> find_closest(const VEC &X, const VECArray &Y, int m, REAL radius)
{
   vector<int> closest;
   // use squared radius to avoid square roots in loop
   REAL r2 = (radius >= infinity)? infinity : radius*radius;

   while (closest.size() < m)
   {
      REAL min_dist=infinity;
      int closest_node=-1;
         
      for (int j=0; j<Y.size(); j++)
      {
         // skip points already chosen 
         if (find(closest.begin(), closest.end(), j) != closest.end()) continue;

         REAL dist = VEC(Y[j] - X).norm2();
         if (dist < min_dist && dist <= r2)
         {
            min_dist = dist;
            closest_node = j;          
         }
      }
      // no more close node ?
      if (closest_node == -1) break;
      
      closest.push_back(closest_node);
   }
   return closest;
}

// Project a point X onto a mesh defined by a list of vertices, triangle normals
// and a list of triangles, which is a continuous list (no separator between triangles) 
// of indices into the vertices list. Therefore, trilist's size must be a multiple of 3.
// The projection is defined as the point on the mesh closest to X.
//
// The returned values is the index(*3) of the triangle that contains P = projection of X, 
// while u,v indicates P's barycentric coordinates in the triangle: (1-u-v).V0 + u.V1 + v.V2 
//
// The accepted template types are containers of VECs : VECArray, vector<VEC>, list<VEC> etc.
//
// The projection is  found by taking the minimum of the:
// i) distance to every vertex of the mesh
// ii) distance to every triangle of the mesh 
// The distance to a triangle is computed as follows: a ray is cast from start point X with
// a direction colinear to the triangle's normal. If a triangle-ray intersection exists, then
// we use the distance from the intersection point to X. Otherwise, it's infinity.

int project_on_mesh_tmp(const VEC &X, const VECArray &vertlist, 
                       const VECArray &normlist, const vector<int> &trilist, 
                       REAL &u, REAL &v, VEC &P)
{
   // extract vertices actually in use (those referenced by the triangles list)
   // start by removing all duplicates from list of vertices indices
   list<int> tmp(trilist.size());
   copy(trilist.begin(), trilist.end(), tmp.begin());
   tmp.sort(); tmp.unique();
   VECArray V(tmp.size());
   int i=0; vector<int> no_dup;
   for (list<int>::const_iterator vtx = tmp.begin(); vtx != tmp.end(); ++vtx)
   {
      V[i++] = vertlist[*vtx];
      no_dup.push_back(*vtx);
   } 

   // find closest vertex to X and its index
   int closest_vtx = find_closest(X, V, 1).front();
   int idx_closest_vtx = no_dup[closest_vtx];
   P = V[closest_vtx];
   REAL min_dist = VEC(P - X).norm2();

   vector<int> closest_tris;  // indices of triangles that contain the closest vertex

   // find closest triangle
   i=0;
   int idx_closest_tri = -1;
   for (vector<int>::const_iterator tri = trilist.begin(); tri != trilist.end(); i++)
   {
      int p[3]; // store three indices of triangle here
      REAL t , u_, v_;
      bool has_closest_vtx = false;

      // loop on each vertex of triangle
      for (int k=0; k<3; k++)
      {
         p[k] = *tri++;
         if ( p[k] == idx_closest_vtx ) has_closest_vtx=true;
      }
      if ( has_closest_vtx ) closest_tris.push_back(i*3);
      
      if ( intersect_triangle (X, normlist[i], vertlist[p[0]], vertlist[p[1]], vertlist[p[2]], 
                               t, u_, v_ ) == true )
      {
         VEC D = t * normlist[i];
         if ( REAL d=D.norm2() < min_dist)
         {
            min_dist = d;
            idx_closest_tri = i*3;
            P = X + D;
            u = u_; v = v_;
         }
      }
   }
   
   // is closest vertex closer than the closest triangle?
   if ( idx_closest_tri == -1)
   {
      // yes, then pick up among all triangles that contain the closest vertex, 
      // the closest triangle (using the barycenters)
      VECArray barycenters(closest_tris.size());
      for (int k=0; k<closest_tris.size(); k++)
      {
         vector<int>::const_iterator tri = trilist.begin() + closest_tris[k];
         barycenters[k] = vertlist[*tri++] + vertlist[*tri++] + vertlist[*tri++];
         barycenters[k] /= 3;
      }
      int closest_barycenter = find_closest(X, barycenters, 1).front();
      idx_closest_tri = closest_tris[closest_barycenter];
      // there remains to compute the parametric u,v
      if (trilist[idx_closest_tri] == idx_closest_vtx) u=v=0;
      else
      { 
         if (trilist[idx_closest_tri+1] == idx_closest_vtx) { u=1; v=0; }
         else { u=0; v=1; }
      }
   }
   return idx_closest_tri;
}


// normals are automatically computed by the function
// The accepted template types are containers of VECs : VECArray, vector<VEC>, list<VEC> etc.
// and containers of ints for the triangle list
int project_on_mesh(const VEC &X, const  VECArray &vertlist, const vector<int> &trilist,
                       REAL &u, REAL &v, VEC &P)
{  
   assert(trilist.size()%3 == 0);

   // compute normals to triangles first
   VECArray normlist(trilist.size());   
   int i=0; 
   VEC V1(dimspace), V2(dimspace);
   for (vector<int>::const_iterator tri = trilist.begin(); tri != trilist.end(); i++)
   {    
      int p0 = *tri++;  
      int p1 = *tri++;  
      int p2 = *tri++;    
      V1 = vertlist[p1] - vertlist[p0];
      V2 = vertlist[p2] - vertlist[p0];
      normlist[i] = cross_prod(V1,V2); // no need to normalize
   }

   return project_on_mesh_tmp(X, vertlist, normlist, trilist, u, v, P);
}

---