---

gunu_raise_degree.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_vector.h"
#include "gunu_raise_degree.h"

namespace gunu {

using gul::rtr;
using gul::Max;
using gul::Min;
using gul::Ptr;
using gul::point;
using gul::hpoint;
using gul::point2;
using gul::hpoint2;
using gul::point1;
using gul::hpoint1;
using gul::set;

/*-------------------------------------------------------------------------*//**
  raise the degree 'p' of a curve to 'p+t'. this function doesn't reserves
  memory for Uh and Q, in the worst case for Q: p + (n-p+1)*(t+1), and for
  Uh: 2*p + (n-p+2)*(t+1) elements are needed, the function returns the new 
  number of controlpoints - 1                                                 */
/*----------------------------------------------------------------------------*/
template< class T, class HP >
int ElevateCurveDegree( int n, int p, const Ptr<T>& U, const Ptr<HP>& Pw, int t,
                        Ptr<T>& Uh, Ptr<HP>& Qw )
{
  int m,i,mpi,j,ph,ph2,mh,kind,r,a,b,cind;
  int mul,oldr,lbz,rbz,k,save,s,first,last,kj,tr;
  T inv,ua,ub,numer,den,bet,alf,gam;
  Ptr< Ptr<T> > bezalfs;
  Ptr<T> alfs;
  Ptr<HP> bpts,ebpts,Nextbpts;

  bezalfs.reserve_pool(p+t+1);
  for( i = 0; i <= p+t; i++ )
    bezalfs[i].reserve_pool(p+1);

  bpts.reserve_pool(p+1);
  ebpts.reserve_pool(p+t+1);

  if( p > 1 )
  {
    Nextbpts.reserve_pool(p-1);
    alfs.reserve_pool(p-1); 
  }
  m = n + p + 1;
  ph = p + t;
  ph2 = ph / 2;

/* --- calc coefficients for raising the degree of the bezier segments ---- */  

  bezalfs[0][0] = (T)1;
  bezalfs[ph][p] = (T)1;
  
  for( i = 1; i <= ph2; i++ )
  {
    inv = (T)1 / rtr<T>::BinCoeff(ph,i);
    mpi = Min( p, i );
    
    for( j = Max( 0,i-t ); j <= mpi; j++ )
      bezalfs[i][j] = inv * rtr<T>::BinCoeff(p,j) * rtr<T>::BinCoeff(t,i-j);
  }    
  for( i = ph2+1; i <= ph-1; i++ )
  {
    mpi = Min( p, i );
    for( j = Max( 0,i-t ); j <= mpi; j++ )
      bezalfs[i][j] = bezalfs[ph-i][p-j];
  }       
  mh = ph;
  kind = ph+1;
  r = -1;
  a = p;
  b = p+1;
  cind = 1;
  ua = U[0];
  Qw[0] = Pw[0];
  
  for( i = 0; i <= ph; i++ )
    Uh[i] = ua;
    
  for( i = 0; i <= p; i++ )
    bpts[i] = Pw[i];  
    
  while( b < m )
  {
    i = b;
    while( (b < m) && (U[b] == U[b+1]) )
      b++;

    mul = b-i+1;
    mh = mh + mul + t;
    ub = U[b];
    oldr = r;
    r = p - mul;
    
    if( oldr > 0)
      lbz = (oldr+2) / 2;
    else
      lbz = 1;

    if( r > 0 )
      rbz = ph - ((r+1) / 2);
    else
      rbz = ph;  

    if( r > 0 )
    {
      numer = ub - ua;
      for( k = p; k > mul; k-- )
        alfs[k-mul-1] = numer / (U[a+k]-ua);
      for( j = 1; j <= r; j++ )  
      {
        save = r - j;
        s = mul + j;            

        for( k = p; k >= s; k-- )
          bpts[k] = ((T)1 - alfs[k-s]) * bpts[k-1] + alfs[k-s] * bpts[k];

        Nextbpts[save] = bpts[p];
      }  
    }
    for( i = lbz; i <= ph; i++ )
    {
      set( ebpts[i], (T)0 );
      mpi = Min( p, i );

      for( j = Max(0,i-t); j <= mpi; j++ )
        ebpts[i] = ebpts[i] + bezalfs[i][j] * bpts[j];
    }

    if( oldr > 1 )
    {
      first = kind-2;
      last = kind;
      den = ub-ua;
      bet = (ub-Uh[kind-1])/den;
      
      for( tr = 1; tr < oldr; tr++ )
      {        
        i = first;
        j = last;
        kj = j-kind+1;
        while( j - i > tr )
        {
          if( i < cind )
          {
            alf = (ub-Uh[i])/(ua-Uh[i]);
            Qw[i] = ((T)1 - alf) * Qw[i-1] + alf * Qw[i];
          }        
          if( j >= lbz )
          {
            if( j - tr <= kind-ph+oldr )
            {  
              gam = (ub-Uh[j-tr])/den;
              ebpts[kj] = ((T)1 - gam) * ebpts[kj+1] + gam * ebpts[kj];
            }
            else
              ebpts[kj] = ((T)1 - bet) * ebpts[kj+1] + bet * ebpts[kj];
          }
          i++;
          j--;
          kj--;
        }      
        
        first--;
        last++;
      }                    
    }    
    if( a != p )
      for( i = 0; i < ph-oldr; i++ )
      {
        Uh[kind] = ua;
        kind++;
      }      
    for( j = lbz; j <= rbz; j++ )
    {
      Qw[cind] = ebpts[j];
      cind++;
    }
    if( b < m )
    {
      for( j = 0; j < r; j++ )
        bpts[j] = Nextbpts[j];
      for( j = r; j <= p; j++ )
        bpts[j] = Pw[b-p+j];
      a = b;
      b++;
      ua = ub;
    }
    else
      for( i = 0; i <= ph; i++ )
        Uh[kind+i] = ub;
  }                      
  return( mh - ph - 1 );
}          
template int ElevateCurveDegree( 
   int n, int p, const Ptr<float>& U, const Ptr< point<float> >& Pw, int t,
   Ptr<float>& Uh, Ptr< point<float> >& Qw );
template int ElevateCurveDegree( 
   int n, int p, const Ptr<float>& U, const Ptr< hpoint<float> >& Pw, int t,
   Ptr<float>& Uh, Ptr< hpoint<float> >& Qw );
template int ElevateCurveDegree( 
   int n, int p, const Ptr<float>& U, const Ptr< point2<float> >& Pw, int t,
   Ptr<float>& Uh, Ptr< point2<float> >& Qw );
template int ElevateCurveDegree( 
   int n, int p, const Ptr<float>& U, const Ptr< hpoint2<float> >& Pw, int t,
   Ptr<float>& Uh, Ptr< hpoint2<float> >& Qw );

template int ElevateCurveDegree( 
   int n, int p, const Ptr<double>& U, const Ptr< point<double> >& Pw, int t,
   Ptr<double>& Uh, Ptr< point<double> >& Qw );
template int ElevateCurveDegree( 
   int n, int p, const Ptr<double>& U, const Ptr< hpoint<double> >& Pw, int t,
   Ptr<double>& Uh, Ptr< hpoint<double> >& Qw );
template int ElevateCurveDegree( 
   int n, int p, const Ptr<double>& U, const Ptr< point2<double> >& Pw, int t,
   Ptr<double>& Uh, Ptr< point2<double> >& Qw );
template int ElevateCurveDegree( 
   int n, int p, const Ptr<double>& U, const Ptr< hpoint2<double> >& Pw, int t,
   Ptr<double>& Uh, Ptr< hpoint2<double> >& Qw );

/*-------------------------------------------------------------------------*//**
  raise the degree 'pu' or 'pv' of a surface by 't'. this function doesn't
  reserves memory for Uh,Vh and Qw, they must be dimensioned in u or v 
  direction like for curves.                                                  */
/*----------------------------------------------------------------------------*/
template< class T, class HP >
void ElevateSurfaceDegreeU( 
              int nu, int pu, const Ptr<T>& U, int nv, int pv, const Ptr<T>& V,
              const Ptr< Ptr<HP> >& Pw, int t,
              int *nhu, Ptr<T>& Uh, int *nhv, Ptr<T>& Vh,
              Ptr< Ptr<HP> >& Qw )
{
  int m,i,mpi,j,ph,ph2,mh,kind,r,a,b,cind,row;
  int mul,oldr,lbz,rbz,k,save,s,first,last,kj,tr;
  T inv,ua,ub,numer,den,bet,alf,gam;
  Ptr< Ptr<T> > bezalfs;
  Ptr<T> alfs;
  Ptr< Ptr< HP > > bpts,ebpts,Nextbpts;

  bezalfs.reserve_pool(pu+t+1);
  for( i = 0; i <= pu+t; i++ )
    bezalfs[i].reserve_pool(pu+1);

  bpts.reserve_pool(nv+1);
  for( i = 0; i <= nv; i++ )
    bpts[i].reserve_pool(pu+1);

  ebpts.reserve_pool(nv+1);
  for( i = 0; i <= nv; i++ )
    ebpts[i].reserve_pool(pu+t+1);

  if( pu > 1 )
  {
    Nextbpts.reserve_pool(nv+1);
    for( i = 0; i <= nv; i++ )
      Nextbpts[i].reserve_pool(pu-1);

    alfs.reserve_pool(pu-1);
  }
  m = nu + pu + 1;
  ph = pu + t;
  ph2 = ph / 2;
  
/* --- calc coefficients for raising the degree of the bezier segments ---- */  

  bezalfs[0][0] = 1.0;
  bezalfs[ph][pu] = 1.0;
  
  for( i = 1; i <= ph2; i++ )
  {
    inv = (T)1 / rtr<T>::BinCoeff( ph, i );
    mpi = Min( pu, i );
    
    for( j = Max( 0,i-t ); j <= mpi; j++ )
      bezalfs[i][j] = inv * rtr<T>::BinCoeff(pu,j) * rtr<T>::BinCoeff(t,i-j);
  }    
  for( i = ph2+1; i <= ph-1; i++ )
  {
    mpi = Min( pu, i );
    for( j = Max( 0,i-t ); j <= mpi; j++ )
    bezalfs[i][j] = bezalfs[ph-i][pu-j];
  }       
  mh = ph;
  kind = ph+1;
  r = -1;
  a = pu;
  b = pu+1;
  cind = 1;
  ua = U[0];
  for( row = 0; row <= nv; row++ )
    Qw[row][0] = Pw[row][0];
  
  for( i = 0; i <= ph; i++ )
    Uh[i] = ua;
    
  for( row = 0; row <= nv; row++ )
    for( i = 0; i <= pu; i++ )
      bpts[row][i] = Pw[row][i];  
    
  while( b < m )
  {
    i = b;
    while( (b < m) && (U[b] == U[b+1]) )
      b++;

    mul = b-i+1;
    mh = mh + mul + t;
    ub = U[b];
    oldr = r;
    r = pu - mul;
    
    if( oldr > 0)
      lbz = (oldr+2) / 2;
    else
      lbz = 1;

    if( r > 0 )
      rbz = ph - ((r+1) / 2);
    else
      rbz = ph;  

    if( r > 0 )
    {
      numer = ub - ua;
      for( k = pu; k > mul; k-- )
        alfs[k-mul-1] = numer / (U[a+k]-ua);
      for( j = 1; j <= r; j++ )  
      {
        save = r - j;
        s = mul + j;            

        for( row = 0; row <= nv; row++ )
        {
          for( k = pu; k >= s; k-- )
          {
            bpts[row][k] = 
               ((T)1 - alfs[k-s]) * bpts[row][k-1] + alfs[k-s] * bpts[row][k];
          }  
          Nextbpts[row][save] = bpts[row][pu];
        }
      }  
    }
    for( row = 0; row <= nv; row++ )
    {
      for( i = lbz; i <= ph; i++ )
      {
        set( ebpts[row][i], (T)0 );
        mpi = Min( pu, i );
        for( j = Max(0,i-t); j <= mpi; j++ )
          ebpts[row][i] = ebpts[row][i] + bezalfs[i][j] * bpts[row][j];      
      }
    }
    if( oldr > 1 )
    {
      first = kind-2;
      last = kind;
      den = ub-ua;
      bet = (ub-Uh[kind-1])/den;
      
      for( tr = 1; tr < oldr; tr++ )
      {        
        i = first;
        j = last;
        kj = j-kind+1;
        while( j - i > tr )
        {
          if( i < cind )
          {
            alf = (ub-Uh[i])/(ua-Uh[i]);
            for( row = 0; row <= nv; row++ )
              Qw[row][i] = ((T)1 - alf) * Qw[row][i-1] + alf * Qw[row][i];
          }        
          if( j >= lbz )
          {
            if( j - tr <= kind-ph+oldr )
            {  
              gam = (ub-Uh[j-tr])/den;
              for( row = 0; row <= nv; row++ )
              {
                ebpts[row][kj] = 
                     ((T)1 - gam) * ebpts[row][kj+1] + gam * ebpts[row][kj];
              }
            }
            else
            {
              for( row = 0; row <= nv; row++ )
              {
                ebpts[row][kj] = 
                     ((T)1 - bet) * ebpts[row][kj+1] + bet * ebpts[row][kj];
              }
            }
          }
          i++;
          j--;
          kj--;
        }      
        
        first--;
        last++;
      }                    
    }    
    if( a != pu )
      for( i = 0; i < ph-oldr; i++ )
      {
        Uh[kind] = ua;
        kind++;
      }
    for( j = lbz; j <= rbz; j++ )
    {
      for( row = 0; row <= nv; row++ )
        Qw[row][cind] = ebpts[row][j];
      cind++;
    }
    if( b < m )
    {
      for( row = 0; row <= nv; row++ )
      {
        for( j = 0; j < r; j++ )
          bpts[row][j] = Nextbpts[row][j];
        for( j = r; j <= pu; j++ )
          bpts[row][j] = Pw[row][b-pu+j];
      }
      a = b;
      b++;
      ua = ub;
    }
    else
      for( i = 0; i <= ph; i++ )
        Uh[kind+i] = ub;
  }                      
  *nhu = mh - ph - 1;

  // knot vector in V direction only needs to be copied
  for( i = 0; i < nv + pv + 2; i++ )
    Vh[i] = V[i];

  *nhv = nv;  
}          
// template instantiation
template void ElevateSurfaceDegreeU( 
      int nu, int pu, const Ptr<float>& U, int nv, int pv, const Ptr<float>& V,
      const Ptr< Ptr< point<float> > >& Pw, int t,
      int *nhu, Ptr<float>& Uh, int *nhv, Ptr<float>& Vh,
      Ptr< Ptr< point<float> > >& Qw );
template void ElevateSurfaceDegreeU( 
      int nu, int pu, const Ptr<float>& U, int nv, int pv, const Ptr<float>& V,
      const Ptr< Ptr< hpoint<float> > >& Pw, int t,
      int *nhu, Ptr<float>& Uh, int *nhv, Ptr<float>& Vh,
      Ptr< Ptr< hpoint<float> > >& Qw );
template void ElevateSurfaceDegreeU( 
      int nu, int pu, const Ptr<float>& U, int nv, int pv, const Ptr<float>& V,
      const Ptr< Ptr< point1<float> > >& Pw, int t,
      int *nhu, Ptr<float>& Uh, int *nhv, Ptr<float>& Vh,
      Ptr< Ptr< point1<float> > >& Qw );
template void ElevateSurfaceDegreeU( 
      int nu, int pu, const Ptr<float>& U, int nv, int pv, const Ptr<float>& V,
      const Ptr< Ptr< hpoint1<float> > >& Pw, int t,
      int *nhu, Ptr<float>& Uh, int *nhv, Ptr<float>& Vh,
      Ptr< Ptr< hpoint1<float> > >& Qw );

template void ElevateSurfaceDegreeU( 
    int nu, int pu, const Ptr<double>& U, int nv, int pv, const Ptr<double>& V,
    const Ptr< Ptr< point<double> > >& Pw, int t,
    int *nhu, Ptr<double>& Uh, int *nhv, Ptr<double>& Vh,
    Ptr< Ptr< point<double> > >& Qw );
template void ElevateSurfaceDegreeU( 
    int nu, int pu, const Ptr<double>& U, int nv, int pv, const Ptr<double>& V,
    const Ptr< Ptr< hpoint<double> > >& Pw, int t,
    int *nhu, Ptr<double>& Uh, int *nhv, Ptr<double>& Vh,
    Ptr< Ptr< hpoint<double> > >& Qw );
template void ElevateSurfaceDegreeU( 
    int nu, int pu, const Ptr<double>& U, int nv, int pv, const Ptr<double>& V,
    const Ptr< Ptr< point1<double> > >& Pw, int t,
    int *nhu, Ptr<double>& Uh, int *nhv, Ptr<double>& Vh,
    Ptr< Ptr< point1<double> > >& Qw );
template void ElevateSurfaceDegreeU( 
    int nu, int pu, const Ptr<double>& U, int nv, int pv, const Ptr<double>& V,
    const Ptr< Ptr< hpoint1<double> > >& Pw, int t,
    int *nhu, Ptr<double>& Uh, int *nhv, Ptr<double>& Vh,
    Ptr< Ptr< hpoint1<double> > >& Qw );


}


        

---