Biomechanical Joint Model  Author: Anderson Maciel

---

blob.cpp

Go to the documentation of this file.

#include <math.h>
#include "blob.h"

#define max(a,b)  (((a) > (b)) ? (a) : (b))
#define min(a,b)  (((a) < (b)) ? (a) : (b))


//namespace MarchingCubes {

#define SQR(x)       ((x) * (x))

static bool ray_intersect_primitive(Ray &ray, const Quadric &quadric, REAL root[], REAL c){

        REAL k2,k1,k0;
        REAL s,t,t1,t2;

        // map ray into model local coordinate system
        VEC X0( quadric.T.fullMult(ray.O) );
        VEC R( quadric.T.mult(ray.D) );

        k2 = c * (SQR(R[0]) / quadric.a2 + 
                  SQR(R[1]) / quadric.b2 + 
                  SQR(R[2]) / quadric.c2);
        k1 = 2*c * (R[0]* X0[0] / quadric.a2 +
                    R[1]* X0[1] / quadric.b2 +
                    R[2]* X0[2] / quadric.c2);
        k0 = c * (SQR(X0[0]) / quadric.a2 +
                  SQR(X0[1]) / quadric.b2 +
                  SQR(X0[2]) / quadric.c2) - 1.;
                
        s = SQR(k1) - 4.*k2*k0;
        root[0] = root[1] = -1.0;
        if (s < 0.) return(false);
        s = sqrt(s);
        t1 = 0.5*(-k1 + s)/k2;
        t2 = 0.5*(-k1 - s)/k2;
        t = max(t1,t2);
        root[0] = min(t1,t2);
        root[1] = t;
        if  (t <0.) return (false);
        ray.t = t;
        return(true);
}

bool Quadric::intersectedBy(Ray &ray, REAL root[]) const{ 
   return ray_intersect_primitive(ray, *this, root, 2.);
}

bool Blob::intersectedBy(Ray &ray, REAL root[]) const{ 
  return ray_intersect_primitive(ray, *this, root, 1.);
}

//                 ELLIPSOIDS                      //
bool Ellipsoids::intersectedBy(Ray &ray) const{
   int n = quadlist.size();
   if ( n <=0 ) return false;

   REAL **interval = new REAL* [n]; // t1 and t2  intersection of ray with quadric 
   for (int i=0; i<n; i++) interval[i] = new REAL [2];

   // 1: find the blobs which have radius of influence to the ray 
   for ( int j=0; j<n; ++j)
      quadlist[j].intersectedBy(ray, interval[j]);
    
    // 2 : find the minimum t
    REAL t_min=infinity;
    for ( int k=0;k<n;k++) 
    {
       if (interval[k][0]>=0  && interval[k][0]<t_min)
          t_min = interval[k][0];

       if (interval[k][1]>=0 && interval[k][1]<t_min)
          t_min = interval[k][1];
    }
    
    ray.t = t_min;
    ray.P = ray.O + t_min*ray.D;
    delete [] interval;
    
    return (ray.segment==true)? (t_min >=0 && t_min <=1) : (t_min >=0 && t_min < infinity);
}

REAL Ellipsoids::getDistance(const VEC &U) const{   

   REAL k, t=-infinity;
   bool outside = true;
   static VEC X(3);

   for (vector<Quadric>::const_iterator quad=quadlist.begin(); quad!=quadlist.end(); ++quad)
   {
      // map point into quadric lcs
      X = quad->T.fullMult(U);
      k =  0.5 - (X[0]*X[0]/quad->a2 + X[1]*X[1]/quad->b2 + X[2]*X[2]/quad->c2);

      // right on the surface?
      if (fabs(k) < epsilon) return 0;

      if (outside && k<0)
      {
         if (k>t) t=k;
      }

      if (k>0) 
      {  
         if (outside == true)
         {          
            outside = false;
            t = k;
         }
         else if (k<t) t=k;
      }
   }
   return t;   
}

