MeasureOfStraightLines.ih

// IMPLEMENTATION of inline methods.

#include <cstdlib>
#include <vector>
#include <cmath>

#include <boost/math/special_functions/asinh.hpp>


// Internal functions                                                        //


inline
DGtal::MeasureOfStraightLines::MeasureOfStraightLines()
{
  myEpsilon = 0.0005;
}


inline
DGtal::MeasureOfStraightLines:: ~MeasureOfStraightLines()
{
}

inline
void 
DGtal::MeasureOfStraightLines::selfDisplay( std::ostream & out ) const
{
  out << "[MeasureOfStraightLines]";
}

inline
bool 
DGtal::MeasureOfStraightLines::isValid() const
{
  return true;
}

inline double 
DGtal::MeasureOfStraightLines::computeMeasureEdge ( double a0,double b0, double a1, double b1 )
{
  //double measure=0;

  //Aligned with the 'b' axis
  if ( ( b0 == 0 ) && ( b1 == 0 ) )
    return 0;
  //Aligned with the 'b' axis
  if ( ( a0 == 0 ) && ( a1 == 0 ) )
    return 0;

  double delta = ( a0*b1 - a1*b0 );

  if ( a0 == 0 )
    return delta * ( 1.0/ ( 1.0 + sqrt ( 1 + pow ( a1,2 ) ) ) );
  if ( a1 == 0 )
    return delta * ( ( -1.0 + sqrt ( 1.0 + pow ( a0,2 ) ) ) /pow ( a0,2 ) );

  if ( a0 != a1 )
    {
      return  delta * ( ( a0 - a1 + sqrt ( 1 + pow ( a0,2 ) ) *a1 - a0*sqrt ( 1 + pow ( a1,2 ) ) ) / ( pow ( a0,2 ) *a1 - a0*pow ( a1,2 ) ) );
    }
  else
    return  delta * a1/ ( a0* ( 1 + pow ( a1,2 ) + sqrt ( 1 + pow ( a1,2 ) ) ) );
}


inline void 
DGtal::MeasureOfStraightLines::setEpsilon(const double aValue)
{
  myEpsilon = aValue;
}

// Implementation of inline methods                                          //

inline double 
DGtal::MeasureOfStraightLines::computeMeasure(const std::vector<double> &a,
                                              const std::vector<double> &b)
{
  double measure = 0;

  ASSERT(a.size() == b.size());
  
  for (unsigned int i=0 ; i < a.size() ; i++)
    {
      measure += computeMeasureEdge ( a[i], b[i], 
                                      a[ ( i+1 ) % a.size()],b[ ( i+1 ) %a.size()] );
    } 

  return measure;
}

inline double 
DGtal::MeasureOfStraightLines::computeCentroidA(const std::vector<double> &a,
                                                const std::vector<double> &b)
{
  double measure = 0;
  double C_a = 0;

  ASSERT(a.size() == b.size());
  
  for (unsigned int i=0 ; i < a.size() ; i++)
    {
      measure += computeMeasureEdge( a[i], b[i], 
                                     a[ ( i+1 ) % a.size()],b[ ( i+1 ) %a.size()] );
      C_a += computeCentroidEdge_a ( a[i], b[i], 
                                     a[ ( i+1 ) % a.size()],b[ ( i+1 ) %a.size()] );
    } 

  return C_a/measure;
}

inline double 
DGtal::MeasureOfStraightLines::computeCentroidB(const std::vector<double> &a,
                                                const std::vector<double> &b)
{
  double measure = 0;
  double C_b = 0;

  ASSERT(a.size() == b.size());
  
  for (unsigned int i=0 ; i < a.size() ; i++)
    {
      measure += computeMeasureEdge( a[i], b[i], 
                                     a[ ( i+1 ) % a.size()],b[ ( i+1 ) %a.size()] );
      C_b += computeCentroidEdge_b ( a[i], b[i], 
                                     a[ ( i+1 ) % a.size()],b[ ( i+1 ) %a.size()] );
    } 

  return C_b/measure;
}


inline
double 
DGtal::MeasureOfStraightLines::__computeCentroidSquare_b ( double x1, double y1, double x2,double y2 )
{
  double val;
  val = ( ( -1.0* ( sqrt ( 1.0 + pow ( x1,2 ) ) *x2 ) + x1*sqrt ( 1 + pow ( x2,2 ) ) ) * ( pow ( y1,2 ) - pow ( y2,2 ) ) ) / ( 2.0*sqrt ( ( 1.0 + pow ( x1,2 ) ) * ( 1.0 + pow ( x2,2 ) ) ) );

  ASSERT ( val>=0 );
  return val;
}


inline
int
DGtal::MeasureOfStraightLines::sign ( const double a )
{
  if ( a>=0 )
    return 1;
  else
    return -1;
}


