spline.C

Go to the documentation of this file.

#ifdef HAVE_CONFIG_H
#include <tumble-conf.h>
#endif /* HAVE_CONFIG_H */

//#define RAW_LS

#include <iostream>
#include <cstddef>

#include "globals.h"
#include "datastore.h"
#include "spline.h"


using namespace std;



/* These are ratios for the constraints on spline movement
   big is the ratio between the differentiability at the end vs. everything else
   constraint_ratio is the ratio from midside constraints to vertex constraints
 */
static const double big=0.0;
static const double constraint_ratio=10.0;



QBSpline::QBSpline(const vector<ControlPoint> &deboor, DataStore *_data_store, bool _closed){
    vector<double> newk;
    initialize(deboor, newk, _data_store, 1.0, _closed);
}

QBSpline::QBSpline(const vector<ControlPoint> &deboor, const vector<double> &newk, DataStore *_data_store, bool _closed){
    initialize(deboor, newk, _data_store, 1.0, _closed);
}

QBSpline::QBSpline(const vector<ControlPoint> &deboor, const vector<double> &newk, DataStore *_data_store, double _restlength, bool _closed)
{
    initialize(deboor, newk, _data_store, _restlength, _closed);
}

void QBSpline::initialize(const vector<ControlPoint> &deboor, const vector<double> &newk, DataStore *_data_store, double _restlength, bool _closed)
{
    unsigned i;
    data_store = _data_store;
    restlength = _restlength;
    closed=_closed;

    if( !data_store )
        FATAL_ERROR("QBSpline::intialize -- data_store is NULL");

    /* closed splines must have the same ControlPoint as the first and last deboor point */
    if( closed && (deboor[0] != deboor[deboor.size()-1]) )
        FATAL_ERROR("QBSpline::intialize -- Closed spline has two distinct endpoints");

    /* closed splines must have at least 5 points */
    if( closed  && deboor.size() < 5 )
        FATAL_ERROR("QBSpline::intialize -- Closed spline has "<<deboor.size()<<" points, but must have at least 5.");

    /* open splines must have at least 3 points */
    if( deboor.size() < 3 )
        FATAL_ERROR("QBSpline::intialize -- Spline has "<<deboor.size()<<" points but must have at least 3.");


    /* Initialize the vector of internal BezierVertexs to NULL
     * There should be |deboor|-3 internal verticies
     */
    vertexs.insert(vertexs.begin(), deboor.size()-3, NULL);

    /* Initialize the knot vector
     * If newk is empty, make our own knots.
     * Otherwise we were given a knot vector, so use it.
     * There should be |deboor|-1 knots.
     */
    if(newk.size() == 0){
        k.push_back(0.0);
        for(i=1; i < deboor.size()-1; i++){
            k.push_back( k.back() + (*deboor[i+1] - *deboor[i-1]).mag());
        }
    } else {
        assert(newk.size() == deboor.size() - 1);
        for(i=0; i < newk.size()-1; i++){
            if( newk[i] >= newk[i+1] ) FATAL_ERROR("QBSpline::intialize -- knot sequence is non-increasing" );
        }
        k=newk;
    }

    /* Populate the b vector there are 2*|deboor|-3 ControlPoints
     * in b.  The new points must be added to the point list of the
     * BoundryComplex.
     */
    Point2D new_point;
    ControlPoint temp;
    b.push_back(deboor[0]);
    b.push_back(deboor[1]);
    for(i=1; i < deboor.size()-2; i++){
        temp=deboor[i+1];
        new_point = *(b.back()) * (dt(i)/(dt(i-1)+dt(i)))
                  + (*temp)   * (dt(i-1)/(dt(i-1)+dt(i)));
        b.push_back(data_store->add_cp(new_point));
        b.push_back(temp);
    }
    b.push_back(deboor.back());
}

QBSpline::QBSpline(const QBSpline &other)
{
    copy(other.b.begin(), other.b.end(), back_inserter(b));
    copy(other.k.begin(), other.k.end(), back_inserter(k));
    copy(other.vertexs.begin(), other.vertexs.end(), back_inserter(vertexs));
    closed = other.closed;
    data_store = other.data_store;
}