BBox Ellipsoids::getBoundingBox()
{
   bbox.makeEmpty();
   VEC X(3), Y(3);

   for (vector<Quadric>::const_iterator quad_it = quadlist.begin(); 
        quad_it != quadlist.end(); ++quad_it)
   {
      BBox bb;
      X[0] = quad_it->a; X[1] = quad_it->b; X[2] = quad_it->c;

      // for each vertex of BBox
      for (int p=0; p<8; p++)
      {
         int sign[3]; 
         sign[0] = p&4; 
         sign[1] = p&2;
         sign[2] = p&1;    
 
         for (int k=0; k<3; k++) 
            Y[k] = (sign[k]==0)? X[k] : -X[k];

         Y = quad_it->M.fullMult(Y);  // point mult
            
         // Y now contains world coordinates of BBox's vertex
         bb.extendBy(Y);
      }
      bbox.extendBy(bb);   
   }
//    bbox.scale(1.05);

   return bbox; 
}

void Ellipsoids::findIntersection(const VEC &P1, const VEC &P2, REAL valp1, REAL valp2, VEC &P, VEC &N) const{
   Ray ray(P1, P2-P1);
   ray.segment = true;

   if (!intersectedBy(ray))
   {
      //cout << "Fatal error: no intersection in Ellipsoids::findIntersection\n";
           printf( "Fatal error: no intersection in Ellipsoids::findIntersection\n" );
      return;
   }
   P = ray.P;

   static VEC X(3);
   for (vector<Quadric>::const_iterator quad=quadlist.begin(); quad!=quadlist.end(); ++quad)
   {
      // map point into quadric lcs
      X = quad->T.fullMult(P);
      REAL k =  0.5 - (X[0]*X[0]/quad->a2 + X[1]*X[1]/quad->b2 + X[2]*X[2]/quad->c2);
      if (fabs(k) < epsilon)
      {
         N = P - quad->M.getTranslation();
         N.normalize();
         break;
      }
   }
   return;
}

//                 METABALLS                       //
BBox Metaballs::getBoundingBox()
{
   bbox.makeEmpty();
   VEC X(3), Y(3);

   for (vector<Blob>::const_iterator blob_it = bloblist.begin(); 
        blob_it != bloblist.end(); ++blob_it)
   {
      BBox bb;

      if (blob_it->w>0)
      {
         X[0] = blob_it->a; X[1] = blob_it->b; X[2] = blob_it->c;
         X *= blob_it->w;

         // for each vertex of BBox
         for (int p=0; p<8; p++)
         {
            int sign[3]; 
            sign[0] = p&4; 
            sign[1] = p&2;
            sign[2] = p&1;         
 
            for (int k=0; k<3; k++) 
               Y[k] = (sign[k]==0)? X[k] : -X[k];

            Y = blob_it->M.fullMult(Y);  // point mult
            
            // Y now contains world coordinates of BBox's vertex
            bb.extendBy(Y);
         }
         bbox.extendBy(bb);
      }
   }
   
//    bbox.scale(1.05);
   return bbox;
}



// comparison routine for quicksort of stdlib
int mycompare(const void *aP, const void *bP) 
{ 
  REAL a=*((float*)aP);
  REAL b=*((float*)bP);

  if (a<b)       return -1;
  else if (a==b) return  0;
  else           return  1;
} 
int (*mycompareP)(const void *, const void *) = mycompare;

#define LESS(a,b)    ((a) <((b)+epsilon))
#define GREAT(a,b)   ((a) >((b)-epsilon))

