---

Geometry.cpp

Go to the documentation of this file.

/*
Copyright (c) 2000-2003, Jelle Kok, University of Amsterdam
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.

3. Neither the name of the University of Amsterdam nor the names of its
contributors may be used to endorse or promote products derived from this
software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

#include "Geometry.h"
#include <stdio.h>    // needed for sprintf

int sign( double d1 )
{
  return (d1>0)?1:-1;
}

double max( double d1, double d2 )
{
  return (d1>d2)?d1:d2;
}

double min( double d1, double d2 )
{
  return (d1<d2)?d1:d2;
}


AngDeg Rad2Deg( AngRad x )
{
  return ( x * 180 / M_PI );
}

AngRad Deg2Rad( AngDeg x )
{
  return ( x * M_PI / 180 );
}

double cosDeg( AngDeg x )
{
  return ( cos( Deg2Rad( x ) ) );
}

double sinDeg( AngDeg x )
{
  return ( sin( Deg2Rad( x ) ) );
}

double tanDeg( AngDeg x )
{
  return ( tan( Deg2Rad( x ) ) );
}

AngDeg atanDeg( double x )
{
  return ( Rad2Deg( atan( x ) ) );
}

double atan2Deg( double x, double y )
{
  if( fabs( x ) < EPSILON && fabs( y ) < EPSILON )
    return ( 0.0 );

  return ( Rad2Deg( atan2( x, y ) ) );
}

AngDeg acosDeg( double x )
{
  if( x >= 1 )
    return ( 0.0 );
  else if( x <= -1 )
    return ( 180.0 );

  return ( Rad2Deg( acos( x ) ) );
}

AngDeg asinDeg( double x )
{
  if( x >= 1 )
    return ( 90.0 );
  else if ( x <= -1 )
    return ( -90.0 );

  return ( Rad2Deg( asin( x ) ) );
}

bool isAngInInterval( AngDeg ang, AngDeg angMin, AngDeg angMax )
{
  // convert all angles to interval 0..360
  if( ( ang    + 360 ) < 360 ) ang    += 360;
  if( ( angMin + 360 ) < 360 ) angMin += 360;
  if( ( angMax + 360 ) < 360 ) angMax += 360;

  if( angMin < angMax ) // 0 ---false-- angMin ---true-----angMax---false--360
    return angMin < ang && ang < angMax ;
  else                  // 0 ---true--- angMax ---false----angMin---true---360
    return !( angMax < ang && ang < angMin );
}

AngDeg getBisectorTwoAngles( AngDeg angMin, AngDeg angMax )
{
  // separate sine and cosine part to circumvent boundary problem
  return VecPosition::normalizeAngle(
            atan2Deg( (sinDeg( angMin) + sinDeg( angMax ) )/2.0,
                      (cosDeg( angMin) + cosDeg( angMax ) )/2.0 ) );
}

/*****************************************************************************/
/*******************   CLASS VECPOSITION   ***********************************/
/*****************************************************************************/

VecPosition::VecPosition( double x, double y, CoordSystemT cs )
{
  setVecPosition( x, y, cs );
}

VecPosition VecPosition::operator - ( )
{
  return ( VecPosition( -m_x, -m_y ) );
}

VecPosition VecPosition::operator + ( const double &d )
{
  return ( VecPosition( m_x + d, m_y + d ) );
}

VecPosition VecPosition::operator + ( const VecPosition &p )
{
  return ( VecPosition( m_x + p.m_x, m_y + p.m_y ) );
}

VecPosition VecPosition::operator - ( const double &d )
{
  return ( VecPosition( m_x - d, m_y - d ) );
}

VecPosition VecPosition::operator - ( const VecPosition &p )
{
  return ( VecPosition( m_x - p.m_x, m_y - p.m_y ) );
}

VecPosition VecPosition::operator * ( const double &d  )
{
  return ( VecPosition( m_x * d, m_y * d  ) );
}

VecPosition VecPosition::operator * ( const VecPosition &p )
{
  return ( VecPosition( m_x * p.m_x, m_y * p.m_y ) );
}

