---

curve.cpp

//$Header: /home/ben/Mapper/c++/RCS/curve.cpp,v 6.6 2002/07/14 17:08:52 ben Exp $
// Copyright Benedict Adamson 2002.
// This file is part of Mapper.

// Mapper is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.

// Mapper is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.

// You should have received a copy of the GNU General Public License
// along with Mapper; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#include <cmath>

#include "curve.h"
#include "func.h"
#include "minimise.h"



static const double pi = M_PI;



curve::curve(void)
{
   //Do nothing
}



curve::~curve(void)
{
   //Do nothing
}



vector<double> curve::shape(
   const double t
   ) const
{
   {//preconditionss:
      assert( 2 <= this->param.size() );
   };
   const unsigned n = this->param.size();
   vector< double > s( n );
   s[0] = 1 - t;
   s[1] = t;
   if( 2 < n ){
      const double pt = pi * t;
      for( unsigned i = 2; i < n; i++ ){
         const double w = sin( pt * double(i-1) );
         s[0] -= w;
         s[i] = w;
      };
   };
   return s;
}



void curve::shape(
   const double t,
   vector< double > &s, //output
   vector< double > &dsdt //output
   ) const
{
   {//preconditionss:
      assert( 2 <= this->param.size() );
   };
   const unsigned n = this->param.size();
   s = vector<double >( n );
   dsdt = vector<double >( n );
   s[0] = 1 - t;
   dsdt[0] = -1;
   s[1] = t;
   dsdt[1] = 1;
   if( 2 < n ){
      const double pt = pi * t;
      for( unsigned i = 2; i < n; i++ ){
         const double f = double(i-1);
         const double a = pt * f;
         const double w = sin( a );
         const double dwdt = pi * f * cos( a );
         s[0] -= w;
         dsdt[0] -= dwdt;
         s[i] = w;
         dsdt[i] = dwdt;
      };
   };
}



vector2 curve::position(
   const double t
   ) const
{
   {//preconditions:
      assert( 2 <= this->param.size() );
      //normally( 0 <= t && t <= 1 );
   };
   vector2 p( 0, 0 );
   const unsigned n = this->param.size();
   vector< double > s( this->shape( t ) );
   assert( s.size() == n );
   for( unsigned i = 0; i < n; i++ ){
     const vector2 *ph = this->param[i];
     assert( ph );
     p += s[i] * (*ph);
   };
   return p;
}



namespace {
   
   
   
   class curve_distance2
      :
   public func::single_dimension
   {
         //Calculates the square of the distance between a given point
         //and points on a given curve.
      public: //constructors and destructor
         
         inline curve_distance2(void);
         //postoncditions:
         //assert( !this->c );
         //assert( !this->p );
         
         inline virtual ~curve_distance2(void);
         
      public: //attributes
         
         virtual double operator() (
            const double t
            ) const;
         //preconditions:
         //assert( this->c );
         //assert( this->p );
         
      public: //relationships

         const curve *c;

         const vector2 *p;
         
   };
   
   
   
   inline curve_distance2::curve_distance2(void)
      :
      c( 0 ),
      p( 0 )
      {
         //Do nothing
         {//postoncditions:
            assert( !this->c );
            assert( !this->p );
         };
      }
   
   
   
   inline curve_distance2::~curve_distance2(void)
      {
         //Do nothing
      }
   
   
   double curve_distance2::operator() (
      const double t
      ) const
      {
         {//preconditions
            assert( c );
            assert( p );
         };
         const vector2 dp( this->c->position( t ) - *(this->p) );
         return dp*dp;
      }
   
};



void nearest_point(
   const curve &cu,
   const vector2 &p,
   const double tol,
   const unsigned max_iter,
   double &t2, //output
   double &d22, //output
   unsigned &iter //output
   )
{
   //Calculates the position on the curve (t on cu
   //which is closest to the given point (p).
   {//preconditions:
      assert( 2 <= cu.param.size() );
      assert( 0.0 < tol );
      assert( 0 < max_iter );
   };
   curve_distance2 cd2; //function object used for minimisation calculation
   cd2.c = &cu;
   cd2.p = &p;
   double t1, t3;
   double d21, d23;
   if( *cu.param[0] != *cu.param[1] ){//curve has length
      {//determine good starting point for bracket
         const unsigned np = cu.param.size();
         if( np <= 2 ){//straight line
            t1 = 0.25;
            t2 = 0.75;
         }else{//curved
            //Must be carefull because there could be multiple minima.
            //Try to find a point (t2) close to the global minimum.
            const unsigned n = (np-2)*2;
            unsigned i2; //the best so far
            for( unsigned i = 0; i <= n; i++ ){
               const double ti = double(i) / double(n);
               const double d2i = cd2( ti );
               if( 0 == i || d2i <= d22 ){//better than the best so far
                  i2 = i;
                  t2 = ti;
                  d22 = d2i;
               };
            };
            assert( i2 <= n );
            assert( 0 <= d22 );
            t1 = double( i2 - 1 ) / double(n);
         };
      };
      unsigned iter_1, iter_2;
      minimise::bracket(
         cd2, max_iter, iter_1,
         t1, t2, t3,
         d21, d22, d23
         );
      if( iter_1 < max_iter ){ //allowed to do better
         minimise::Brent(
            cd2,
            tol, max_iter - iter_1,
            t1, t2, t3,
            d21, d22, d23,
            iter_2
            );
      }else{ //failed to find bracket
         iter_2 = 0;
      };
      iter = iter_1 + iter_2;
      const double t2o = t2;
      t2 = min<double>( t2, 1 );
      t2 = max<double>( t2, 0 );
      if( t2o != t2 )
         d22 = cd2( t2 );
      ;
   }else{ //zero-length curve
      t2 = 0.5;
      iter = 0;
      d22 = cd2( t2 );
   };
   {//postconditions:
      assert( 0 <= t2 && t2 <= 1 );
      assert( iter <= max_iter );
      assert( 0 <= d22 );
   };
}



void curve::position(
   const double t,
   vector2 &p,
   vector2 &dpdt
   ) const
{
   {//preconditions:
      assert( 2 <= this->param.size() );
      //normally( 0 <= t && t <= 1 );
   };
   {//initialise
      p = vector2( 0, 0 );
      dpdt = vector2( 0, 0 );
   };
   vector< double > s, dsdt;
   this->shape( t, s, dsdt );
   for( unsigned i = 0; i < this->param.size(); i++ ){
      const vector2 *ph = this->param[i];
      assert( ph );
      const vector2 v( *ph );
      p += s[i] * v;
      dpdt += dsdt[i] * v;
   };
}

---