triangle3d.hpp

// Copyright (C) 2002-2012 Nikolaus Gebhardt
// This file is part of the "Irrlicht Engine".
// For conditions of distribution and use, see copyright notice in nirtcpp/nirtcpp.hpp
 
#ifndef NIRT_TRIANGLE_3D_HPP_INCLUDED
#define NIRT_TRIANGLE_3D_HPP_INCLUDED
 
#include <nirtcpp/core/engine/vector3d.hpp>
#include <nirtcpp/core/engine/line3d.hpp>
#include <nirtcpp/core/engine/plane3d.hpp>
#include <nirtcpp/core/engine/aabbox3d.hpp>
 
namespace nirt
{
namespace core
{
 
    template <class T>
    class triangle3d
    {
    public:
 
        triangle3d() {}
        triangle3d(const vector3d<T>& v1, const vector3d<T>& v2, const vector3d<T>& v3) : pointA(v1), pointB(v2), pointC(v3) {}
 
        bool operator==(const triangle3d<T>& other) const
        {
            return other.pointA==pointA && other.pointB==pointB && other.pointC==pointC;
        }
 
        bool operator!=(const triangle3d<T>& other) const
        {
            return !(*this==other);
        }
 
 
        bool isTotalInsideBox(const aabbox3d<T>& box) const
        {
            return (box.isPointInside(pointA) &&
                box.isPointInside(pointB) &&
                box.isPointInside(pointC));
        }
 
 
        bool isTotalOutsideBox(const aabbox3d<T>& box) const
        {
            return ((pointA.X > box.MaxEdge.X && pointB.X > box.MaxEdge.X && pointC.X > box.MaxEdge.X) ||
 
                (pointA.Y > box.MaxEdge.Y && pointB.Y > box.MaxEdge.Y && pointC.Y > box.MaxEdge.Y) ||
                (pointA.Z > box.MaxEdge.Z && pointB.Z > box.MaxEdge.Z && pointC.Z > box.MaxEdge.Z) ||
                (pointA.X < box.MinEdge.X && pointB.X < box.MinEdge.X && pointC.X < box.MinEdge.X) ||
                (pointA.Y < box.MinEdge.Y && pointB.Y < box.MinEdge.Y && pointC.Y < box.MinEdge.Y) ||
                (pointA.Z < box.MinEdge.Z && pointB.Z < box.MinEdge.Z && pointC.Z < box.MinEdge.Z));
        }
 
 
        core::vector3d<T> closestPointOnTriangle(const core::vector3d<T>& p) const
        {
            const core::vector3d<T> rab = line3d<T>(pointA, pointB).getClosestPoint(p);
            const core::vector3d<T> rbc = line3d<T>(pointB, pointC).getClosestPoint(p);
            const core::vector3d<T> rca = line3d<T>(pointC, pointA).getClosestPoint(p);
 
            const T d1 = rab.getDistanceFrom(p);
            const T d2 = rbc.getDistanceFrom(p);
            const T d3 = rca.getDistanceFrom(p);
 
            if (d1 < d2)
                return d1 < d3 ? rab : rca;
 
            return d2 < d3 ? rbc : rca;
        }
 
        /*
        \param p Point to test. Assumes that this point is already
        on the plane of the triangle.
        \return True if the point is inside the triangle, otherwise false. */
        bool isPointInside(const vector3d<T>& p) const
        {
            vector3d<f64> af64((f64)pointA.X, (f64)pointA.Y, (f64)pointA.Z);
            vector3d<f64> bf64((f64)pointB.X, (f64)pointB.Y, (f64)pointB.Z);
            vector3d<f64> cf64((f64)pointC.X, (f64)pointC.Y, (f64)pointC.Z);
            vector3d<f64> pf64((f64)p.X, (f64)p.Y, (f64)p.Z);
            return (isOnSameSide(pf64, af64, bf64, cf64) &&
                isOnSameSide(pf64, bf64, af64, cf64) &&
                isOnSameSide(pf64, cf64, af64, bf64));
        }
 
 
        bool isPointInsideFast(const vector3d<T>& p) const
        {
            const vector3d<T> a = pointC - pointA;
            const vector3d<T> b = pointB - pointA;
            const vector3d<T> c = p - pointA;
 
            const f64 dotAA = a.dotProduct( a);
            const f64 dotAB = a.dotProduct( b);
            const f64 dotAC = a.dotProduct( c);
            const f64 dotBB = b.dotProduct( b);
            const f64 dotBC = b.dotProduct( c);
 
            // get coordinates in barycentric coordinate system
            const f64 invDenom =  1/(dotAA * dotBB - dotAB * dotAB);
            const f64 u = (dotBB * dotAC - dotAB * dotBC) * invDenom;
            const f64 v = (dotAA * dotBC - dotAB * dotAC ) * invDenom;
 
            // We count border-points as inside to keep downward compatibility.
            // Rounding-error also needed for some test-cases.
            return (u > -ROUNDING_ERROR_f32) && (v >= 0) && (u + v < 1+ROUNDING_ERROR_f32);
 
        }
 
 
 