QBSpline& QBSpline::operator=(const QBSpline &other)
{
    copy(other.b.begin(), other.b.end(), back_inserter(b));
    copy(other.k.begin(), other.k.end(), back_inserter(k));
    copy(other.vertexs.begin(), other.vertexs.end(), back_inserter(vertexs));
    closed = other.closed;
    data_store = other.data_store;
    return *this;
}

ControlPoint QBSpline::get_segment_center(double u0, double u1) const
{
    int u0_idx, u1_idx;
    double temp;
    if(u0 > u1){
        temp = u0;
        u0 = u1;
        u1 = temp;
    }

    u0_idx = get_closest_knot(u0);
    u1_idx = get_closest_knot(u1);
    //cout<<"u0: "<<u0<<" u0-idx: "<<u0_idx<<"  u1: "<<u1<<" u1_idx: "<<u1_idx<<endl;
    assert( u1_idx - u0_idx == 1 );
    assert( u0_idx < (int)k.size()-1);
    assert( u1_idx < (int)k.size() );
    assert( u0_idx >= 0);
    assert( u1_idx >= 1);
    return b[u0_idx*2+1];
}

ControlPoint QBSpline::get_deboor(int i) const
{
  int d_size=((b.size()+1)/2) + 1;

  if( i <= 0 )        {return b.front();}
  if( i >= d_size-1 ) {return b.back();}
  return b[2*i-1];
}

void QBSpline::getip(double t, int *i, double *p) const
{
  *i=-1;
  double den;
  for(unsigned j=0; j < k.size()-1 ;j++){
    if((k[j]<=t)&&(k[j+1]>t)){
      *i=j;
      break;
    }
  }
  if(k.back()==t){
    //this is possibly wrong
    if(closed){
      *i=0;
      *p=0.0;
      return;
    }
    else{
      *i=k.size()-2;
    }
  }

  if(*i==-1){
    cout<<"QBSpline::getIP: parameter t is out of range"<<endl;
    *i=0;
    *p=0.0;
    return;
  }
  else{
    den=dt(*i);
    if(den!=0.0){
      *p=(t-k[*i])/den;
    }
    else{
      *p=0.0;
    }
  }
}

int QBSpline::get_closest_knot(double u) const
{
    double best_dist= fabs(k[0] - u);
    double dist;
    int best=0;
    int i;
    for(i=1; i < (int)k.size(); i++){
        dist= fabs(k[i] - u);
        if(dist < best_dist){
            best_dist=dist;
            best=i;
        }
    }
    return best;
}

int QBSpline::get_closest_bezier(double u) const
{
    return get_closest_knot(u) * 2 ;
}

int QBSpline::get_edge_cps(double u0, double u1, ControlPoint &cp0, ControlPoint &cp1, ControlPoint &cp2, int &cp0_is_vert, int &cp2_is_vert) const
{
    int cp0_idx, cp1_idx, cp2_idx;
    cp0_is_vert=-1;
    cp2_is_vert=-1;
    //if(u1 <= u0) return FALSE;
    cp0_idx = get_closest_bezier(u0);
    cp2_idx = get_closest_bezier(u1);
    //printf("u0: %f  u1: %f  cp0: %i  cp2: %i spline: %x spline_max: %i\n",(float)u0,(float)u1,cp0_idx,cp2_idx,(unsigned)this,b.size()-1);
    if(cp0_idx < 0 || cp2_idx < 0 ){
        printf("QBSpline::get_edge_cps error bad u value\n");
        return FALSE;
    }
    if(((cp0_idx==0) && (cp2_idx == (int)b.size()-2))||
        ((cp2_idx==0) && (cp2_idx == (int)b.size()-2))){
        cp1_idx = b.size()-1;
    }
    else if(cp2_idx-cp0_idx==2){
        cp1_idx = cp0_idx+1;
    }
    else if(cp0_idx-cp2_idx==2){
        cp1_idx = cp0_idx-1;
    }
    else{
        printf("QBSpline::get_edge_cps error, incorrect bezier points returned\n");
        return FALSE;
    }
    cp0 = b[cp0_idx];
    cp1 = b[cp1_idx];
    cp2 = b[cp2_idx];

    /* Determine if cp0 or cp2 are BoundaryVertexs.  Return -1 for no, 0 for vertex0, and 1 for vertex1 */
    if(cp0_idx == 0) cp0_is_vert=0;
    if(cp0_idx == (int)b.size()-1) cp0_is_vert=1;
    if(cp2_idx == 0) cp2_is_vert=0;
    if(cp2_idx == (int)b.size()-1) cp2_is_vert=1;
    return TRUE;
}

