PointTriangleNearest.h

Go to the documentation of this file.

/*  Copyright (C) 2003-2021 Free Electron Organization
    Any use of this software requires a license.  If a valid license
    was not distributed with this file, visit freeelectron.org. */

/** @file */

#ifndef __geometry_PointTriangleNearest_h__
#define __geometry_PointTriangleNearest_h__

#define FE_PTN_DEBUG    FALSE

namespace fe
{
namespace ext
{

/**************************************************************************//**
    @brief Find point nearest to a triangle

    @ingroup geometry
*//***************************************************************************/
template <typename T>
class PointTriangleNearest
{
    public:
static  T   solve(const Vector<3,T>& vert0,
                    const Vector<3,T>& vert1,
                    const Vector<3,T>& vert2,
                    const Vector<3,T>& origin,
                    Vector<3,T>& direction,
                    Vector<3,T>& intersection,
                    Barycenter<T>& barycenter);

};

template <typename T>
inline T PointTriangleNearest<T>::solve(const Vector<3,T>& vert0,
        const Vector<3,T>& vert1,
        const Vector<3,T>& vert2,
        const Vector<3,T>& origin,
        Vector<3,T>& direction,
        Vector<3,T>& intersection,
        Barycenter<T>& barycenter)
{
#if FE_PTN_DEBUG
    feLog("\nPointTriangleNearest<T>\n");
    feLog("vert0 %s\n",c_print(vert0));
    feLog("vert1 %s\n",c_print(vert1));
    feLog("vert2 %s\n",c_print(vert2));
    feLog("origin %s\n",c_print(origin));
#endif
    const Vector<3,T> perp=cross(vert1-vert0,vert2-vert0);
    if(magnitudeSquared(perp)<fe::tol*fe::tol)
    {
#if FE_PTN_DEBUG
        feLog("cross mag %.6G below tolerance %.6G\n",magnitude(perp),fe::tol);
#endif
        return T(-1);
    }
    const Vector<3,T> normal=unit(cross(vert1-vert0,vert2-vert0));
    const Vector<3,T> to_origin=origin-vert0;
    const T displacement=dot(normal,to_origin);
    intersection=origin-normal*displacement;

    barycenter.solve(vert0,vert1,vert2,intersection);

#if FE_PTN_DEBUG
    feLog("normal %s\n",c_print(normal));
    feLog("intersection %s\n",c_print(intersection));
    feLog("barycenter %s\n",c_print(barycenter));
#endif

    //* clamp to triangle
    barycenter.clamp(vert0,vert1,vert2);

    intersection=location(barycenter,vert0,vert1,vert2);

    direction=intersection-origin;
    T range=magnitude(direction);
    if(range>T(0))
    {
        direction*=T(1)/range;
    }

#if FE_PTN_DEBUG
    feLog("intersection %s\n",c_print(intersection));
    feLog("direction %s\n",c_print(direction));
    feLog("barycenter %s\n",c_print(barycenter));
    feLog("range %.6G\n",range);
#endif

    return range;
}

} /* namespace ext */
} /* namespace fe */

#endif /* __geometry_PointTriangleNearest_h__ */