VecPosition VecPosition::operator / ( const double &d )
{
  return ( VecPosition( m_x / d, m_y / d  ) );
}

VecPosition VecPosition::operator / ( const VecPosition &p )
{
  return ( VecPosition( m_x / p.m_x, m_y / p.m_y ) );
}

void VecPosition::operator = ( const double &d )
{
  m_x = d;
  m_y = d;
}

void VecPosition::operator +=( const VecPosition &p )
{
  m_x += p.m_x;
  m_y += p.m_y;
}

void VecPosition::operator += ( const double &d )
{
  m_x += d;
  m_y += d;
}

void VecPosition::operator -=( const VecPosition &p )
{
  m_x -= p.m_x;
  m_y -= p.m_y;
}

void VecPosition::operator -=( const double &d )
{
  m_x -= d;
  m_y -= d;
}

void VecPosition::operator *=( const VecPosition &p )
{
  m_x *= p.m_x;
  m_y *= p.m_y;
}

void VecPosition::operator *=( const double &d )
{
  m_x *= d;
  m_y *= d;
}

void VecPosition::operator /=( const VecPosition &p )
{
  m_x /= p.m_x;
  m_y /= p.m_y;
}

void VecPosition::operator /=( const double &d )
{
  m_x /= d;
  m_y /= d;
}

bool VecPosition::operator !=( const VecPosition &p )
{
  return ( ( m_x != p.m_x ) || ( m_y != p.m_y ) );
}

bool VecPosition::operator !=( const double &d )
{
  return ( ( m_x != d ) || ( m_y != d ) );
}

bool VecPosition::operator ==( const VecPosition &p )
{
  return ( ( m_x == p.m_x ) && ( m_y == p.m_y ) );
}

bool VecPosition::operator ==( const double &d )
{
  return ( ( m_x == d ) && ( m_y == d ) );
}

ostream& operator <<( ostream &os, VecPosition v )
{
  return ( os << "( " << v.m_x << ", " << v.m_y << " )" );
}

void VecPosition::show( CoordSystemT cs )
{
  if( cs == CARTESIAN )
    cout << *this << endl;
  else
    cout << "( r: " << getMagnitude( ) << ", phi: " << getDirection( ) << "  )";
}

string VecPosition::str( CoordSystemT cs )
{
  char buf[ 1024 ];

  if( cs == CARTESIAN )
    sprintf( buf, "( %f, %f )", getX( ), getY( ) );
  else
    sprintf( buf, "( r: %f, phi: %f )", getMagnitude( ), getDirection( ) );

  string str( buf );
  return ( str );
}

bool VecPosition::setX( double dX )
{
  m_x = dX;
  return ( true );
}

double VecPosition::getX( ) const
{
  return ( m_x );
}

bool VecPosition::setY( double dY )
{
  m_y = dY;
  return ( true );
}

double VecPosition::getY( ) const
{
  return ( m_y );
}

void VecPosition::setVecPosition( double dX, double dY, CoordSystemT cs)
{
  if( cs == CARTESIAN )
  {
    m_x = dX;
    m_y = dY;
  }
  else
    *this = getVecPositionFromPolar( dX, dY );
}

double VecPosition::getDistanceTo( const VecPosition p )
{
  return ( ( *this - p ).getMagnitude( ) );
}

VecPosition VecPosition::setMagnitude( double d )
{
  if( getMagnitude( ) > EPSILON )
     ( *this ) *= ( d / getMagnitude( ) );

  return ( *this );
}

double VecPosition::getMagnitude( ) const
{
  return ( sqrt( m_x * m_x + m_y * m_y ) );
}

AngDeg VecPosition::getDirection( ) const
{
  return ( atan2Deg( m_y, m_x ) );
}

bool VecPosition::isInFrontOf( const VecPosition &p )
{
  return ( ( m_x > p.getX( ) ) ? true : false );
}

bool VecPosition::isInFrontOf( const double &d )
{
  return ( ( m_x > d ) ? true : false );
}

