src/rcrt/primitives/ConvexQuad.cpp

Go to the documentation of this file.
00001 #include "ConvexQuad.h"
00002 #include <limits>
00003 
00004 using namespace std;
00005 
00006 namespace rcrt
00007 {
00008 
00009 ConvexQuad::ConvexQuad(SolidObject* parent, const Point3D& pa, const Point3D& pb, const Point3D& pc,const Point3D& pd)
00010 :Primitive(parent), va(pa), vb(pb), vc(pc), vd(pd), center((pa+(pb+(pc+pd)))/4), normal((pb-pa).crossP(pd-pa)), box(pa, pb)
00011 {
00012         box.extend(vc);
00013         box.extend(vd);
00014         normal.normalize();
00015         //normal = Vec3D(1,1,1);
00016 }
00017 
00018 ConvexQuad::~ConvexQuad()
00019 {
00020 }
00021 
00022 bool ConvexQuad::isInside(Intersection& is) const
00023 {
00024         char    LargestComponent;             /* of the normal vector         */
00025     Vec2D   A, B, C, D;                   /* Projected vertices           */
00026     Vec2D   M;                            /* Projected intersection point */
00027     Vec2D   AB, BC, CD, AD, AM, AE;       /* Miscellanous 3D-vectors      */
00028     float       u, v;                         /* Parametric coordinates       */
00029     float       a, b, c, SqrtDelta;           /* Quadratic equation           */
00030     bool       isValid = false;         /* isValid flag            */
00031     Vec2D   Vector;                       /* Temporary 2D-vector          */
00032     
00033 
00034     /*
00035     ** Projection on the plane that is most parallel to the facet
00036     */
00037     Vec3D absNormal = absNormal.abs();
00038     
00039     LargestComponent = (fabs(normal.x()) > fabs(normal.y()) ? (fabs(normal.x()) > fabs(normal.z()) ? 'x' : 'z') : (fabs(normal.y()) > fabs(normal.z()) ? 'y' : 'z'));
00040 
00041     if (LargestComponent == 'x') {
00042         A[0] = va[1]; B[0] = vb[1]; C[0] = vc[1]; D[0] = vd[1];
00043         M[0] = is.getPosition()[1];
00044     }
00045     else {
00046         A[0] = va[0]; B[0] = vb[0]; C[0] = vc[0]; D[0] = vd[0];
00047         M[0] = is.getPosition()[0];
00048     }
00049 
00050     if (LargestComponent == 'z') {
00051         A[1] = va[1]; B[1] = vb[1]; C[1] = vc[1]; D[1] = vd[1];
00052         M[1] = is.getPosition()[1];
00053     }
00054     else {
00055         A[1] = va[2]; B[1] = vb[2]; C[1] = vc[2]; D[1] = vd[2];
00056         M[1] = is.getPosition()[2];
00057     }
00058 
00059     AB = B - A; BC = C - B;
00060     CD = D - C; AD = D - A;
00061     AE = CD + AB; AE = AE * -1; AM = M - A;
00062 
00063     if (fabs(AB.det(CD)) < numeric_limits<float>::epsilon())             /* case AB // CD  - ZERO_TOL (DELTA_VEC2(AB, CD), MY_TOL) */
00064     {
00065         Vector = AB - CD;
00066         v = AM.det(Vector) / AD.det(Vector);
00067         if ((v >= 0.0) && (v <= 1.0)) {
00068             b = AB.det(AD) - AM.det(AE);
00069             c = AM.det(AD);
00070             u = fabs(b) < numeric_limits<float>::epsilon() ? -1.0 : c/b;
00071             isValid = ((u >= 0.0) && (u <= 1.0));
00072         }
00073     }
00074     else if (fabs(BC.det(AD)) < numeric_limits<float>::epsilon())         /* case AD // BC */
00075     {
00076         Vector = AD + BC;
00077         u = AM.det(Vector) / AB.det(Vector);
00078         if ((u >= 0.0) && (u <= 1.0)) {
00079             b = AD.det(AB) - AM.det(AE);
00080             c = AM.det(AB);
00081             v = fabs(b) < numeric_limits<float>::epsilon() ? -1.0 : c/b;
00082             isValid = ((v >= 0.0) && (v <= 1.0));
00083         }
00084     }
00085     else                                                    /* general case */
00086     {
00087         a = AB.det(AE); c = - AM.det(AD);
00088         b = AB.det(AD) - AM.det(AE);
00089         a = -0.5/a; b *= a; c *= (a + a); SqrtDelta = b*b + c;
00090         if (SqrtDelta >= 0.0) {
00091             SqrtDelta = sqrt(SqrtDelta);
00092             u = b - SqrtDelta;
00093             if ((u < 0.0) || (u > 1.0))        /* we want u between 0 and 1 */
00094                 u = b + SqrtDelta;
00095             if ((u >= 0.0) && (u <= 1.0)) {
00096                 v = AD[0] + u * AE[0];
00097                 if (fabs(v) < numeric_limits<float>::epsilon())
00098                     v = (AM[1] - u * AB[1]) / (AD[1] + u * AE[1]);
00099                 else
00100                     v = (AM[0] - u * AB[0]) / v;
00101                 isValid = ((v >= 0.0) && (v <= 1.0));
00102             }
00103         }
00104     }
00105 
00106     if (isValid) {
00107         is.setUV(Point2D(u,v));
00108     }
00109     return (isValid);
00110 }
00111 
00112 const AABB& ConvexQuad::getBoundingBox() const
00113 {
00114         return box;
00115 }
00116 
00117 const Point3D& ConvexQuad::getCentroid() const
00118 {
00119         return center;
00120 }
00121 
00122 Intersection ConvexQuad::intersect(Ray& r) const
00123 {
00124     /* if the ray is parallel to the quadrangle, there is no intersection */
00125     float dist = r.dir() * normal;
00126     if (fabs(dist) < numeric_limits<float>::epsilon())
00127         return Intersection();
00128 
00129     /* compute ray intersection with the plane of the quadrangle */
00130     Vec3D v = va - r.org();
00131     dist = (v*normal) / dist;
00132     
00133     Intersection is(dist, this, r.atDistance(dist), 0, 0);
00134     /* is the point in the quadrangle ? */
00135     if (isInside(is)) {
00136         is.setSNormalL(normal.abs());
00137         return is;
00138     }
00139     return Intersection();
00140 }
00141 
00142 Vec3D ConvexQuad::getSNormal(float a, float b) const
00143 {
00144         return normal;
00145 }
00146 
00147 Vec3D ConvexQuad::getGNormal(float a, float b) const
00148 {
00149         return normal;
00150 }
00151 
00152 Point2D ConvexQuad::getUV(float a, float b) const
00153 {
00154         
00155     return Point2D(a,b);
00156 }
00157 }

Generated on Thu Jan 31 19:26:19 2008 for RenderingCompetitionRayTracer by  doxygen 1.5.3