Triangle.cpp

Go to the documentation of this file.

#include <Triangle.h>


Triangle::Triangle(){}
Triangle::Triangle(const V_3D& x0,const V_3D& x1,const V_3D& x2)
{x[0]=x0;x[1]=x1;x[2]=x2;}
Triangle::Triangle(const Triangle& t)
{for(int k_dim=0;k_dim<3;k_dim++){x[k_dim]=t.x[k_dim];}}
Triangle::~Triangle(){destroy();}

int Triangle::destroy(){for(int k_dim=0;k_dim<3;k_dim++)x[k_dim].set(0.0,0.0,0.0);return 0;}

V_3D  Triangle::operator()(int index) const
{
  if(index<0 || index>2){printf("Error in operator(%d) in Triangle, index must be in [0,3]\n",index);exit(-1);}
  return x[index];
}
V_3D& Triangle::operator()(int index)
{
  if(index<0 || index>2){printf("Error in operator(%d) in Triangle, index must be in [0,3]\n",index);exit(-1);}
  return x[index];
}
V_3D  Triangle::operator[](int index) const
{
  if(index<0 || index>2){printf("Error in operator[%d] in Triangle, index must be in [0,3]\n",index);exit(-1);}
  return x[index];
}
V_3D  Triangle::operator[](int index)
{
  if(index<0 || index>2){printf("Error in operator[%d] in Triangle, index must be in [0,3]\n",index);exit(-1);}
  return x[index];
}
Triangle& Triangle::operator=(const Triangle& t)
{for(int k_dim=0;k_dim<3;k_dim++){x[k_dim]=t.x[k_dim];} return *this;}


int Triangle::barycentric_coordinates(const V_3D& _x,double *alpha,double *beta,double *gamma) const
{
  // check if _x is in the plane
  double epsilon=0.0001;
  double L = (_x-x[0]).dot(normal());
  if(L>epsilon)//not in the plane
    return -1;


  
  // calculate the barycentric coordinates
  double area_1=0.0,area_2=0.0,area_3=0.0,area_tot=0.0;
  area_1 = (x[1]-x[0]).area(_x-x[0]);
  area_2 = (x[2]-x[1]).area(_x-x[1]);
  area_3 = (x[0]-x[2]).area(_x-x[2]);
  area_tot = area();

  *alpha = area_1/area_tot;
  *beta  = area_2/area_tot;
  *gamma = area_3/area_tot;

  return 0;
}

int Triangle::is_vertex_inside(const V_3D _x) const
{
  double alpha=-1.0,beta=-1.0,gamma=-1.0;
  if(barycentric_coordinates(_x,&alpha,&beta,&gamma)==-1) //not in plane
    return 0;

  double epsilon=0.0001;
  if(alpha>1.0 || beta>1.0 || gamma>1.0 || (alpha+beta+gamma-1.0)>epsilon) // outside the polygon
    return 0;

  // in the polygon
  return 1;
}

V_3D Triangle::normal() const
{return (x[2]-x[0]).vector_prod(x[1]-x[0]).normalized();}

double Triangle::area() const
{return (x[2]-x[0]).area(x[1]-x[0]);}

V_3D Triangle::closest_point(const V_3D& _x,int *type) const
{
  double epsilon=0.00001;

  // First project and look if the projection is inside the triangle (barycentric coordinates)
  V_3D n=normal();
  V_3D proj = _x-((_x-x[0]).dot(n))*n;

  // if the proj is inside the triangle, this is the closest point
  if(is_vertex_inside(proj)==1)
    {*type=0;return proj;}

      
  // if not, take the closest line
  *type=1;
  Segment s0=get_segment(0),s1=get_segment(1),s2=get_segment(2);
  V_3D y0=s0.closest_point(_x),y1=s1.closest_point(_x),y2=s2.closest_point(_x);
  double dist_1=(y0-_x).norm(),dist_2=(y1-_x).norm(),dist_3=(y2-_x).norm();
  if(dist_1<=dist_2 && dist_1<=dist_3)
    {
      //check if its a vertex point
      for(int k=0;k<3;k++){if((y0-x[k]).norm()<epsilon) *type=2;}
      return y0;
    }
  else if(dist_2<=dist_1 && dist_2<=dist_3)
    {
      for(int k=0;k<3;k++){if((y1-x[k]).norm()<epsilon) *type=2;}
      return y1;
    }
  else 
    {
      for(int k=0;k<3;k++){if((y2-x[k]).norm()<epsilon) *type=2;}
      return y2;
    }

}

Segment Triangle::get_segment(int k_edge) const
{
  if(k_edge<0 || k_edge>2)
    {printf("Error in get_segment(%d) in Triangle, k_edge has to be between 0 and 2\n",k_edge);exit(-1);}

  Segment s(x[k_edge],x[(k_edge+1)%3]);
  return s;
}

std::vector <V_3D> Triangle::plane_intersection(const V_3D& n,const V_3D& x0,int *type) const
{int type_0=1,type_1=-1,type_2=1;return plane_intersection(n,x0,type,&type_0,&type_1,&type_2);}
std::vector <V_3D> Triangle::plane_intersection(const V_3D& n,const V_3D& x0,int *type,int *type_0,int *type_1,int *type_2) const
{
  Segment s0=get_segment(0),s1=get_segment(1),s2=get_segment(2);
  *type_0=-1,*type_1=-1,*type_2=-1;
  V_3D inter_0(-1.0,-1.0,-1.0),inter_1(-1.0,-1.0,-1.0),inter_2(-1.0,-1.0,-1.0);

  // take the intersection with the three segments of the polygon
  inter_0=s0.plane_intersection(n,x0,type_0);
  inter_1=s1.plane_intersection(n,x0,type_1);
  inter_2=s2.plane_intersection(n,x0,type_2);

  std::vector <V_3D> intersections;

  //1 - no intersection
  if(*type_0==0 && *type_1==0 && *type_2==0)
    {*type=0;return intersections;}


  //2 - parrallel
  if(*type_0==3){*type=3;intersections.push_back(x[0]),intersections.push_back(x[1]);return intersections;}
  if(*type_1==3){*type=3;intersections.push_back(x[1]),intersections.push_back(x[2]);return intersections;}
  if(*type_2==3){*type=3;intersections.push_back(x[2]),intersections.push_back(x[0]);return intersections;}

  //3 - vertex intersection (and only with a vertex)
  if(*type_0==2 && *type_1!=1 && *type_2!=1){*type=2;intersections.push_back(x[0]);return intersections;}
  if(*type_1==2 && *type_0!=1 && *type_2!=1){*type=2;intersections.push_back(x[1]);return intersections;}
  if(*type_2==2 && *type_0!=1 && *type_1!=1){*type=2;intersections.push_back(x[2]);return intersections;}

  //4 - two distinct points
  if( (*type_0==1 || *type_0==2) && (*type_1==1 || *type_1==2) )
    {*type=1;intersections.push_back(inter_0);intersections.push_back(inter_1);return intersections;}
  if( (*type_0==1 || *type_0==2) && (*type_2==1 || *type_2==2) )
    {*type=1;intersections.push_back(inter_0);intersections.push_back(inter_2);return intersections;}
  if( (*type_1==1 || *type_1==2) && (*type_2==1 || *type_2==2) )
    {*type=1;intersections.push_back(inter_1);intersections.push_back(inter_2);return intersections;}

  
  // else unknown configuration
  cout<<"Error in plane_intersection("<<n<<","<<x0<<") ";
  printf("Unknown configuration type=(%d,%d,%d)\n",*type_0,*type_1,*type_2);
  exit(-1);
}