---

gunu_derivatives.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_float.h"
#include "gul_vector.h"
#include "guar_bincoeff.h"
#include "gunu_basics.h"
#include "gunu_derivatives.h"


namespace gunu {

using gul::rtr;
using gul::Ptr;
using gul::point;
using gul::point1;
using gul::point2;
using gul::hpoint;
using gul::hpoint1;
using gul::hpoint2;
using gul::set;
using gul::euclid;
using gul::ortho;
using gul::Min;
using gul::Max;
using gul::cross_product;
using gul::normalize;
using gul::rel_equal;

/*------------------------------------------------------------------------
  Calculate NURBS basis functions and their first 'n' derivatives.
  Only the non-vanishing basis functions are calculated (see "The NURBS Book")
------------------------------------------------------------------------*/  
template< class T >
void CalcDersBasisFuns( 
            const T u, const int i, const int p, const int n,
            const Ptr< T >& U, Ptr< Ptr< T > >& ders )
{
  T saved,temp,d;
  int j,k,r,s1,s2,rk,pk,j1,j2;

  Ptr< Ptr< T > > ndu,a;
  Ptr< T > left,right;
  
  ndu.reserve_place( reserve_stack(Ptr< T >,p+1), p+1 );
  for( j = 0; j < p+1; j++ )
    ndu[j].reserve_place( reserve_stack(T,p+1), p+1 );
 
  a.reserve_place( reserve_stack(Ptr< T >,2), 2 );
  for( j = 0; j < 2; j++ )
    a[j].reserve_place( reserve_stack(T,p+1), p+1 );

  left.reserve_place( reserve_stack(T,p+1), p+1 );
  right.reserve_place( reserve_stack(T,p+1), p+1 );
  
  ndu[0][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++ )
    {
      ndu[j][r] = right[r+1] + left[j-r];
      temp = ndu[r][j-1]/ndu[j][r];
      
      ndu[r][j] = saved + right[r+1]*temp;
      saved = left[j-r]*temp;
    }
    ndu[j][j] = saved;
  }

  
  for( j = 0; j <= p; j++ )
  {
    ders[0][j] = ndu[j][p];        /* Basis Funktionen = 0.Ableitung */
  }

  for( r = 0; r <= p; r++ )
  {
    s1 = 0;    s2 = 1;
    a[0][0] = 1.0;

    for( k = 1; k <= n; k++ )
    {
      d = 0.0;
      rk = r-k;  pk = p-k;
      
      if( r >= k )
      {
        a[s2][0] = a[s1][0] / ndu[pk+1][rk];
        d = a[s2][0] * ndu[rk][pk];
      }  
      if( rk >= -1 )
        j1 = 1;
      else
        j1 = -rk;  
      if( r-1 <= pk )
        j2 = k-1;
      else
        j2 = p-r;  

      for( j = j1; j <= j2; j++ )
      {
        a[s2][j] = (a[s1][j] - a[s1][j-1]) / ndu[pk+1][rk+j];
        d += a[s2][j] * ndu[rk+j][pk];
      }  
      if( r <= pk )
      {
        a[s2][k] = -a[s1][k-1] / ndu[pk+1][r];
        d += a[s2][k] * ndu[r][pk];
      }  
      ders[k][r] = d;
      j = s1; s1 = s2; s2 = j;
    }        
  }
  
  r = p;
  for( k = 1; k <= n; k++ )
  {
    for( j = 0; j <= p; j++ )
      ders[k][j] *= (T)r;
    r *= p-k;
  }    
}
// template instantiation
template void CalcDersBasisFuns( 
          const float u, const int i, const int p, const int n, 
          const Ptr<float>& U, Ptr< Ptr<float> >& ders );
template void CalcDersBasisFuns( 
          const double u, const int i, const int p, const int n,
          const Ptr<double>& U, Ptr< Ptr<double> >& ders );

