emptyspline.C

Go to the documentation of this file.

#ifndef _EMPTYSPLINE_C
#define _EMPTYSPLINE_C

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif /* HAVE_CONFIG_H */


#include <iostream>
#include <cstddef>

#include "boundarymesh.h"
#include "emptyspline.h"

EmptySpline::EmptySpline(const vector<ControlPoint> &deboor, BoundaryMesh *_bdry_mesh, bool _closed){
    vector<double> newk;
    initialize(deboor, newk, _bdry_mesh, 1.0, _closed);
}

EmptySpline::EmptySpline(const vector<ControlPoint> &deboor, const vector<double> &newk, BoundaryMesh *_bdry_mesh, bool _closed){
    initialize(deboor, newk, _bdry_mesh, 1.0, _closed);
}

EmptySpline::EmptySpline(const vector<ControlPoint> &deboor, const vector<double> &newk, BoundaryMesh *_bdry_mesh, double _restlength, bool _closed)
{
    initialize(deboor, newk,_bdry_mesh, _restlength, _closed);
}

void EmptySpline::initialize(const vector<ControlPoint> &deboor, const vector<double> &newk, BoundaryMesh *_bdry_mesh, double _restlength, bool _closed)
{
  unsigned i;
  bdry_mesh = _bdry_mesh;
  restlength = _restlength;

  /* A closed spline is a spline where the first and last deboor point
   * is the same.  Open splines must have at least 3 deboor points,
   * closed splines must have at least 5 deboor points
   */
   closed=_closed;
   
   if( closed  && deboor.size() < 5 ) {
     cout<<"Error: closed spline has "<<deboor.size()<<" points but must have at least 5."<<endl;
     exit(1);
   }
    
   if( deboor.size() < 3 ) {
     cout<<"Error: spline has "<<deboor.size()<<" points but must have at least 3."<<endl;
     exit(1);
   }
  
  /* Initialize the vertexs vector */
  vertexs.insert(vertexs.begin(), deboor.size()-3, NULL); 

  /* if newk is empty, make our own 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());
    }
  }    

  /* Store knot vector.  There should be d.size()-1 knots. */
  else{ 
    if(newk.size() != deboor.size()-1){
    printf("Error: spline has %i deboor points but %i knots (sould be deboor-1)\n",deboor.size(),newk.size());
    exit(1);
    }
    k=newk;
  }

  /* Populate the b vector there are 2*deboor.size()-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( bdry_mesh->add_cp(new_point) );
    b.push_back(temp);
  }
  b.push_back(deboor.back());
}

EmptySpline::EmptySpline(const EmptySpline &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;
    bdry_mesh = other.bdry_mesh;    
}

EmptySpline& EmptySpline::operator=(const EmptySpline &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;
    bdry_mesh = other.bdry_mesh;
    return *this;    
}

void EmptySpline::compute_knot_vector()
{
  k[0]=0.0;
  if(b.size() == 3){ 
    k[1] = (*(b[2]) - *(b[0])).mag();
  } else {
    k[1]=( *(b[3]) - *(b[0])).mag();
    for(unsigned i=2; i < k.size()-1; i++){
      k[i]=k[i-1] + (*(b[2*i+1]) - *(b[2*i-3])).mag();
    }
    k[k.size()-1] = k[k.size()-2]+ ( *(b[ b.size()-4 ])- *(b[ b.size()-1 ]) ).mag();
  }
}

void EmptySpline::compute_bezier_points()
{
  for(unsigned i=1; 2*i < 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)));
  }
}

ControlPoint EmptySpline::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 EmptySpline::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 EmptySpline::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<<"EmptySpline::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 EmptySpline::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 EmptySpline::get_closest_bezier(double u) const
{    
    return get_closest_knot(u) * 2 ;
}

/* This is code for a reverse lookup of the control points on a spline given the u values into it.
 * It is necessary in order to rebuild a mesh from the ruby output file format
 * Also, we require u0 < u1 for this to work
 * Furthermoore, this function will check if cp0 or cp2 are at the ends of the spline.  This data will
 * be returned via the cp0_is_vert and cp2_is_vert variables.  
 * We let -1=not a vertex
 *         0=vertex0
 *         1=vertex1
 *
 * This way we can determine where each BezierVertex is in respect to the boundary and can thus set the
 * boundary data correctly with the BezierVertex::set_bdry() functions. 
 *
 * This whole function is only used when reading in from the ruby file format.  This is the conssession
 * we must make for not having cached control point values.  Since it is only called during the read
 * operation the inefficientcies here are not important. 
 */
int EmptySpline::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("EmptySpline::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("EmptySpline::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 EmptySpline::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 EmptySpline::evaluate(double u) const
{
  int i;
  double p;
  getip(u,&i,&p);
  return evaluate_segment(i,p);
}

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


void EmptySpline::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 EmptySpline::add_knot(double u,ControlPoint &cp0, ControlPoint &cp1, ControlPoint &cp2)
{
  Point2D p0,p1,p2;
  int i;
  if(u<=0 || u > k.back()){return;}
  evaluate_full(u,&p1,&p0,&p2,&i);
  
  cp0 = b[2*i+1];
  *cp0 = p0; 
  cp1 = bdry_mesh->add_cp(p1);
  cp2 = bdry_mesh->add_cp(p2);
  
  b.insert(b.begin()+(2*i+2),2,cp2);
  b[2*i+2] = cp1;

  k.insert(k.begin()+(i+1),u);  
  vertexs.insert(vertexs.begin() + (i),NULL);
  compute_bezier_points();
}

/* This function must be called AFTER the Bezier Edges and vertexs 
 * associated with the knot to be removed have been destroyed.  This is
 * because the data points are owned by the bezier Edges and vertexs and
 * so must be deleted before the control points are deleted by this function.  
 */
bool EmptySpline::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<<"EmptySpline::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.
     */
    bdry_mesh->rem_cp( b[2*i]  );
    bdry_mesh->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] =   *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]   =   *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; 
}

