---

gunu_raise_degree2.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 '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 ElevateSurfaceDegreeV( 
             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,col;
  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(pv+t+1);
  for( i = 0; i <= pv+t; i++ ) 
    bezalfs[i].reserve_pool(pv+1);

  bpts.reserve_pool(pv+1);
  for( i = 0; i <= pv; i++ )
    bpts[i].reserve_pool(nu+1);

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

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

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

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

    mul = b-i+1;
    mh = mh + mul + t;
    ub = V[b];
    oldr = r;
    r = pv - 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 = pv; k > mul; k-- )
        alfs[k-mul-1] = numer / (V[a+k]-ua);
      for( j = 1; j <= r; j++ )  
      {
        save = r - j;
        s = mul + j;            

        for( k = pv; k >= s; k-- )
        {
          for( col = 0; col <= nu; col++ )
          {
            bpts[k][col] = 
              ((T)1 - alfs[k-s]) * bpts[k-1][col] + alfs[k-s] * bpts[k][col];
          }  
        }
        for( col = 0; col <= nu; col++ )
          Nextbpts[save][col] = bpts[pv][col];
      }  
    }
    for( i = lbz; i <= ph; i++ )
    {
      for( col = 0; col <= nu; col++ )
        set( ebpts[i][col], (T)0 );

      mpi = Min( pv, i );
      for( j = Max(0,i-t); j <= mpi; j++ )
      {
        for( col = 0; col <= nu; col++ )
          ebpts[i][col] = ebpts[i][col] + bezalfs[i][j] * bpts[j][col];
      }      
    }
    if( oldr > 1 )
    {
      first = kind-2;
      last = kind;
      den = ub-ua;
      bet = (ub-Vh[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-Vh[i])/(ua-Vh[i]);
            for( col = 0; col <= nu; col++ )
              Qw[i][col] = ((T)1 - alf) * Qw[i-1][col] + alf * Qw[i][col];
          }        
          if( j >= lbz )
          {
            if( j - tr <= kind-ph+oldr )
            {  
              gam = (ub-Vh[j-tr])/den;
              for( col = 0; col <= nu; col++ )
              {
                ebpts[kj][col] = 
                      ((T)1 - gam) * ebpts[kj+1][col] + gam * ebpts[kj][col];
              }
            }
            else
            {
              for( col = 0; col <= nu; col++ )
              {
                ebpts[kj][col] =
                      ((T)1 - bet) * ebpts[kj+1][col] + bet * ebpts[kj][col];
              }
            }
          }
          i++;
          j--;
          kj--;
        }      
        first--;
        last++;
      }                    
    }    
    if( a != pv )
      for( i = 0; i < ph-oldr; i++ )
      {
        Vh[kind] = ua;
        kind++;
      }
    for( j = lbz; j <= rbz; j++ )
    {
      for( col = 0; col <= nu; col++ )
        Qw[cind][col] = ebpts[j][col];
      cind++;
    }
    if( b < m )
    {
      for( j = 0; j < r; j++ )
        for( col = 0; col <= nu; col++ )
          bpts[j][col] = Nextbpts[j][col];
      for( j = r; j <= pv; j++ )
        for( col = 0; col <= nu; col++ )
          bpts[j][col] = Pw[b-pv+j][col];

      a = b;
      b++;
      ua = ub;
    }
    else
      for( i = 0; i <= ph; i++ )
        Vh[kind+i] = ub;
  }                  
    
  *nhv = mh - ph - 1;

  // just copy knot vector for U-direction
  for( i = 0; i < nu + pu + 2; i++ )
    Uh[i] = U[i];

  *nhu = nu;  
}          
// template instantiation
template void ElevateSurfaceDegreeV( 
      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 ElevateSurfaceDegreeV( 
      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 ElevateSurfaceDegreeV( 
      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 ElevateSurfaceDegreeV( 
      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 ElevateSurfaceDegreeV( 
    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 ElevateSurfaceDegreeV( 
    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 ElevateSurfaceDegreeV( 
    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 ElevateSurfaceDegreeV( 
    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 );



}


        

---