Biomechanical Joint Model  Author: Anderson Maciel

---

primitive.cpp

Go to the documentation of this file.

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

 FILE               : primitive.C
 CREATION DATE      : Tue May  9 10:10:40 MDT 2000
 CREATION AUTHOR(S) : aubel 
 ------------------------------------------------------------------------------
 KEYWORD(S)         : marching cubes, voxelisation, volume
 DEFINED TYPE(S)    : 
 RELATED TYPE(S)    : 

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

//namespace MarchingCubes{  
 
bool
ImplicitPrimitive::intersectedBy(Ray &ray) const
{
   assert(ray.segment == true); // intersection w.r.t. infinite ray not implemented yet 

   fprintf(stderr,"ImplicitPrimitive::intersectedBy\n");
   REAL val_ref = getDistance(ray.O);
   if (val_ref < 0)
   {      
      VEC P(3), Pprev(3);
      REAL ini = 1E-4;
      REAL step = ini;
      int j=0;

      for (REAL i=1; step < 0.9; i+=0.1, j++)
      {
        P = ray.O + step*ray.D;
        REAL val = getDistance(P);
#if 0
        printf("P=[%lf,%lf,%lf]->%lf, ray.O=[%lf,%lf,%lf]->%lf, step*ray.D=[%lf,%lf,%lf]\n",P[0],P[1],P[2],\
               val, ray.O[0],ray.O[1],ray.O[2], val_ref, step*ray.D[0],step*ray.D[1],step*ray.D[2]);
        getchar();
#endif

         if (val >= -epsilon)
         {
           //ray.P = dichotomy(ray.O, P, val_ref, val);
           ray.P = dichotomy(Pprev, P, val_ref, val);
           ray.t = VEC(ray.P-ray.O).norm() / ray.D.norm();
           return true;
         }
         
         step = pow(ini, 1/i);
         Pprev=P; val_ref=val;
//       cout << step<< " ";
//       if (j%16==0) cout << endl;
      }
   }
   return false;
}

REAL 
ImplicitPrimitive::getDistance(const VEC &X) const 
{ 
   REAL dist = (!bbox.isEmpty() && bbox.contains(X))? density(X)-iso_value : -iso_value;
//   cout << "distance = " << dist << endl;

   return (!bbox.isEmpty() && bbox.contains(X))? density(X)-iso_value : -iso_value; 
}

// return point at iso-surface between P1 and P2
// if no such point exists, return a 0-length vector
/*
VEC ImplicitPrimitive::dichotomy(const VEC &P1, const VEC &P2, REAL valp1, REAL valp2) const
{
#if 0
   printf("dichotomy:P1[%lf,%lf,%lf]<->P2[%lf,%lf,%lf], valp1=%lf, valp2=%lf\n",P1[0],P1[1],P1[2],P2[0],P2[1],P2[2],valp1,valp2);
   getchar();
#endif
   if (abs(valp1) < epsilon) return P1;   
   if (abs(valp2) < epsilon) return P2;   

   VEC P = P1 + 0.5*(P2-P1);
   if (abs(valp1-valp2) < epsilon*10e3) return P;

   REAL val = density(P) - iso_value;
   if (valp1*val < 0) return dichotomy(P1, P, valp1, val);
   if (valp2*val < 0) return dichotomy(P, P2, val, valp2);
   
   return VEC();
}
*/


VEC ImplicitPrimitive::dichotomy(const VEC &P1, const VEC &P2, REAL valp1, REAL valp2) const{

  VEC P(3), Pt1=P1, Pt2=P2;
  REAL val;
  const double EPS = 10e-5;

  while(fabs(valp1-valp2)>EPS){
    P=Pt1+0.5*(Pt2-Pt1);
    val=density(P)-iso_value;
    if(val<0.0){
      Pt2=P;
      valp2=val;
    }
    else{
      valp1=val;
      Pt1=P;
    }
  }
  return P;

}



void ImplicitPrimitive::findIntersection(const VEC &P1, const VEC &P2, REAL valp1, REAL valp2, 
                                    VEC &P, VEC &N) const
{
   //   Linearly interpolate the position where an isosurface cuts
   //   an edge between P1 and P2, each with their own scalar value
   if (fabs(valp1) < epsilon) P = P1;
   else if (fabs(valp2) < epsilon) P = P2;
   else if (fabs(valp1-valp2) < epsilon) P = P1;
   else 
   {
      REAL mu = fabs(valp1 / (valp2 - valp1));
      P[0] = P1[0] + mu*(P2[0]-P1[0]);
      P[1] = P1[1] + mu*(P2[1]-P1[1]);
      P[2] = P1[2] + mu*(P2[2]-P1[2]);
   }

   N = getNormal(P);
   N.normalize();
}
//}

---