bool VecPosition::isBehindOf( const VecPosition &p )
{
  return ( ( m_x < p.getX( ) ) ? true : false );
}

bool VecPosition::isBehindOf( const double &d )
{
  return ( ( m_x < d ) ? true : false );
}

bool VecPosition::isLeftOf( const VecPosition &p )
{
  return ( ( m_y < p.getY( ) ) ? true : false );
}

bool VecPosition::isLeftOf( const double &d )
{
  return ( ( m_y < d ) ? true : false );
}

bool VecPosition::isRightOf( const VecPosition &p )
{
  return ( ( m_y > p.getY( ) ) ? true : false );
}

bool VecPosition::isRightOf( const double &d )
{
  return ( ( m_y > d ) ? true : false );
}

bool VecPosition::isBetweenX( const VecPosition &p1, const VecPosition &p2 )
{
  return ( ( isInFrontOf( p1 ) && isBehindOf( p2 ) ) ? true : false );
}

bool VecPosition::isBetweenX( const double &d1, const double &d2 )
{
  return ( ( isInFrontOf( d1 ) && isBehindOf( d2 ) ) ? true : false );
}

bool VecPosition::isBetweenY( const VecPosition &p1, const VecPosition &p2 )
{
  return ( ( isRightOf( p1 ) && isLeftOf( p2 ) ) ? true : false );
}

bool VecPosition::isBetweenY( const double &d1, const double &d2 )
{
  return ( ( isRightOf( d1 ) && isLeftOf( d2 ) ) ? true : false );
}

VecPosition VecPosition::normalize( )
{
  return ( setMagnitude( 1.0 ) );
}

VecPosition VecPosition::rotate( AngDeg angle )
{
  // determine the polar representation of the current VecPosition
  double dMag    = this->getMagnitude( );
  double dNewDir = this->getDirection( ) + angle;  // add rotation angle to phi
  setVecPosition( dMag, dNewDir, POLAR );          // convert back to Cartesian
  return ( *this );
}

VecPosition VecPosition::globalToRelative( VecPosition origin, AngDeg ang )
{
  // convert global coordinates into relative coordinates by aligning
  // relative frame and world frame. First perform translation to make
  // origins of both frames coincide. Then perform rotation to make
  // axes of both frames coincide (use negative angle since you rotate
  // relative frame to world frame).
  *this -= origin;
  return ( rotate( -ang ) );
}

VecPosition VecPosition::relativeToGlobal( VecPosition origin, AngDeg ang )
{
  // convert relative coordinates into global coordinates by aligning
  // world frame and relative frame. First perform rotation to make
  // axes of both frames coincide (use positive angle since you rotate
  // world frame to relative frame). Then perform translation to make
  // origins of both frames coincide.
  rotate( ang );
  *this += origin;
  return ( *this );
}

VecPosition VecPosition::getVecPositionOnLineFraction( VecPosition &p,
                                                       double      dFrac )
{
  // determine point on line that lies at fraction dFrac of whole line
  // example: this --- 0.25 ---------  p
  // formula: this + dFrac * ( p - this ) = this - dFrac * this + dFrac * p =
  //          ( 1 - dFrac ) * this + dFrac * p
  return ( ( *this ) * ( 1.0 - dFrac ) + ( p * dFrac ) );
}

VecPosition VecPosition::getVecPositionFromPolar( double dMag, AngDeg ang )
{
  // cos(phi) = x/r <=> x = r*cos(phi); sin(phi) = y/r <=> y = r*sin(phi)
  return ( VecPosition( dMag * cosDeg( ang ), dMag * sinDeg( ang ) ) );
}

AngDeg VecPosition::normalizeAngle( AngDeg angle )
{
  while( angle > 180.0  ) angle -= 360.0;
  while( angle < -180.0 ) angle += 360.0;

  return ( angle );
}


/*****************************************************************************/
/********************** CLASS GEOMETRY ***************************************/
/*****************************************************************************/

