CircleFrom3Points.ih

#include <cstdlib>

// IMPLEMENTATION of inline methods.

// ----------------------- Standard services ------------------------------

template <typename TPoint>
inline
DGtal::CircleFrom3Points<TPoint>::~CircleFrom3Points()
{
}

template <typename TPoint>
inline
DGtal::CircleFrom3Points<TPoint>::CircleFrom3Points()
{
}

template <typename TPoint>
inline
void
DGtal::CircleFrom3Points<TPoint>::init(
        const Point& aFirstPoint, 
        const Point& aSecondPoint,
        const Point& aThirdPoint)
{
  myP = aFirstPoint;
  myQ = aSecondPoint;
  myR = aThirdPoint;
}


template <typename TPoint>
inline
DGtal::CircleFrom3Points<TPoint>::CircleFrom3Points(
        const Point& aFirstPoint, 
        const Point& aSecondPoint,
        const Point& aThirdPoint):
        myP(aFirstPoint),
        myQ(aSecondPoint),
        myR(aThirdPoint)
{
}


template <typename TPoint>
inline
DGtal::CircleFrom3Points<TPoint>::CircleFrom3Points(
        const CircleFrom3Points & other):
        myP(other.myP),
        myQ(other.myQ),
        myR(other.myR)
{
}


template <typename TPoint>
inline
DGtal::CircleFrom3Points<TPoint>&
DGtal::CircleFrom3Points<TPoint>::operator=(
        const CircleFrom3Points & other)
{
  myP = other.myP;
  myQ = other.myQ;
  myR = other.myR;
  return *this;
}



template <typename TPoint>
inline
typename DGtal::CircleFrom3Points<TPoint>::Distance
DGtal::CircleFrom3Points<TPoint>::signedDistance(const Point& aP) const
{
  Vector u( (myP[0]-aP[0])*(myR[1]-aP[1])-(myR[0]-aP[0])*(myP[1]-aP[1]), 
                  (myR[0]-aP[0])*(myR[0]-myP[0])+(myR[1]-aP[1])*(myR[1]-myP[1]) );
  Vector v( (myP[0]-aP[0])*(myQ[1]-aP[1])-(myQ[0]-aP[0])*(myP[1]-aP[1]), 
                  (myQ[0]-aP[0])*(myQ[0]-myP[0])+(myQ[1]-aP[1])*(myQ[1]-myP[1]) );

  return -( (u[0] * v[1]) - (u[1] * v[0]) );
}

template <typename TPoint>
inline
void
DGtal::CircleFrom3Points<TPoint>::getParameters(double& cx, double& cy, double& rr) const
{
  Vector u(myP - myQ);
  Vector v(myR - myP);

  Coordinate area = (u[0] * v[1]) - (u[1] * v[0]);

  if (area == 0) {
    std::cerr << "Infinite radius detected in method getParameters() of CircleFrom3Points" << std::endl;
    throw InfiniteNumberException();
  } else {

    //center
    double areaD = NumberTraits<Coordinate>::castToDouble(area)*2.0; 
    double pz, qz, rz;
    pz = std::pow( NumberTraits<Coordinate>::castToDouble(myP[0]),2.0 )/areaD 
        + std::pow( NumberTraits<Coordinate>::castToDouble(myP[1]),2.0 )/areaD;
    qz = std::pow( NumberTraits<Coordinate>::castToDouble(myQ[0]),2.0 )/areaD 
        + std::pow( NumberTraits<Coordinate>::castToDouble(myQ[1]),2.0 )/areaD;
    rz = std::pow( NumberTraits<Coordinate>::castToDouble(myR[0]),2.0 )/areaD 
        + std::pow( NumberTraits<Coordinate>::castToDouble(myR[1]),2.0 )/areaD;

    cx = NumberTraits<Coordinate>::castToDouble(myP[1])*(qz-rz)
        -NumberTraits<Coordinate>::castToDouble(myQ[1])*(pz-rz)
        +NumberTraits<Coordinate>::castToDouble(myR[1])*(pz-qz);
    cy = -NumberTraits<Coordinate>::castToDouble(myP[0])*(qz-rz)
         +NumberTraits<Coordinate>::castToDouble(myQ[0])*(pz-rz)
         -NumberTraits<Coordinate>::castToDouble(myR[0])*(pz-qz);

    //radius
    double a, b, c; 
    a = std::pow( NumberTraits<Coordinate>::castToDouble(myQ[0]-myP[0]),2.0 )
       +std::pow( NumberTraits<Coordinate>::castToDouble(myQ[1]-myP[1]),2.0 );
    b = std::pow( NumberTraits<Coordinate>::castToDouble(myR[0]-myQ[0]),2.0 )
       +std::pow( NumberTraits<Coordinate>::castToDouble(myR[1]-myQ[1]),2.0 );
    c = std::pow( NumberTraits<Coordinate>::castToDouble(myR[0]-myP[0]),2.0 )
       +std::pow( NumberTraits<Coordinate>::castToDouble(myR[1]-myP[1]),2.0 );

    rr = ( ( std::sqrt(a)/std::abs(areaD) )*std::sqrt(b)*std::sqrt(c) );

  }
}

template <typename TPoint>
inline
double
DGtal::CircleFrom3Points<TPoint>::getCurvature() const
{
  Vector u(myP - myQ);
  Vector v(myR - myP);

  //4 times the area of the triangle 
  Coordinate area = (u[0] * v[1]) - (u[1] * v[0]);
  double areaD = NumberTraits<Coordinate>::castToDouble(area)*2.0; 
  
  //length of the sides of the triangle 
  double a, b, c; 
  a = std::pow( NumberTraits<Coordinate>::castToDouble(myQ[0]-myP[0]),2.0 )
     +std::pow( NumberTraits<Coordinate>::castToDouble(myQ[1]-myP[1]),2.0 );
  b = std::pow( NumberTraits<Coordinate>::castToDouble(myR[0]-myQ[0]),2.0 )
     +std::pow( NumberTraits<Coordinate>::castToDouble(myR[1]-myQ[1]),2.0 );
  c = std::pow( NumberTraits<Coordinate>::castToDouble(myR[0]-myP[0]),2.0 )
     +std::pow( NumberTraits<Coordinate>::castToDouble(myR[1]-myP[1]),2.0 );

  //curvature
  return -( areaD/ (std::sqrt(a)*std::sqrt(b)*std::sqrt(c)) );
}


// Interface - public :

template <typename TPoint>
inline
std::string
DGtal::CircleFrom3Points<TPoint>::className() const
{
  return "CircleFrom3Points";
}

template <typename TPoint>
inline
void
DGtal::CircleFrom3Points<TPoint>::selfDisplay ( std::ostream & out ) const
{
  out << "[CircleFrom3Points] passing through:\n";
  out << myP << myQ << myR;
}

template <typename TPoint>
inline
bool
DGtal::CircleFrom3Points<TPoint>::isValid() const
{
    return true;
}