V_3D.cpp

Go to the documentation of this file.

#include <V_3D.h>




V_3D::V_3D()
{
  for(int k_dim=0;k_dim<3;k_dim++)
    V[k_dim]=0;
}

V_3D::~V_3D()
{}

V_3D::V_3D(const V_3D &v_3d)
{
  for(int k_dim=0;k_dim<3;k_dim++)
    V[k_dim] = v_3d.V[k_dim];
}

V_3D::V_3D(double x,double y,double z)
{
  V[0]=x;
  V[1]=y;
  V[2]=z;
}

int V_3D::set(double *_v)
{
  for(int k_dim=0;k_dim<3;k_dim++)
    V[k_dim]=_v[k_dim];
  return 0;
}
int V_3D::set(int k_dim,double value)
{
  V[k_dim]=value;
  return 0;
}

int V_3D::set(V_3D _V)
{
  for(int k_dim=0;k_dim<3;k_dim++)
    V[k_dim]=_V.get(k_dim);
  return 0;
}

int V_3D::set(double x,double y,double z)
{
  V[0] = x;
  V[1] = y;
  V[2] = z;
  return 0;
}
  
double V_3D::get(int k_dim) const
{
  return V[k_dim];
}

int V_3D::get(double *_v) const
{
  for(int k_dim=0;k_dim<3;k_dim++)
    _v[k_dim]=V[k_dim];
  return 0;
}

V_3D V_3D::get() const
{
  return *this;
}

double *V_3D::get_pointer()
{
  return V;
}

V_3D operator+(const V_3D& v1,const V_3D& v2)
{
  V_3D Z;
  int k_dim=0;
  for(k_dim=0;k_dim<3;k_dim++)
    Z.set(k_dim,v1.V[k_dim]+v2.V[k_dim]);
  return Z;
}

V_3D operator-(const V_3D& v1,const V_3D& v2)
{
  V_3D Z;
  int k_dim=0;
  for(k_dim=0;k_dim<3;k_dim++)
    Z.set(k_dim,v1.V[k_dim]-v2.V[k_dim]);
  return Z;
}

V_3D operator*(const double& alpha,const V_3D& v1)
{
  V_3D Z;
  int k_dim=0;
  for(k_dim=0;k_dim<3;k_dim++)
    Z.set(k_dim,alpha*v1.V[k_dim]);
  return Z;
}

V_3D operator*(const V_3D& v1,const double& alpha)
{
  V_3D Z;
  int k_dim=0;
  for(k_dim=0;k_dim<3;k_dim++)
    Z.set(k_dim,alpha*v1.V[k_dim]);
  return Z;
}


V_3D operator/(const V_3D& v1,const double& alpha)
{
  V_3D Z;
  int k_dim=0;
  for(k_dim=0;k_dim<3;k_dim++)
    Z.set(k_dim,v1.V[k_dim]/alpha);
  return Z;
}


double V_3D::dot(V_3D _v) const
{
  double res=0.0;
  int k_dim=0;
  for(k_dim=0;k_dim<3;k_dim++)
    res += V[k_dim]*_v.get(k_dim);
  return res;
}

double V_3D::norm() const
{
  double res=0.0;
  int k_dim=0;
  for(k_dim=0;k_dim<3;k_dim++)
    res += V[k_dim]*V[k_dim];
  
  return powf(res,0.5);
}

V_3D& V_3D::operator=(const V_3D& _v)
{
  for(int k_dim=0;k_dim<3;k_dim++)
    V[k_dim] = _v.V[k_dim];
  return *this;
}

V_3D V_3D::vector_prod(V_3D _v) const
{
  V_3D res;

  res.set(0,V[1]*_v.get(2)-V[2]*_v.get(1));
  res.set(1,V[2]*_v.get(0)-V[0]*_v.get(2));
  res.set(2,V[0]*_v.get(1)-V[1]*_v.get(0));
  
  return res;
}

double V_3D::area(V_3D _v) const
{
  return vector_prod(_v).norm()/2;
}


ostream& operator << (ostream& flux, V_3D _v)
{
  flux<<"("<<_v.get(0)<<","<<_v.get(1)<<","<<_v.get(2)<<")";
  return flux;
}

bool operator==(const V_3D& v1,const V_3D& v2)
{
  double epsilon=0.0001;
  if(fabs(v1.V[0]-v2.V[0])<epsilon
     && fabs(v1.V[1]-v2.V[1])<epsilon
     && fabs(v1.V[2]-v2.V[2])<epsilon)
    return true;
  else
    return false;
}

bool operator!=(const V_3D& v1,const V_3D& v2)
{
  double epsilon=0.0001;
  if(fabs(v1.V[0]-v2.V[0])>=epsilon
     || fabs(v1.V[1]-v2.V[1])>=epsilon
     || fabs(v1.V[2]-v2.V[2])>=epsilon)
    return true;
  else
    return false;
}



double& V_3D::operator[](int index)
{
  if(index<0 || index>2)
    {printf("error in operator [] in V_3d, index [%d] not correct\n",index);exit(-1);}
  return V[index];
}
double V_3D::operator[](int index) const
{
  if(index<0 || index>2)
    {printf("error in operator [] in V_3d, index [%d] not correct\n",index);exit(-1);}
  return V[index];
}

V_3D& V_3D::operator+=(const V_3D& _v)
{
  V[0] += _v[0];
  V[1] += _v[1];
  V[2] += _v[2];
  return *this;
}
V_3D& V_3D::operator/=(const double& alpha)
{
  for(int k_dim=0;k_dim<3;k_dim++)
    V[k_dim] /= alpha;
  return *this;
}