//This function does NOT move the endpoint(s)
//Although it does respect their apparent motion as given in dp
void EmptySpline::move(Point2D *dp, Point2D *dq){
    Point2D *new_d;
    int i;

    /* Print out dp and dq for debugging purposes
    cout << "-------------------------"<<endl;
    cout << "dp: ";
    for (i=0;i<k.size();i++)
      cout << dp[i] << " - ";
    cout << endl << "=====" << endl<<"dq: ";
    for (i=0;i<k.size()-1;i++)
      cout << dq[i] << " - ";
    cout <<endl;
    */

    //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];
    }
    

    /*
    cout << "Printing new vectors" << endl;
    for (int i=0; i< k.size()-1;i++)
      {
    cout << new_d[i] << " - ";
      }
    cout << endl;
    */

    // 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];

    for(i=1; i < get_num_deboor()-1; i++){
      *get_deboor(i) = *b[2*i-1] + new_d[i-1];
    }
    compute_bezier_points();
    delete[] new_d;


    /*
    //debugging print errors on differentiability condition
    if (closed)
      {
    Point2D err = *b[1]*0.5 + *b[b.size()-2]*0.5;
    err -= *b[0];
    cout << "Diff error is: "<<err<<endl;
      }
    */
}

Point2D* EmptySpline::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;
}
/*
 this function takes constraint matrix A and array of vectors to solve,
 bs. A and bs are not modified.  The function then computes A^T * A and
 stores it as a Matrix5 (5 banded matrix).  Then Crout LU decomposition
 and a LUsolver [both are in Matrix5] are used to solve for each vector.
 Since A^T*A is actually a cyclic banded matrix the Sherman-Morrison
 Formula is used to correct for the errors caused by the non-zero 
 off-diagonal elements.  Finally the function will return an array of
 solution vectors.
*/
double** EmptySpline::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** EmptySpline::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 EmptySpline::douglas_peucker(double tol, hash_set<BezierVertex*, CellHasher<BezierVertex>, CellEqual<BezierVertex> > &keep) const{
    int nk = k.size(); 
    if(nk <= 2) return;
    
    /* we create a vector of bezier points who coorespond to
     * a BezierVertex.  This has size =  k.size().
     */
    bool *mark = new bool[nk];
    unsigned i;

    /* 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 ednpoints
     * as they are already added 
     */
    for(i=1; i< (unsigned)nk-1; i++){
        if(!mark[i]){
            keep.insert(vertexs[i-1]);
        }
    }
    delete[] mark;
}

void EmptySpline::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 EmptySpline::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]; }
  }
}




ostream& operator<<(ostream& stream,const EmptySpline& qbs)
{
    unsigned i;
    unsigned size_b=qbs.b.size();
    unsigned size_k=qbs.k.size();
  
    stream.precision(3);
    stream<<"EmptySpline:= { 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;
}

#endif

---