Point2D QBSpline::evaluate_segment(int s, double u) const
{
  s*=2;
  double v=1-u;
  return ( (*b[s])*v*v )+( (*b[s+1])*2*u*v )+( (*b[s+2])*u*u );
}


Point2D QBSpline::evaluate(double u) const
{
  int i;
  double p;
  getip(u,&i,&p);
  return evaluate_segment(i,p);
}

void QBSpline::print() const
{
  cout<<*this;
}

void QBSpline::evaluate_full(double u,Point2D *t,Point2D *b1,Point2D *b2,int *i) const
{
  double p;
  getip(u,i,&p);
  int j = *i;
  *b1 = (*(b[2*j]))*(1-p) + (*(b[2*j+1]))*p;
  *b2 = (*(b[2*j+1]))*(1-p) + (*(b[2*j+2]))*p;
  *t = *b1 * (1-p) + *b2 * p;
}

void QBSpline::split(int i, double u, ControlPoint& cp0,
                     ControlPoint& cp1, ControlPoint& cp2)
{
    assert(i >= 0); assert((unsigned)i < k.size());
    assert(u >= 0.0); assert(u <= 1.0);

    add_knot(k[i] + u*(k[i+1]-k[i]), cp0, cp1, cp2);
}

void QBSpline::add_knot(double u,ControlPoint &cp0, ControlPoint &cp1, ControlPoint &cp2)
{
    Point2D p0,p1,p2;
    int i;
    if(u<=0 || u > k.back())
        FATAL_ERROR("QBSpline::add_knot -- knot value: "<<u<<" is out of range");

    /* Run Decasteluj's Algorithm to get new points */
    evaluate_full(u,&p1,&p0,&p2,&i);

    cp0  = b[2*i+1];
    cp0.access() = p0;

    /* allocate two more points and insert them in 'b' */
    cp1 = data_store->add_cp(p1);
    cp2 = data_store->add_cp(p2);
    b.insert(b.begin()+(2*i+2),2,cp2);
    b[2*i+2] = cp1;
    /* at this point we have:
     * b[2*i+1] = cp0 == deboor(i+1)
     * b[2*i+2] = cp1
     * b[2*i+3] = cp2 == deboor(i+2)
     * segment(i) = [2*i, 2*i+1, 2*i+2]
     * segment(i+1) = [2*i+2, 2*i+3, 2*i+4]
     */

    /* add the new knot value */
    k.insert(k.begin()+(i+1),u);

    /* make room for the new internal vertex */
    vertexs.insert(vertexs.begin() + (i),NULL);
}

bool QBSpline::remove_knot(double u, ControlPoint &cp, double &u0, double &u1){
    int i;
    double p;
    if(u<=0 || u >= k.back()) return false;
    getip(u,&i,&p);
    if(i<=0 || (unsigned)i == k.size()-1){
        cout<<"QBSpline::remove_knot: cannot remove end point"<<endl;
        return false;
    }

    /* Now the two segments on either side of the knot will become one.
     * We will remove two bezier points and the knot value.  Also, we
     * will remove the control points from the mesh.
     */
    data_store->rem_cp( b[2*i]  );
    data_store->rem_cp( b[2*i+1] );
    b.erase(b.begin() + (2*i), b.begin() + (2*i+2) );
    k.erase(k.begin()+i);
    vertexs.erase(vertexs.begin() + (i-1));

    /* update preceeding bezier point unless it is the first point of the spline
     * in which case it is a deboor point
     */
    if(i > 1){
        b[2*i-2].access()
      = *b[2*i-3] * (dt(i-1) / ( dt(i-2) + dt(i-1) ))
      + *b[2*i-1] * (dt(i-2) / ( dt(i-2) + dt(i-1) ));
    }

    /* update proceeding bezier point unless it is the last point of the spline
     * in which case it is a deboor point
     */
    if((unsigned)2*i < b.size() - 1 ){
        b[2*i].access()
          = *b[2*i-1] * (dt(i)   / ( dt(i-1) + dt(i)   ))
      + *b[2*i+1] * (dt(i-1) / ( dt(i-1) + dt(i)   ));
    }

    /* Now send back the new center control point of the final segment */
    cp = b[2*i-1];

    /* Also send back the u values for the new bezier_edge to be created
     * We can't rely on the BezierVertex's for this info, as the D0 vertexs
     * don't have any u values
     */
    //compute_bezier_points();
    u0 = k[i-1];
    u1 = k[i];
    return true;
}


