A point within the tetrahedron

출처 : http://steve.hollasch.net/cgindex/geometry/ptintet.html

 

어떤 포인트가 tetrahedron(사면체)안에 있는지 확인하는 방법.

사면체를 이루는 정점 v1, v2, v3, v4과 입력된 포인트 v에 대해서 적절하게 4x4 determinant(행렬식)들을 구해 그 관계를 따져보면 된다.

 

4개의 정점에 대해서 행렬식은 세개의 정점이 이루는 평면을 기준으로 네 번째 점이 평면의 위에있는지, 아래에있는지, 동일평면상(coplanar)인지를 판단하는 잦대가 된다. 예를들어 네개의 정점에 대한 다음 행렬식은,

 

 

V0, V1, V2가 이루는 CCW(Counter-Clock Wise) front face 평면에 대하여 V3의 위치가 행렬식의 값이 음수라면 앞면, 양수라면 뒷면, 0이라면 동일평면상에 존재한다는것을 나타낸다.

(평면을 형성할 기준정점을 무엇으로 선택하냐에 따라 행렬식에 의한 앞, 뒤 위치판단은 달라질 수 있다. 예로 V1, V2, V3을 기준평면 정점으로 잡는다면 행렬식의 값이 양수일때 앞면, 음수일때 뒷면, 0일때 동일평면상을 나타낸다.)

이러한 사실을 바탕으로 5개의 정점(사면체 4개의 정점 + 테스트할 입력정점)에서 나올 수 있는 4가지 정점의 모든 조합에 대한 행렬식들을 구하며 이들의 관계에 따라 입력된 정점이 사면체와 어떤 위치관계를 갖는지를 파악하게된다.

