MRInSphere.h

Go to the documentation of this file.

#pragma once
 
#include "MRPrecisePredicates3.h"
#include "MRHighPrecision.h"
#include "MRPch/MRBindingMacros.h"
#include <cassert>
#include <type_traits>
 
namespace MR
{
 


enum class InSphereResult
{
    NoSphere,
    Outside,
    OnSphere,
    Inside
};

[[nodiscard]] MRMESH_API InSphereResult inSphere( const Vector3i & a, const Vector3i & b, const Vector3i & c,
    const Vector3i & d, std::int64_t rSq );

[[nodiscard]] MRMESH_API InSphereResult inSphere( const std::array<PreciseVertCoords, 4> & vs, std::int64_t rSq );

template <typename T>
class InSphereTester
{
    static_assert( std::is_floating_point_v<T> );
public:
    bool reset( const Vector3<T> & va, const Vector3<T> & vb, const Vector3<T> & vc, T sqRadius )
    {
        a = va;
        rSq = sqRadius;
        E = -1;
        u = vb - va;
        v = vc - va;

        const T uu = u.lengthSq();
        const T vv = v.lengthSq();
        if ( uu > 4 * rSq || vv > 4 * rSq || ( v - u ).lengthSq() > 4 * rSq )
            return false;
 
        w = cross( u, v );
        W = w.lengthSq();
        if ( W <= 0 )
            return false; 

        const T uv = dot( u, v );
        M = ( vv * ( uu - uv ) ) * u + ( uu * ( vv - uv ) ) * v;

        E = 4 * rSq * W * W - M.lengthSq();
        return E >= 0;
    }

    [[nodiscard]] InSphereResult operator()( const Vector3<T> & d ) const
    {
        assert( E >= 0 ); 
        const auto q = d - a;

        const T qq = q.lengthSq();
        if ( qq > 4 * rSq )
            return InSphereResult::Outside;
 
        const T A = W * qq - dot( q, M ); 
        const T t = dot( q, w ); 

        if ( A < 0 && t >= 0 )
            return InSphereResult::Inside;
        if ( A >= 0 && t <= 0 )
            return ( A == 0 && ( t == 0 || E == 0 ) ) ? InSphereResult::OnSphere : InSphereResult::Outside;
        const T lhs = A * A * W;
        const T rhs = E * t * t;
        if ( lhs == rhs )
            return InSphereResult::OnSphere;
        return ( A < 0 ) == ( lhs > rhs ) ? InSphereResult::Inside : InSphereResult::Outside;
    }
 
private:
    Vector3<T> a;    
    Vector3<T> u, v; 
    Vector3<T> w;    
    T W = 0;         
    Vector3<T> M;    
    T E = -1;        
    T rSq = 0;       
};
 

template <>
class InSphereTester<int>
{
public:
    MRMESH_API bool reset( const Vector3i & a, const Vector3i & b, const Vector3i & c, std::int64_t rSq );

    [[nodiscard]] MRMESH_API InSphereResult operator()( const Vector3i & d ) const;

    [[nodiscard]] MRMESH_API InSphereResult operator()( const std::array<PreciseVertCoords, 4> & vs ) const;
 
private:
    Vector3i a;           
    Vector3i64 u, v;      
    Vector3i64 w;         
    Int256 W;             
    Vector3i256 M;        
    Int512 E = -1;        
    std::int64_t rSq = 0; 
};
 
using InSphereTesterf = InSphereTester<float>;
using InSphereTesterd = InSphereTester<double>;
using InSphereTesteri = InSphereTester<int>;

template <typename T>
[[nodiscard]] std::enable_if_t<std::is_floating_point_v<T>, InSphereResult> inSphere( const Vector3<T> & a, const Vector3<T> & b, const Vector3<T> & c,
    const Vector3<T> & d, T rSq )
{
    InSphereTester<T> tester;
    if ( !tester.reset( a, b, c, rSq ) )
        return InSphereResult::NoSphere;
    return tester( d );
}
 
MR_BIND_TEMPLATE( InSphereResult inSphere( const Vector3f & a, const Vector3f & b, const Vector3f & c, const Vector3f & d, float rSq ) );
MR_BIND_TEMPLATE( InSphereResult inSphere( const Vector3d & a, const Vector3d & b, const Vector3d & c, const Vector3d & d, double rSq ) );

 
}