//
// ray tracing MUST be used for Shen's density function as it is zero
// after a certain distance threshold.
//
// Substitute the ray equation into the blob quadric equations,thus t is broken
// into intervals where the ellipsoids of various radii intersect the ray. 
// We are interested only int the outmost intersection point. Suppose 
// [tmin,tmax] is the interval of outmost quadric, we only test the quadrics
// whose interval overlaps this outmost interval. The intervals of those blobs
// are superimposed and broken into new intervals with [tmin,tmax] each of 
// which has a sum of quadratic densities associated with it. We test 
// sebsequently from right to left of new intervals until we find a solution.
// ##Update 1(July 15, 94): k2,k1,k0, s change from float to double type.
// first compute the only the outmost metaball and record its t value.
// then compute additive t value. If the additive t-value < single t-value,
// return the single t-value;
// 
bool 
Metaballs_Shen::intersectedBy(Ray &ray) const
{
    int n = bloblist.size();
    if (n<=0) return false;
    REAL    s, t, t1, t2, tmin, tmax;

    REAL *t_value = new REAL [2*n];
    REAL **interval = new REAL* [n]; // t1 and t2  intersction of ray with quardic 
    for (int i=0; i<n; i++) interval[i] = new REAL [2];
    
    REAL k2,k1,k0;
    int k;
    VEC R,X0; // equivalent ray base and direction vector in blob's lcs

    // 1: find the blobs which have radius of influence to the ray 
    for (int j=0; j<n; ++j){
       bloblist[j].intersectedBy(ray, interval[j]);
    }
    
    // 2 : find the maximum t
    k =0;
    for (int ii=0;ii< n;++ii) {
        t_value[k] = interval[ii][0];
        t_value[k+1] = interval[ii][1];
        k +=2;
    }
    qsort(t_value, k, sizeof(REAL), mycompareP);

    if (t_value[2*n-1] <0.) 
    { 
       delete [] t_value; 
       delete [] interval; 
       return(false); 
    } 
 
    bool found = false;
    ray.t = -1;

    for (k=1; k<2*n && t_value[k]<0; ++k);
    // 3: solve the equation f = iso_value in each interval 
    // we exit the loop as soon as an intersection point is found 
    // start with the smallest interval to ensure we retrieve the closest point 
    for (; k<2*n; ++k){
        k2 = k1 = k0 = 0.;
        tmin = t_value[k-1]; tmax = t_value[k];// loop through interval
//      if (tmax < 0.) break;  
        if (tmin <0.) tmin =0; // intersection point must be along direction
        if (fabs(tmax-tmin) < epsilon) continue;
        if (ray.segment==true && tmax >1.) break;

        // search the blobs that has density contribution in the interval
        int kk=0; 
        for (vector<Blob>::const_iterator blob=bloblist.begin();blob!=bloblist.end();blob++,++kk){
            if ((LESS(interval[kk][0],tmin))
                &&(GREAT(interval[kk][1],tmax)))
            {
                // map ray into blob lcs
                X0 = blob->T.fullMult(ray.O);
                R = blob->T.mult(ray.D);

                k2 += blob->w*(R[0]*R[0]/blob->a2 +
                              R[1]*R[1]/blob->b2 +R[2]*R[2]/blob->c2);
                k1 += 2.0*blob->w*(X0[0]*R[0]/blob->a2+
                                  X0[1]*R[1]/blob->b2 +X0[2]*R[2]/blob->c2);
                k0 +=  blob->w*(X0[0]*X0[0]/blob->a2 +
                               X0[1]*X0[1]/blob->b2+X0[2]*X0[2]/blob->c2)-blob->w;
            }
        } // end of for (i=...)
        k0 += iso_value;
        s = k1*k1 - 4.*k2*k0;
        if (s < 0.) continue;
        s = sqrt(s);
        t1 = 0.5*(-k1 + s)/k2;
        t2 = 0.5*(-k1 - s)/k2;
        t = max(t1,t2);
        if  ( GREAT(t,tmin) && LESS(t,tmax))
        {
           found = true;
           ray.t =t;
           break;
        }
        else {
            t = min(t1,t2);
            if  ( GREAT(t,tmin) && LESS(t,tmax))
            {
               found = true;
               ray.t =t;
               break;
            }
        }
    } // end of for (k= ...) 

    delete [] t_value; 
    delete [] interval;

    if (found) ray.P = ray.O + ray.t*ray.D;
    return found;
}


// return field value at location U 
REAL
Metaballs_Shen::density(const VEC &U) const
{
  REAL k, t=0;
  static VEC X(3);
  
  for (vector<Blob>::const_iterator blob=bloblist.begin(); blob!=bloblist.end(); ++blob)
  {
     // map point into blob lcs
     X = blob->T.fullMult(U);

     k =  1.0 - (X[0]*X[0]/blob->a2 + X[1]*X[1]/blob->b2 + X[2]*X[2]/blob->c2);
     if ( k>0 ) t += blob->w * k;
  }
  return t;
}

/*
// Normally we should use the following code to compute the gradient (which yields the normal)
// But Shen's density function is only C0, so discontinuities exist in the gradient
// As a result, normals will appear "wrong" at places where metaballs blend
VEC
Metaballs_Shen::gradient(const VEC &U) const
{
  REAL k;
  static VEC X(3);
  VEC W(3);

  for (vector<Blob>::const_iterator blob=bloblist.begin(); blob!=bloblist.end(); blob++)
  {
      X = blob->T.fullMult(U);
      k =  1 - (X[0]*X[0]/blob->a2 + X[1]*X[1]/blob->b2 + X[2]*X[2]/blob->c2);
      if ( k>0 )
      {
          W[0] += 2*blob->w*X[0]/blob->a2;
          W[1] += 2*blob->w*X[1]/blob->b2;
          W[2] += 2*blob->w*X[2]/blob->c2;
      }
  }  
  return W;
}
*/