void QBSpline::move(Point2D *dp, Point2D *dq, unsigned xCoord, unsigned yCoord){

#ifdef RAW_LS
  move_raw(dp,dq);
  return;

#endif
#ifndef RAW_LS
      Point2D *new_d;
      int i;

    //get new d points
    if(k.size() >= 3){
      //Solve for the displacements of the deBoor points
      new_d = compute_new_points(dp,dq);
    }
    else{
        //too small to build the matrix
        new_d = new Point2D[1];
        new_d[0] = dq[0];
    }

    // These two lines will set the endpoints appropriately, however we aren't move the endpoints
    //*get_deboor(0) = *b[0] + dp[0];
    //*get_deboor(get_num_deboor()-1) = *b[b.size()-1] + dp[k.size()-1];

    Point2D newPos;

    //    cout << "spline:move attempting w/ xC, yC = "<<xCoord<<", "<<yCoord<<endl;

    for(i=1; i < get_num_deboor()-1; i++){
      const ControlPoint cp = get_deboor(i);
      DataPoint dp = data_store -> get_data(cp);

      //      cout << "cp = "<<cp<<" :::: dp = "<<dp<<endl;

      newPos = *b[2*i-1] + new_d[i-1];
      dp->set(newPos.x(), xCoord);
      dp->set(newPos.y(), yCoord);

      //      cout << "dp is "<<dp<<endl;
    }

    //    cout << "spline:move ->  Finished moving deboor points " << endl;

    for(i=1; 2*i < ((int)b.size())-1; i++){
      const ControlPoint cp = b[2*i];
      DataPoint dp = data_store -> get_data(cp);

      //      cout << "cp = "<<cp<<" :::: dp = "<<dp<<endl;

      newPos =  *b[2*i-1] * ( dt(i)   / (dt(i-1) + dt(i)) )
    + *b[2*i+1] * ( dt(i-1) / (dt(i-1) + dt(i)) );
      dp->set(newPos.x(), xCoord);
      dp->set(newPos.y(), yCoord);

      //      cout << "dp is "<<dp<<endl;
    }

    delete[] new_d;

#endif
}

/* 
   Interpolate the motion of the DeBoor points exactly. 
   Let the Bezier Points float.
   This is designed for a solver that plans to respect differentiability
   This won't work for closed splines
 */
#ifdef RAW_LS
void QBSpline::move_raw(Point2D *dp, Point2D *dq){

  // These two lines will set the endpoints appropriately, however we aren't move the endpoints
  //*get_deboor(0) = *b[0] + dp[0];
  //*get_deboor(get_num_deboor()-1) = *b[b.size()-1] + dp[k.size()-1];

  int i;
  // Set the deboor points exactly
  for (i=1; i< get_num_deboor()-1; i++){
    *get_deboor(i) = *b[2*i-1] + dq[i-1]*2.0 - (dp[i-1] + dp[i])*0.5;
  }

  // Float the knot points to their appropriate location
  for(i=1; 2*i < ((int)b.size())-1; i++){
    *b[2*i] = *b[2*i-1] * ( dt(i)   / (dt(i-1) + dt(i)) )
      + *b[2*i+1] * ( dt(i-1) / (dt(i-1) + dt(i)) );
  }

}
#endif

