📅 2010-Oct-31 ⬩ ✍️ Ashwin Nanjappa ⬩ 🏷️ cpp, orientation ⬩ 📚 Archive
Orientation is a fundamental concept in computational geometry. It is the most important predicate too. Calculating the orientation of a triangle (2D), tetrahedron (3D) or pentachoron (4D) reduces to the calculation of the sign of a determinant. The following code is in C++, but can easily be adapted to any language.
First, it is helpful to define an orientation type:
enum ORIENT_TYPE
{
ORIENT_NEG = -1,
ORIENT_ZERO = +0,
ORIENT_POS = +1,
};ORIENT_ZERO indicates a degeneracy in the simplex. For example, in 2D this means that the 3 points are collinear. This is typically not preferred and symbolic perturbation techniques like Simulation of Simplicity might have to be used to handle such cases.
A helper function to convert a determinant value to orientation:
ORIENT_TYPE detToOrient( const int det )
{
return ( det > 0 ) ? ORIENT_POS : ( ( det < 0 ) ? ORIENT_NEG : ORIENT_ZERO );
}2D Orientation:
ORIENT_TYPE orientation2
(
const int* p0,
const int* p1,
const int* p2
)
{
// | m00 m01 |
// | m10 m11 |
const int m00 = p1[0] - p0[0];
const int m01 = p1[1] - p0[1];
const int m10 = p2[0] - p0[0];
const int m11 = p2[1] - p0[1];
const int det = ( m00 * m11 ) - ( m10 * m01 );
return detToOrient( det );
}3D Orientation:
ORIENT_TYPE orientation3
(
const int* p0,
const int* p1,
const int* p2,
const int* p3
)
{
// | m00 m01 m02 |
// | m10 m11 m12 |
// | m20 m21 m22 |
const int m00 = p1[0] - p0[0];
const int m01 = p1[1] - p0[1];
const int m02 = p1[2] - p0[2];
const int m10 = p2[0] - p0[0];
const int m11 = p2[1] - p0[1];
const int m12 = p2[2] - p0[2];
const int m20 = p3[0] - p0[0];
const int m21 = p3[1] - p0[1];
const int m22 = p3[2] - p0[2];
const int m00Det = ( m11 * m22 ) - ( m21 * m12 );
const int m01Det = ( m10 * m22 ) - ( m20 * m12 );
const int m02Det = ( m10 * m21 ) - ( m20 * m11 );
const int det = ( m00 * m00Det ) - ( m01 * m01Det ) + ( m02 * m02Det );
return detToOrient( det );
}4D Orientation:
ORIENT_TYPE orientation4
(
const int* p0,
const int* p1,
const int* p2,
const int* p3,
const int* p4
)
{
// | m00 m01 m02 m03 |
// | m10 m11 m12 m13 |
// | m20 m21 m22 m23 |
// | m30 m31 m32 m33 |
const int m00 = p1[0] - p0[0];
const int m01 = p2[0] - p0[0];
const int m02 = p3[0] - p0[0];
const int m03 = p4[0] - p0[0];
const int m10 = p1[1] - p0[1];
const int m11 = p2[1] - p0[1];
const int m12 = p3[1] - p0[1];
const int m13 = p4[1] - p0[1];
const int m20 = p1[2] - p0[2];
const int m21 = p2[2] - p0[2];
const int m22 = p3[2] - p0[2];
const int m23 = p4[2] - p0[2];
const int m30 = p1[3] - p0[3];
const int m31 = p2[3] - p0[3];
const int m32 = p3[3] - p0[3];
const int m33 = p4[3] - p0[3];
const int m00m11Det = ( m22 * m33 ) - ( m32 * m23 );
const int m00m12Det = ( m21 * m33 ) - ( m31 * m23 );
const int m00m13Det = ( m21 * m32 ) - ( m31 * m22 );
const int m01m10Det = m00m11Det;
const int m01m12Det = ( m20 * m33 ) - ( m30 * m23 );
const int m01m13Det = ( m20 * m32 ) - ( m30 * m22 );
const int m02m10Det = m00m12Det;
const int m02m11Det = m01m12Det;
const int m02m13Det = ( m20 * m31 ) - ( m30 * m21 );
const int m03m10Det = m00m13Det;
const int m03m11Det = m01m13Det;
const int m03m12Det = m02m13Det;
const int m00Det =
+ ( m11 * m00m11Det )
- ( m12 * m00m12Det )
+ ( m13 * m00m13Det );
const int m01Det =
+ ( m10 * m01m10Det )
- ( m12 * m01m12Det )
+ ( m13 * m01m13Det );
const int m02Det =
+ ( m10 * m02m10Det )
- ( m11 * m02m11Det )
+ ( m13 * m02m13Det );
const int m03Det =
+ ( m10 * m03m10Det )
- ( m11 * m03m11Det )
+ ( m12 * m03m12Det );
const int det =
+ ( m00 * m00Det )
- ( m01 * m01Det )
+ ( m02 * m02Det )
- ( m03 * m03Det );
return detToOrient( det );
}