// We do not use Shen's density function to compute the gradient as it's dicontinuous
// we use instead an exponential function (closely resembles Shen's density function)
// return blob's gradient(U) = normal at point U
VEC Metaballs_Shen::gradient(const VEC &U) const{
   REAL k;
   static VEC X(3), Y(3);
   VEC W(3);
   REAL coeff=2;

   for (vector<Blob>::const_iterator blob=bloblist.begin(); blob!=bloblist.end(); ++blob)
   {
      X = blob->T.fullMult(U);

      k = -2 * blob->w * coeff *
         exp(-coeff * (X[0]*X[0]/blob->a2 + X[1]*X[1]/blob->b2 + X[2]*X[2]/blob->c2));

      Y[0] = X[0]*k/blob->a2;
      Y[1] = X[1]*k/blob->b2;
      Y[2] = X[2]*k/blob->c2;      
      
      // map back vector into g.c.s.
      W += blob->M.mult(Y);
   }
   return W;
}


//  Metaballs with an exponential density function //
REAL Metaballs_exp::density(const VEC &U) const{
   REAL t=0;
   static VEC X(3);

   for (vector<Blob>::const_iterator blob=bloblist.begin(); blob!=bloblist.end(); ++blob)
   {
      // map point into blob lcs
      X = blob->T.fullMult(U);
      t += blob->w *  
         exp(-coeff * (X[0]*X[0]/blob->a2 + X[1]*X[1]/blob->b2 + X[2]*X[2]/blob->c2));
   }
   return t;
}


// return blob's gradient(U) = normal at point U
VEC Metaballs_exp::gradient(const VEC &U) const
{
   REAL k;
   static VEC X(3), Y(3);
   VEC W(3);

   for (vector<Blob>::const_iterator blob=bloblist.begin(); blob!=bloblist.end(); ++blob)
   {
      X = blob->T.fullMult(U);

      k = -2 * blob->w * coeff *
         exp(-coeff * (X[0]*X[0]/blob->a2 + X[1]*X[1]/blob->b2 + X[2]*X[2]/blob->c2));

      Y[0] = X[0]*k/blob->a2;
      Y[1] = X[1]*k/blob->b2;
      Y[2] = X[2]*k/blob->c2;      
      
      // map back vector into g.c.s.
      W += blob->M.mult(Y);
   }
   return W;
}


REAL  MetaballsSphere::density(const VEC &U) const{
  REAL r, den, t=0;
#if 0
  fprintf(stderr,"MetaballsSphere::density\n");
#endif
  for(vector<Sphere>::const_iterator sph=spherelist.begin(); sph!=spherelist.end(); sph++){
    r = (U - sph->C).norm();
    den=0;
    if(r<sph->R)
      if(r<sph->ro)
        den=sph->k*(sph->ro-r)+1;
      else
        den=0.25*((sph->k*(r-sph->ro))-2.0)*((sph->k*(r-sph->ro))-2.0);

    t+=den;
  }
  return t;
}

VEC   MetaballsSphere::gradient(const VEC& U) const{
 REAL r, aux;
 VEC grad(3), g(3);

 for(vector<Sphere>::const_iterator sph=spherelist.begin(); sph!=spherelist.end(); sph++){
   r = (U - sph->C).norm();
   aux = 0;
   if (r<sph->R){
     if (r<sph->ro){
       aux=-sph->k;
     }
     else{
       aux=sph->k*((sph->k*(r-sph->ro))-2.0)/2.0;
     }
   }

   grad=U-sph->C;
   grad*=aux/r;
   g+=grad;
 }
 return g;
}


BBox  MetaballsSphere::getBoundingBox(){
  VEC UR(3), LL(3);
  bbox.makeEmpty();
  for(vector<Sphere>::iterator sph=spherelist.begin(); sph!=spherelist.end(); ++sph){
     UR[0] = sph->C[0] + sph->R;
     UR[1] = sph->C[1] + sph->R;
     UR[2] = sph->C[2] + sph->R;

     LL[0] = sph->C[0] - sph->R;
     LL[1] = sph->C[1] - sph->R;
     LL[2] = sph->C[2] - sph->R; 
     
     bbox.extendBy(BBox(LL,UR));
  }

  return bbox;

}


