RayCylinderIntersect.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 __solve_RayCylinderIntersect_h__
#define __solve_RayCylinderIntersect_h__

namespace fe
{
namespace ext
{

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

    @ingroup solve

    "Intersecting a Ray with a Cylinder"
    by Joseph M. Cychosz and Warren N. Waggenspack, Jr.,
    (3ksnn64@ecn.purdue.edu, mewagg@mewnw.dnet.lsu.edu)
    in "Graphics Gems IV", Academic Press, 1994
*//***************************************************************************/
template <typename T>
class RayCylinderIntersect
{
    public:
static  T       solve(const Vector<3,T>& base,const Vector<3,T> axis,
                        const T radius,
                        const Vector<3,T>& origin,
                        const Vector<3,T>& direction,
                        T& along,
                        Vector<3,T>& intersection);
static  void    resolveContact(const Vector<3,T>& base,const Vector<3,T> axis,
                        const T radius,
                        const Vector<3,T>& origin,
                        const Vector<3,T>& direction,
                        const T range,const T along,
                        const Vector<3,T>& intersection,
                        Vector<3,T>& normal);
};

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

    const Vector<3,T> RC=origin-base;
    Vector<3,T> n;
    cross3(n,direction,axis);
    const T nlength=magnitude(n);
    if(nlength<fe::tol)
    {
        return T(-1);
    }

    n*=T(1)/nlength;
    const T d=fabs(dot(RC,n));  // distance from axis
    if(d>radius)
    {
        return T(-1);
    }

    Vector<3,T> tmp;
    cross3(tmp,RC,axis);
    T t= -dot(tmp,n)/nlength;
    cross3(tmp,n,axis);
    normalize(tmp);
    T range=t-fabs(sqrtf(radius*radius-d*d)/dot(direction,tmp));

    const T axislength2=magnitudeSquared(axis);
    intersection=origin+range*direction;
    along=dot(axis,intersection-base)/axislength2;
//  feLog("intersection %s along %.6G range %.6G\n",
//          c_print(intersection),along,range);

    return (along>0 && along<1.0f)? range: T(-1);
}

// direction must be normalized
template <typename T>
inline void RayCylinderIntersect<T>::resolveContact(
        const Vector<3,T>& base,const Vector<3,T> axis,T radius,
        const Vector<3,T>& origin,const Vector<3,T>& direction,
        const T range,const T along,
        const Vector<3,T>& intersection,Vector<3,T>& normal)
{
    Vector<3,T> centerAlong=base+axis*along;
    normal=(intersection-centerAlong)/radius;
}

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

#endif /* __solve_RayCylinderIntersect_h__ */