/* ------------------------------------------------------------------------- *\   QUAD.C :   by Christophe Schlick and Gilles Subrenat (15 May 1994)   "Ray Intersection of Tessellated Surfaces : Quadrangles versus Triangles"   in Graphics Gems V (edited by A. Paeth), Academic Press\* ------------------------------------------------------------------------- */#include "quad.h"/*** Macro definitions*/#define MY_TOL                ((real) 0.0001)#define LARGEST_COMPONENT(A)  (ABS((A).x) > ABS((A).y) ? \                              (ABS((A).x) > ABS((A).z) ? 'x' : 'z') : \                              (ABS((A).y) > ABS((A).z) ? 'y' : 'z'))/*** Check if the point is in the quadrangle*/static bool point_in_quad (QUAD *Quad, HIT *Hit){    char       LargestComponent;             /* of the normal vector         */    realvec2   A, B, C, D;                   /* Projected vertices           */    realvec2   M;                            /* Projected intersection point */    realvec2   AB, BC, CD, AD, AM, AE;       /* Miscellanous 3D-vectors      */    real       u, v;                         /* Parametric coordinates       */    real       a, b, c, SqrtDelta;           /* Quadratic equation           */    bool       Intersection = FALSE;         /* Intersection flag            */    realvec2   Vector;                       /* Temporary 2D-vector          */    /*    ** Projection on the plane that is most parallel to the facet    */    LargestComponent = LARGEST_COMPONENT(Quad->Normal);    if (LargestComponent == 'x') {        A.x = Quad->A.y; B.x = Quad->B.y; C.x = Quad->C.y; D.x = Quad->D.y;        M.x = Hit->Point.y;    }    else {        A.x = Quad->A.x; B.x = Quad->B.x; C.x = Quad->C.x; D.x = Quad->D.x;        M.x = Hit->Point.x;    }    if (LargestComponent == 'z') {        A.y = Quad->A.y; B.y = Quad->B.y; C.y = Quad->C.y; D.y = Quad->D.y;        M.y = Hit->Point.y;    }    else {        A.y = Quad->A.z; B.y = Quad->B.z; C.y = Quad->C.z; D.y = Quad->D.z;        M.y = Hit->Point.z;    }    SUB_VEC2 (AB, B, A); SUB_VEC2 (BC, C, B);    SUB_VEC2 (CD, D, C); SUB_VEC2 (AD, D, A);    ADD_VEC2 (AE, CD, AB); NEG_VEC2 (AE, AE); SUB_VEC2 (AM, M, A);    if (ZERO_TOL (DELTA_VEC2(AB, CD), MY_TOL))             /* case AB // CD */    {        SUB_VEC2 (Vector, AB, CD);        v = DELTA_VEC2(AM, Vector) / DELTA_VEC2(AD, Vector);        if ((v >= 0.0) && (v <= 1.0)) {            b = DELTA_VEC2(AB, AD) - DELTA_VEC2(AM, AE);            c = DELTA_VEC2 (AM, AD);            u = ZERO_TOL(b, MY_TOL) ? -1.0 : c/b;            Intersection = ((u >= 0.0) && (u <= 1.0));        }    }    else if (ZERO_TOL(DELTA_VEC2(BC, AD), MY_TOL))         /* case AD // BC */    {        ADD_VEC2 (Vector, AD, BC);        u = DELTA_VEC2(AM, Vector) / DELTA_VEC2(AB, Vector);        if ((u >= 0.0) && (u <= 1.0)) {            b = DELTA_VEC2(AD, AB) - DELTA_VEC2(AM, AE);            c = DELTA_VEC2 (AM, AB);            v = ZERO_TOL(b, MY_TOL) ? -1.0 : c/b;            Intersection = ((v >= 0.0) && (v <= 1.0));        }    }    else                                                    /* general case */    {        a = DELTA_VEC2(AB, AE); c = - DELTA_VEC2 (AM,AD);        b = DELTA_VEC2(AB, AD) - DELTA_VEC2(AM, AE);        a = -0.5/a; b *= a; c *= (a + a); SqrtDelta = b*b + c;        if (SqrtDelta >= 0.0) {            SqrtDelta = sqrt(SqrtDelta);            u = b - SqrtDelta;            if ((u < 0.0) || (u > 1.0))        /* we want u between 0 and 1 */                u = b + SqrtDelta;            if ((u >= 0.0) && (u <= 1.0)) {                v = AD.x + u * AE.x;                if (ZERO_TOL(v, MY_TOL))                    v = (AM.y - u * AB.y) / (AD.y + u * AE.y);                else                    v = (AM.x - u * AB.x) / v;                Intersection = ((v >= 0.0) && (v <= 1.0));            }        }    }    if (Intersection) {        Hit->u = u;        Hit->v = v;    }    return (Intersection);}/*** Search for an intersection between a quadrangle and a ray*/bool hit_ray_quad (RAY *Ray, QUAD *Quad, HIT *Hit){    realvec3     Point;    /* if the ray is parallel to the quadrangle, there is no intersection */    Hit->Distance = DOT_VEC3 (Ray->Vector, Quad->Normal);    if (ZERO_TOL(Hit->Distance, MY_TOL)) return (FALSE);    /* compute ray intersection with the plane of the quadrangle */    SUB_VEC3 (Point, Quad->A, Ray->Point);    Hit->Distance = DOT_VEC3 (Point, Quad->Normal) / Hit->Distance;    MULS_VEC3 (Hit->Point, Ray->Vector, Hit->Distance);    INC_VEC3 (Hit->Point, Ray->Point);    /* is the point in the quadrangle ? */    return (point_in_quad(Quad, Hit));}/* ------------------------------------------------------------------------- */