double Geometry::getLengthGeomSeries( double dFirst, double dRatio, double dSum )
{
  if( dRatio < 0 )
    cerr << "(Geometry:getLengthGeomSeries): negative ratio" << endl;

  // s = a + ar + ar^2 + .. + ar^n-1 and thus sr = ar + ar^2 + .. + ar^n
  // subtract: sr - s = - a + ar^n) =>  s(1-r)/a + 1 = r^n = temp
  // log r^n / n = n log r / log r = n = length
  double temp = (dSum * ( dRatio - 1 ) / dFirst) + 1;
  if( temp <= 0 )
    return -1.0;
  return log( temp ) / log( dRatio ) ;
}

double Geometry::getSumGeomSeries( double dFirst, double dRatio, double dLength)
{
  // s = a + ar + ar^2 + .. + ar^n-1 and thus sr = ar + ar^2 + .. + ar^n
  // subtract: s - sr = a - ar^n) =>  s = a(1-r^n)/(1-r)
  return dFirst * ( 1 - pow( dRatio, dLength ) ) / ( 1 - dRatio ) ;
}

double Geometry::getSumInfGeomSeries( double dFirst, double dRatio )
{
  if( dRatio > 1 )
    cerr << "(Geometry:CalcLengthGeomSeries): series does not converge" <<endl;

  // s = a(1-r^n)/(1-r) with n->inf and 0<r<1 => r^n = 0
  return dFirst / ( 1 - dRatio );
}

double Geometry::getFirstGeomSeries( double dSum, double dRatio, double dLength)
{
  // s = a + ar + ar^2 + .. + ar^n-1 and thus sr = ar + ar^2 + .. + ar^n
  // subtract: s - sr = a - ar^n) =>  s = a(1-r^n)/(1-r) => a = s*(1-r)/(1-r^n)
  return dSum *  ( 1 - dRatio )/( 1 - pow( dRatio, dLength ) ) ;
}

double Geometry::getFirstInfGeomSeries( double dSum, double dRatio )
{
  if( dRatio > 1 )
    cerr << "(Geometry:getFirstInfGeomSeries):series does not converge" <<endl;

  // s = a(1-r^n)/(1-r) with r->inf and 0<r<1 => r^n = 0 => a = s ( 1 - r)
  return dSum * ( 1 - dRatio );
}

int Geometry::abcFormula(double a, double b, double c, double *s1, double *s2)
{
  double dDiscr = b*b - 4*a*c;       // discriminant is b^2 - 4*a*c
  if (fabs(dDiscr) < EPSILON )       // if discriminant = 0
  {
    *s1 = -b / (2 * a);              //  only one solution
    return 1;
  }
  else if (dDiscr < 0)               // if discriminant < 0
    return 0;                        //  no solutions
  else                               // if discriminant > 0
  {
    dDiscr = sqrt(dDiscr);           //  two solutions
    *s1 = (-b + dDiscr ) / (2 * a);
    *s2 = (-b - dDiscr ) / (2 * a);
    return 2;
  }
}

/*****************************************************************************/
/********************* CLASS CIRCLE ******************************************/
/*****************************************************************************/

Circle::Circle( VecPosition pos, double dR )
{
  setCircle( pos, dR );
}

Circle::Circle( )
{
  setCircle( VecPosition(-1000.0,-1000.0), 0);
}

void Circle::show( ostream& os)
{
  os << "c:" << m_posCenter << ", r:" << m_dRadius;
}

bool Circle::setCircle( VecPosition pos, double dR )
{
  setCenter( pos );
  return setRadius( dR  );
}
bool Circle::setRadius( double dR )
{
  if( dR > 0 )
  {
    m_dRadius = dR;
    return true;
  }
  else
  {
    m_dRadius = 0.0;
    return false;
  }
}

double Circle::getRadius()
{
  return m_dRadius;
}

bool Circle::setCenter( VecPosition pos )
{
  m_posCenter = pos;
  return true;
}

VecPosition Circle::getCenter()
{
  return m_posCenter;
}

