src/rcrt/primitives/Triangle.cpp

Go to the documentation of this file.

#include "Triangle.h"

using namespace std;

namespace rcrt
{


float Triangle::volume(const Point3D& o, const Point3D& e,
                        const Point3D& a, const Point3D& b) const
{
        return (a-o).crossP((b-o)) * (e-o)/6;
}


bool Triangle::sameSign(const float& a, const float& b) const
{
        return (a < 0 == b < 0) || a == 0 || b == 0;
}

Triangle::Triangle(const Point3D& pa,const Point3D& pb,const Point3D& pc,
                        const Vec3D& fNo, SolidObject* parent):Primitive(parent),
                        a(0),b(0),c(0),center((pa+pb.getPosVec3D()+pc.getPosVec3D())/3.0f),
                        faceNormal(fNo.normalized()),smooth(false)
{
        a = new Vertex(pa,faceNormal);
        b = new Vertex(pb,faceNormal);
        c = new Vertex(pc,faceNormal);
        ownVertices = true;
        updateBox();
}

Triangle::Triangle(Vertex* av,Vertex* bv,Vertex* cv, SolidObject* parent, bool smo):Primitive(parent),
a(av),b(bv),c(cv), center((a->pos()+b->pos().getPosVec3D()+c->pos().getPosVec3D())/3.0f),
faceNormal((a->normal()+b->normal()+c->normal())/3.0f),smooth(smo)
{
        updateBox();
        faceNormal.normalize();
        Vec3D fNormal = (b->pos()-a->pos()).crossP(c->pos()-a->pos());
        fNormal.normalize();
        if(fNormal*faceNormal < 0)
                fNormal = fNormal * -1;
        faceNormal = fNormal;
        ownVertices = false;    
        
        
        //calculation of vTangent
        
        Vec2D v1 = b->uv()-a->uv();
        Vec2D v2 = c->uv()-a->uv();
        
//      std::cout << "v1: (" << v1[0] << "," << v1[1] << ")" << std::endl;
//      std::cout << "v2: " << v2[0] << "," << v2[1] << ")"<< std::endl;
        
//      float e = 1 / (v2.y()-(v2.x()/v1.x())*v1.y());
//      float d = (-e*v2.x())/v1.x();
        float x = v2.det(v1);
        float d = -(v2[0]*(a->uv()[1]-1)-v2[1]*a->uv()[0])/x;
        float e =  (v1[0]*(a->uv()[1]-1)-v1[1]*a->uv()[0])/x;
        
//      std::cout << "e: " << e << std::endl;
//      std::cout << "d: " << e << std::endl;
        
        Vec3D edge1 = b->pos()-a->pos();
        Vec3D edge2 = c->pos()-a->pos();
        
//      std::cout << "edge1: " << edge1 << std::endl;
//      std::cout << "edge2: " << edge2 << std::endl;
        
        tangent = edge1 * d + edge2 * e;
//      std::cout << "tangent: " << tangent << std::endl;
}


Triangle::~Triangle()
{
        if(ownVertices){
                if(a) delete a;
                if(b) delete b;
                if(c) delete c;
        }
}

Vertex* Triangle::operator[] (int i) const
{
        switch(i){
        case 0: return a;
        case 1: return b;
        case 2: return c;
        default: throw exception();
        }
}

Vertex*& Triangle::operator[] (int i)
{
        switch(i){
        case 0: return a;
        case 1: return b;
        case 2: return c;
        default: throw exception();
        }
}

/*
Intersection Triangle::intersect(Ray& ray) const
{
        ray.tris++;
        const Point3D end = ray.org() + ray.dir();
        float V1 = volume(ray.org(), end, a->pos(), b->pos());
        //if(V1 == 0)
                //std::cout << "[TRI-Intersect] V1 == 0" << std::endl;
        float V2 = volume(ray.org(), end, b->pos(), c->pos());
        //if(V2 == 0)
                //std::cout << "[TRI-Intersect] V2 == 0" << std::endl;
        float V3 = volume(ray.org(), end, c->pos(), a->pos());
        //if(V3 == 0)
                //std::cout << "[TRI-Intersect] V3 == 0" << std::endl;
                
        if(!sameSign(V1,V2) || !sameSign(V2,V3) || !sameSign(V1,V3)){
                std::cout << "[TRI-Intersect] V1 ="<< V1 << std::endl;
                std::cout << "[TRI-Intersect] V2 ="<< V2 << std::endl;
                std::cout << "[TRI-Intersect] V3 ="<< V3 << std::endl;
                std::cout << "[TRI-Intersect] signs:"<<((V1 < 0) ? "-":"+") <<((V2 < 0) ? ",-":",+")<<((V3 < 0) ? ",-":",+") << std::endl;
                float Va = volume(a->pos(), b->pos(), ray.org(), c->pos());
                float V = V1+V2+V3;
                float t = Va/V;
                std::cout << "[TRI-Intersect] => no intersection dist="<<t << std::endl;
                return Intersection();
        }
        float Va = volume(a->pos(), b->pos(), ray.org(), c->pos());
        float V = V1+V2+V3;
        float t = Va/V;
        if(t < numeric_limits<float>::epsilon()){
                //std::cout << "[TRI-Intersect] t < "<< numeric_limits<float>::epsilon() << std::endl;
                //std::cout << "[TRI-Intersect] => no intersection" << std::endl;
                return Intersection();
        }
        //barycentric coords:
        float alpha = V2/V;
        float beta = V3/V;
        float gamma = 1-alpha-beta;
        
        bool flip = ray.dir()*faceNormal >= numeric_limits<float>::epsilon();
        //std::cout << "[TRI-Intersect] t = " << t << " backside=" << flip << " distToOrg=" << (ray.org()-getCentroid()).norm() << " angle="<<ray.dir()*normal << std::endl;
        
        Intersection is(t, this, ray.atDistance(t), beta, gamma, flip);
        
        Point2D uv(alpha * a->uv().x() + beta * b->uv().x() + gamma * c->uv().x(),
                        alpha * a->uv().y() + beta * b->uv().y() + gamma * c->uv().y());
        
        is.setUV(uv);
        
        return is;
}*/

Intersection Triangle::intersect(Ray& ray) const
{
        ray.tris++;
        const Vec3D edge1 = b->pos()-a->pos();
    const Vec3D edge2 = c->pos()-a->pos();

    const Vec3D pvec = ray.dir().crossP(edge2);

    const float det = edge1 * pvec;
    if (fabs(det) < numeric_limits<float>::epsilon()) return Intersection();

    const float inv_det = 1.0f / det;

    const Vec3D tvec = ray.org()-a->pos();
    float lambda = tvec * pvec;
    lambda *= inv_det;

    if (lambda < 0.0f || lambda > 1.0f) return Intersection();

    const Vec3D qvec = tvec.crossP(edge1);
    float mue = ray.dir() * qvec;
    mue *= inv_det;

    if (mue < 0.0f || mue+lambda > 1.0f) return Intersection();

    float f = edge2 * qvec;
    f = f * inv_det - numeric_limits<float>::epsilon();
    if (ray.maxDist() <= f || f <  numeric_limits<float>::epsilon()  ) return Intersection();

    //Vec3D normal = lambda * b->normal() + mue * c->normal() + (1-lambda-mue) * a->normal();
    //normal.normalize();
    bool flip = ray.dir()*faceNormal >= numeric_limits<float>::epsilon();
    //std::cout << "[TRI-Intersect] t = " << f << " backside=" << flip << " distToOrg=" << (ray.org()-getCentroid()).norm() << " angle="<<ray.dir()*normal << std::endl;
    
    Intersection is = Intersection(f, this, ray.atDistance(f), lambda, mue, flip);
    
    is.setUV(Point2D(lambda,mue));
    
    return is;
}

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

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

Vec3D Triangle::getSNormal(float beta, float gamma) const
{
        if(smooth)
                return (beta * b->normal() + gamma * c->normal() + (1-beta-gamma) * a->normal()).normalized();
        else
                return faceNormal;
}

Vec3D Triangle::getGNormal(float beta, float gamma) const
{
        return faceNormal;
}

Point2D Triangle::getUV(float beta, float gamma) const
{
        float alpha = 1-beta-gamma;
        return Point2D(alpha * a->uv().x() + beta * b->uv().x() + gamma * c->uv().x(),
                                alpha * a->uv().y() + beta * b->uv().y() + gamma * c->uv().y());
}

Point3D Triangle::getPoint(float beta, float gamma) const
{
        float alpha = 1-beta-gamma;
        return  (c->pos()*gamma)+((a->pos()*alpha) + (b->pos()*beta));
}

Point2D Triangle::getUV(const Point2D& p) const
{
        float beta = p[0];
        float gamma = p[1];
        return getUV(beta,gamma);
}

Point3D Triangle::getPoint(const Point2D& p) const
{
        float beta = p[0];
        float gamma = p[1];
        return  getPoint(beta,gamma);
}

Vec3D Triangle::getTangent(float beta, float gamma) const
{
        return tangent; 
}

void Triangle::updateBox()
{
        Point3D maxp = Point3D(max(a->pos().x(),max(b->pos().x(),c->pos().x())),
                        max(a->pos().y(),max(b->pos().y(),c->pos().y())),
                        max(a->pos().z(),max(b->pos().z(),c->pos().z())));
        //fixing triangles with zero extension box
        //the BIH doesn't seem to like them
        maxp[0] += 0.001;
        maxp[1] += 0.001;
        maxp[2] += 0.001;
        Point3D minp = Point3D(
                        min(a->pos().x(),min(b->pos().x(),c->pos().x())),
                        min(a->pos().y(),min(b->pos().y(),c->pos().y())),
                        min(a->pos().z(),min(b->pos().z(),c->pos().z())));
        minp[0] -= 0.001;
        minp[1] -= 0.001;
        minp[2] -= 0.001;
        box.extend(maxp);
        box.extend(minp);
        
}

bool Triangle::isSmooth() const
{
        return smooth;
}

void Triangle::clipPlane(Axis axis, float plane, AABB& lBox, AABB& rBox) const {
        int max=0;
        Vec3D rays[2][2];
        Point3D p0(0);
        
        if (a->pos()[axis] <= plane){
                lBox.extend(a->pos());
                if (b->pos()[axis] > plane){
                        rays[max][0] = a->pos()-p0;
                        rays[max++][1]= b->pos()-a->pos();
                }
                if (c->pos()[axis] > plane){
                        rays[max][0]=a->pos()-p0;
                        rays[max++][1]=c->pos()-a->pos();
                }
        }
        else
                rBox.extend(a->pos());

        if (b->pos()[axis] <= plane){
                lBox.extend(b->pos());
                if (a->pos()[axis] > plane){
                        rays[max][0]=b->pos()-p0;
                        rays[max++][1]=a->pos()-b->pos();
                }
                if (c->pos()[axis] > plane){
                        rays[max][0]=b->pos()-p0;
                        rays[max++][1]=c->pos()-b->pos();
                }
        }
        else
                rBox.extend(b->pos());

        if (c->pos()[axis] <= plane){
                lBox.extend(c->pos());
                if (a->pos()[axis] > plane){
                        rays[max][0]=c->pos()-p0;
                        rays[max++][1]=a->pos()-c->pos();
                }
                if (b->pos()[axis] > plane){
                        rays[max][0]=c->pos()-p0;
                        rays[max++][1]=b->pos()-c->pos();
                }
        }
        else
                rBox.extend(c->pos());

        // calculate new Bounds of the Triangle (left and right)
        Point3D point;
        for (int i=0; i<max; i++)
        {
                (rays[i][1]).normalize();
                float t=(plane - rays[i][0][axis])/rays[i][1][axis];

                point=Point3D(rays[i][0]+t*rays[i][1]);
                lBox.extend(point);
                rBox.extend(point);
        }
}

}

---