/*-------------------------------------------------------------------------
  Calculates the first 'd' derivatives of a curve at a point 'u'. The results
  are homogeneous points (like the control points), returned in 'CK[]'
  (see "The NURBS Book")
--------------------------------------------------------------------------- */  
template< class T, class EP >
GULAPI void BSPCurveDerivatives( 
           const T u, const int n, const int p, const gul::Ptr< T >& U,
           const gul::Ptr< EP >& Pw, const int d, gul::Ptr< EP >& CK )
{
  int du,k,j,span;
  EP v1;
  Ptr< Ptr< T > > ders;
    
  du = Min( d, p );

  ders.reserve_place( reserve_stack(Ptr< T >,du+1), du+1 );
  for( j = 0; j <du+1; j++ )
    ders[j].reserve_place( reserve_stack(T,p+1), p+1 );

  for( k = p+1; k <= d; k++ )
    set( CK[k], (T)0.0 );   /* Ableitungen > Grad = 0 */

  span = FindSpan( u, n, p, U );
  CalcDersBasisFuns( u, span, p, du, U, ders );
  
  for( k = 0; k <= du; k++ )
  {
    set( CK[k], (T)0.0 );
    for( j = 0; j <= p; j++ )
    {
      v1 = ders[k][j] * Pw[span-p+j];
      CK[k] = CK[k] + v1;
    }
  }
}  
// template instantiation
template GULAPI void BSPCurveDerivatives( 
     const float u, const int n, const int p, const Ptr< float >& U,
     const Ptr< point<float> >& Pw, const int d, Ptr< point<float> >& CK );
template GULAPI void BSPCurveDerivatives( 
      const double u, const int n, const int p, const Ptr< double >& U,
      const Ptr< point<double> >& Pw, const int d, Ptr< point<double> >& CK );
template GULAPI void BSPCurveDerivatives( 
       const float u, const int n, const int p, const Ptr< float >& U,
       const Ptr< point2<float> >& Pw, const int d, Ptr< point2<float> >& CK );
template GULAPI void BSPCurveDerivatives( 
      const double u, const int n, const int p, const Ptr< double >& U,
      const Ptr< point2<double> >& Pw, const int d, Ptr< point2<double> >& CK );

/*---------------------------------------------------------------------------
  Calculates mixed partial derivatives of a NURBS surface. Let 'k' and 'l' be
  the number of derivatives in 'u' and 'v' direction, then all mixed ders with
  'k+l <= d' are calculated. They are returned in 'SKL[k][l]'.  
  (see "The NURBS Book")
------------------------------------------------------------------------- */  
template< class T, class EP >
void BSPSurfaceDerivatives(
                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< EP > >& Pw,
                const int d, Ptr< Ptr< EP > >& SKL )
{                
  int du,dv,k,l,uspan,vspan,s,r,dd;
  EP v1;
  Ptr< EP > temp;
  Ptr< Ptr< T > > Nu,Nv;

  du = Min( d, pu );
  dv = Min( d, pv );

  Nu.reserve_place( reserve_stack(Ptr< T >,du+1), du+1 );
  for( k = 0; k < pu+1; k++ )
    Nu[k].reserve_place( reserve_stack(T,pu+1), pu+1 );

  Nv.reserve_place( reserve_stack(Ptr< T >,dv+1), dv+1 );
  for( k = 0; k < pv+1; k++ )
    Nv[k].reserve_place( reserve_stack(T,pv+1), pv+1 );

  temp.reserve_place( reserve_stack(EP,pv+1), pv+1 );

  for( k = pu+1; k <= d; k++ )
    for( l = 0; l <= d - k; l++ )
      set( SKL[k][l], (T)0.0 );
      
  for( l = pv+1; l <= d; l++ )
    for( k = 0; k <= d - l; k++ )
      set( SKL[k][l], (T)0.0 );
      
  uspan = FindSpan( u, nu, pu, U );
  CalcDersBasisFuns( u, uspan, pu, du, U, Nu );
 
  vspan = FindSpan( v, nv, pv, V );
  CalcDersBasisFuns( v, vspan, pv, dv, V, Nv );
  
  for( k = 0; k <= du; k++ )
  {
    for( s = 0; s <= pv; s++ )
    {
      set( temp[s], (T)0.0 );
      for( r = 0; r <= pu; r++ )
      {
        v1 = Nu[k][r] * Pw[vspan-pv+s][uspan-pu+r];
        temp[s] = temp[s] + v1;
      }  
    }
    dd = Min( d - k, dv );
    for( l = 0; l <= dd; l++ )
    {
      set( SKL[k][l], (T)0.0 );
      for( s = 0; s <= pv; s++ )
      {
        v1 = Nv[l][s] * temp[s];
        SKL[k][l] = SKL[k][l] + v1;
      }
    }
  }
}  
// template instantiation
template void BSPSurfaceDerivatives(
                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,
                const int d, Ptr< Ptr< point<float> > >& SKL );           