V_3D V_3D::get_min_scalar(const V_3D& _v) const
{
  V_3D temp;
  for(int k_dim=0;k_dim<3;k_dim++)
    temp[k_dim]=(V[k_dim]<_v[k_dim]?V[k_dim]:_v[k_dim]);
  return temp;
}
V_3D V_3D::get_max_scalar(const V_3D& _v) const
{
  V_3D temp;
  for(int k_dim=0;k_dim<3;k_dim++)
    temp[k_dim]=(V[k_dim]>_v[k_dim]?V[k_dim]:_v[k_dim]);
  return temp;
}

double V_3D::get_max_of_coeff() const
{return max(max(V[0],V[1]),V[2]);}
double V_3D::get_min_of_coeff() const
{return min(min(V[0],V[1]),V[2]);}


int V_3D::equal(const V_3D& v,double epsilon) const
{
  if(powf((V[0]-v[0])*(V[0]-v[0])+(V[1]-v[1])*(V[1]-v[1])+(V[2]-v[2])*(V[2]-v[2]),0.5)<=epsilon)
    return 1;
  else
    return 0;
}

bool V_3D::operator>=(const V_3D &v) const
{
  if(V[0]>=v[0] && V[1]>=v[1] && V[2]>=v[2])
    return 1;
  return 0;
}
bool V_3D::operator>(const V_3D &v) const
{
  if(V[0]>v[0] && V[1]>v[1] && V[2]>v[2])
    return 1;
  return 0;
}
bool V_3D::operator<=(const V_3D &v) const
{
  if(V[0]<=v[0] && V[1]<=v[1] && V[2]<=v[2])
    return 1;
  return 0;
}
bool V_3D::operator<(const V_3D &v) const
{
  if(V[0]<v[0] && V[1]<v[1] && V[2]<v[2])
    return 1;
  return 0;
}

double V_3D::tetrahedra_volume(const V_3D& p2,const V_3D& p3,const V_3D& p4) const
{
  //compute
  // V = 1/2*det(A,B,C,D) = 1/2* (ABxAC).AD

  V_3D AB=p2-*this;
  V_3D AC=p3-*this;
  V_3D AD=p4-*this;

  double Vol=0.0;
  Vol = 0.5 * (AB.vector_prod(AC)).dot(AD);
  return Vol;
}

int V_3D::get_barycenter_coordinates_tetrahedra(const V_3D& p1,const V_3D& p2,const V_3D& p3,const V_3D& p4,double *alpha,double *beta,double *gamma,double *delta) const
{
  double Va=(*this).tetrahedra_volume(p2,p3,p4);
  double Vb=p1.tetrahedra_volume(*this,p3,p4);
  double Vc=p1.tetrahedra_volume(p2,*this,p4);
  double Vd=p1.tetrahedra_volume(p2,p3,*this);
  double Vtot=p1.tetrahedra_volume(p2,p3,p4);

  *alpha=Va/Vtot;
  *beta =Vb/Vtot;
  *gamma=Vc/Vtot;
  *delta=Vd/Vtot;

  return 0;
}

int V_3D::is_inside_volume(const V_3D& p1,const V_3D& p2,const V_3D& p3,const V_3D& p4) const
{
  double alpha=0.0,beta=0.0,gamma=0.0,delta=0.0;
  get_barycenter_coordinates_tetrahedra(p1,p2,p3,p4,&alpha,&beta,&gamma,&delta);

  if(alpha>0 && beta>0 && gamma>0 && delta>0)
    return 1;
  return 0;
}

double V_3D::distance_to_oriented_plane(const V_3D& n,const V_3D& x0) const
{
  //normalize
  double epsilon=0.000001;
  double norm_n=n.norm();
  if(norm_n<epsilon)
    {printf("Error n is null in distance to oriented_plane in V_3D...\n");exit(-1);}
  V_3D n_normalized=n/norm_n;

  double proj = ((*this)-x0).dot(n_normalized);
  return proj;
}

double V_3D::cos_angle(const V_3D& v) const
{
  V_3D current = (*this)/(*this).norm();
  V_3D vn = v/v.norm();

  return current.dot(vn);
}

double V_3D::angle(const V_3D& v)const
{return acos(cos_angle(v));}

double V_3D::cos_angle_to_plane_with_normal(const V_3D& normal_to_plane) const
{
  double norm_normal=normal_to_plane.norm();
  if(norm_normal<=0.000001)
    {printf("Error in angle_to_plane_with_normal in V_3D, normal has norm=0\n");exit(-1);}
  double current_norm=norm();
  V_3D this_normed = (*this)/current_norm;
  V_3D projected = this_normed - (this_normed.dot(normal_to_plane))*normal_to_plane;
  double projected_norm= this_normed.norm();
  if(projected_norm<0.00001)
    {printf("Error in angle_to_plane_with_normal in V_3D, projected_normal has norm=0\n");exit(-1);}
  projected = projected/projected_norm;
  double dot_prod=this_normed.dot(projected);

  return dot_prod;
}

V_3D V_3D::normalized() const
{
  double epsilon=0.00001;
  double current_norm=norm();
   if(current_norm<epsilon)
     {
       //printf("Error norm is zero in [V_3D].normalized() function\n");
       return V_3D(0,0,1);
     }
  return (*this)/current_norm;
}