Point2D* QBSpline::compute_new_points(Point2D *dp,Point2D *dq) const{
    int i,rows,cols;
    int size_dp = k.size();
    int size_dq = k.size() - 1;

    cols = size_dp - 1;
    rows = size_dp - 2 + size_dq;

    // Closed Splines have one extra constraint for differentiability at the d0 vertex
    if (closed) rows++;

    Matrix A(rows,cols);
    Point2D *B = new Point2D[rows];


    //construct matrix
    // add constraints for each p affectng the deBoor point to it's right and left
    int p_affecting;
    int col_affected1;
    int col_affected2;
    // i is constraint(row number)
    for(i=0; i < size_dp-2; i++)
      {
    p_affecting = i+1;
    col_affected1 = i;
    col_affected2 = i+1;
    A.m[i][col_affected1] = ratio(p_affecting,p_affecting);
    A.m[i][col_affected2] = ratio(p_affecting+1,p_affecting);
    B[i] = dp[p_affecting] ;
      }

    //add constraint for each dq affecting the deBoor point above it, and to it's right and left
    int q_affecting;
    int col_affected3;
    double conr = constraint_ratio;
    // i is constraint(row number)
    for(i= size_dp-2; i < size_dp- 2 + size_dq;i++)
      {
    q_affecting = i-size_dp+2;
    col_affected1 = q_affecting - 1;
    col_affected2 = q_affecting;
    col_affected3 = q_affecting + 1;
    if (q_affecting > 0)
      if (q_affecting < size_dq - 1)
        {
          // normal interior case
          A.m[i][col_affected1] = (ratio(q_affecting, q_affecting) * 0.25) * conr;
          A.m[i][col_affected2] = ((ratio(q_affecting+1,q_affecting) * 0.25) + 0.5 + (ratio(q_affecting+1, q_affecting+1) * 0.25)) * conr;
          A.m[i][col_affected3] = (ratio(q_affecting+2,q_affecting+1) * 0.25) * conr;
          B[i] = dq[q_affecting] * conr;
        }
      else
        {  //last dq case, only affects two columns
          A.m[i][col_affected1] = (ratio(q_affecting, q_affecting) * 0.25) * conr;
          A.m[i][col_affected2] = ((ratio(q_affecting+1,q_affecting) * 0.25) + 0.5)*conr ;
          //affect of the end point goes on the right hand side
          B[i] = (dq[q_affecting] - dp[size_dp-1]*0.25) * conr ;
        }
    else
      {
        // first dq case, only affects two columns
        A.m[i][col_affected2] = (0.5 + (ratio(q_affecting+1, q_affecting+1) * 0.25)) * conr;
        A.m[i][col_affected3] = (ratio(q_affecting+2,q_affecting+1) * 0.25) * conr;
        //affect of the end point goes on the right hand side
        B[i] = (dq[q_affecting] - dp[0]*0.25) * conr ;
      }
      }

    // Add the extra constraint for differentiability at d0 vertex of closed splines
    if (closed)
      {
    // Note that these 0.5 are arbitrary, otherwise the problem is nonlinear
    // This paradigm is equivalent to that of fixing the knot sequence everywhere else
    A.m[rows-1][0] = 0.5;
    A.m[rows-1][cols-1] = 0.5;
    // Note this constraint is not solely in terms of displacements
    B[rows-1] = dp[0] + *b[0] - *b[1]*0.5 - *b[b.size()-2]*0.5;

    // This is a numerical trick to make this constraint big times as important as the others
    A.m[rows-1][0] = A.m[rows-1][0] * big;
    A.m[rows-1][cols-1] = A.m[rows-1][cols-1] * big;
    B[rows-1] = B[rows-1] * big;
      }

        /* Print out matrix and vector for debugging purposes
    if (closed)
    for (int x=0;x<rows;x++)
      {
    for (int y=0;y<cols;y++)
      cout << A.m[x][y] << "\t";
    cout << B[x];
    cout << endl;
      }
       */

    /* this solver system is set up so that multiple vectors can
     * be solved using the same matrix.  Right now only x and
     * y coordinates are being solved for, but if there is more
     * data in the future this will be trivial to upgrade and still
     * maintain optimal efficiency (the matrix is factored only once
     * and used to solve many systems)  I can change this to be
     * faster if we don't need this robustness
     */

    //construct array of vectors to be solved
    int num_vectors = 2;
    double **bs = new double*[num_vectors];
    for(i=0; i < num_vectors; i++){
        bs[i] = new double[rows];
    }

    for(i=0; i < rows; i++){
        bs[0][i] = B[i][0];
        bs[1][i] = B[i][1];
    }

    double **xs;
    if(k.size() > 4){
        //send matrix and array of vectors to solver
        xs = solve_system(A,bs,cols,num_vectors);
    } else {
        //send to matrix inverse solver since it is too small for the
        //band diagonal solver
        xs = solve_small_system(A,bs,num_vectors);
    }

    //destroy array of vectors to be solved;
    for(i=0; i < num_vectors; i++) delete[] bs[i];
    delete[] bs;

    //create new points
    Point2D *nd=new Point2D[cols];
    for(i=0;i<cols;i++){
        nd[i].assign(xs[0][i], xs[1][i]);
    }

    //destroy array of solution vectors
    for(i=0; i < num_vectors; i++) delete[] xs[i];
    delete[] xs;
    delete[] B;

    return nd;
}