double Circle::getCircumference()
{
  return 2.0*M_PI*getRadius();
}

double Circle::getArea()
{
  return M_PI*getRadius()*getRadius();
}

bool Circle::isInside( VecPosition pos )
{
  return m_posCenter.getDistanceTo( pos ) < getRadius() ;
}

int Circle::getIntersectionPoints( Circle c, VecPosition *p1, VecPosition *p2)
{
    double x0, y0, r0;
    double x1, y1, r1;

    x0 = getCenter( ).getX();
    y0 = getCenter( ).getY();
    r0 = getRadius( );
    x1 = c.getCenter( ).getX();
    y1 = c.getCenter( ).getY();
    r1 = c.getRadius( );

    double      d, dx, dy, h, a, x, y, p2_x, p2_y;

    // first calculate distance between two centers circles P0 and P1.
    dx = x1 - x0;
    dy = y1 - y0;
    d = sqrt(dx*dx + dy*dy);

    // normalize differences
    dx /= d; dy /= d;

    // a is distance between p0 and point that is the intersection point P2
    // that intersects P0-P1 and the line that crosses the two intersection
    // points P3 and P4.
    // Define two triangles: P0,P2,P3 and P1,P2,P3.
    // with distances a, h, r0 and b, h, r1 with d = a + b
    // We know a^2 + h^2 = r0^2 and b^2 + h^2 = r1^2 which then gives
    // a^2 + r1^2 - b^2 = r0^2 with d = a + b ==> a^2 + r1^2 - (d-a)^2 = r0^2
    // ==> r0^2 + d^2 - r1^2 / 2*d
    a = (r0*r0 + d*d - r1*r1) / (2.0 * d);

    // h is then a^2 + h^2 = r0^2 ==> h = sqrt( r0^2 - a^2 )
    double      arg = r0*r0 - a*a;
    h = (arg > 0.0) ? sqrt(arg) : 0.0;

    // First calculate P2
    p2_x = x0 + a * dx;
    p2_y = y0 + a * dy;

    // and finally the two intersection points
    x =  p2_x - h * dy;
    y =  p2_y + h * dx;
    p1->setVecPosition( x, y );
    x =  p2_x + h * dy;
    y =  p2_y - h * dx;
    p2->setVecPosition( x, y );

    return (arg < 0.0) ? 0 : ((arg == 0.0 ) ? 1 :  2);
}

double Circle::getIntersectionArea( Circle c )
{
  VecPosition pos1, pos2, pos3;
  double d, h, dArea;
  AngDeg ang;

  d = getCenter().getDistanceTo( c.getCenter() ); // dist between two centers
  if( d > c.getRadius() + getRadius() )           // larger than sum radii
    return 0.0;                                   // circles do not intersect
  if( d <= fabs(c.getRadius() - getRadius() ) )   // one totally in the other
  {
    double dR = min( c.getRadius(), getRadius() );// return area smallest circ
    return M_PI*dR*dR;
  }

  int iNrSol = getIntersectionPoints( c, &pos1, &pos2 );
  if( iNrSol != 2 )
    return 0.0;

  // the intersection area of two circles can be divided into two segments:
  // left and right of the line between the two intersection points p1 and p2.
  // The outside area of each segment can be calculated by taking the part
  // of the circle pie excluding the triangle from the center to the
  // two intersection points.
  // The pie equals pi*r^2 * rad(2*ang) / 2*pi = 0.5*rad(2*ang)*r^2 with ang
  // the angle between the center c of the circle and one of the two
  // intersection points. Thus the angle between c and p1 and c and p3 where
  // p3 is the point that lies halfway between p1 and p2.
  // This can be calculated using ang = asin( d / r ) with d the distance
  // between p1 and p3 and r the radius of the circle.
  // The area of the triangle is 2*0.5*h*d.

  pos3 = pos1.getVecPositionOnLineFraction( pos2, 0.5 );
  d = pos1.getDistanceTo( pos3 );
  h = pos3.getDistanceTo( getCenter() );
  ang = asin( d / getRadius() );

  dArea = ang*getRadius()*getRadius();
  dArea = dArea - d*h;

  // and now for the other segment the same story
  h = pos3.getDistanceTo( c.getCenter() );
  ang = asin( d / c.getRadius() );
  dArea = dArea + ang*c.getRadius()*c.getRadius();
  dArea = dArea - d*h;

  return dArea;
}


