Barycenter.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 __math_Barycenter_h__
#define __math_Barycenter_h__

namespace fe
{

/**************************************************************************//**
    @brief Barycentric coordinates for a triangle

    @ingroup math

    b2 is not stored, but implied as (1 - b0 - b1)

    Usually, b0 >= 0, b1 >= 0, and b0 + b1 <= 1,
    but this is not a restriction.

    coincident point = vert0 * b0 + vert1 * b1 + vert2 * b2
*//***************************************************************************/
template <typename T>
class Barycenter: public Vector<2,T>
{
    public:
                    Barycenter(void)                                        {}

                    Barycenter(T a_b0, T a_b1)
                    {   set(*this,a_b0,a_b1); }

                    Barycenter(T a_b0, T a_b1, T a_b2)
                    {   set(*this,a_b0,a_b1,a_b2); }

                    Barycenter(const Barycenter<T>& a_rOther):
                        Vector<2,T>()
                    {   Vector<2,T>::operator=(a_rOther); }

                    Barycenter(const Vector<3,T>& a_rOther):
                        Vector<2,T>()
                    {   set(*this,a_rOther[0],a_rOther[1]); }

        Barycenter& operator=(const Barycenter<T>& a_rOther)
                    {   Vector<2,T>::operator=(a_rOther);
                        return *this; }

        Barycenter& operator=(const Vector<3,T>& a_rOther)
                    {   Vector<2,T>::operator=(a_rOther);
                        return *this; }

                    //* TODO verify using RayTriangleIntersect
                    //* I think their U and V are normalized weights
                    //* toward vert1 and vert2 from vert0.
        Barycenter& setUV(T a_u,T a_v)
                    {   set(*this,T(1)-a_u-a_v,a_u);
                        return *this; }

        operator Vector<3,T>(void) const
                    {   return Vector<3,T>((*this)[0],(*this)[1],
                                T(1)-(*this)[0]-(*this)[1]); }

        void        solve(
                        const Vector<3,T>& a_vert0,
                        const Vector<3,T>& a_vert1,
                        const Vector<3,T>& a_vert2,
                        const Vector<3,T>& a_point);

                    //* clamp to triangle
        void        clamp(
                        const Vector<3,T>& a_vert0,
                        const Vector<3,T>& a_vert1,
                        const Vector<3,T>& a_vert2);

    private:

static  T           along( const Vector<3,T>& start,
                            const Vector<3,T>& end,
                            const Vector<3,T>& point);
};

template<class T, typename U, typename V, typename W>
inline Barycenter<T>& set(Barycenter<T> &a_lhs, U a_b0, V a_b1, W a_b2)
{
    FEASSERT(fabs(a_b0+a_b1+a_b2-T(1))<T(0.001));
    set(a_lhs,a_b0,a_b1);
    return a_lhs;
}

template<int N,class T, typename U, typename V, typename W>
inline Vector<N,T> location(const Barycenter<T>& a_barycenter,
    const Vector<N,U>& a_vert0,
    const Vector<N,V>& a_vert1,
    const Vector<N,W>& a_vert2)
{
    return a_vert0*a_barycenter[0]+a_vert1*a_barycenter[1]+
            a_vert2*(T(1)-a_barycenter[0]-a_barycenter[1]);
}

template<class T, typename U>
inline U blend(const Barycenter<T>& a_barycenter,
    U a_value0, U a_value1, U a_value2)
{
    return a_value0*a_barycenter[0]+a_value1*a_barycenter[1]+
            a_value2*(T(1)-a_barycenter[0]-a_barycenter[1]);
}

template <typename T>
inline void Barycenter<T>::solve(
    const Vector<3,T>& a_vert0,
    const Vector<3,T>& a_vert1,
    const Vector<3,T>& a_vert2,
    const Vector<3,T>& a_point)
{
    const T A=a_vert0[0]-a_vert2[0];
    const T B=a_vert1[0]-a_vert2[0];
    const T C=a_vert2[0]-a_point[0];

    const T D=a_vert0[1]-a_vert2[1];
    const T E=a_vert1[1]-a_vert2[1];
    const T F=a_vert2[1]-a_point[1];

    const T G=a_vert0[2]-a_vert2[2];
    const T H=a_vert1[2]-a_vert2[2];
    const T I=a_vert2[2]-a_point[2];

    const T d0=(A*(E+H)-B*(D+G));
    const T d1=(B*(D+G)-A*(E+H));

    const T b0=(d0==0.0)? 0.0: (B*(F+I)-C*(E+H))/d0;
    const T b1=(d1==0.0)? 0.0: (A*(F+I)-C*(D+G))/d1;

    set(*this,b0,b1);
}

//* TODO can this be done without verts simply in unit space?
template <typename T>
inline void Barycenter<T>::clamp(
    const Vector<3,T>& a_vert0,
    const Vector<3,T>& a_vert1,
    const Vector<3,T>& a_vert2)
{
    const Vector<3,T> point=location(*this,a_vert0,a_vert1,a_vert2);

    if((*this)[0]<T(0))
    {
        (*this)[0]=T(0);
        (*this)[1]=along(a_vert2,a_vert1,point);
    }
    else if((*this)[1]<T(0))
    {
        (*this)[0]=along(a_vert2,a_vert0,point);
        (*this)[1]=T(0);
    }
    else
    {
        const T sum=(*this)[0]+(*this)[1];
        if(sum>T(1))
        {
            (*this)[0]=along(a_vert1,a_vert0,point);
            (*this)[1]=T(1)-(*this)[0];
        }
    }
}

template <typename T>
inline T Barycenter<T>::along(const Vector<3,T>& start,
        const Vector<3,T>& end,
        const Vector<3,T>& point)
{
    const Vector<3,T> edge=end-start;
    const T mag2=magnitudeSquared(edge);
    if(mag2==T(0))
    {
        return T(0);
    }

    T result=dot(edge,point-start)/mag2;
    if(result<T(0))
    {
        return T(0);
    }
    if(result>T(1))
    {
        return T(1);
    }
    return result;
}

typedef Barycenter<Real>    SpatialBary;

} /* namespace */

#endif /* __math_Barycenter_h__ */