src/rcrt/primitives/ConvexQuad.cpp

Go to the documentation of this file.

#include "ConvexQuad.h"
#include <limits>

using namespace std;

namespace rcrt
{

ConvexQuad::ConvexQuad(SolidObject* parent, const Point3D& pa, const Point3D& pb, const Point3D& pc,const Point3D& pd)
: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)
{
        box.extend(vc);
        box.extend(vd);
        normal.normalize();
        //normal = Vec3D(1,1,1);
}

ConvexQuad::~ConvexQuad()
{
}

bool ConvexQuad::isInside(Intersection& is) const
{
        char    LargestComponent;             /* of the normal vector         */
    Vec2D   A, B, C, D;                   /* Projected vertices           */
    Vec2D   M;                            /* Projected intersection point */
    Vec2D   AB, BC, CD, AD, AM, AE;       /* Miscellanous 3D-vectors      */
    float       u, v;                         /* Parametric coordinates       */
    float       a, b, c, SqrtDelta;           /* Quadratic equation           */
    bool       isValid = false;         /* isValid flag            */
    Vec2D   Vector;                       /* Temporary 2D-vector          */
    

    /*
    ** Projection on the plane that is most parallel to the facet
    */
    Vec3D absNormal = absNormal.abs();
    
    LargestComponent = (fabs(normal.x()) > fabs(normal.y()) ? (fabs(normal.x()) > fabs(normal.z()) ? 'x' : 'z') : (fabs(normal.y()) > fabs(normal.z()) ? 'y' : 'z'));

    if (LargestComponent == 'x') {
        A[0] = va[1]; B[0] = vb[1]; C[0] = vc[1]; D[0] = vd[1];
        M[0] = is.getPosition()[1];
    }
    else {
        A[0] = va[0]; B[0] = vb[0]; C[0] = vc[0]; D[0] = vd[0];
        M[0] = is.getPosition()[0];
    }

    if (LargestComponent == 'z') {
        A[1] = va[1]; B[1] = vb[1]; C[1] = vc[1]; D[1] = vd[1];
        M[1] = is.getPosition()[1];
    }
    else {
        A[1] = va[2]; B[1] = vb[2]; C[1] = vc[2]; D[1] = vd[2];
        M[1] = is.getPosition()[2];
    }

    AB = B - A; BC = C - B;
    CD = D - C; AD = D - A;
    AE = CD + AB; AE = AE * -1; AM = M - A;

    if (fabs(AB.det(CD)) < numeric_limits<float>::epsilon())             /* case AB // CD  - ZERO_TOL (DELTA_VEC2(AB, CD), MY_TOL) */
    {
        Vector = AB - CD;
        v = AM.det(Vector) / AD.det(Vector);
        if ((v >= 0.0) && (v <= 1.0)) {
            b = AB.det(AD) - AM.det(AE);
            c = AM.det(AD);
            u = fabs(b) < numeric_limits<float>::epsilon() ? -1.0 : c/b;
            isValid = ((u >= 0.0) && (u <= 1.0));
        }
    }
    else if (fabs(BC.det(AD)) < numeric_limits<float>::epsilon())         /* case AD // BC */
    {
        Vector = AD + BC;
        u = AM.det(Vector) / AB.det(Vector);
        if ((u >= 0.0) && (u <= 1.0)) {
            b = AD.det(AB) - AM.det(AE);
            c = AM.det(AB);
            v = fabs(b) < numeric_limits<float>::epsilon() ? -1.0 : c/b;
            isValid = ((v >= 0.0) && (v <= 1.0));
        }
    }
    else                                                    /* general case */
    {
        a = AB.det(AE); c = - AM.det(AD);
        b = AB.det(AD) - AM.det(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[0] + u * AE[0];
                if (fabs(v) < numeric_limits<float>::epsilon())
                    v = (AM[1] - u * AB[1]) / (AD[1] + u * AE[1]);
                else
                    v = (AM[0] - u * AB[0]) / v;
                isValid = ((v >= 0.0) && (v <= 1.0));
            }
        }
    }

    if (isValid) {
        is.setUV(Point2D(u,v));
    }
    return (isValid);
}

const AABB& ConvexQuad::getBoundingBox() const
{
        return box;
}

const Point3D& ConvexQuad::getCentroid() const
{
        return center;
}

Intersection ConvexQuad::intersect(Ray& r) const
{
    /* if the ray is parallel to the quadrangle, there is no intersection */
    float dist = r.dir() * normal;
    if (fabs(dist) < numeric_limits<float>::epsilon())
        return Intersection();

    /* compute ray intersection with the plane of the quadrangle */
    Vec3D v = va - r.org();
    dist = (v*normal) / dist;
    
    Intersection is(dist, this, r.atDistance(dist), 0, 0);
    /* is the point in the quadrangle ? */
    if (isInside(is)) {
        is.setSNormalL(normal.abs());
        return is;
    }
    return Intersection();
}

Vec3D ConvexQuad::getSNormal(float a, float b) const
{
        return normal;
}

Vec3D ConvexQuad::getGNormal(float a, float b) const
{
        return normal;
}

Point2D ConvexQuad::getUV(float a, float b) const
{
        
    return Point2D(a,b);
}
}

---