        bool getIntersectionWithLimitedLine(const line3d<T>& line,
            vector3d<T>& outIntersection) const
        {
            return getIntersectionWithLine(line.start,
                line.getVector(), outIntersection) &&
                outIntersection.isBetweenPoints(line.start, line.end);
        }
 
 
 
        bool getIntersectionWithLine(const vector3d<T>& linePoint,
            const vector3d<T>& lineVect, vector3d<T>& outIntersection) const
        {
            if (getIntersectionOfPlaneWithLine(linePoint, lineVect, outIntersection))
                return isPointInside(outIntersection);
 
            return false;
        }
 
 
 
        bool getIntersectionOfPlaneWithLine(const vector3d<T>& linePoint,
            const vector3d<T>& lineVect, vector3d<T>& outIntersection) const
        {
            // Work with f64 to get more precise results (makes enough difference to be worth the casts).
            const vector3d<f64> linePointf64(linePoint.X, linePoint.Y, linePoint.Z);
            const vector3d<f64> lineVectf64(lineVect.X, lineVect.Y, lineVect.Z);
            vector3d<f64> outIntersectionf64;
 
            core::triangle3d<nirt::f64> trianglef64(vector3d<f64>((f64)pointA.X, (f64)pointA.Y, (f64)pointA.Z)
                                        ,vector3d<f64>((f64)pointB.X, (f64)pointB.Y, (f64)pointB.Z)
                                        , vector3d<f64>((f64)pointC.X, (f64)pointC.Y, (f64)pointC.Z));
            const vector3d<nirt::f64> normalf64 = trianglef64.getNormal().normalize();
            f64 t2;
 
            if ( core::iszero ( t2 = normalf64.dotProduct(lineVectf64) ) )
                return false;
 
            f64 d = trianglef64.pointA.dotProduct(normalf64);
            f64 t = -(normalf64.dotProduct(linePointf64) - d) / t2;
            outIntersectionf64 = linePointf64 + (lineVectf64 * t);
 
            outIntersection.X = (T)outIntersectionf64.X;
            outIntersection.Y = (T)outIntersectionf64.Y;
            outIntersection.Z = (T)outIntersectionf64.Z;
            return true;
        }
 
 
 
        vector3d<T> getNormal() const
        {
            return (pointB - pointA).crossProduct(pointC - pointA);
        }
 
 
        bool isFrontFacing(const vector3d<T>& lookDirection) const
        {
            const vector3d<T> n = getNormal().normalize();
            const f32 d = (f32)n.dotProduct(lookDirection);
            return F32_LOWER_EQUAL_0(d);
        }
 
        plane3d<T> getPlane() const
        {
            return plane3d<T>(pointA, pointB, pointC);
        }
 
        T getArea() const
        {
            return (pointB - pointA).crossProduct(pointC - pointA).getLength() * 0.5f;
 
        }
 
        void set(const core::vector3d<T>& a, const core::vector3d<T>& b, const core::vector3d<T>& c)
        {
            pointA = a;
            pointB = b;
            pointC = c;
        }
 
        vector3d<T> pointA;
        vector3d<T> pointB;
        vector3d<T> pointC;
 
    private:
        // Using f64 instead of <T> to avoid integer overflows when T=int (maybe also less floating point troubles).
        bool isOnSameSide(const vector3d<f64>& p1, const vector3d<f64>& p2,
            const vector3d<f64>& a, const vector3d<f64>& b) const
        {
            vector3d<f64> bminusa = b - a;
            vector3d<f64> cp1 = bminusa.crossProduct(p1 - a);
            vector3d<f64> cp2 = bminusa.crossProduct(p2 - a);
            f64 res = cp1.dotProduct(cp2);
            if ( res < 0 )
            {
                // This catches some floating point troubles.
                // Unfortunately slightly expensive and we don't really know the best epsilon for iszero.
                vector3d<f64> cp1n = bminusa.normalize().crossProduct((p1 - a).normalize());
                if (core::iszero(cp1n.X, (f64)ROUNDING_ERROR_f32)
                    && core::iszero(cp1n.Y, (f64)ROUNDING_ERROR_f32)
                    && core::iszero(cp1n.Z, (f64)ROUNDING_ERROR_f32) )
                {
                    res = 0.f;
                }
            }
            return (res >= 0.0f);
        }
    };
 
 
    using triangle3df = triangle3d<f32>;
 
    using triangle3di = triangle3d<s32>;
 
} // end namespace core
} // end namespace nirt
 
#endif