---

gunu_basics.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_basics.h"


namespace gunu {

using gul::point1;
using gul::point2;
using gul::hpoint;
using gul::hpoint1;
using gul::hpoint2;
using gul::set;
using gul::euclid;

/*************************************************************************

  BASICAL NURBS FUNCTIONS

**************************************************************************/

/*-------------------------------------------------------------------------
  Search for the knotspan containing 'u' in a knot vector 'U' (clamped,degree
  'p','n'+1 knots (see "The NURBS book"))
-------------------------------------------------------------------------*/
template< class T >
int FindSpan( const T u, const int n, const int p, 
              const Ptr< T >& U )
{
  int low,mid,high;
  
/* --- search knotspan (from right):   */

  if( u == U[n+1] )
  {
    if( u == U[n+p+1] ) // clamped
      return(n);
                      
    return( n+1 );   // unclamped, knot span = n+1
  }
     
  low = p;
  high = n+1;
  mid = (low+high) / 2 ;  

  while( (u < U[mid]) || (u >= U[mid+1]) ) 
  {
    if( u < U[mid] )
      high = mid;
    else
      low = mid;
      
    mid = (low+high) / 2;
  }
  
  return(mid);     
}           
// template instantiation
template int FindSpan( const float u, const int n, const int p, 
                       const Ptr< float >& U );
template int FindSpan( const double u, const int n, const int p, 
                       const Ptr< double >& U );

/*-------------------------------------------------------------------------
  Finds the knotspan containing 'u' in a knot vector 'U' (clamped,degree
  'p','n'+1 knots (see "The NURBS book")), and counts the multiplicity
  of 'u'
-------------------------------------------------------------------------*/
template< class T >
int FindSpanMultip( const T u, const int n,  const int p, 
                    const Ptr< T >& U, int *s )
{
  int l,low,mid,high;
  
/* --- Knotenspanne suchen (von rechts):   */

  if( u == U[n+1] )
  {
    if( u == U[n+p+1] ) // clamped
    {
      *s = p + 1;
      return(n);
    }
                      
    l = n;           // if unclamped
    while( l >= 0 )  //   count multiplicity
    {
      if( U[l] != u ) break;
      l--;
    }
    *s = n+1 - l;
    return( n+1 );   // knot span = n+1
  }
     
  low = p;
  high = n+1;
  mid = (low+high) / 2 ;  

  while( (u < U[mid]) || (u >= U[mid+1]) ) 
  {
    if( u < U[mid] )
      high = mid;
    else
      low = mid;
      
    mid = (low+high) / 2;
  }
  
  l = mid;
  while( l >= 0 )
  {
    if( U[l] != u ) break;
    l--;
  }
  *s = mid - l;

  return(mid);     
}              
// template instantiation
template int FindSpanMultip( const float u, const int n,  const int p, 
                             const Ptr< float >& U, int *s );
template int FindSpanMultip( const double u, const int n,  const int p, 
                             const Ptr< double >& U, int *s );


/* --------------------------------------------------------------------
  Checks whether a knot vector is valid, clamped and normalized
---------------------------------------------------------------------*/
template< class T >
GULAPI bool ValidateKnotVec( const int n, const int p, const Ptr<T>& knt,
                             bool& is_clamped, bool& is_normalized )
{
  int m, multi, i;
  bool clamped0=false, clamped1=false;
  T u;

  m =  knt.nElems() - 1;
  if( m != n+p+1 ) return false;

  u = knt[0];
  multi = 1;
  for( i = 1; i <= m; i++ )
  {
    if( knt[i] < knt[i-1] ) 
      return false;    // knots must be increasing

    if( knt[i] == u )
      multi++;
    else
    {
      u = knt[i];
      multi = 1;
    }

    if( multi > p )
    {
      if( multi == p+1 )
      {
        if( i == p )
        {
          clamped0 = true;
          continue;
        }
        if( i == m )
        {
          clamped1 = true;
          continue;
        }
      }
      return false;
    }
  }

  is_clamped = clamped0 && clamped1; 
  is_normalized = ((knt[p] == (T)0.0) && (knt[n+1] == (T)1.0));

  return true;
}
// template instantiation
template GULAPI bool ValidateKnotVec( 
                          const int n, const int p, const Ptr<float>& knt,
                          bool& is_clamped, bool& is_normalized );
template GULAPI bool ValidateKnotVec( 
                          const int n, const int p, const Ptr<double>& knt,
                          bool& is_clamped, bool& is_normalized );

/*------------------------------------------------------------------
  Calculates a single basis function
------------------------------------------------------------------*/
template< class T >
GULAPI T CalcBasisFunction( int p, int i, T u, int n, const Ptr<T>& U ) 
{
  Ptr<T> N;
  int j,k;
  T saved,Uleft,Uright,temp;

  N.reserve_pool(p+1);

  if( ((i == 0) && (u == U[0])) ||
      ((i == n) && (u == U[n+1])) )  // since clamped knot vector
    return (T)1;

  if( (u < U[i]) || (u >= U[i+p+1]) )
    return 0;

  // zeroth-degree basis functions
  for( j = 0; j <= p; j++ )
  {
    if( (u >= U[i+j]) && (u < U[i+j+1]) )
      N[j] = (T)1; 
    else
      N[j] = (T)0;
  }

  // compute triangular table
  for( k = 1; k <= p; k++ )
  {
    if( N[0] == (T)0 )
      saved = (T)0;
    else
      saved = ((u-U[i])*N[0])/(U[i+k]-U[i]);

    for( j = 0; j < p-k+1; j++ )
    {
      Uleft = U[i+j+1];
      Uright = U[i+j+k+1];

      if( N[j+1] == (T)0 )
      {
        N[j] = saved;
        saved = (T)0;
      }
      else
      {
        temp = N[j+1]/(Uright-Uleft);
        N[j] = saved + (Uright-u)*temp;
        saved = (u-Uleft)*temp;
      }
    }
  }

  return N[0];
}
// template instantiation
template GULAPI float CalcBasisFunction( 
                 int p, int i, float u, int n, const Ptr<float>& U );
template GULAPI double CalcBasisFunction( 
                 int p, int i, double u, int n, const Ptr<double>& U );


/*------------------------------------------------------------------
  Calculates the non vanishing NURBS basis functions (see "The NURBS book")
-------------------------------------------------------------------*/
template< class T >
void CalcBasisFunctions( const T u, const int i, const int p, const Ptr< T >& U,
                         Ptr< T >& N )
{
  Ptr< T > left,right;
  T saved,temp;
  int j,r;
  
  left.reserve_place( reserve_stack(T,p+1), p+1 );
  right.reserve_place( reserve_stack(T,p+1), p+1 );
  
  N[0] = 1.0;
  
  for( j = 1; j <= p; j++ )
  {
    left[j] = u - U[i+1-j];
    right[j] = U[i+j] - u;
    saved = 0.0;
    
    for( r = 0; r < j; r++ )
    {
      temp = N[r] / (right[r+1] + left[j-r]);
      N[r] = saved + (right[r+1] * temp);
      saved = left[j-r] * temp;
    }
    
    N[j] = saved;
  }  
}
// template instantiation
template void CalcBasisFunctions( 
       const float u, const int i, const int p, const Ptr<float>& U,
       Ptr<float>& N );
template void CalcBasisFunctions( 
       const double u, const int i, const int p, const Ptr<double>& U,
       Ptr<double>& N );


/*--------------------------------------------------------------------------
  Constructs a uniform clamped knot vector ( u_{i+1}-u_{i}) is equal
  for all inner nodes
--------------------------------------------------------------------------*/ 
template< class T >
GULAPI void UniformKnotVector( const int n, const int p, Ptr< T >& U )
{
  int segs,i;
  T seglen;
  
  for( i = 0; i <= p; i++ )
  {
    U[i] = 0.;
    U[n+1+i] = 1.;
  }

  segs = n-p+1;
  seglen = (T)1. / ((T)segs);
  for( i = p+1; i <= n; i++ )
  {
    // U[i] = U[i-1] + seglen;
    U[i] = U[p] + seglen*(T)(i-p);
  }
}
// template instantiation
template GULAPI void UniformKnotVector( const int n, const int p, Ptr< float >& U );
template GULAPI void UniformKnotVector( const int n, const int p, Ptr< double >& U );

/*------------------------------------------------------------------------
  Calculates a point on a NURBS curve (see "The NURBS Book")
----------------------------------------------------------------------- */
template< class T, class HP, class EP >
GULAPI void CurvePoint( 
                 const T u, const int n, const int p, const Ptr< T >& U,
                 const Ptr< HP >& Pw, EP *C )
{
  Ptr< T > N;
  HP v1,h;
  int span,i;

  N.reserve_place( reserve_stack(T,p+1), p+1 );
  
  span = FindSpan( u, n, p, U );
  CalcBasisFunctions<T>( u, span, p, U, N );

  set( h, (T)0.0 );
  
  for( i = 0; i <= p; i++ )
  {
    v1 = N[i] * Pw[span-p+i];
    h = h + v1;
  } 
  *C = euclid( h );
}
// template instantiation
template GULAPI void CurvePoint( 
          const float u, const int n, const int p, const Ptr< float >& U,
          const Ptr< hpoint<float> >& Pw, point<float> *C );
template GULAPI void CurvePoint( 
          const double u, const int n, const int p, const Ptr< double >& U,
          const Ptr< hpoint<double> >& Pw, point<double> *C );         
template GULAPI void CurvePoint( 
          const float u, const int n, const int p, const Ptr< float >& U,
          const Ptr< hpoint2<float> >& Pw, point2<float> *C );
template GULAPI void CurvePoint( 
          const double u, const int n, const int p, const Ptr< double >& U,
          const Ptr< hpoint2<double> >& Pw, point2<double> *C );

template GULAPI void CurvePoint( 
          const float u, const int n, const int p, const Ptr< float >& U,
          const Ptr< point<float> >& Pw, point<float> *C );
template GULAPI void CurvePoint( 
          const double u, const int n, const int p, const Ptr< double >& U,
          const Ptr< point<double> >& Pw, point<double> *C );         
template GULAPI void CurvePoint( 
          const float u, const int n, const int p, const Ptr< float >& U,
          const Ptr< point2<float> >& Pw, point2<float> *C );
template GULAPI void CurvePoint( 
          const double u, const int n, const int p, const Ptr< double >& U,
          const Ptr< point2<double> >& Pw, point2<double> *C );



/*------------------------------------------------------------------------
  Calculates a point on a NURBS surface (see "The NURBS Book")
----------------------------------------------------------------------- */
template< class T, class HP, class EP > 
GULAPI
void SurfacePoint( 
                const T u, const T v, 
                const int nu, const int pu, const Ptr< T >& U,
                const int nv, const int pv, const Ptr< T >& V,                    
                const Ptr< Ptr< HP > >& Pw, EP *retS )
{                          
  Ptr< T > Nu,Nv;
  int uspan,vspan,uind,vind,l,k;
  HP S,temp,v1;

  Nu.reserve_place( reserve_stack(T,pu+1), pu+1 );
  Nv.reserve_place( reserve_stack(T,pv+1), pv+1 );
 
  uspan = FindSpan( u, nu, pu, U );
  CalcBasisFunctions( u, uspan, pu, U, Nu );
  
  vspan = FindSpan( v, nv, pv, V );
  CalcBasisFunctions( v, vspan, pv, V, Nv );

  uind = uspan-pu;
  set( S, (T)0.0 );
  for( l = 0; l <= pv; l++ )
  {
    set( temp, (T)0.0 );
    vind = vspan - pv + l;
    for( k = 0; k <= pu; k++ )
    {
      v1 = Nu[k] * Pw[vind][uind+k]; /* meine Kontrollpunktmatrix ist falschrum,
                                        muss daher transponiert werden */
      temp = temp + v1;
    }
    v1 = Nv[l] * temp;
    S = S + v1;
  }   
  *retS = euclid( S );
}          
// template instantiation
template GULAPI void SurfacePoint( 
                const float u, const float v, 
                const int nu, const int pu, const Ptr< float >& U,
                const int nv, const int pv, const Ptr< float >& V,                    
                const Ptr< Ptr< hpoint<float> > >& Pw, point<float> *retS );
template GULAPI void SurfacePoint( 
                const double u, const double v, 
                const int nu, const int pu, const Ptr< double >& U,
                const int nv, const int pv, const Ptr< double >& V,                    
                const Ptr< Ptr< hpoint<double> > >& Pw, point<double> *retS );
template GULAPI void SurfacePoint( 
                const float u, const float v, 
                const int nu, const int pu, const Ptr< float >& U,
                const int nv, const int pv, const Ptr< float >& V,                    
                const Ptr< Ptr< hpoint1<float> > >& Pw, point1<float> *retS );
template GULAPI void SurfacePoint( 
                const double u, const double v, 
                const int nu, const int pu, const Ptr< double >& U,
                const int nv, const int pv, const Ptr< double >& V,                    
                const Ptr< Ptr< hpoint1<double> > >& Pw, point1<double> *retS );

template GULAPI void SurfacePoint( 
                const float u, const float v, 
                const int nu, const int pu, const Ptr< float >& U,
                const int nv, const int pv, const Ptr< float >& V,                    
                const Ptr< Ptr< point<float> > >& Pw, point<float> *retS );
template GULAPI void SurfacePoint( 
                const double u, const double v, 
                const int nu, const int pu, const Ptr< double >& U,
                const int nv, const int pv, const Ptr< double >& V,                    
                const Ptr< Ptr< point<double> > >& Pw, point<double> *retS );
template GULAPI void SurfacePoint( 
                const float u, const float v, 
                const int nu, const int pu, const Ptr< float >& U,
                const int nv, const int pv, const Ptr< float >& V,                    
                const Ptr< Ptr< point1<float> > >& Pw, point1<float> *retS );
template GULAPI void SurfacePoint( 
                const double u, const double v, 
                const int nu, const int pu, const Ptr< double >& U,
                const int nv, const int pv, const Ptr< double >& V,                    
                const Ptr< Ptr< point1<double> > >& Pw, point1<double> *retS );


/*------------------------------------------------------------------------
  Calculates a point on a Bezier surface (for Bezier surfaces there is
  no need to search the knot span)
----------------------------------------------------------------------- */

template< class T, class HP, class EP > 
GULAPI
void BezierSurfacePoint( 
                const T u, const T v, 
                const int pu, const Ptr< T >& U,
                const int pv, const Ptr< T >& V,                    
                const Ptr< Ptr< HP > >& Pw, EP *retS )
{                          
  Ptr< T > Nu,Nv;
  int uspan,vspan,uind,vind,l,k;
  HP S,temp,v1;

  Nu.reserve_place( reserve_stack(T,pu+1), pu+1 );
  Nv.reserve_place( reserve_stack(T,pv+1), pv+1 );
 
  uspan = pu;
  CalcBasisFunctions( u, uspan, pu, U, Nu );
  
  vspan = pv;
  CalcBasisFunctions( v, vspan, pv, V, Nv );

  uind = uspan-pu;
  set( S, (T)0.0 );
  for( l = 0; l <= pv; l++ )
  {
    set( temp, (T)0.0 );
    vind = vspan - pv + l;
    for( k = 0; k <= pu; k++ )
    {
      v1 = Nu[k] * Pw[vind][uind+k]; /* meine Kontrollpunktmatrix ist falschrum,
                                        muss daher transponiert werden */
      temp = temp + v1;
    }
    v1 = Nv[l] * temp;
    S = S + v1;
  }   
  *retS = euclid( S );
}          
// template instantiation
template GULAPI void BezierSurfacePoint( 
                const float u, const float v, 
                const int pu, const Ptr< float >& U,
                const int pv, const Ptr< float >& V,                    
                const Ptr< Ptr< hpoint<float> > >& Pw, point<float> *retS );
template GULAPI void BezierSurfacePoint( 
                const double u, const double v, 
                const int pu, const Ptr< double >& U,
                const int pv, const Ptr< double >& V,                    
                const Ptr< Ptr< hpoint<double> > >& Pw, point<double> *retS );
template GULAPI void BezierSurfacePoint( 
                const float u, const float v, 
                const int pu, const Ptr< float >& U,
                const int pv, const Ptr< float >& V,                    
                const Ptr< Ptr< hpoint1<float> > >& Pw, point1<float> *retS );
template GULAPI void BezierSurfacePoint( 
                const double u, const double v, 
                const int pu, const Ptr< double >& U,
                const int pv, const Ptr< double >& V,                    
                const Ptr< Ptr< hpoint1<double> > >& Pw, point1<double> *retS );

template GULAPI void BezierSurfacePoint( 
                const float u, const float v, 
                const int pu, const Ptr< float >& U,
                const int pv, const Ptr< float >& V,                    
                const Ptr< Ptr< point<float> > >& Pw, point<float> *retS );
template GULAPI void BezierSurfacePoint( 
                const double u, const double v, 
                const int pu, const Ptr< double >& U,
                const int pv, const Ptr< double >& V,                    
                const Ptr< Ptr< point<double> > >& Pw, point<double> *retS );
template GULAPI void BezierSurfacePoint( 
                const float u, const float v, 
                const int pu, const Ptr< float >& U,
                const int pv, const Ptr< float >& V,                    
                const Ptr< Ptr< point1<float> > >& Pw, point1<float> *retS );
template GULAPI void BezierSurfacePoint( 
                const double u, const double v, 
                const int pu, const Ptr< double >& U,
                const int pv, const Ptr< double >& V,                    
                const Ptr< Ptr< point1<double> > >& Pw, point1<double> *retS );




}

---