template void BSPSurfaceDerivatives(
                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,
                const int d, Ptr< Ptr< point<double> > >& SKL );

template void BSPSurfaceDerivatives(
                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,
                const int d, Ptr< Ptr< hpoint<float> > >& SKL );           
template void BSPSurfaceDerivatives(
                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,
                const int d, Ptr< Ptr< hpoint<double> > >& SKL );

/*-------------------------------------------------------------------------
  Calculates the first 'd' partial derivatives of the projection of a NURBS
  curve into euclidian space (3-dimensions), so the calculated ders are points
  in 3-dimensional space   
  (see "The NURBS Book")
------------------------------------------------------------------------- */  
template< class T, class HP, class EP >
GULAPI void CurveDerivatives(
                    const T u, const int n, const int p,
                    const Ptr< T >& U, const Ptr< HP >& Pw, const int d,
                    Ptr< EP >& CK )             
{
  Ptr< HP > CKh;
  EP v,v1;
  T w0;
  int k,i;

  CKh.reserve_place( reserve_stack(HP,d+1), d+1 );
     
  BSPCurveDerivatives<T,HP>( u, n, p, U, Pw, d, CKh );

  w0 = CKh[0].w;
  
  for( k = 0; k <= d; k++ )
  {
    v = ortho( CKh[k] );
    
    for( i = 1; i <= k; i++ )
    {
      v1 = (rtr<T>::BinCoeff(k,i) * CKh[i].w) * CK[k-i];
      v = v - v1;
    }
    CK[k] = ((T)1.0 / w0) * v;       
  }
}
// template instantiation
template GULAPI void CurveDerivatives(
          const float u, const int n, const int p,
          const Ptr< float >& U, const Ptr< hpoint<float> >& Pw, const int d,
          Ptr< point<float> >& CK );
template GULAPI void CurveDerivatives(
          const double u, const int n, const int p,
          const Ptr< double >& U, const Ptr< hpoint<double> >& Pw, const int d,
          Ptr< point<double> >& CK );

template GULAPI void CurveDerivatives(
         const float u, const int n, const int p,
         const Ptr< float >& U, const Ptr< hpoint2<float> >& Pw, const int d,
         Ptr< point2<float> >& CK );
template GULAPI void CurveDerivatives(
         const double u, const int n, const int p,
         const Ptr< double >& U, const Ptr< hpoint2<double> >& Pw, const int d,
         Ptr< point2<double> >& CK );