REAL  MetaballsSphere_exp::density(const VEC &U) const{
  REAL r, F;
#if 0
  fprintf(stderr,"MetaballsSphere_exp::density\n");
#endif
  F=0;
  for(vector<Sphere>::const_iterator sph=spherelist.begin(); sph!=spherelist.end(); ++sph){
    r = (U - sph->C).norm();
    F += exp(-coeff*r);
  }
  return F;
}

VEC   MetaballsSphere_exp::gradient(const VEC& U) const{
 VEC df_da(3), dF_da(3), r;

 dF_da[0]=dF_da[1]=dF_da[2]=0.0;

 for(vector<Sphere>::const_iterator sph=spherelist.begin(); sph!=spherelist.end(); ++sph){
   r = (U - sph->C);
   df_da= -2*coeff*exp(-coeff*r.norm())*r;
   dF_da+=df_da;
 }
 return dF_da;
}

double MetaballsTriangle::EuclidDistFromTriangle(const VEC &P1, const VEC &P2, const VEC &P3,
                                                 const VEC &Pt, const double *Plane, VEC &N) const{

  VEC na(3), nb(3), nc(3); // normals of the planes pependicular to the Facet and passing throught the Facet edges
  VEC a(3), b(3), c(3);
  double Pt_aplane, Pt_bplane, Pt_cplane, dist=0.0, D;
  VEC n(3);
  n[0]=Plane[0], n[1]=Plane[1], n[2]=Plane[2], D=Plane[3];

  a = P3-P2;
  b = P1-P2;
  c = P3-P1;
  
  na =  cross_prod(n,a);
  nb = -cross_prod(n,b);
  nc = -cross_prod(n,c);

  // chech where is the point
  Pt_aplane = dot_prod(Pt-P2,na);
  Pt_bplane = dot_prod(Pt-P2,nb);
  Pt_cplane = dot_prod(Pt-P1,nc);

#if 0
  cout<<"a="<<a<<endl;
  cout<<"b="<<b<<endl;
  cout<<"c="<<c<<endl;
  
  cout<<"na="<<na<<endl;
  cout<<"nb="<<nb<<endl;
  cout<<"nc="<<nc<<endl;  

  cout<<"Pt in respect to plane na ="<<Pt_aplane<<endl;
  cout<<"Pt in respect to plane na ="<<Pt_bplane<<endl;
  cout<<"Pt in respect to plane nc ="<<Pt_cplane<<endl;
#endif

  if(Pt_aplane>0.0 && Pt_bplane>0.0 && Pt_cplane>0.0){ // point projects on the Facet
    dist = (Pt[0]*n[0] + Pt[1]*n[1] + Pt[2]*n[2] - D)/n.norm();
        std::cout<<"Zone 1"<<endl;
  }
  else if(Pt_aplane<0.0 && Pt_bplane>0.0 && Pt_cplane>0.0){ // point project on the side of P2P3 edge
    dist = cross_prod((Pt-P2),a).norm()/a.norm();
    VEC Pt_prj = Pt + (D-dot_prod(Pt,n)/n.norm2())*n;
    double q = (Pt_prj-P2).norm(), cos_alfa = dot_prod(Pt_prj-P2,a)/((Pt_prj-P2).norm()*a.norm());
    double p = q*cos_alfa;
    VEC P2P = p*a/a.norm(), N = (Pt-P2) - P2P;
    std::cout<<"Zone 2"<<endl;
  }
  else if(Pt_aplane>0.0 && Pt_bplane>0.0 && Pt_cplane<0.0){ // point project on the side of P1P3 edge
    dist = cross_prod((Pt-P1),c).norm()/c.norm();
    VEC Pt_prj = Pt + (D-dot_prod(Pt,n)/n.norm2())*n;
    double q = (Pt_prj-P1).norm(), cos_alfa = dot_prod(Pt_prj-P1,c)/((Pt_prj-P1).norm()*c.norm());
    double p = q*cos_alfa;
    VEC P1P = p*c/c.norm(), N = (Pt-P1) - P1P;
    std::cout<<"Zone 3"<<endl;
  }
  else if(Pt_aplane>0.0 && Pt_bplane<0.0 && Pt_cplane>0.0){ // point project on the side of P1P2 edge
    dist = cross_prod((Pt-P2),b).norm()/b.norm();
    VEC Pt_prj = Pt + (D-dot_prod(Pt,n)/n.norm2())*n;
    double q = (Pt_prj-P2).norm(), cos_alfa = dot_prod(Pt_prj-P2,b)/((Pt_prj-P2).norm()*b.norm());
    double p = q*cos_alfa;
    VEC P2P = p*b/b.norm(), N = (Pt-P2) - P2P;
    std::cout<<"Zone 4"<<endl;
  }
  else if(Pt_aplane<0.0 && Pt_bplane<0.0 && Pt_cplane>0.0){ // point projects on the side of P2 point
    dist = (Pt-P2).norm();
    std::cout<<"Zone 5"<<endl;
  }
  else if(Pt_aplane<0.0 && Pt_bplane>0.0 && Pt_cplane<0.0){ // point projects on the side of P3 point
    dist = (Pt-P3).norm();
    std::cout<<"Zone 6"<<endl;
  }
  else if(Pt_aplane>0.0 && Pt_bplane<0.0 && Pt_cplane<0.0){ // point projects on the side of P1 point
    dist = (Pt-P1).norm();
    std::cout<<"Zone 7"<<endl;
  } 
  return dist;
}

