#include #include "GraphicsGems.h" /* ---- intqdr - Intersect a ray with a quadric surface. -------------- */ /* */ /* */ /* Description: */ /* Intqdr determines the intersection of a ray with a quadric */ /* surface in matrix form. */ /* */ /* On entry: */ /* raybase = The coordinate defining the base of the */ /* intersecting ray. */ /* raycos = The direction cosines of the above ray. */ /* base = The location of the base of the quadric. */ /* axis = The axis of symmetry for the quadric. */ /* radius = The bounding radius from the axis of symmetry. */ /* Q = A 4x4 coefficient matrix representing the quadric */ /* surface. */ /* */ /* On return: */ /* rin = The entering distance of the intersection. */ /* rout = The leaving distance of the intersection. */ /* */ /* */ /* Returns: True if the ray intersects the quadric surface. */ /* */ /* -------------------------------------------------------------------- */ boolean intqdr (raybase,raycos,base,axis,radius,Q,rin,rout) Point3 raybase; /* Base of the intersection ray */ Vector3 raycos; /* Direction cosines of the ray */ Vector3 base; /* Base point for the quadric */ Vector3 axis; /* Axis of symmetry */ double radius; /* Bounding radius */ Matrix4 Q; /* Quadric coefficient matrix */ double *rin; /* Entering distance */ double *rout; /* Leaving distance */ { boolean hit; /* True if ray intersects quad. */ double qa, qb, qc, qd, qe, /* Coefficient matrix terms */ qf, qg, qh, qi, qj; double k2, k1, k0; /* Quadric coefficients */ double disc; Vector3 D; /* Ray base to quadric base */ double d; /* Shortest distance between */ /* ray and axis of symmetry */ Vector3 N; /* Compute the perpendicular distance between the symmetry axis */ /* and the ray. */ V3Sub (&base,&raybase,&D); /* compute shortest distance */ V3Cross (&raycos,&axis,&N); if ((d=V3Length(&N)) == 0) { /* if parallel */ V3Cross (&raycos,&D,&N); d = fabs (V3Length(&N)); } else { V3Normalize (&N); d = fabs (V3Dot(&D,&N)/d); } hit = (d <= radius); if (!hit) return (hit); /* If outside of cylinder */ qa = Q.element[0][0]; qb = Q.element[0][1] + Q.element[1][0]; qc = Q.element[0][2] + Q.element[2][0]; qd = Q.element[0][3] + Q.element[3][0]; qe = Q.element[1][1]; qf = Q.element[1][2] + Q.element[2][1]; qg = Q.element[1][3] + Q.element[3][1]; qh = Q.element[2][2]; qi = Q.element[2][3] + Q.element[3][2]; qj = Q.element[3][3]; k2 = raycos.x * (raycos.x*qa + raycos.y*qb) + raycos.y * (raycos.y*qe + raycos.z*qf) + raycos.z * (raycos.z*qh + raycos.z*qc); k1 = raycos.x * ((raybase.x*2.*qa + raybase.y*qb + raybase.z*qc) + qd) + raycos.y * ((raybase.x*qb + raybase.y*2.*qe + raybase.z*qf) + qg) + raycos.z * ((raybase.x*qc + raybase.y*qg + raybase.z*2.*qh) + qi); k0 = raybase.x * (raybase.x*qa + raybase.y*qb + raybase.z*qc + qd) + raybase.y * ( raybase.y*qe + raybase.z*qf + qg) + raybase.z * ( raybase.z*qh + qi) + qj; disc = k1*k1 - 4.*k2*k0; hit = (disc >= 0.0); if (hit) { /* If ray hits quadric */ disc = sqrt(disc); *rin = (-k1 - disc) / (2.*k2); /* entering distance */ *rout = (-k1 + disc) / (2.*k2); /* leaving distance */ } return (hit); }