double** QBSpline::solve_system(Matrix& A,double **bs, int cols,int num_vectors) const{
    int i,j;

    double *B;
    double fact;
    double **xs=new double*[num_vectors];
    double *u=new double[cols];

    Matrix *ATA  = A.transpose_times_self();
    Matrix5 M(*ATA);

    double alpha = ATA->m[cols-1][0];
    double beta  = ATA->m[0][cols-1];
    double gamma = -M.d3[0];

    /* if closed we need to solve the cyclic banded case */
    if(closed){
        M.d3[0]      -= gamma;
        M.d3[cols-1] -= alpha * beta/gamma;
        for(i=0; i < num_vectors; i++){
            xs[i] = new double[cols];

            u[0] = gamma;
            u[cols-1] = alpha;
            for(j=1; j < cols-1; j++) { u[j] = 0.0; }

            B = A.transpose_times_vector(bs[i]);
            B = M.LUsolve(B);
            u = M.LUsolve(u);

            fact = (B[0] + beta * (B[cols-1]/gamma) ) / (1.0 + u[0] + beta * (u[cols-1]/gamma));
            for(j=0; j < cols; j++){
                xs[i][j] = B[j] - (fact * u[j]);
            }
            delete[] B;
        }
    }

    /* otherwise it is a simple banded case to solve */
    else {
        for(i=0; i < num_vectors; i++){
            B = A.transpose_times_vector(bs[i]);
            B = M.LUsolve(B);
            xs[i] = B; /* use vectors returned by transpose_times_vector for xs (will be deleted later)*/
        }
    }

    delete ATA;
    delete[] u;
    return xs;
}

double** QBSpline::solve_small_system(Matrix &A,double **bs, int num_vectors) const{
    int i;
    double **xs=new double*[num_vectors];
    double *B;

    Matrix *ATA = A.transpose_times_self();
    Matrix *ATAinv = ATA->inverse();

    for(i=0;i<num_vectors;i++){
        B = A.transpose_times_vector(bs[i]);
        xs[i] = ATAinv->times_vector(B);
        delete[] B;
    }
    delete ATA;
    delete ATAinv;

    return xs;
}

void QBSpline::douglas_peucker(double tol, hashers::hash_set<BezierVertex*> &keep) const
{
    unsigned i;
    int nk = k.size();
    if(nk <= 2) return; /* no internal knots */

    /* if no tolerance, keep all knots */
    if(tol== 0.0){
        for(i=0; i< (unsigned)nk-2; i++) keep.insert( vertexs[i] );
        return;
    }

    double top, bot, left, right;
    bbox(&top, &bot, &left, &right);
    double size = MAX(fabs(top-bot), fabs(left-right));
    tol = size*tol;

    /* we create a vector of bezier points who coorespond to
     * a BezierVertex.  This has size =  k.size().
     */
    bool mark[nk];

    /* true=throw away, false=keep*/
    for(i=0; i < k.size(); i++) mark[i]=true;
    simplify(tol, mark, 0, nk-1);

    /*never remove first or last vertex */
    mark[0] = false;
    mark[nk-1] = false;

    /* add all keeper vertexs, we can skip endpoints
     * as they are already added
     */
    for(i=1; i< (unsigned)nk-1; i++){
        if(!mark[i]) keep.insert(vertexs[i-1]);
    }
}