/*****************************************************************************/
/***********************  CLASS LINE *****************************************/
/*****************************************************************************/

Line::Line( double dA, double dB, double dC )
{
  m_a = dA;
  m_b = dB;
  m_c = dC;
}

ostream& operator <<(ostream & os, Line l)
{
  double a = l.getACoefficient();
  double b = l.getBCoefficient();
  double c = l.getCCoefficient();

  // ay + bx + c = 0 -> y = -b/a x - c/a
  if( a == 0 )
    os << "x = " << -c/b;
  else
  {
    os << "y = ";
    if( b != 0 )
      os << -b/a << "x ";
    if( c > 0 )
       os << "- " <<  fabs(c/a);
    else if( c < 0 )
       os << "+ " <<  fabs(c/a);
  }
  return os;
}

void Line::show( ostream& os)
{
  os << *this;
}

VecPosition Line::getIntersection( Line line )
{
  VecPosition pos;
  double x, y;
  if( ( m_a / m_b ) ==  (line.getACoefficient() / line.getBCoefficient() ))
    return pos; // lines are parallel, no intersection
  if( m_a == 0 )            // bx + c = 0 and a2*y + b2*x + c2 = 0 ==> x = -c/b
  {                          // calculate x using the current line
    x = -m_c/m_b;                // and calculate the y using the second line
    y = line.getYGivenX(x);
  }
  else if( line.getACoefficient() == 0 )
  {                         // ay + bx + c = 0 and b2*x + c2 = 0 ==> x = -c2/b2
   x = -line.getCCoefficient()/line.getBCoefficient(); // calculate x using
   y = getYGivenX(x);       // 2nd line and calculate y using current line
  }
  // ay + bx + c = 0 and a2y + b2*x + c2 = 0
  // y = (-b2/a2)x - c2/a2
  // bx = -a*y - c =>  bx = -a*(-b2/a2)x -a*(-c2/a2) - c ==>
  // ==> a2*bx = a*b2*x + a*c2 - a2*c ==> x = (a*c2 - a2*c)/(a2*b - a*b2)
  // calculate x using the above formula and the y using the current line
  else
  {
    x = (m_a*line.getCCoefficient() - line.getACoefficient()*m_c)/
                    (line.getACoefficient()*m_b - m_a*line.getBCoefficient());
    y = getYGivenX(x);
  }

  return VecPosition( x, y );
}


int Line::getCircleIntersectionPoints( Circle circle,
              VecPosition *posSolution1, VecPosition *posSolution2 )
{
  int    iSol;
  double dSol1=0, dSol2=0;
  double h = circle.getCenter().getX();
  double k = circle.getCenter().getY();

  // line:   x = -c/b (if a = 0)
  // circle: (x-h)^2 + (y-k)^2 = r^2, with h = center.x and k = center.y
  // fill in:(-c/b-h)^2 + y^2 -2ky + k^2 - r^2 = 0
  //         y^2 -2ky + (-c/b-h)^2 + k^2 - r^2 = 0
  // and determine solutions for y using abc-formula
  if( fabs(m_a) < EPSILON )
  {
    iSol = Geometry::abcFormula( 1, -2*k, ((-m_c/m_b) - h)*((-m_c/m_b) - h)
              + k*k - circle.getRadius()*circle.getRadius(), &dSol1, &dSol2);
    posSolution1->setVecPosition( (-m_c/m_b), dSol1 );
    posSolution2->setVecPosition( (-m_c/m_b), dSol2 );
    return iSol;
  }

  // ay + bx + c = 0 => y = -b/a x - c/a, with da = -b/a and db = -c/a
  // circle: (x-h)^2 + (y-k)^2 = r^2, with h = center.x and k = center.y
  // fill in:x^2 -2hx + h^2 + (da*x-db)^2 -2k(da*x-db) + k^2 - r^2 = 0
  //         x^2 -2hx + h^2 + da^2*x^2 + 2da*db*x + db^2 -2k*da*x -2k*db
  //                                                         + k^2 - r^2 = 0
  //       (1+da^2)*x^2 + 2(da*db-h-k*da)*x + h2 + db^2  -2k*db + k^2 - r^2 = 0
  // and determine solutions for x using abc-formula
  // fill in x in original line equation to get y coordinate
  double da = -m_b/m_a;
  double db = -m_c/m_a;

  double dA = 1 + da*da;
  double dB = 2*( da*db - h - k*da );
  double dC = h*h + db*db-2*k*db + k*k - circle.getRadius()*circle.getRadius();

  iSol = Geometry::abcFormula( dA, dB, dC, &dSol1, &dSol2 );

  posSolution1->setVecPosition( dSol1, da*dSol1 + db );
  posSolution2->setVecPosition( dSol2, da*dSol2 + db );
  return iSol;

}