inline
double 
DGtal::MeasureOfStraightLines::__computeCentroidEdgeApprox_b ( double a0, double b0,double a1,double b1 )
{
  double measure,y2;
  unsigned int nb_step;
  if ( a1 < a0 )
    {
      double tmp = a0;
      a0 = a1;
      a1 = tmp;
      tmp = b0;
      b0 = b1;
      b1 = tmp;
    }

  measure = 0;
  nb_step = ( unsigned int ) floor ( fabs ( a1-a0 ) /myEpsilon );
  double slope = ( b1-b0 ) / ( a1-a0 );
  double decal = b1 - ( b1-b0 ) / ( a1-a0 ) *a1;

  if ( slope == 0 )
    return __computeCentroidSquare_b ( a0, 0, a1, b0 );

  for ( unsigned int i=0; i < nb_step ; i++ )
    {
      y2 = ( b1-b0 ) / ( a1-a0 ) * ( a0+ ( i+1 ) *myEpsilon ) + decal;
      if ( y2>0 )
        measure += __computeCentroidSquare_b ( a0+i*myEpsilon, 0, a0+ ( i+1 ) *myEpsilon, y2 );
      else
        measure += __computeCentroidSquare_b ( a0+i*myEpsilon, y2, a0+ ( i+1 ) *myEpsilon, 0 );
    }
  return measure;
}

inline
double 
DGtal::MeasureOfStraightLines::__computeCentroidTriApprox_b ( double a0, double b0,double a1,double b1 )
{
  int signe = sign ( a0*b1 - a1*b0 );
  double measure = 0;

  double m_A0A1=0, m_0A0=0, m_0A1=0;
  if ( ( b0 != 0 ) && ( a0 != 0 ) )
    m_0A0 = __computeCentroidEdgeApprox_b ( 0,0,a0,b0 );
  if ( ( b1 != 0 )  && ( a1 != 0 ) )
    m_0A1 = __computeCentroidEdgeApprox_b ( 0,0,a1,b1 );
  if ( ( a1 != a0 ) )
    m_A0A1 = __computeCentroidEdgeApprox_b ( a0,b0,a1,b1 );

  ASSERT ( m_0A0>=0 );
  ASSERT ( m_0A1>=0 );
  ASSERT ( m_A0A1>=0 );

  if ( a0 < a1 )
    {
      double det = ( a1 * ( b0 - b1 ) - b1* ( a0 - a1 ) );
      if ( det >0 )
        measure = m_0A0 + m_A0A1 - m_0A1;
      else
        measure = m_0A1 - m_0A0 - m_A0A1;
    }
  else
    {
      double det = ( a0 * ( b1 - b0 ) - b0* ( a1 - a0 ) );
      if ( det >0 )
        measure = m_0A1 + m_A0A1 - m_0A0;
      else
        measure = m_0A0 - m_0A1 - m_A0A1;
    }
  ASSERT ( measure>=0 );
  return signe*measure;
}

inline
double 
DGtal::MeasureOfStraightLines::computeCentroidEdge_b ( double a0,double b0, double a1, double b1 )
{
  //double measure=0;
  double delta= ( a0*b1 - a1*b0 );

  if ( ( b0 == 0 ) && ( b1 == 0 ) )
    return 0;
  if ( ( a0 == 0 ) && ( a1 == 0 ) )
    return 0;

  if ( a0 == 0 )
    return delta* ( a1* ( ( -2 + sqrt ( 1 + a1*a1 ) ) *b0 + 2*b1 ) + ( b0 - 2*b1 ) *boost::math::asinh ( a1 ) ) / ( 2*a1*a1*a1 );
  if ( a1 == 0 )
    return delta* ( a0* ( 2*b0 + ( -2 + sqrt ( 1 + a0*a0 ) ) *b1 ) + ( -2*b0 + b1 ) *boost::math::asinh ( a0 ) ) / ( 2*a0*a0*a0 );

  if ( a0 == a1 )
    return delta* ( ( - ( ( a1* ( b0 + a1* ( -a1 + sqrt ( 1 + pow ( a1,-2 ) ) *sqrt ( 1 + pow ( a1,2 ) ) ) *
                                  b1 ) ) /sqrt ( 1.0 + pow ( a1,2 ) ) ) + ( b0 + b1 ) *boost::math::asinh ( a1 ) ) /
                    ( 2.*pow ( a1,3 ) ) );

  return __computeCentroidTriApprox_b ( a0,b0,a1,b1 );

}

inline
double 
DGtal::MeasureOfStraightLines::computeCentroidEdge_a ( double a0,double b0, double a1, double b1 )
{
  //double measure=0;
  double delta= ( a0*b1 - a1*b0 );

  if ( ( b0 == 0 ) && ( b1 == 0 ) )
    return 0;
  if ( ( a0 == 0 ) && ( a1 == 0 ) )
    return 0;

  if ( a0 == 0 )
    return delta* ( a1 - boost::math::asinh ( a1 ) ) / ( a1*a1 );
  if ( a1 == 0 )
    return delta* ( a0 - boost::math::asinh ( a0 ) ) / ( a0*a0 );

  if ( a0 == a1 )
    return delta* ( ( -sqrt ( 1 + pow ( a1,-2 ) ) - 1.0/ ( a1*sqrt ( 1 + pow ( a1,2 ) ) ) + a1/sqrt ( 1 + pow ( a1,2 ) ) +
                      ( 2*boost::math::asinh ( a1 ) )  /pow ( a1,2 ) ) /2. );

  return delta* ( ( - ( a1*boost::math::asinh ( a0 ) ) + a0*boost::math::asinh ( a1 ) ) / ( a0* ( a0 - a1 ) *a1 ) );
}




// Implementation of inline functions and external operators                 //

inline
std::ostream&
DGtal::operator<<( std::ostream & out, 
                   const DGtal::MeasureOfStraightLines & object )
{
  object.selfDisplay( out );
  return out;
}

//                                                                           //