void QBSpline::simplify(double tol,bool *mark, int j,int z) const
{
     if(z == j + 1) return;
     int maxi=0;
     double maxd=0.0;
     double dv;
     Point2D pj = *b[2*j];
     Point2D pz = *b[2*z];

    for(int i=j+1; i < z; i++){
        dv =  b[2*i]->dist_from_line(pj,pz);
        if(dv >= maxd){
            maxd = dv;
            maxi = i;
        }
    }
    if( maxd > tol ){
        mark[maxi] = false;
        simplify(tol, mark, j, maxi);
        simplify(tol, mark, maxi, z);
    }
}

void QBSpline::bbox(double *top,double *left,double *bot,double *right) const
{
    *left=*right=(*(b[0]))[0];
    *top=*bot=(*(b[0]))[1];
    for(unsigned i=1; i<b.size(); i++){
        if( (*(b[i]))[0] < *left ) { *left = (*(b[i]))[0]; }
        if( (*(b[i]))[0] > *right ) { *right = (*(b[i]))[0]; }
        if( (*(b[i]))[1] > *top ) { *top = ( *(b[i]) )[1]; }
        if( (*(b[i]))[1] < *bot ) { *bot = ( *(b[i]) )[1]; }
    }
}

int QBSpline::get_num_deboor() const
{
    return ((b.size()+1)/2)+1;
}

int QBSpline::get_num_segments() const
{
    return k.size() - 1;
}

void QBSpline::set_bezier_vertex(double u, BezierVertex *v)
{
    vertexs[get_closest_knot(u)-1]=v;
}

void QBSpline::set_bezier_vertex(int idx, BezierVertex *v)
{
    //exclude first and last knots
    assert( (idx >= 1) && (idx <= (int)k.size()-2));
    vertexs[idx-1]=v;
}

unsigned QBSpline::get_bezier_vertex_idx(const BezierVertex *v) const
{
    unsigned j;
    std::vector<BezierVertex*>::const_iterator i;
    for(j=0, i=vertexs.begin(); i!= vertexs.end(); ++i, j++){
        if( v == *i) return j+1; //add one for compatability with set_bezier_vertex
    }
    return 0;
}

BezierVertex *QBSpline::get_bezier_vertex(unsigned idx) const
{
  return vertexs[idx - 1];
}

ControlPoint QBSpline::get_control_point_at_idx(unsigned idx) const
{
    assert( idx < k.size() );
    return b[idx*2];
}

ControlPoint QBSpline::get_edge_center_cp(unsigned u0_idx, unsigned u1_idx) const
{
    unsigned idx;
    assert( abs( (int) u0_idx - (int) u1_idx) == 1 );
    assert( u0_idx < k.size());
    assert( u1_idx < k.size());
    idx = u0_idx + u1_idx;
    return b[idx];
}

double QBSpline::get_u_at_idx(unsigned idx) const
{
    assert(idx < k.size());
    return k[idx];
}

double QBSpline::dt(int i) const
{
  return k[i+1]-k[i];
}

double QBSpline::ratio(int a,int b) const
{
  return (k[a]-k[a-1])/(k[b+1]-k[b-1]);
}

ostream& operator<<(ostream& stream,const QBSpline& qbs)
{
    unsigned i;
    unsigned size_b=qbs.b.size();
    unsigned size_k=qbs.k.size();

    stream.precision(3);
    stream<<"QBSpline:= { closed="<<qbs.closed<<", restlength="<<qbs.restlength<<","<<endl;
    stream<<"size_b="<< size_b <<", b=[";
    for(i=0; i<size_b-1; i++) { stream<<*(qbs.b[i])<<", "; }
    stream<<*(qbs.b[size_b-1])<<"],"<<endl;

    stream<<"size_k="<<size_k <<", k=[";
    for(i=0; i<size_k-1; i++){ stream<<qbs.k[i]<<", "; }
    stream<<qbs.k[size_k-1]<<"]"<<endl;
    stream<<"Interior Vertexs(num="<<size_k-2<<"): [End-Point, ";
    for(i=0; i<size_k-2; i++){ stream<<qbs.vertexs[i]<<", "; }
    stream<<"End-Point]}"<<endl;
    return stream;
}

---