---

gunu_knot_removal.cpp

Go to the documentation of this file.

/* LIBGUL - Geometry Utility Library
 * Copyright (C) 1998-1999 Norbert Irmer
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library 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
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */

#include "stdafx.h"

#include <iostream>

#include "gul_types.h"
#include "gul_float.h"
#include "gul_vector.h"
#include "guma_sorting.h"
#include "gunu_basics.h"
#include "gunu_parametrize.h"
#include "gunu_global_approximate.h"
#include "gunu_derivatives.h"

#include "gunu_knot_removal.h"

namespace gunu {

using gul::Ptr;
using gul::hpoint2;
using gul::point2;
using gul::hpoint;
using gul::point;
using gul::rel_equal;
using gul::abs_equal;
using gul::euclidian_distance;
using gul::length;
using gul::rtr;
using guma::_BinarySearch;

/*------------------------------------------------------------------------
  Tries to remove the knot U[r] of multiplicity 's' (U[r]!=U[r+1]) from
  'U' 'num' times. It returns how often the knot actually was removed,
  the tolerance 'tol' is used to decide if a knot is removable.
  The new knot vector and control points overwrite the old ones.
------------------------------------------------------------------------*/

template< class T, class HP >
int RemoveCurveKnot( int num, int r, int s, T tol, 
                     int n, int p, Ptr<T>& U, Ptr<HP>& Pw )

{
  T u,alfi,alfj;
  int m,ord,fout,last,first,t,off,i,j,ii,jj,k;
  bool remflag;
  Ptr<HP> temp;

  temp.reserve_pool(2*p+1);
  u = U[r];
  
  m = n+p+1;
  ord = p+1;
  
  fout = (2*r - s - p)/2;
  last = r-s;
  first = r-p;
  
  for( t = 0; t < num; t++ )
  {
    off = first-1;
    temp[0] = Pw[off];
    temp[last+1-off] = Pw[last+1];

    i = first;
    j = last;
    ii = 1;
    jj = last - off;
    remflag = false;
    
    while( j-i > t )
    {
      alfi = (u-U[i])/(U[i+ord+t]-U[i]);
      alfj = (u-U[j-t])/(U[j+ord]-U[j-t]);
      temp[ii] = ((T)1/alfi)*(Pw[i]-((T)1-alfi)*temp[ii-1]);
      temp[jj] = ((T)1/((T)1-alfj))*(Pw[j]-alfj*temp[jj+1]);
      
      i = i+1;
      j = j-1;
      ii = ii+1;
      jj = jj-1;
    }

    if( j-i < t )
    {
      if( abs_equal(temp[ii-1], temp[jj+1], tol) )
        remflag = true;
    }
    else
    {
      alfi = (u-U[i])/(U[i+ord+t]-U[i]);
      if( abs_equal( Pw[i], alfi*temp[ii+t+1] + ((T)1-alfi)*temp[ii-1], tol) )
        remflag = true;
    }
    
    if( !remflag )
      break;
    
    i = first;
    j = last;
    
    while( j-i > t )
    {
      Pw[i] = temp[i-off];
      Pw[j] = temp[j-off];

      i = i+1;
      j = j-1;
    }
      
    first = first-1;
    last = last+1;    
  }
  
  if( t == 0 )
    return 0;
  
  for( k = r+1; k <= m; k++ )
    U[k-t] = U[k];
    
  i = j = fout;
  
  for( k = 1; k < t; k++ )
  {
    if( k%2 )
      i = i+1;
    else
      j = j-1;
  }
  
  for( k = i+1; k <= n; k++, j++ )
  {
    Pw[j] = Pw[k];
  }  
    
  return t;
}

// template instantiation
template int RemoveCurveKnot( 
                int num, int r, int s, float tol, 
                int n, int p, Ptr<float>& U, Ptr< hpoint<float> >& Pw );
template int RemoveCurveKnot( 
                int num, int r, int s, float tol, 
                int n, int p, Ptr<float>& U, Ptr< hpoint2<float> >& Pw );

template int RemoveCurveKnot( 
                int num, int r, int s, double tol, 
                int n, int p, Ptr<double>& U, Ptr< point<double> >& Pw );
template int RemoveCurveKnot( 
                int num, int r, int s, double tol, 
                int n, int p, Ptr<double>& U, Ptr< point2<double> >& Pw );


/*------------------------------------------------------------------------
  Removes the knot U[r] of multiplicity 's' (U[r]!=U[r+1]) from
  'U' 'num' times (without tolerance check). 
  The new knot vector and control points overwrite the old ones.
------------------------------------------------------------------------*/

template< class T, class HP >
void RemoveCurveKnot( int num, int r, int s, 
                     int n, int p, Ptr<T>& U, Ptr<HP>& Pw )

{
  T u,alfi,alfj;
  int m,ord,fout,last,first,t,off,i,j,ii,jj,k;
  Ptr<HP> temp;

  if( num == 0 )  // nothing to do
    return;

  temp.reserve_pool(2*p+1);
  u = U[r];
  
  m = n+p+1;
  ord = p+1;
  
  fout = (2*r - s - p)/2;
  last = r-s;
  first = r-p;
  
  for( t = 0; t < num; t++ )
  {
    off = first-1;
    temp[0] = Pw[off];
    temp[last+1-off] = Pw[last+1];

    i = first;
    j = last;
    ii = 1;
    jj = last - off;
      
    while( j-i > t )
    {
      alfi = (u-U[i])/(U[i+ord+t]-U[i]);
      alfj = (u-U[j-t])/(U[j+ord]-U[j-t]);
      temp[ii] = ((T)1/alfi)*(Pw[i]-((T)1-alfi)*temp[ii-1]);
      temp[jj] = ((T)1/((T)1-alfj))*(Pw[j]-alfj*temp[jj+1]);
      
      i = i+1;
      j = j-1;
      ii = ii+1;
      jj = jj-1;
    }

    i = first;
    j = last;
    
    while( j-i > t )
    {
      Pw[i] = temp[i-off];
      Pw[j] = temp[j-off];

      i = i+1;
      j = j-1;
    }
      
    first = first-1;
    last = last+1;    
  }
    
  for( k = r+1; k <= m; k++ )
    U[k-t] = U[k];
    
  i = j = fout;
  
  for( k = 1; k < t; k++ )
  {
    if( k%2 )
      i = i+1;
    else
      j = j-1;
  }
  
  for( k = i+1; k <= n; k++, j++ )
  {
    Pw[j] = Pw[k];
  }  
    
  return;
}

// template instantiation
template void RemoveCurveKnot( 
                int num, int r, int s,
                int n, int p, Ptr<float>& U, Ptr< hpoint<float> >& Pw );
template void RemoveCurveKnot( 
                int num, int r, int s,
                int n, int p, Ptr<float>& U, Ptr< hpoint2<float> >& Pw );

template void RemoveCurveKnot( 
                int num, int r, int s,
                int n, int p, Ptr<double>& U, Ptr< point<double> >& Pw );
template void RemoveCurveKnot( 
                int num, int r, int s,
                int n, int p, Ptr<double>& U, Ptr< point2<double> >& Pw );


/*-----------------------------------------------------------------------------
  Calculates a knot removal error factor
-----------------------------------------------------------------------------*/
template< class T, class HP >
T RemovalError( int r, int s, int n, int p, const Ptr<T>& U, const Ptr<HP>& Pw )

{
  T u,alfi,alfj,Br;
  int ord,last,first,off,i,j,ii,jj;
  Ptr<HP> temp;

  temp.reserve_pool(2*p+1);
  u = U[r];
  
  ord = p+1;
  last = r-s;
  first = r-p;

  off = first-1;

  temp[0] = Pw[off];
  temp[last+1-off] = Pw[last+1];

  i = first;
  j = last;
  ii = 1;
  jj = last - off;
  
  while( j-i > 0 )
  {
    alfi = (u-U[i])/(U[i+ord]-U[i]);
    alfj = (u-U[j])/(U[j+ord]-U[j]);

    temp[ii] = ((T)1/alfi)*(Pw[i]-((T)1-alfi)*temp[ii-1]);
    temp[jj] = ((T)1/((T)1-alfj))*(Pw[j]-alfj*temp[jj+1]);
      
    i = i+1;
    j = j-1;
    ii = ii+1;
    jj = jj-1;
  }

  if( j-i < 0 )
  {
    Br = euclidian_distance( temp[ii-1], temp[jj+1] );
  }
  else
  {
    alfi = (u-U[i])/(U[i+ord]-U[i]);
    Br = euclidian_distance( Pw[i], alfi*temp[ii+1] + ((T)1-alfi)*temp[ii-1] );
  }
    
  return Br; 
}

// template instantiation
template float RemovalError( 
      int r, int s, int n, int p, const Ptr<float>& U, 
          const Ptr< point<float> >& Pw );

/*----------------------------------------------------------------------
 Calculates the number of distinct knots of a knot vector, and their 
 multiplicities
----------------------------------------------------------------------*/

template< class T >
int GetMultiplicities( int n, int p, const Ptr<T>& U, Ptr<int>& S, Ptr<int>& R )
{
  int i,k,s;

  k = 0;
  i = p+1;

  while( i <= n )  
  {
    s = 1;

    while( U[i+s] == U[i] )
      s++;

    S[k] = s;
    R[k] = i + s - 1;

    i += s;
    k++;
  }

  return k;
}

// template instantiation
template int GetMultiplicities( 
             int n, int p, const Ptr<float>& U, Ptr<int>& S, Ptr<int>& R );



/*-------------------------------------------------------------------------
  calc range of parameter indices which fall into the domain of each basis
  function
--------------------------------------------------------------------------*/

template< class T >
GULAPI void CalcParameterRanges(
       int           n,
       int           p,
       const Ptr<T>& U,
       int           nP,
       const Ptr<T>& P,
       Ptr<int>&     R,
       Ptr<int>&     S )

{
  int   i,ilo,ihi;
  bool  finished_early;

  ilo = 0;
  ihi = 0;
  finished_early = false;

  for( i = 0; i <= n; i++ )
  {
    while( P[ilo] < U[i] ) // search first param. >= ulo
    {
      ilo++;

      // if all remaining parameter ranges are empty
      if( ilo >= nP )
      {
        finished_early = true;
        break;
      }
    }

    if( finished_early )
      break;

    R[i] = ilo;

    while( U[i+p+1] >= P[ihi] ) // search last param. <= uhi
    {
      ihi++;
      if( ihi == nP )
        break;
    }
    ihi--;                  

    // if parameter range is empty, start next search with 0
    // again
    if( ihi == -1 )
    {
      ihi = 0;
      S[i] = 0;
    }
    else
    {
      S[i] = ihi-ilo+1;

      // if not the last interval, skip last point, if it is on
      // the interval end

      if( (i+p+1 <= n) && (P[ihi] == U[i+p+1]) ) 
      {  
        S[i]--;
        ihi--;
        if( ihi < 0 )
          ihi = 0;
      }
    }
  }

  if( finished_early )
  {
    for( ; i <= n; i++ )
    {
      R[i] = nP;
      S[i] = 0;
    }
  }
}
// template instantiation
template GULAPI void CalcParameterRanges(
       int               n,
       int               p,
       const Ptr<float>& U,
       int               nP,
       const Ptr<float>& P,
       Ptr<int>&         R,
       Ptr<int>&         S );


/*-------------------------------------------------------------------------
  calculates the maximum number of indices which fall into one kotspan
  (this function is a hack, like the function above :)
--------------------------------------------------------------------------*/

template< class T >
GULAPI int CalcMaxPointsPerSpan(
       int           n,
       int           p,
       const Ptr<T>& U,
       int           nP,
       const Ptr<T>& P )
{
  int i1,i2,i;
  T u1,u2;
  int maxpps;

  i = 1;
  i1 = 0;
  i2 = 0;
  maxpps = 0;

  for(;;)
  {
    // search next knotspan [u1,u2]
    while( (i <= n+p+1) && (U[i] == U[i-1]) ) i++;
    if( i == n+p+2 ) break; // not found
    u1 = U[i-1];
    u2 = U[i];
    i++;

    // search first parameter >= u1
    while( (i1 < nP) && (P[i1] < u1) ) i1++;
    if( (i1 == nP) || (P[i1] < u1) ) // not found
      break;
  
    // search last parameter <= u2
    while( (i2 < nP) && (P[i2] <= u2) ) i2++;
    i2--;
    if( (i2 < 0) || (P[i2] > u2) ) // not found
      continue;
 
    if( i2-i1+1 > maxpps )
      maxpps = i2-i1+1;

    i1 = i2;
  }

  return maxpps;
}

/*----------------------------------------------------------------------
  Removes (roughly) as many knots as possible from a curve. it returns the
  index of the last inner knot of the new knot vector (new 'n'). U,Pw,
  and QE are input _and_ output arrays.
-----------------------------------------------------------------------*/  
template< class T, class EP >
GULAPI int RemoveCurveKnots(
        int             n,
        int             p,
        Ptr<T>&         U,
        Ptr<EP>&        Pw,  
        int             nQ,
        const Ptr<T>&   QU,
        T               tol,
        Ptr<T>&         QE )
{
  Ptr<int> R,S,QR,QS;
  Ptr<T>   B,newerr;
  int      i,mini,nK,r,s,k,ilo,ihi,dummi;
  T        minB,N,alfa;
  bool     removable;

  S.reserve_pool(n-p);
  R.reserve_pool(n-p);
  B.reserve_pool(n-p);
  QR.reserve_pool(n+1);
  QS.reserve_pool(n+1);
  newerr.reserve_pool(nQ);

  nK = GetMultiplicities(n,p,U,S,R);
  if( nK == 0 ) return n;  // nothing to remove

  // calc removal error factors for each distinct interior knot

  minB = B[0] = RemovalError(R[0],S[0],n,p,U,Pw);
  mini = 0;

  for( i = 1; i < nK; i++ )
  {
    B[i] = RemovalError(R[i],S[i],n,p,U,Pw);
    if( B[i] <= minB )
    {
      minB = B[i];
      mini = i;
    }
  }

  CalcParameterRanges(n,p,U,nQ,QU,QR,QS);

  for(;;)
  {
    if( minB == rtr<T>::maximum() ) break;

    r = R[mini];
    s = S[mini];

    removable = true;

    if( ((p+s)%2) == 0 )
    {
      k = (p+s)/2; 

      for( i = QR[r-k]; i < QR[r-k]+QS[r-k]; i++ ) 
      {
        N = CalcBasisFunction(p,r-k,QU[i],n,U);
        newerr[i] = N*minB + QE[i]; 
        if( newerr[i] > tol )
        {
          removable = false;
          B[mini] = rtr<T>::maximum();
          break;
        }
      }
      if( removable )
      {
        for( i = QR[r-k]; i < QR[r-k]+QS[r-k]; i++ ) 
          QE[i] = newerr[i];
      }
    }
    else
    {
      k = (p+s+1)/2;
      alfa = (U[r]-U[r-k+1])/(U[r-k+p+2]-U[r-k+1]);
      minB *= (T)1 - alfa;

      for( i = QR[r-k+1]; i < QR[r-k+1]+QS[r-k+1]; i++ ) 
      {
        N = CalcBasisFunction(p,r-k+1,QU[i],n,U);
        newerr[i] = N*minB + QE[i]; 
        if( newerr[i] > tol )
        {
          removable = false;
          B[mini] = rtr<T>::maximum();
          break;
        }
      }
      if( removable )
      {
        for( i = QR[r-k+1]; i < QR[r-k+1]+QS[r-k+1]; i++ ) 
          QE[i] = newerr[i];
      }
    }
  
    if( removable )
    {
      RemoveCurveKnot(1,r,s,n,p,U,Pw);

      // U now contains the new knot vector
      n--;

      // remove entries for removed knot/basis function

      if( s == 1 )
      {
        nK--;
        if( nK == 0 ) break;  // no more internal knots 
      
        for( i = mini; i < nK; i++ )
        {
          B[i] = B[i+1];
          R[i] = R[i+1]-1;
          S[i] = S[i+1];
        }
      }
      else
      {
        S[mini]--;

        for( i = mini; i < nK; i++ )
          R[i]--;
      } 

      // calculate changed 'Br' values:
      
      // for these inner knots new span index = old span index:

      i = mini-1;
      while( (i >= 0) && (R[i] >= r-p) )
      {
        B[i] = RemovalError(R[i],S[i],n,p,U,Pw);
        i--;
      }

      // for these inner knots new span index = old span index - 1:

      i = mini;
      while( (i < nK) && (R[i] <= r+p-s) )
      {
        B[i] = RemovalError(R[i],S[i],n,p,U,Pw);
        i++;
      }

      // update parameter index ranges
      for( i = r-p-1; i <= r-s; i++ )
      {
        if( QS[i] == 0 ) // if parameter range was empty for this basis func.
        {
          _BinarySearch(U[i],nQ,QU,&dummi,&ilo);
          QR[i] = ihi = ilo;
        }
        else
          ilo = ihi = QR[i];

        if( ilo < 0 ) continue; // still no points in this domain
  
        while( QU[ihi] <= U[i+p+1] )
        {
          ihi++;
          if( ihi >= nQ ) break;
        }
        ihi--;

        // if not the last interval, skip last point if it is on
        // the interval end
        if( (i+p+1 <= n)  && (QU[ihi] == U[i+p+1]) )
          ihi--;        

        QS[i] = ihi-ilo+1;
      }

      // remove parameter range entries for removed basis function

      for( i = r; i <= n; i++ )
      {
        QR[i] = QR[i+1];
        QS[i] = QS[i+1];
      }
    }      

    // search knot with smallest 'Br'
    mini = 0;
    minB = B[0];
    for( i = 1; i < nK; i++ )
    {
      if( B[i] <= minB )
      {
        mini = i;
        minB = B[i];
      }
    }

  }

  return n;
}

// template instantiation
template GULAPI int RemoveCurveKnots(
        int                         n,
        int                         p,
        Ptr<float>&                 U,
        Ptr< point<float> >&        Pw,  
        int                         nQ,
        const Ptr<float>&           QU,
        float                       tol,
        Ptr<float>&                 QE );



/*-------------------------------------------------------------------
  newton iteration
--------------------------------------------------------------------*/

template< class T, class EP >
bool PTCCheck( const EP& P, const Ptr<EP>& D, T eps, T *err )
{
  EP C,V,A;
  T num,den;

  C = D[0];
  V = D[1];  
  A = D[2];  
  
  // check point coincidence
  if( rel_equal(C,P,eps) )
    return true;

  // check zero cosine
  num = V*(C-P);
  den = rtr<T>::sqrt(V*V) * rtr<T>::sqrt((C-P)*(C-P));
  if( rel_equal(num+den,den,eps) )
    return true;

  *err = rtr<T>::fabs(num/den);

  return false;
}

template< class T, class HP, class EP >
bool ProjectToCurve( const EP& P, T u, T eps, int n, int p, 
                     const Ptr<T>& U, const Ptr<HP>& Pw,
                     T *ret_u, Ptr<EP>& D )
{
  bool found;
  EP C,V,A;
  T u1,err,err1,a,b;
  int t,s,sav,i,j;

  CurveDerivatives(u,n,p,U,Pw,2,D);
  found = PTCCheck( P, D, eps, &err );
  if( found ) 
  {
    *ret_u = u;
    return true;  
  }

  // initialize parameter range interval [a,b]
  // (C2 continuity required for changing knot spans)

  t = sav = FindSpan(u,n,p,U);
  for(;;)  // find left boundary
  {
    s = 1;
    while( (s <= p) && (U[t-s] == U[t]) ) s++;
    if( p-s < 2 )
    {
      a = U[t];
      break;
    }
    t = t-s;
  }

  t = sav + 1;
  for(;;)  // find right boundary
  {
    s = 1;
    while( (s <= p) && (U[t+s] == U[t]) ) s++;
    if( p-s < 2 )
    {
      b = U[t];
      break;
    }
    t = t+s;
  }

  // main loop
  for( j = 0; j < 100; j++ )  // max iterations = 100 ?
  {
    // calc next value for u
    C = D[0];
    V = D[1];  
    A = D[2];  
    u1 = u - (V*(C-P))/(A*(C-P) + V*V);

    // clip with search interval
    if( u1 < a ) u1 = a;
    if( u1 >= b ) u1 = b - b*rtr<T>::epsilon();

    // check that parameter changes enough
    if( rel_equal(u,u1,eps) )
    {
      *ret_u = u;
      return true;
    }

    // check step width (error should decrease when doing a newton step)
    for( i = 0; i < 10; i++ )
    {
      CurveDerivatives(u1,n,p,U,Pw,2,D);
      found = PTCCheck( P, D, eps, &err1 );
      if( found ) 
      {
        *ret_u = u1;
        return true;  
      }

      if( err1 < err )
        break;

      u1 = u + (T)0.5*(u1-u);
    }
    // step width adjustment not successful
    if( err1 > err )
      return false;

    u = u1;
    err = err1;
  }

  return false;
}

/*-----------------------------------------------------------------------
  increase multiplicity of all knots in a knot vector by one,
  the new knot vector must be reserved by the caller, it returns
  the index of the last knot in the new knot vector
------------------------------------------------------------------------*/
template< class T >
inline int IncMultiplicities( int n, int p, const Ptr<T>& U, Ptr<T>& Uh )
{
  int i,j;

  Uh[0] = U[0];
  j = 1; 

  for( i = 1; i <= n+p+1; i++ )
  {
    if( U[i] != U[i-1] )
    {
      Uh[j] = U[i-1];
      j++;
    }
    Uh[j] = U[i];
    j++;
  }

  Uh[j] = U[n+p+1];

  return j;
}

/*-----------------------------------------------------------------------
  global curve approximation with error bound
------------------------------------------------------------------------*/
template< class T, class EP >
GULAPI bool GlobalCurveApproximationE( 
        int            nQ, 
        const Ptr<EP>& Q,
        T              tol,   // a bound for the absolute error of the curve
        T              eps,   // a bound for the relative error (needed internally
                              // for measuring point coincidence and zero cosine)
        int            degree,
        int           *ret_n,
        Ptr<T>        *ret_U,
        Ptr<EP>       *ret_Pw,
        Ptr<T>        *ret_QE,  // error at each data point
        Ptr<T>        *ret_QU ) // parameter value of each data point 
{
  Ptr<T> d,QU,QE,U;
  Ptr<EP> Pw;
  T sum,u;
  int i,n,nn,p;
  int j,nh,ph,maxpps,deg;
  Ptr<T> Uh;
  Ptr<EP> Ph,Pwh;
  bool result;
  Ptr<EP> D;

  d.reserve_pool(nQ);
  QU.reserve_pool(nQ);
  QE.reserve_pool(nQ);
  Pw.reserve_pool(nQ);
  U.reserve_pool(nQ+2);
  D.reserve_pool(3);

  sum = ChordLength( nQ, Q, d );
  if( sum == (T)0 ) return false;
  ChordLengthParameters( nQ, sum, d, QU );

  // build a linear curve interpolating the data points
  U[0] = (T)0;
  U[1] = (T)0;
  for( i = 1; i < nQ-1; i++ )
    U[i+1] = QU[i];
  U[nQ] = (T)1;
  U[nQ+1] = (T)1;

  for( i = 0; i < nQ; i++ )
    Pw[i] = Q[i];

  // initialize errors to zero
  for( i = 0; i < nQ; i++ )
    QE[i] = (T)0;

  n = nQ-1;
  p = 1;

  for( deg = 1; deg <= degree; deg++ )
  { 
    nn = RemoveCurveKnots( n, p, U, Pw, nQ, QU, tol, QE );

    // build a knot vector with degree p + 1 (increases multiplicities by 1)
    if( p < degree )
    {
      n = nn;

      Uh.reserve_pool(2*(n+p+2));
      j = IncMultiplicities(n,p,U,Uh);

      ph = p+1;
      nh = j-ph-1;
      Pwh.reserve_pool(nh+1);  
    }
    // when the final degree has been reached, do a final fit with a curve of
    // the same degree, if the knot removal was successful
    else
    {
      if( nn == n )
        break; 

      n = nn;

      Uh = U;
      ph = p;
      nh = n;
      Pwh.reserve_pool(nh+1);  
    }
    maxpps = CalcMaxPointsPerSpan(nh,ph,Uh,nQ,QU);
    result = DoGlobalCurveApproximation(nQ,Q,QU,nh,ph,Uh,maxpps,Pwh);
    if( !result ) return false;

    for( i = 0; i < nQ; i++ )
    {
      result = ProjectToCurve<T,EP,EP>( Q[i], QU[i], eps, 
                                        nh, ph, Uh, Pwh, &u, D );
      if( !result ) return false;

      if( (i > 0) && (u <= QU[i-1]) ) // for the time being
        return false;
  
      QU[i] = u;   
      QE[i] = length(D[0]-Q[i]);
    }

    Pw = Pwh;
    U = Uh;
    n = nh;
    p = ph;
  }
 
  *ret_n = n;
  (*ret_U) = U;
  (*ret_Pw) = Pw;
  (*ret_QE) = QE;
  (*ret_QU) = QU;

  return true;
}

// template instantiation
template GULAPI bool GlobalCurveApproximationE< float, point<float> >( 
        int                 nQ, 
        const               Ptr<point<float> >& Q,
        float               tol,   // a bound for the absolute error of the curve
        float               eps,   // a bound for the relative error (needed internally
                                   // for measuring point coincidence and zero cosine)
        int                 degree,
        int                *ret_n,
        Ptr<float>         *ret_U,
        Ptr<point<float> > *ret_Pw, 
        Ptr<float>         *ret_QE,   // error at each data point
        Ptr<float>         *ret_QU ); // parameter value of each data point

}

---