Code Yarns ‍👨‍💻
Tech BlogPersonal Blog

Orientation

📅 2010-Oct-31 ⬩ ✍️ Ashwin Nanjappa ⬩ 📚 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 );
}