MRInSphere.h

Go to the documentation of this file.

#pragma once
 
#include "MRPrecisePredicates3.h"
#include "MRFastInt.h"
#include "MRPch/MRBindingMacros.h"
#include <array>
#include <cassert>
#include <type_traits>
#include <utility>
 
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;
    }

    void flip()
    {
        assert( E >= 0 ); 
        std::swap( u, v );
        w = -w;
    }

    [[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 );

    void flip()
    {
        assert( E >= 0 ); 
        std::swap( u, v );
        w = -w;
    }

    [[nodiscard]] MRMESH_API InSphereResult operator()( const Vector3i & d ) const;
 
protected:
    Vector3i a;           
    Vector3i64 u, v;      
    Vector3i64 w;         
    FastInt<192> W;       
    std::array<FastInt<192>, 3> M; 
    FastInt<384> E = -1;  
    std::int64_t rSq = 0; 
};
 
using InSphereTesterf = InSphereTester<float>;
using InSphereTesterd = InSphereTester<double>;
using InSphereTesteri = InSphereTester<int>;

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

    void flip()
    {
        if ( degenerateTriangle_ )
        {
            std::swap( u, v ); 
            pairSigma_ = -pairSigma_; 
        }
        else
            InSphereTester<int>::flip();
        std::swap( vb_, vc_ );
    }

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

    [[nodiscard]] bool degenerateTriangle() const { return degenerateTriangle_; }


    [[nodiscard]] MR_BIND_IGNORE const Vector3i64 & normal() const { return w; }

    [[nodiscard]] MR_BIND_IGNORE const FastInt<192> & normalSq() const { return W; }

    [[nodiscard]] MR_BIND_IGNORE const FastInt<384> & heightSq() const { return E; }
 
protected:
    Vector3i64 pairV_;    
    Vector3i pairPt_;     
    VertId va_, vb_, vc_; 
    int pairSigma_ = 0;   
    bool degenerateTriangle_ = false; 
};

class FastInSphereTesterSoS : public InSphereTesterSoS
{
public:
    MRMESH_API bool reset( const PreciseVertCoords & a, const PreciseVertCoords & b, const PreciseVertCoords & c, std::int64_t rSq );

    void flip()
    {
        InSphereTesterSoS::flip();
        hn_ = -hn_;
    }

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

    [[nodiscard]] MRMESH_API bool outsideBothSpheres( const Vector3i & d ) const;
 
private:
    Vector3d cc_, hn_;

    double tol_ = 0;
};

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 ) );

 
}