Ellipse2D.ih

#include <cstdlib>

// IMPLEMENTATION of inline methods.

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

template <typename T>
inline
DGtal::Ellipse2D<T>::~Ellipse2D()
{
}

template <typename T>
inline
DGtal::Ellipse2D<T>::Ellipse2D(const double x0, const double y0,
             const double a0, const double a1, const double theta)
  : myCenter(x0,y0), myAxis1(a0),myAxis2(a1),myTheta(theta)
{}


template <typename T>
inline
DGtal::Ellipse2D<T>::Ellipse2D(const RealPoint2D &aPoint,
             const double a0, const double a1, const double theta)
  : myCenter(aPoint), myAxis1(a0),myAxis2(a1),myTheta(theta)
{}

template <typename T>
inline
DGtal::Ellipse2D<T>::Ellipse2D(const Point &aPoint,
             const double a0, const double a1, const double theta) 
  : myAxis1(a0),myAxis2(a1),myTheta(theta)
{
  myCenter = aPoint;
}

// ------------- Implementation of 'StarShaped' services ------------------

template <typename T>
inline
double
DGtal::Ellipse2D<T>::parameter( const RealPoint2D & pp ) const
{
  RealPoint2D v2d( pp );
  v2d -= myCenter;
  
  double angle; 

  if ( v2d[0] == 0.0 ) 
    {
      if ( v2d[1] >0 )
        angle = M_PI/2.0; 
      else 
        angle = 1.5*M_PI; 
    }
  else if (  ( v2d[0] > 0.0 ) && (   v2d[1] >= 0.0 ) )
    angle = atan(v2d[1]/v2d[0]);
  else if (  ( v2d[0] > 0.0 ) && (   v2d[1] <= 0.0 ) )
    angle = 2*M_PI + atan(v2d[1]/v2d[0]);
  else if (  ( v2d[0] < 0.0 ) && (   v2d[1] >= 0.0 ) )
    angle = atan(v2d[1]/v2d[0]) + M_PI;
  else // (  ( v2d[0] < 0.0 ) && (   v2d[1] <= 0.0 ) )
    angle = atan(v2d[1]/v2d[0]) + M_PI;

   return angle;
}

template <typename T>
inline
typename DGtal::Ellipse2D<T>::RealPoint2D 
DGtal::Ellipse2D<T>::x( double t ) const
{
  double a2 = myAxis1*myAxis1;
  double b2 = myAxis2*myAxis2;
  double costth = cos( t - myTheta );
  //  double sintth = sin( t - myTheta );
  double cost = cos( t );
  double sint = sin( t );
  double rho = myAxis2 / sqrt( 1.0 - ((a2-b2)/a2)*costth*costth);
    //myAxis2*myAxis1 / sqrt( a2 - (a2-b2)*costth*costth);
  RealPoint2D v( rho*cost, 
     rho*sint );
  v += myCenter;
  return v;
}


template <typename T>
inline
typename DGtal::Ellipse2D<T>::RealVector2D 
DGtal::Ellipse2D<T>::xp( const double t ) const
{
  double a2 = myAxis1*myAxis1;
  double b2 = myAxis2*myAxis2;
  double costth = cos( t - myTheta );
  double sintth = sin( t - myTheta );
  double cost = cos( t );
  double sint = sin( t );
  double rho = myAxis2 / sqrt( 1.0 - ((a2-b2)/a2)*costth*costth);
  double a = myAxis1; 
  double b = myAxis2;
  double rhod= a*b*(b2-a2)*sintth*costth
    / std::pow( a2*sintth*sintth + b2*costth*costth, 1.5 );
  RealPoint2D v( rhod*cost - rho*sint, rhod*sint + rho*cost );
  
  return v;
}

template <typename T>
inline
typename DGtal::Ellipse2D<T>::RealVector2D
DGtal::Ellipse2D<T>::xpp( const double t ) const
{
  double a2 = myAxis1*myAxis1;
  double b2 = myAxis2*myAxis2;
  double costth = cos( t - myTheta );
  double sintth = sin( t - myTheta );
  double cost = cos( t );
  double sint = sin( t );
  double rho = myAxis2 / sqrt( 1.0 - ((a2-b2)/a2)*costth*costth);

  double a = myAxis1; 
  double b = myAxis2;
  double rhod = a*b*(b2-a2)*sintth*costth
    / std::pow( a2*sintth*sintth + b2*costth*costth, 1.5 );
  double rhodd = a*b*(b2-a2)
    / std::pow( a2*sintth*sintth + b2*costth*costth, 2.5 )
    * ( (costth*costth - sintth*sintth) * (a2*sintth*sintth + b2*costth*costth)
  + 3.0*(b2-a2)*sintth*sintth*costth*costth );
  
  RealPoint2D v( rhodd*cost - 2.0*rhod*sint - rho*cost, 
     rhodd*sint + 2.0*rhod*cost - rho*sint );
  return v;
}


// Interface - public :

template <typename T>
inline
void
DGtal::Ellipse2D<T>::selfDisplay ( std::ostream & out ) const
{
  out << "[Ellipse2D] center= "<<myCenter<<" big axis="<<myAxis1
      <<" small axis="<<myAxis2<<" phase="<<myTheta;
}

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



// Implementation of inline functions                                        //

template <typename T>
inline
std::ostream&
DGtal::operator<< ( std::ostream & out, 
      const Ellipse2D<T> & object )
{
  object.selfDisplay( out );
  return out;
}

//                                                                           //