( 행렬식 구하는 방법이 궁금한사람은 바로 전 포스트를 보면 된다. http://alleysark.tistory.com/entry/matrix-determinant  )

 

5가지 행렬식은 아래와 같다.

 

 

D0 == 0 이라면 사면체를 이루는 네 정점은 동일평면상(coplanar)에 존재하는것이다.

Di == 0 이라면 입력정점 V는 V와 Vi를 제외한 세 정점이 이루는 평면과 동일한 평면상에 존재한다.

  동일평면상에 존재하는 경우라면 세 점이 이루는 삼각형의 밖에있는지, 안에있는지, 모서리에 존재하는지 다시 확인해줘야한다.

  3차원 공간에서 삼각형과 포인트의 위치관계는 노말벡터를 이용했다 (행렬식을 이용할수있지 않을까..?)

어떤 Di 중 하나라도 D0와 부호가 다르면 V는 사면체 밖에 존재한다.

모든 Di 와 D0의 부호가 같다면 V는 사면체 않에 존재한다.

 

구현이야 별건 없지만..  일단은 이런식으로 사용하면 된다.

 

class Triangle{
public:
 ...

 enum OptionalConstant{
  INVALID_TRIANGLE=-1, OUTSIDE, INSIDE,
  EDGE_01, EDGE_02, EDGE_12
 };

 

 // ccw front normal calculation
 static Vector CalculateNormal(const Vertex& va, const Vertex& vb, const Vertex& vc){
  Vector vab = vb - va;
  Vector vac = vc - va;

  Vector norm = vab.Cross(vac);
  norm.Normalize();
  return norm;
 }

 

 static OptionalConstant PointOnFace(const Vertex& va, const Vertex& vb, const Vertex& vc, const Vertex& vp){
  Vector n1 = CalculateNormal(va, vb, vp);
  if( n1.MagnitudeSq() == 0 )
   return EDGE_01;

  Vector n2 = CalculateNormal(vb, vc, vp);
  if( n2.MagnitudeSq() == 0 )
   return EDGE_12;

  if( n1.Dot(n2) < 0 )
   return OUTSIDE;

  Vector n3 = CalculateNormal(vc, va, vp);
  if( n3.MagnitudeSq() == 0 )
   return EDGE_12;

  if( n1.Dot(n3) < 0 )
   return OUTSIDE;
  return INSIDE;
 }
};

 

 

class Tetrahedron{
public:
 const Vertex* v[4];

 Tetrahedron(){
  memset(v, NULL, sizeof(v));
 }

 Tetrahedron(const Vertex& va, const Vertex& vb, const Vertex& vc, const Vertex& vd){
  v[0] = &va; v[1] = &vb; v[2] = &vc; v[3] = &vd;
 }
 
 enum OptionalConstant{
  INVALID_TETRAHEDRON=-1, OUTSIDE, INSIDE,
  FACE_012, FACE_013, FACE_023, FACE_123,
  EDGE_01, EDGE_02, EDGE_03, EDGE_12, EDGE_13, EDGE_23,
  VTX_0, VTX_1, VTX_2, VTX_3
 };

 
 OptionalConstant IsWithin(const Vertex& vtx){
  for(unsigned i=0; i<4; i++)
   if( (*v[i]) == vtx )
    return OptionalConstant(VTX_0 + i);
  
  int d0sign = sign( det4x4(
   v[0]->x, v[0]->y, v[0]->z, 1,
   v[1]->x, v[1]->y, v[1]->z, 1,
   v[2]->x, v[2]->y, v[2]->z, 1,
   v[3]->x, v[3]->y, v[3]->z, 1 ) );
  if( d0sign == 0 )
   return INVALID_TETRAHEDRON; //all vertices of tetrahedron are coplanar

  int d1sign = sign( det4x4(
   vtx.x, vtx.y, vtx.z, 1,
   v[1]->x, v[1]->y, v[1]->z, 1,
   v[2]->x, v[2]->y, v[2]->z, 1,
   v[3]->x, v[3]->y, v[3]->z, 1 ) );
  if( d1sign == 0 ){
   Triangle::OptionalConstant onWhere = Triangle::PointOnFace(*v[1], *v[2], *v[3], vtx);
   switch(onWhere){
   case Triangle::EDGE_01:
    return EDGE_12;
   case Triangle::EDGE_02:
    return EDGE_13;
   case Triangle::EDGE_12:
    return EDGE_23;
   case Triangle::INSIDE:
    return FACE_123;
   }
   return OUTSIDE;
  }
  if( d0sign != d1sign )
   return OUTSIDE;

  int d2sign = sign( det4x4(
   v[0]->x, v[0]->y, v[0]->z, 1,
   vtx.x, vtx.y, vtx.z, 1,
   v[2]->x, v[2]->y, v[2]->z, 1,
   v[3]->x, v[3]->y, v[3]->z, 1 ) );
  if( d2sign == 0 ){
   Triangle::OptionalConstant onWhere = Triangle::PointOnFace(*v[0], *v[2], *v[3], vtx);
   switch(onWhere){
   case Triangle::EDGE_01:
    return EDGE_02;
   case Triangle::EDGE_02:
    return EDGE_03;
   case Triangle::EDGE_12:
    return EDGE_23;
   case Triangle::INSIDE:
    return FACE_023;
   }
   return OUTSIDE;
  }
  if( d0sign != d2sign )
   return OUTSIDE;

  int d3sign = sign( det4x4(
   v[0]->x, v[0]->y, v[0]->z, 1,
   v[1]->x, v[1]->y, v[1]->z, 1,
   vtx.x, vtx.y, vtx.z, 1,
   v[3]->x, v[3]->y, v[3]->z, 1 ) );
  if( d3sign == 0 ){
   Triangle::OptionalConstant onWhere = Triangle::PointOnFace(*v[0], *v[1], *v[3], vtx);
   switch(onWhere){
   case Triangle::EDGE_01:
    return EDGE_01;
   case Triangle::EDGE_02:
    return EDGE_03;
   case Triangle::EDGE_12:
    return EDGE_13;
   case Triangle::INSIDE:
    return FACE_013;
   }
   return OUTSIDE;
  }
  if( d0sign != d3sign )
   return OUTSIDE;

  int d4sign = sign( det4x4(
   v[0]->x, v[0]->y, v[0]->z, 1,
   v[1]->x, v[1]->y, v[1]->z, 1,
   v[2]->x, v[2]->y, v[2]->z, 1,
   vtx.x, vtx.y, vtx.z, 1) );
  if( d4sign == 0 ){
   Triangle::OptionalConstant onWhere = Triangle::PointOnFace(*v[0], *v[1], *v[2], vtx);
   switch(onWhere){
   case Triangle::EDGE_01:
    return EDGE_01;
   case Triangle::EDGE_02:
    return EDGE_02;
   case Triangle::EDGE_12:
    return EDGE_12;
   case Triangle::INSIDE:
    return FACE_012;
   }
   return OUTSIDE;
  }
  if( d0sign != d4sign )
   return OUTSIDE;

  return INSIDE;
 }

};