// VC currently (VC7) has problems with the ordering of function templates,
// so i have to declare the following as specializations
#ifdef _MSC_VER
template<> GULAPI void CurveDerivatives(
     const float u, const int n, const int p,
     const gul::Ptr<float>& U, const gul::Ptr<gul::point2<float> >& Pw, 
     const int d, gul::Ptr<gul::point2<float> >& CK )  
{
  BSPCurveDerivatives(u,n,p,U,Pw,d,CK);
}
template<> GULAPI void CurveDerivatives(
     const double u, const int n, const int p,
     const gul::Ptr<double>& U, const gul::Ptr<gul::point2<double> >& Pw, 
     const int d, gul::Ptr<gul::point2<double> >& CK )  
{
  BSPCurveDerivatives(u,n,p,U,Pw,d,CK);
}
template<> GULAPI void CurveDerivatives(
     const float u, const int n, const int p,
     const gul::Ptr<float>& U, const gul::Ptr<gul::point<float> >& Pw, 
     const int d, gul::Ptr<gul::point<float> >& CK )  
{
  BSPCurveDerivatives(u,n,p,U,Pw,d,CK);
}
template<> GULAPI void CurveDerivatives(
     const double u, const int n, const int p,
     const gul::Ptr<double>& U, const gul::Ptr<gul::point<double> >& Pw, 
     const int d, gul::Ptr<gul::point<double> >& CK )  
{
  BSPCurveDerivatives(u,n,p,U,Pw,d,CK);
}
#endif

/*-------------------------------------------------------------------------
  Calculates the first 'd' mixed partial derivatives of the projection of a
  NURBS surface into euclidian space (3-dimensions), so the calculated ders
  are points in 3-dimensional space. They are returned in 'SKL[k][l]'.  
  (see "The NURBS Book")
------------------------------------------------------------------------- */
template< class T, class HP, class EP >          
void SurfaceDerivatives(
                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,
                const int d, Ptr< Ptr < EP > >& SKL )
{
  int i,j,k,l;
  EP v1,v2,vh;
  T w00;
  Ptr< Ptr < HP > > SKLh;
  
  SKLh.reserve_place( reserve_stack(Ptr < HP >,d+1), d+1 );
  for( i = 0; i < d+1; i++ )
    SKLh[i].reserve_place( reserve_stack(HP,d+1), d+1 );
    
  BSPSurfaceDerivatives<T,HP>( u, v, nu, pu, U, nv, pv, V, Pw, d, SKLh );

  w00 = SKLh[0][0].w;
  
  for( k = 0; k <= d; k++ )
  {
    for( l = 0; l <= d-k; l++ )
    {
      v1 = ortho( SKLh[k][l] );

      for( j = 1; j <= l; j++ )
      {
        vh = (rtr<T>::BinCoeff(l,j) * SKLh[0][j].w) * SKL[k][l-j]; 
        v1 = v1 - vh;
      }  
      for( i = 1; i <= k; i++ )
      {
        vh = (rtr<T>::BinCoeff(k,i) * SKLh[i][0].w) * SKL[k-i][l];
        v1 = v1 - vh;

        set( v2, (T)0.0 );
        for( j = 1; j <= l; j++ )
        {
          vh = (rtr<T>::BinCoeff(l,j) * SKLh[i][j].w) * SKL[k-i][l-j];
          v2 = v2 + vh;
        }
        vh =  rtr<T>::BinCoeff(k,i) * v2;
        v1 = v1 - vh;
      }
      SKL[k][l] = ((T)1.0 / w00) * v1;
    }    
  }
}
// template instantiation
template void SurfaceDerivatives(
                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,
                const int d, Ptr< Ptr < point<float> > >& SKL );
template void SurfaceDerivatives(
                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,
                const int d, Ptr< Ptr < point<double> > >& SKL );

template void SurfaceDerivatives(
                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,
                const int d, Ptr< Ptr < point1<float> > >& SKL );
template void SurfaceDerivatives(
                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,
                const int d, Ptr< Ptr < point1<double> > >& SKL );

/*---------------------------------------------------------------------
  Calculates the normal vector of a NURBS surface at u=0,v=0
---------------------------------------------------------------------*/