REAL  MetaballsTriangle::density(const VEC &U) const{
  REAL r, den, t=0;
  VEC N(3);
#if 0
  fprintf(stderr,"MetaballsTriangle::density\n");
#endif
   for(vector<TrianglePrimitive>::const_iterator trg=trianglelist.begin(); trg!=trianglelist.end(); ++trg){
    r = EuclidDistFromTriangle(trg->P1, trg->P2, trg->P3, U, trg->n, N);
    den=0;
    if(r<trg->R)
      if(r<trg->ro)
        den=trg->k*(trg->ro-r)+1;
      else
        den=0.25*((trg->k*(r-trg->ro))-2.0)*((trg->k*(r-trg->ro))-2.0);
    t+=den;
  }
  return t;
}

VEC   MetaballsTriangle::gradient(const VEC& U) const{
 REAL r, aux;
 VEC grad(3), g(3);

 for(vector<TrianglePrimitive>::const_iterator trg=trianglelist.begin(); trg!=trianglelist.end(); ++trg){
   r = EuclidDistFromTriangle(trg->P1, trg->P2, trg->P3, U, trg->n, grad);
   aux = 0;
   if (r<trg->R){
     if (r<trg->ro){
       aux=-trg->k;
     }
     else{
       aux=trg->k*((trg->k*(r-trg->ro))-2.0)/2.0;
     }
   }
   grad*=aux/r;
   g+=grad;
 }
 return g;
}

BBox  MetaballsTriangle::getBoundingBox(){
  VEC UR(3), LL(3);
  int i=0;
  bbox.makeEmpty();
  for(vector<TrianglePrimitive>::iterator trg=trianglelist.begin(); trg!=trianglelist.end(); ++trg, i++){
#if 0
    cout<<"Facet "<<i<<"  :["<<trg->P1[0]<<","<<trg->P1[1]<<","<<trg->P1[2]<<"]"<<endl;
    cout<<"         :["<<trg->P2[0]<<","<<trg->P2[1]<<","<<trg->P2[2]<<"]"<<endl;
    cout<<"         :["<<trg->P3[0]<<","<<trg->P3[1]<<","<<trg->P3[2]<<"]"<<endl;
#endif
    UR[0] = max(max(trg->P1[0],trg->P2[0]),trg->P3[0]) + trg->R;
    UR[1] = max(max(trg->P1[1],trg->P2[1]),trg->P3[1]) + trg->R;
    UR[2] = max(max(trg->P1[2],trg->P2[2]),trg->P3[2]) + trg->R;

    LL[0] = min(min(trg->P1[0],trg->P2[0]),trg->P3[0]) - trg->R;
    LL[1] = min(min(trg->P1[1],trg->P2[1]),trg->P3[1]) - trg->R;
    LL[2] = min(min(trg->P1[1],trg->P2[1]),trg->P3[1]) - trg->R;
     
    bbox.extendBy(BBox(LL,UR));
  }
#if 0
  cout<<"Bbox is:"<<endl;
  cout<<"LL :"<<LL<<endl;
  cout<<"UR :"<<UR<<endl;
#endif
  return bbox;
}


//} // end of the namespace 

---