Line Line::getTangentLine( VecPosition pos )
{
  // ay + bx + c = 0 -> y = (-b/a)x + (-c/a)
  // tangent: y = (a/b)*x + C1 -> by - ax + C2 = 0 => C2 = ax - by
  // with pos.y = y, pos.x = x
  return Line( m_b, -m_a, m_a*pos.getX() - m_b*pos.getY() );
}

VecPosition Line::getPointOnLineClosestTo( VecPosition pos )
{
  Line l2 = getTangentLine( pos );  // get tangent line
  return getIntersection( l2 );     // and intersection between the two lines
}

double Line::getDistanceWithPoint( VecPosition pos )
{
  return pos.getDistanceTo( getPointOnLineClosestTo( pos ) );
}

bool Line::isInBetween( VecPosition pos, VecPosition point1,VecPosition point2)
{
  pos          = getPointOnLineClosestTo( pos ); // get closest point
  double dDist = point1.getDistanceTo( point2 ); // get distance between 2 pos

  // if the distance from both points to the projection is smaller than this
  // dist, the pos lies in between.
  return pos.getDistanceTo( point1 ) <= dDist &&
         pos.getDistanceTo( point2 ) <= dDist;
}

double Line::getYGivenX( double x )
{
 if( m_a == 0 )
 {
   cerr << "(Line::getYGivenX) Cannot calculate Y coordinate: " ;
   show( cerr );
   cerr << endl;
   return 0;
 }
  // ay + bx + c = 0 ==> ay = -(b*x + c)/a
  return -(m_b*x+m_c)/m_a;
}

double Line::getXGivenY( double y )
{
 if( m_b == 0 )
 {
   cerr << "(Line::getXGivenY) Cannot calculate X coordinate\n" ;
   return 0;
 }
  // ay + bx + c = 0 ==> bx = -(a*y + c)/a
  return -(m_a*y+m_c)/m_b;
}

Line Line::makeLineFromTwoPoints( VecPosition pos1, VecPosition pos2 )
{
  // 1*y + bx + c = 0 => y = -bx - c
  // with -b the direction coefficient (or slope)
  // and c = - y - bx
  double dA, dB, dC;
  double dTemp = pos2.getX() - pos1.getX(); // determine the slope
  if( fabs(dTemp) < EPSILON )
  {
    // ay + bx + c = 0 with vertical slope=> a = 0, b = 1
    dA = 0.0;
    dB = 1.0;
  }
  else
  {
    // y = (-b)x -c with -b the slope of the line
    dA = 1.0;
    dB = -(pos2.getY() - pos1.getY())/dTemp;
  }
  // ay + bx + c = 0 ==> c = -a*y - b*x
  dC =  - dA*pos2.getY()  - dB * pos2.getX();
  return Line( dA, dB, dC );
}