template< class T, class HP >
void CalcNormal_U0V0(
          const int nu, const int pu,
          const int nv, const int pv,
          const Ptr<Ptr<HP> >& Pw, point<T> *Normal )
{
  int i,j,k,m;
  point<T> v1,v2,a,b;

  set(a,(T)0);
  set(b,(T)0);

  m = Max(pu,pv);
  
  for( k=0; k<m; k++ )
  {
    i=0;
    j=k+1;
    
    while( i <= k )
    {
      if( (i <= pv) && (j <= pu) )
      {
        v1 = euclid( Pw[i][j] );
        v2 = euclid( Pw[i][0] );
        a = v1 - v2;
        if( !rel_equal(v1,v2,rtr<T>::zero_tol()) )
        {
          k = m;     // Vektor in Tangetialebene gefunden, fertig
          break;
        }
      }
      i++;
      j--;  
    }
  }

  for( k=0; k<m; k++ )
  {
    i=0;
    j=k+1;
    
    while( i <= k )
    {
      if( (j <= pv) && (i <= pu) )
      {
        v1 = euclid( Pw[j][i] );
        v2 = euclid( Pw[0][i] );
        b = v1 - v2;
        if( !rel_equal(v1,v2,rtr<T>::zero_tol()) )
        {
          k = m;     // Vektor in Tangetialebene gefunden, fertig
          break;
        }
      }
      i++;
      j--;  
    }
  }

  v1 = cross_product( a, b );
  normalize( Normal, v1 );
}
// template instantiation
template void CalcNormal_U0V0( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< hpoint<float> > >& Pw, point<float> *Normal );
template void CalcNormal_U0V0( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< point<float> > >& Pw, point<float> *Normal );

template void CalcNormal_U0V0( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< hpoint<double> > >& Pw, point<double> *Normal );
template void CalcNormal_U0V0( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< point<double> > >& Pw, point<double> *Normal );

/*---------------------------------------------------------------------
  Calculates the normal vector of a NURBS surface at u=1,v=0
---------------------------------------------------------------------*/
template< class T, class HP >
void CalcNormal_U1V0( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< HP > >& Pw, point<T> *Normal  )
{
  int i,j,k,m;
  point<T> v1,v2,a,b;

  set(a,(T)0);
  set(b,(T)0);

  m = Max(pu,pv);
  
  for( k=0; k<m; k++ )
  {
    i=0;
    j=k+1;
    
    while( i <= k )
    {
      if( (i <= pv) && (j <= pu) )
      {
        v1 = euclid( Pw[i][nu] );
        v2 = euclid( Pw[i][nu-j] );
        a =  v1 - v2;
        if( !rel_equal(v1,v2,rtr<T>::zero_tol()) )
        {
          k = m;     /* Vektor in Tangetialebene gefunden, fertig */
          break;
        }
      }
      i++;
      j--;  
    }
  }

  for( k=0; k<m; k++ )
  {
    i=0;
    j=k+1;
    
    while( i <= k )
    {
      if( (j <= pv) && (i <= pu) )
      {
        v1 = euclid( Pw[j][nu-i] );
        v2 = euclid( Pw[0][nu-i] );
        b = v1 - v2;
        if( !rel_equal(v1,v2,rtr<T>::zero_tol()) )
        {
          k = m;     /* Vektor in Tangetialebene gefunden, fertig */
          break;
        }
      }
      i++;
      j--;  
    }
  }

  v1 = cross_product( a, b );
  normalize( Normal, v1 );
}
// template instantiation
template void CalcNormal_U1V0( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< hpoint<float> > >& Pw, point<float> *Normal );
template void CalcNormal_U1V0( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< point<float> > >& Pw, point<float> *Normal );

template void CalcNormal_U1V0( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< hpoint<double> > >& Pw, point<double> *Normal );
template void CalcNormal_U1V0( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< point<double> > >& Pw, point<double> *Normal );


