SphericalTriangle.h

 
#pragma once
 
#if defined(SphericalTriangle_RECURSES)
#error Recursive header files inclusion detected in SphericalTriangle.h
#else // defined(SphericalTriangle_RECURSES)
#define SphericalTriangle_RECURSES
 
#if !defined SphericalTriangle_h
#define SphericalTriangle_h
 
// Inclusions
#include <iostream>
#include <map>
#include "DGtal/base/Common.h"
 
 
namespace DGtal
{
 
  // class SphericalTriangle
  template <typename TSpace>
  class SphericalTriangle
  {
    BOOST_CONCEPT_ASSERT(( concepts::CSpace<TSpace> ));
    typedef TSpace                         Space;
    typedef SphericalTriangle<Space>       Self;
    typedef typename Space::RealPoint      RealPoint;
    typedef typename Space::RealVector     RealVector;
    typedef typename RealVector::Component Scalar;
    
    // Checks that dimension is 3.
    BOOST_STATIC_ASSERT(( Space::dimension == 3 ));
    
    // ----------------------- Standard services ------------------------------
  public:
    
    ~SphericalTriangle() {}
    
    SphericalTriangle( const RealVector& va, const RealVector& vb, const RealVector& vc,
                       bool normalize = true )
    {
      setA( va, normalize );
      setB( vb, normalize );
      setC( vc, normalize );
    }
    
    SphericalTriangle ( const SphericalTriangle & other ) = default;
    
    SphericalTriangle & operator= ( const SphericalTriangle & other ) = default;
 
    const RealVector& A() const { return myA; }
    const RealVector& B() const { return myB; }
    const RealVector& C() const { return myC; }
 
    void setA( const RealVector& va, bool normalize = true )
    {
      myA = va;
      if ( normalize )
        {
          Scalar n = myA.norm();
          if ( fabs( n ) > 1e-8 ) myA /= n;
          else myA = RealVector::zero;
        }
    }
 
    void setB( const RealVector& vb, bool normalize = true )
    {
      myB = vb;
      if ( normalize )
        {
          Scalar n = myB.norm();
          if ( fabs( n ) > 1e-8 ) myB /= n;
          else myB = RealVector::zero;
        }
    }
    void setC( const RealVector& vc, bool normalize = true )
    {
      myC = vc;
      if ( normalize )
        {
          Scalar n = myC.norm();
          if ( fabs( n ) > 1e-8 ) myC /= n;
          else myC = RealVector::zero;
        }
    }
 
    bool isDegenerate() const
    {
      Scalar d[ 3 ] = { ( myA - myB ).norm(),
                        ( myA - myC ).norm(),
                        ( myB - myC ).norm() };
      // Checks that the spherical triangle is small or thin.
      if ( ( d[ 0 ] < 1e-8 ) || ( d[ 1 ] < 1e-8 ) || ( d[ 2 ] < 1e-8 ) )
        return true;
      // Checks that the spherical triangle is flat.
      Dimension m = 0;
      if ( d[ 1 ] > d[ m ] ) m = 1;
      if ( d[ 2 ] > d[ m ] ) m = 2;
      return ( fabs( d[ m ] - d[ (m+1)%3 ] - d[ (m+2)%3 ] ) < 1e-8 );
    }
    
    Self polarTriangle() const
    {
      auto Ap = myB.crossProduct(myC);
      auto Bp = myC.crossProduct(myA);
      auto Cp = myA.crossProduct(myB);
      // Reorient points.
      if ( Ap.dot( myA ) < 0.0 ) Ap = -Ap;
      if ( Bp.dot( myB ) < 0.0 ) Bp = -Bp;
      if ( Cp.dot( myC ) < 0.0 ) Cp = -Cp;
      return Self( Ap, Bp, Cp, true );
    }
 
    void interiorAngles( Scalar& alpha, Scalar& beta, Scalar& gamma ) const
    {
      Self    T = polarTriangle();
      if ( T.A() == RealVector::zero || T.B() == RealVector::zero || T.C() == RealVector::zero )
        alpha = beta = gamma = 0.0;
      else
        {
          Scalar ca = std::max( -1.0, std::min( 1.0, T.B().dot( T.C() ) ) );
          Scalar cb = std::max( -1.0, std::min( 1.0, T.C().dot( T.A() ) ) );
          Scalar cc = std::max( -1.0, std::min( 1.0, T.A().dot( T.B() ) ) );
          alpha     = acos( ca );
          beta      = acos( cb );
          gamma     = acos( cc );
        }
    }
 
    Scalar area() const
    {
      Scalar alpha, beta, gamma;
      if ( isDegenerate() ) return 0.0;
      interiorAngles( alpha, beta, gamma );
      return ( (alpha == 0.0) || (beta == 0.0) || (gamma == 0.0) )
               ? 0.0 : 2.0*M_PI - alpha - beta - gamma;
    }
 
    Scalar algebraicArea() const
    {
      Scalar     S = area();
      RealVector M = myA + myB + myC;
      RealVector X = ( myB - myA ).crossProduct( myC - myA );
      if ( M.norm1() <= 1e-8 || X.norm1() <= 1e-8 ) return 0.0;
      return M.dot( X ) < 0.0 ? -S : S;
    }
    
    // ------------------------- Private Datas --------------------------------
  private:
    // ------------------------- Hidden services ------------------------------
  protected:
    RealVector myA;
    RealVector myB;
    RealVector myC;
    
    // ------------------------- Internals ------------------------------------
  private:
    
  }; // end of class SphericalTriangle
  
 
} // namespace DGtal
 
 
// Includes inline functions.
 
//                                                                           //
 
#endif // !defined SphericalTriangle_h
 
#undef SphericalTriangle_RECURSES
#endif // else defined(SphericalTriangle_RECURSES)