Line Line::makeLineFromPositionAndAngle( VecPosition vec, AngDeg angle )
{
  // calculate point somewhat further in direction 'angle' and make
  // line from these two points.
  return makeLineFromTwoPoints( vec, vec+VecPosition(1,angle,POLAR));
}

double Line::getACoefficient() const
{
  return m_a;
}

double Line::getBCoefficient() const
{
 return m_b;
}

double Line::getCCoefficient() const
{
 return m_c;
}

/*****************************************************************************/
/********************* CLASS RECTANGLE ***************************************/
/*****************************************************************************/

Rect::Rect( VecPosition pos, VecPosition pos2 )
{
  setRectanglePoints( pos, pos2 );
}

void Rect::setRectanglePoints( VecPosition pos1, VecPosition pos2 )
{
  m_posLeftTop.setX    ( max( pos1.getX(), pos2.getX() ) );
  m_posLeftTop.setY    ( min( pos1.getY(), pos2.getY() ) );
  m_posRightBottom.setX( min( pos1.getX(), pos2.getX() ) );
  m_posRightBottom.setY( max( pos1.getY(), pos2.getY() ) );
}

void Rect::show( ostream& os )
{
  os << "rect(" << m_posLeftTop << " " << m_posRightBottom << ")";
}

bool Rect::isInside( VecPosition pos )
{
  return pos.isBetweenX( m_posRightBottom.getX(), m_posLeftTop.getX() ) &&
         pos.isBetweenY( m_posLeftTop.getY(),     m_posRightBottom.getY() );

}

bool Rect::setPosLeftTop( VecPosition pos )
{
  m_posLeftTop = pos;
  return true;
}

VecPosition Rect::getPosLeftTop(  )
{
  return m_posLeftTop;
}

bool Rect::setPosRightBottom( VecPosition pos )
{
  m_posRightBottom = pos;
  return true;
}

VecPosition Rect::getPosRightBottom(  )
{
  return m_posRightBottom;
}

/*****************************************************************************/
/********************** TESTING PURPOSES *************************************/
/*****************************************************************************/

/*
#include<iostream.h>

int main( void )
{
  double dFirst = 1.0;
  double dRatio = 2.5;
  double dSum   = 63.4375;
  double dLength = 4.0;

  printf( "sum: %f\n", Geometry::getSumGeomSeries( dFirst, dRatio, dLength));
  printf( "length: %f\n", Geometry::getLengthGeomSeries( dFirst, dRatio, dSum));
}

int main( void )
{
  Line l1(1,-1,3 );
  Line l2(1,-0.2,10 );
 Line l3 = Line::makeLineFromTwoPoints( VecPosition(1,-1), VecPosition(2,-2) );
 l3.show();
 cout << endl;
 l1.show();
 l2.show();
  l1.getIntersection( l2 ).show();
}


int main( void )
{
  Line l( 1, -1, 0 );
  VecPosition s1, s2;
  int i = l.getCircleIntersectionPoints( Circle( VecPosition(1,1),1) &s1,&s2 );
  printf( "number of solutions: %d\n", i );
  if( i == 2 )
  {
    cout << s1 << " " << s2 ;
  }
  else if( i == 1 )
  {
    cout << s1;
  }
  cout << "line: " << l;
}

int main( void )
{
  Circle c11( VecPosition( 10, 0 ), 10);
  Circle c12( VecPosition( 40, 3 ), 40 );
  Circle c21( VecPosition(  0,0 ), 5);
  Circle c22( VecPosition(  3,0 ), 40 );

  VecPosition p1, p2;

  cout << c11.getIntersectionArea( c21 ) << endl;
  cout << c12.getIntersectionArea( c21 ) << endl;
  cout << c22.getIntersectionArea( c11 ) << endl;
  cout << c12.getIntersectionArea( c22 ) << endl;
  return 0;
}

int main( void )
{
  cout << getBisectorTwoAngles( -155.3, 179.0 ) << endl;
  cout << getBisectorTwoAngles( -179.3, 179.0 ) << endl;
}
*/

---