/*---------------------------------------------------------------------
  Calculates the normal vector of a NURBS surface at u=0,v=1
---------------------------------------------------------------------*/
template< class T, class HP >
void CalcNormal_U0V1( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< HP > >& Pw, point<T> *Normal  )
{
  int i,j,k,m;
  point<T> v1,v2,a,b;

  set(a,(T)0);
  set(b,(T)0);
  
  m = Max(pu,pv);
  
  for( k=0; k<m; k++ )
  {
    i=0;
    j=k+1;
    
    while( i <= k )
    {
      if( (i <= pv) && (j <= pu) )
      {
        v1 = euclid( Pw[nv-i][j] );
        v2 = euclid( Pw[nv-i][0] );
        a = v1 - v2;
        if( !rel_equal(v1,v2,rtr<T>::zero_tol()) )
        {
          k = m;     /* Vektor in Tangetialebene gefunden, fertig */
          break;
        }       
      }
      i++;
      j--;  
    }
  }

  for( k=0; k<m; k++ )
  {
    i=0;
    j=k+1;
    
    while( i <= k )
    {
      if( (j <= pv) && (i <= pu) )
      {
        v1 = euclid( Pw[nv][i] );
        v2 = euclid( Pw[nv-j][i] );
        b =  v1 - v2;
        if( !rel_equal(v1,v2,rtr<T>::zero_tol()) )
        {
          k = m;     /* Vektor in Tangetialebene gefunden, fertig */
          break;
        }
      }
      i++;
      j--;  
    }
  }

  v1 = cross_product( a, b );
  normalize( Normal, v1 );
}
// template instantiation
template void CalcNormal_U0V1( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< hpoint<float> > >& Pw, point<float> *Normal );
template void CalcNormal_U0V1( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< point<float> > >& Pw, point<float> *Normal );

template void CalcNormal_U0V1( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< hpoint<double> > >& Pw, point<double> *Normal );
template void CalcNormal_U0V1( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< point<double> > >& Pw, point<double> *Normal );


/*---------------------------------------------------------------------
  Calculates the normal vector of a NURBS surface at u=1,v=1
---------------------------------------------------------------------*/
template< class T, class HP >
void CalcNormal_U1V1( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< HP > >& Pw, point<T> *Normal  )
{
  int i,j,k,m;
  point<T> v1,v2,a,b;

  set(a,(T)0);
  set(b,(T)0);
  
  m = Max(pu,pv);
  
  for( k=0; k<m; k++ )
  {
    i=0;
    j=k+1;
    
    while( i <= k )
    {
      if( (i <= pv) && (j <= pu) )
      {
        v1 = euclid( Pw[nv-i][nu] );
        v2 = euclid( Pw[nv-i][nu-j] );
        a = v1 - v2;
        if( !rel_equal(v1,v2,rtr<T>::zero_tol()) )
        {
          k = m;     /* Vektor in Tangetialebene gefunden, fertig */
          break;
        }
      }
      i++;
      j--;  
    }
  }

  for( k=0; k<m; k++ )
  {
    i=0;
    j=k+1;
    
    while( i <= k )
    {
      if( (j <= pv) && (i <= pu) )
      {
        v1 = euclid( Pw[nv][nu-i] );
        v2 = euclid( Pw[nv-j][nu-i] );
        b = v1 - v2;
        if( !rel_equal(v1,v2,rtr<T>::zero_tol()) )
        {
          k = m;     /* Vektor in Tangetialebene gefunden, fertig */
          break;
        }
      }
      i++;
      j--;  
    }
  }

  v1 = cross_product( a, b );
  normalize( Normal, v1 );
}
// template instantiation
template void CalcNormal_U1V1( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< hpoint<float> > >& Pw, point<float> *Normal );
template void CalcNormal_U1V1( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< point<float> > >& Pw, point<float> *Normal );

template void CalcNormal_U1V1( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< hpoint<double> > >& Pw, point<double> *Normal );
template void CalcNormal_U1V1( 
           const int nu, const int pu, 
           const int nv, const int pv,
           const Ptr< Ptr< point<double> > >& Pw, point<double> *Normal );

}

---