RaySphereIntersect.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_RaySphereIntersect_h__
#define __geometry_RaySphereIntersect_h__

namespace fe
{
namespace ext
{

/**************************************************************************//**
    @brief Find intersection between ray and sphere

    @ingroup geometry
*//***************************************************************************/
template <typename T>
class RaySphereIntersect
{
    public:
static  T       solve(const Vector<3,T>& center,const T radius,
                        const Vector<3,T>& origin,
                        const Vector<3,T>& direction);
static  void    resolveContact(const Vector<3,T>& center,const T radius,
                        const Vector<3,T>& origin,
                        const Vector<3,T>& direction,
                        const T range,
                        Vector<3,T>& intersection,
                        Vector<3,T>& normal);
};

// direction must be normalized
template <typename T>
inline T RaySphereIntersect<T>::solve(const Vector<3,T>& center,T radius,
        const Vector<3,T>& origin, const Vector<3,T>& direction)
{
//  feLog("center %s  radius %.6G\n",c_print(center),radius);
//  feLog("origin: %s\n",c_print(origin));
//  feLog("direction: %s\n",c_print(direction));

    const Vector<3,T> to_origin=origin-center;
    const T b=T(2)*dot(direction,to_origin);
    const T c=dot(to_origin,to_origin)-radius*radius;
    const T i=b*b-T(4)*c;

//  feLog("b %.6G c %.6G i %.6G\n",b,c,i);

    if(i<T(0))
    {
        return T(-1);
    }
#if FALSE
    const T sq=sqrtf(i);
    const T t1=(-b+sq);
    const T t2=(-b-sq);

//  feLog("t1 %.6G t2 %.6G\n",t1,t2);

    T range=T(0.5)*((t1<t2)? t1: t2);
#else
    T range=T(-0.5)*(b+sqrtf(i));
    if(range<0.0)
    {
        range=T(-0.5)*(b-sqrtf(i));
    }
#endif

    return range;
}

// direction must be normalized
template <typename T>
inline void RaySphereIntersect<T>::resolveContact(
        const Vector<3,T>& center,const T radius,
        const Vector<3,T>& origin,
        const Vector<3,T>& direction,
        const T range,
        Vector<3,T>& intersection,
        Vector<3,T>& normal)
{
    FEASSERT(radius>0.0);
    intersection=origin+direction*range;
    normal=(intersection-center)/radius;
}

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

#endif /* __geometry_RaySphereIntersect_h__ */