---

nurbsS.cpp

/*=============================================================================
        File: nurbsS.cpp
     Purpose:       
    Revision: $Id: nurbsS.cpp,v 1.4 2003/01/13 19:42:01 philosophil Exp $
  Created by: Philippe Lavoie          (3 Oct, 1996)
 Modified by: 

 Copyright notice:
          Copyright (C) 1996-1997 Philippe Lavoie
 
          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., 675 Mass Ave, Cambridge, MA 02139, USA.
=============================================================================*/

#ifndef PLIB_NURBS_NURBSS_SOURCE
#define PLIB_NURBS_NURBSS_SOURCE


#include <nurbs.h>
#include <string.h>
#include <matrixRT.h>
#include <math.h>
#include <nurbsS.h>
#include "integrate.h"

#ifdef USING_VCC
#include <malloc.h>
#endif

namespace PLib {

template <class T, int N>
NurbsSurface<T,N>::NurbsSurface(): 
  ParaSurface<T,N>(), U(1),V(1),P(1,1),degU(0),degV(0)
{
}

template <class T, int N>
NurbsSurface<T,N>::NurbsSurface(const NurbsSurface<T,N>& s): 
  ParaSurface<T,N>(), U(s.U), V(s.V), P(s.P), degU(s.degU), degV(s.degV)
{
}


template <class T, int N>
NurbsSurface<T,N>::NurbsSurface(int DegU, int DegV, const Vector<T>& Uk, const Vector<T>& Vk, const Matrix< HPoint_nD<T,N> >& Cp) : 
  ParaSurface<T,N>(), U(Uk),V(Vk),P(Cp),degU(DegU),degV(DegV)
{
  int bad = 0 ;

  if(U.n() != P.rows()+degU+1){
#ifdef USE_EXCEPTION
    throw NurbsSizeError(P.rows(),U.n(),degU) ;
#else
    Error err("NurbsSurface<T,N> constructor") ;
    err << "The U knot vector and the number of rows of the control points are incompatible\n" ;
    err << "P.rows() = " << P.rows() << ", U.n() = " << U.n() << endl ;
    bad = 1 ;
    err.fatal() ;
#endif
  }
  if(V.n() != P.cols()+degV+1){
#ifdef USE_EXCEPTION
    throw NurbsSizeError(P.cols(),V.n(),degV);
#else
    Error err("NurbsSurface<T,N> constructor") ;    
    err << "The V knot vector and the number of columns of the control points are incompatible\n" ;
    err << "P.cols() = " << P.cols() << ", V.n() = " << V.n() << endl ;
    bad = 1 ;
    err.fatal() ;
#endif
  }
  if(bad)
    exit(-1) ;
}

template <class T, int N>
NurbsSurface<T,N>::NurbsSurface(int DegU, int DegV, Vector<T>& Uk, Vector<T>& Vk, Matrix< Point_nD<T,N> >& Cp, Matrix<T>& W) : 
  ParaSurface<T,N>(), U(Uk),V(Vk),P(Cp.rows(),Cp.cols()),degU(DegU),degV(DegV)
{
  int bad = 0 ;

  if(U.n() != Cp.rows()+degU+1){
#ifdef USE_EXCEPTION
    throw NurbsSizeError(P.rows(),U.n(),degU);
#else
    Error err("NurbsSurface<T,N> constructor") ;
    err << "The U knot vector and the number of rows of the control points are incompatible\n" ;
    err << "P.rows() = " << P.rows() << ", U.n() = " << U.n() << endl ;
    bad = 1 ;
    err.fatal() ;
#endif
  }
  if(V.n() != Cp.cols()+degV+1){
#ifdef USE_EXCEPTION
    throw NurbsSizeError(P.cols(),V.n(),degV);
#else
    Error err("NurbsSurface<T,N> constructor") ;
    err << "The V knot vector and the number of columns of the control points are incompatible\n" ;
    err << "P.cols() = " << P.cols() << ", V.n() = " << V.n() << endl ;
    bad = 1 ;
    err.fatal() ;
#endif
  }
  if(W.rows() != Cp.rows()){
#ifdef USE_EXCEPTION
    throw NurbsInputError(W.rows(),Cp.rows()) ;
#else
    Error err("NurbsSurface<T,N> constructor") ;
    err << "The dimension of the weights are incompatible with the dimension of the control points\n" ;
    err << "W( " << W.rows() << ", " << W.cols() << ") and P( " << P.rows() << ", " << P.cols() << ") \n" ;
    bad = 1 ;
    err.fatal() ;
#endif
  }
  if(W.cols() != Cp.cols()){
#ifdef USE_EXCEPTION
    throw NurbsInputError(W.cols(),Cp.cols()) ;
#else
    Error err("NurbsSurface<T,N> constructor") ;
    err << "The dimension of the weights are incompatible with the dimension of the control poitns\n" ;
    err << "W( " << W.rows() << ", " << W.cols() << ") and P( " << P.rows() << ", " << P.cols() << ") \n" ;
    bad = 1 ;
    err.fatal() ;
#endif
  }

  for(int i=0;i<Cp.rows();i++)
    if(bad)
      break ;
    else{
      for(int j=0;j<Cp.cols();j++)
        if(W(i,j)==0.0){
#ifdef USE_EXCEPTION
          throw NurbsInputError();
#else
          Error err("NurbsSurface<T,N> constructor") ;
          err << "Error initializing a NurbsSurface.\n" ;
          err << "The weight W(" << i << ", " << j << ") is equal to 0.\n" ;
          bad = 1 ;
          err.fatal() ;
          break ;
#endif
        }
        else {
          P(i,j) = Cp(i,j) ;
          P(i,j) *= W(i,j) ;
        }
    }

  if(bad)
    exit(-1) ;
}

template <class T, int N>
NurbsSurface<T,N>& NurbsSurface<T,N>::operator=(const NurbsSurface<T,N>& nS){
  P = nS.P ;
  U = nS.U ;
  V = nS.V ;
  degU = nS.degU ;
  degV = nS.degV ;
  return *this ;
}

template <class T, int N>
struct AreaData {
  T v;
  T eps;
  T knotUi;
  T knotUii;
  T knotVj;
  T knotVjj;
  const NurbsSurface<T,N>& s ; 
  const BasicArray<T>  w ; 
  AreaData(const NurbsSurface<T,N>& surf,T e,
           const BasicArray<T>& ww): 
    s(surf),v(T(0)),eps(e),w(ww),
    knotUi(T(0)), knotUii(T(1)) {}
};

template <class T, int N> 
struct OpAreaAuxFcn : public ClassPOvoid<T> {
  T operator()(T u, void* data){
    AreaData<T,N>* p = (AreaData<T,N>*)data ;
    return (p->s).areaF(u,p->v) ; 
  }
};

template <class T, int N>
struct OpAreaFcn : public ClassPOvoid<T> {
  T operator()(T v, void* data){
    static Vector<T> w;
    T err;
    OpAreaAuxFcn<T,N>  f;
    AreaData<T,N>* Data = (AreaData<T,N>*)data ;
    Data->v = v ;            
    return intcc2((ClassPOvoid<T>*)&f,Data,
                  Data->knotUi,Data->knotUii,
                  Data->eps,Data->w,err);  
  }
};

template <class T, int N>
T NurbsSurface<T,N>::area(T eps,int n) const {
  T a = T(0.0) ;
  T err ; 
  static Vector<T> bufFcn(0) ;

  if(bufFcn.n() != n){
    bufFcn.resize(n) ;
    intccini(bufFcn) ;
  }

  AreaData<T,N> data(*this,eps,bufFcn) ; 
  OpAreaFcn<T,N> op;
  

  for(int i=degU;i<P.rows();++i){
    if(U[i] >= U[i+1] || U[i]>=T(1.0))
      continue ;
    data.knotUi  = U[i] ; 
    data.knotUii = U[i+1] ; 
    for(int j=degV;j<P.cols();++j){
      if(V[j] >= V[j+1] || V[j]>=T(1.0))
        continue ;
      data.knotVj  = V[j] ; 
      data.knotVjj = V[j+1] ; 
      a += intcc2((ClassPOvoid<T>*)&op,(void*)&data,
                  data.knotVj,data.knotVjj,
                  eps,bufFcn,err) ;
    }
  }
  return a ; 
}

template <class T, int N>
T NurbsSurface<T,N>::areaIn(T us, T ue, T vs, T ve, T eps, int n) const {
  T l = T() ;
  T err ; 
  T a ;
  bool bLastU = false;
  bool bLastV = false;
  static Vector<T> bufFcn ;

  if(bufFcn.n() != n){
    bufFcn.resize(n) ;
    intccini(bufFcn) ;
  }
  
  AreaData<T,N> data(*this,eps,bufFcn) ; 
  OpAreaFcn<T,N> op;
  for(int i=degU;i<P.rows();++i){
    if(U[i] >= U[i+1] || U[i] >= T(1))
      continue;
    if(i<findSpanU(us))
      continue ; 
    if(us>=U[i] && ue<=U[i+findMultU(i)]){
      data.knotUi  = us;
      data.knotUii = ue;
      bLastU = true;
      goto Integrate_I;
    }
    if(us>=U[i]){
      data.knotUi  = us;
      data.knotUii = U[i+findMultU(i)];
      bLastU = false;
      goto Integrate_I;
    }
    if(ue<=U[i+1]){
      data.knotUi  = U[i];
      data.knotUii = ue;
      bLastU = true;
      goto Integrate_I;
    }
    data.knotUi  = U[i] ; 
    data.knotUii = U[i+findMultU(i)] ; 
  Integrate_I:
    for(int j=degV;j<P.cols();++j){
      if(V[j] >= V[j+1] || V[j]>=T(1))
        continue ;
      if(j<findSpanV(vs))
        continue ; 
      if(vs>=V[j] && ve<=V[j+findMultV(j)]){
        data.knotVj  = vs;
        data.knotVjj = ve;
        bLastV = true;
        goto Integrate_II;
      }
      if(vs>=V[j]){
        data.knotVj  = vs;
        data.knotVjj = V[j+findMultV(j)];
        bLastV = false;
        goto Integrate_II;
      }
      if(ve<=V[j+1]){
        data.knotVj  = V[j];
        data.knotVjj = ve;
        bLastV = true;
        goto Integrate_II;
      }
      data.knotVj  = V[j] ; 
      data.knotVjj = V[j+findMultV(j)] ; 
    Integrate_II:
      a += intcc2((ClassPOvoid<T>*)&op,(void*)&data,data.knotVj,data.knotVjj,eps,bufFcn,err) ;
      if (bLastU && bLastV)
        return a;
      if (bLastV)
        break;
    }
  }
  return a ; 
}

// the definitions are in f_nurbs.cpp and d_nurbs.cpp


template <class T, int N>
  T NurbsSurface<T,N>::areaF(T u, T v) const {
  
  Matrix<Point_nD<T,N> > Skl(2,2) ; 
  deriveAt(u,v,1,Skl);
  T tmp = norm(crossProduct(Skl(1,0),Skl(0,1)));
  return tmp ; 
}
 
template <class T, int N>
int NurbsSurface<T,N>::ok() {
  if(P.rows() <= degU)
    return 0 ;
  if(P.cols() <= degV)
    return 0 ;
  if(P.rows() != U.n()+degU+1)
    return 0 ;
  if(P.cols() != V.n()+degV+1)
    return 0 ;
  return 1 ;
}
template <class T, int N>
void NurbsSurface<T,N>::reset(const Matrix< HPoint_nD<T,N> >& Pts, const Vector<T> &U1, const Vector<T> &V1){
  P = Pts;
  U = U1;
  V = V1;
  degU = U.size()-P.rows()-1;
  degV = V.size()-P.cols()-1;
}

template <class T, int N>
void NurbsSurface<T,N>::resize(int Pu, int Pv, int DegU, int DegV){
  P.resize(Pu,Pv) ;
  degU = DegU ;
  degV = DegV ;
  U.resize(Pu+DegU+1) ;
  V.resize(Pv+DegV+1) ;
}

template <class T, int N>
void NurbsSurface<T,N>::resizeKeep(int Pu, int Pv, int DegU, int DegV){
  P.resizeKeep(Pu,Pv) ;
  degU = DegU ;
  degV = DegV ;
  U.resize(Pu+DegU+1) ;
  V.resize(Pv+DegV+1) ;
}

template <class T, int D>
int NurbsSurface<T,D>::skinV(NurbsCurveArray<T,D>& ca, int dV) {
  Vector<T> vk(ca.n()) ;
  //Vector<T> d(ca.n()) ;
  T* d ;
  int i,k,N,K ;

  if(ca.n()<dV){
#ifdef USE_EXCEPTION
    throw NurbsInputError();
#else
    Error err("NurbsSurface<T,D> skinV") ;
    err << "The number of curves are insufficient for the degree of interpolation specified.\n" ;
    err << "n of curves = " << ca.n() << ", degree of interpolation = " << dV << endl ;
    err.warning() ;
    return 0 ;
#endif
  }

  generateCompatibleCurves(ca) ;

  K = ca.n() ;
  N = ca[0].ctrlPnts().n() ;

  resize(N,K,ca[0].degree(),dV) ;

  //d.resize(ca[0].ctrlPnts().n()) ;
  d = new T[ca[0].ctrlPnts().n()] ; 
  for(i=0;i<N;i++){
    d[i] = 0 ;
    for(k=1;k<K;k++){
      d[i] += norm(ca[k].ctrlPnts(i)-ca[k-1].ctrlPnts(i)) ;
    }
  }

  vk[0] = 0 ;
  for(k=1;k<K;k++){
    vk[k] = 0 ;
    for(i=0;i<N;i++)
      vk[k] += norm(ca[k].ctrlPnts(i)-ca[k-1].ctrlPnts(i))/d[i] ;
    vk[k] /= (T)N ;//+ 1.0 ;
    vk[k] += vk[k-1] ;
  }
  vk[vk.n()-1] = 1.0 ;

  for(i=1;i<K-degV;i++){
    V[i+degV] = 0 ;
    for(k=i;k<i+degV;k++)
      V[i+degV] += vk[k] ;
    V[i+degV] /= (T)degV ;
  }
  for(i=0;i<=degV;i++) V[i] = 0.0 ;
  for(i=V.n()-degV-1;i<V.n();i++) V[i] = 1.0 ;


  Vector< HPoint_nD<T,D> > Q(K) ;
  NurbsCurve<T,D> i_curve ;


  for(i=0;i<N;i++){
    for(k=0;k<K;k++)
      Q[k] = ca[k].ctrlPnts(i) ;
    i_curve.globalInterpH(Q,vk,V,degV);
    for(k=0;k<K;k++)
      P(i,k) = i_curve.ctrlPnts(k) ;
  }

  U = ca[0].knot() ;

  delete []d ; 
  return 1 ;
}

template <class T, int nD>
int NurbsSurface<T,nD>::skinU(NurbsCurveArray<T,nD>& ca, int dU) {
  Vector<T> uk(ca.n());
  //Vector<T> d(ca.n()) ;
  T* d; 
  int i,k,N,K ;

  if(ca.n()<dU){
#ifdef USE_EXCEPTION
    throw NurbsInputError();
#else
    Error err("NurbsSurface<T,N> skinU") ;
    err << "The number of curves are insufficient for the degree of interpolation specified.\n" ;
    err << "n of curves = " << ca.n() << ", degree of interpolation = " << dU << endl ;
    err.warning() ;
    return 0 ;
#endif
  }

  generateCompatibleCurves(ca) ;

  K = ca.n() ;
  N = ca[0].ctrlPnts().n() ;

  resize(K,N,dU,ca[0].degree()) ;

  //d.resize(ca[0].ctrlPnts().n()) ;
  d = new T[ca[0].ctrlPnts().n()] ; 
  for(i=0;i<N;i++){
    d[i] = 0 ;
    for(k=1;k<K;k++){
      d[i] += norm(ca[k].ctrlPnts(i)-ca[k-1].ctrlPnts(i)) ;
    }
  }

  uk[0] = 0 ;
  for(k=1;k<K;k++){
    uk[k] = 0 ;
    for(i=0;i<N;i++)
      uk[k] += norm(ca[k].ctrlPnts(i)-ca[k-1].ctrlPnts(i))/d[i] ;
    uk[k] /= (T)N ;//+ 1.0 ;
    uk[k] += uk[k-1] ;
  }
  uk[uk.n()-1] = 1.0 ;

  for(i=1;i<K-degU;i++){
    U[i+degU] = 0 ;
    for(k=i;k<i+degU;k++)
      U[i+degU] += uk[k] ;
    U[i+degU] /= (T)degU ;
  }
  for(i=0;i<=degU;i++) U[i] = 0.0 ;
  for(i=U.n()-degU-1;i<U.n();i++) U[i] = 1.0 ;


  Vector< HPoint_nD<T,nD> > Q(K) ;
  NurbsCurve<T,nD> i_curve ;


  for(i=0;i<N;i++){
    for(k=0;k<K;k++)
      Q[k] = ca[k].ctrlPnts(i) ;
    i_curve.globalInterpH(Q,uk,U,degU) ;
    for(k=0;k<K;k++)
      P(k,i) = i_curve.ctrlPnts(k) ;
  }

  V = ca[0].knot() ;

  delete []d ; 

  return 1 ;
}

template <class T, int N>
int surfMeshParams(const Matrix< Point_nD<T,N> >& Q, Vector<T>& uk, Vector<T>& vl){
  int n,m,k,l,num ;
  double d,total ;
  //Vector<T> cds(Q.rows()) ;
  T* cds = new T[maximum(Q.rows(),Q.cols())] ; // alloca might not have enough space 

  n = Q.rows() ;
  m = Q.cols() ;
  uk.resize(n) ;
  vl.resize(m) ;
  num = m ;
  

  // Compute the uk
  uk.reset(0) ;

  for(l=0;l<m;l++){
    total = 0.0 ;
    for(k=1;k<n;k++){
      cds[k] = norm(Q(k,l)-Q(k-1,l)) ;
      total += cds[k] ;
    }
    if(total==0.0) 
      num-- ;
    else {
      d = 0.0 ;
      for(k=1;k<n;k++){
        d += cds[k] ;
        uk[k] += d/total ;
      }
    }
  }

  if(num==0) {
    delete []cds ; 
    return 0 ;
  }
  for(k=1;k<n-1;k++)
    uk[k] /= num ;
  uk[n-1] = 1.0 ;

  // Compute the vl
  vl.reset(0) ;

  //cds.resize(m) ; // this line removed since the maximum is allocated at the beginning

  num = n ;

  for(k=0;k<n;k++){
    total = 0.0 ;
    for(l=1;l<m;l++){
      cds[l] = norm(Q(k,l)-Q(k,l-1)) ;
      total += cds[l] ;
    }
    if(total==0.0) 
      num-- ;
    else {
      d = 0.0 ;
      for(l=1;l<m;l++){
        d += cds[l] ;
        vl[l] += d/total ;
      }
    }
  }

  delete []cds ; 

  if(num==0) 
    return 0 ;
  for(l=1;l<m-1;l++)
    vl[l] /= num ;
  vl[m-1] = 1.0 ;


  return 1 ;
}

template <class T, int N>
int surfMeshParamsH(const Matrix< HPoint_nD<T,N> >& Q, Vector<T>& uk, Vector<T>& vl){
  int n,m,k,l,num ;
  double d,total ;
  //Vector<T> cds(Q.rows()) ;
  T* cds = new T[maximum(Q.rows(),Q.cols())] ;

  n = Q.rows() ;
  m = Q.cols() ;
  uk.resize(n) ;
  vl.resize(m) ;
  num = m ;
  

  // Compute the uk
  uk.reset(0) ;

  for(l=0;l<m;l++){
    total = 0.0 ;
    for(k=1;k<n;k++){
      cds[k] = norm(Q(k,l)-Q(k-1,l)) ;
      total += cds[k] ;
    }
    if(total==0.0) 
      num-- ;
    else {
      d = 0.0 ;
      for(k=1;k<n;k++){
        d += cds[k] ;
        uk[k] += d/total ;
      }
    }
  }

  if(num==0) {
    delete []cds ; 
    return 0 ;
  }
  for(k=1;k<n-1;k++)
    uk[k] /= num ;
  uk[n-1] = 1.0 ;

  // Compute the vl
  vl.reset(0) ;
  //cds.resize(m) ; // taking the maximum so this is not needed

  num = n ;

  for(k=0;k<n;k++){
    total = 0.0 ;
    for(l=1;l<m;l++){
      cds[l] = norm(Q(k,l)-Q(k,l-1)) ;
      total += cds[l] ;
    }
    if(total==0.0) 
      num-- ;
    else {
      d = 0.0 ;
      for(l=1;l<m;l++){
        d += cds[l] ;
        vl[l] += d/total ;
      }
    }
  }

  delete []cds ; 

  if(num==0) 
    return 0 ;
  for(l=1;l<m-1;l++)
    vl[l] /= num ;
  vl[m-1] = 1.0 ;
  return 1 ;
}

template <class T, int N>
void NurbsSurface<T,N>::globalInterp(const Matrix< Point_nD<T,N> >& Q, int pU, int pV){
  Vector<T> vk,uk ;

  resize(Q.rows(),Q.cols(),pU,pV) ;

  surfMeshParams(Q,uk,vk) ;
  knotAveraging(uk,pU,U) ;
  knotAveraging(vk,pV,V) ;
  

  Vector< HPoint_nD<T,N> > Pts(Q.rows()) ;
  NurbsCurve<T,N> R ;
  
  int i,j ;

  for(j=0;j<Q.cols();j++){
    for(i=0;i<Q.rows();i++)
      Pts[i] = Q(i,j) ;
    R.globalInterpH(Pts,uk,U,pU);
    for(i=0;i<Q.rows();i++)
      P(i,j) = R.ctrlPnts(i) ;
  }

  Pts.resize(Q.cols()) ;
  for(i=0;i<Q.rows();i++){
    for(j=0;j<Q.cols();j++)
      Pts[j] = P(i,j) ;
    R.globalInterpH(Pts,vk,V,pV) ;
    for(j=0;j<Q.cols();j++)
      P(i,j) = R.ctrlPnts(j) ;
  }
}

template <class T, int N>
void NurbsSurface<T,N>::globalInterpH(const Matrix< HPoint_nD<T,N> >& Q, int pU, int pV){
  Vector<T> vk,uk ;

  resize(Q.rows(),Q.cols(),pU,pV) ;

  surfMeshParamsH(Q,uk,vk) ;
  knotAveraging(uk,pU,U) ;
  knotAveraging(vk,pV,V) ;
  

  Vector< HPoint_nD<T,N> > Pts(Q.rows()) ;
  NurbsCurve<T,N> R ;
  
  int i,j ;

  for(j=0;j<Q.cols();j++){
    for(i=0;i<Q.rows();i++)
      Pts[i] = Q(i,j) ;
    R.globalInterpH(Pts,uk,U,pU);
    for(i=0;i<Q.rows();i++)
      P(i,j) = R.ctrlPnts(i) ;
  }

  Pts.resize(Q.cols()) ;
  for(i=0;i<Q.rows();i++){
    for(j=0;j<Q.cols();j++)
      Pts[j] = P(i,j) ;
    R.globalInterpH(Pts,vk,V,pV) ;
    for(j=0;j<Q.cols();j++)
      P(i,j) = R.ctrlPnts(j) ;
  }
}

template <class T, int N>
void NurbsSurface<T,N>::leastSquares(const Matrix< Point_nD<T,N> >& Q, int pU, int pV, int nU, int nV){
  Vector<T> vk,uk ;

  resize(nU,nV,pU,pV) ;

  surfMeshParams(Q,uk,vk) ;

  Vector< HPoint_nD<T,N> > Pts(Q.rows()) ;
  NurbsCurve<T,N> R ;
  
  int i,j ;

  Matrix< HPoint_nD<T,N> > P2 ;

  P2.resize(nU,Q.cols()) ;

  for(j=0;j<Q.cols();j++){
    for(i=0;i<Q.rows();i++)
      Pts[i] = Q(i,j) ;
    R.leastSquaresH(Pts,pU,nU,uk);
    for(i=0;i<P.rows();i++)
      P2(i,j) = R.ctrlPnts(i) ;
    if(j==0)
      U = R.knot() ;
  }

  Pts.resize(Q.cols()) ;
  for(i=0;i<P.rows();i++){
    for(j=0;j<Q.cols();j++)
      Pts[j] = P2(i,j) ;
    R.leastSquaresH(Pts,pV,nV,vk) ;
    for(j=0;j<P.cols();j++)
      P(i,j) = R.ctrlPnts(j) ;
    if(i==0)
      V = R.knot() ;
  }
}

template <class T, int N>
void globalSurfInterpXY(const Matrix< Point_nD<T,N> >& Q, int pU, int pV, NurbsSurface<T,N>& S) {
  Vector<T> uk,vk ;
  T um,uM ;
  T vm,vM ;

  um = Q(0,0).y() ;
  vm = Q(0,0).x() ;
  uM = Q(Q.rows()-1,0).y() ;
  vM = Q(0,Q.cols()-1).x() ;

  uk.resize(Q.rows()) ;
  vk.resize(Q.cols()) ;

  uk[0] = 0.0 ;
  vk[0] = 0.0 ;
  uk[uk.n()-1] = 1.0 ;
  vk[vk.n()-1] = 1.0 ;

  T dU = uM-um ;
  T dV = vM-vm ;

  int i ;
  for(i=1;i<uk.n()-1;++i){
    uk[i] = Q(i,0).y()/dU ;
  }
  
  for(i=1;i<vk.n()-1;++i){
    vk[i] = Q(0,i).x()/dV ;
  }
  
  globalSurfInterpXY(Q,pU,pV,S,uk,vk) ;
}

template <class T, int N>
void globalSurfInterpXY(const Matrix< Point_nD<T,N> >& Q, int pU, int pV, NurbsSurface<T,N>& S, const Vector<T>& uk, const Vector<T>& vk){
  Vector<T> V,U ;
  int i,j ;


  knotAveraging(uk,pU,U) ;
  knotAveraging(vk,pV,V) ;
  
  Vector< HPoint_nD<T,N> > P(Q.rows()) ;
  NurbsCurve<T,N> R ;
  
  S.resize(Q.rows(),Q.cols(),pU,pV) ;
  S.U = U ;
  S.V = V ;

  for(j=0;j<Q.cols();j++){
    for(i=0;i<Q.rows();i++)
      P[i] = Q(i,j) ;
    R.globalInterpH(P,uk,U,pU) ;
    for(i=0;i<Q.rows();i++)
      S.P(i,j) = R.ctrlPnts(i) ;
  }

  P.resize(Q.cols()) ;
  for(i=0;i<Q.rows();i++){
    for(j=0;j<Q.cols();j++)
      P[j] = S.P(i,j) ;
    R.globalInterpH(P,vk,V,pV) ;
    for(j=0;j<Q.cols();j++)
      S.P(i,j) = R.ctrlPnts(j) ;
  }
}


template <class T, int N>
void globalSurfApprox(const Matrix< Point_nD<T,N> >& Q, int pU, int pV, NurbsSurface<T,N>& S, double error){

  NurbsCurveArray<T,N> R ; 
  Vector< Point_nD<T,N> > P ;

  Matrix< HPoint_nD<T,N> > St ;
  Vector<T> Ut,Vt ;

  Vector<T> vk,uk ;

  surfMeshParams(Q,uk,vk) ;
  
  R.resize(Q.cols()) ;
  P.resize(Q.rows()) ;

  int i,j ;

  for(j=0;j<Q.cols();j++){
    for(i=0;i<Q.rows();i++)
      P[i] = Q(i,j) ;
    R[j].globalApproxErrBnd(P,uk,pU,error) ;
  }


  generateCompatibleCurves(R) ;

  Ut.resize(R[0].knot().n()) ;
  Ut = R[0].knot() ;

  St.resize(R[0].ctrlPnts().n(),R.n()) ;

  for(i=0;i<R[0].ctrlPnts().n();i++){
    for(j=0;j<R.n();j++)
      St(i,j) = R[j].ctrlPnts(i) ;
  }

  P.resize(St.cols()) ;
  R.resize(St.rows()) ;
  for(i=0;i<St.rows();i++){
    for(j=0;j<St.cols();j++)
      P[j] = project(St(i,j)) ;
    R[i].globalApproxErrBnd(P,vk,pV,error) ;
  }

  generateCompatibleCurves(R) ;

  Vt.resize(R[0].knot().n()) ;
  Vt = R[0].knot() ;
  
  S.resize(St.rows(),R[0].ctrlPnts().n(),pU,pV) ;
  for(i=0;i<S.ctrlPnts().rows();i++)
    for(j=0;j<S.ctrlPnts().cols();j++){
      S.P(i,j) = R[i].ctrlPnts(j) ;
    }
  S.U = Ut ;
  S.V = Vt ;
}


template <class T, int N>
void NurbsSurface<T,N>::degreeElevate(int tU, int tV) {
  degreeElevateU(tU) ;
  degreeElevateV(tV) ;
}

template <class T, int N>
void NurbsSurface<T,N>::degreeElevateU(int t) {
  if(t<=0){
    return ;
  }

  NurbsSurface<T,N> S(*this) ;
  
  int i,j,k ;
  int n = S.ctrlPnts().rows()-1;
  int p = S.degU ;
  int m = n+p+1;
  int ph = p+t ;
  int ph2 = ph/2 ;
  Matrix<T> bezalfs(p+t+1,p+1) ; // coefficients for degree elevating the Bezier segment
  Matrix< HPoint_nD<T,N> > bpts(p+1,S.P.cols()) ; // pth-degree Bezier control points of the current segment
  Matrix< HPoint_nD<T,N> > ebpts(p+t+1,S.P.cols()) ; // (p+t)th-degree Bezier control points of the  current segment
  Matrix< HPoint_nD<T,N> > Nextbpts(p-1,S.P.cols()) ; // leftmost control points of the next Bezier segment
  Vector<T> alphas(p-1) ; // knot instertion alphas.
  
  // Compute the Binomial coefficients
  Matrix<T> Bin(ph+1,ph2+1) ;
  binomialCoef(Bin) ;
  
  // Compute Bezier degree elevation coefficients
  T inv,mpi ;
  bezalfs(0,0) = bezalfs(ph,p) = 1.0 ;
  for(i=1;i<=ph2;i++){
    inv= 1.0/Bin(ph,i) ;
    mpi = minimum(p,i) ;
    for(j=maximum(0,i-t); j<=mpi; j++){
      bezalfs(i,j) = inv*Bin(p,j)*Bin(t,i-j) ;
    }
  }
  
  for(i=ph2+1;i<ph ; i++){
    mpi = minimum(p,i) ;
    for(j=maximum(0,i-t); j<=mpi ; j++)
      bezalfs(i,j) = bezalfs(ph-i,p-j) ;
  }
  
  // Allocate more control points than necessary
  resize(S.P.rows()+S.P.rows()*t,S.P.cols(),ph,S.degV) ; 
  
  int colJ ;
  int mh = ph ;
  int kind = ph+1 ;
  T ua = S.U[0] ;
  T ub = 0.0 ;
  int r=-1 ; 
  int oldr ;
  int a = p ;
  int b = p+1 ; 
  int cind = 1 ;
  int rbz,lbz = 1 ; 
  int mul,save,s;
  T alf ;
  int first, last, kj ;
  T den,bet,gam,numer ;
  
  for(i=0;i<S.P.cols();i++)
    P(0,i) = S.P(0,i) ;
  for(i=0; i <= ph ; i++){
    U[i] = ua ;
  }
  
  // Initialize the first Bezier segment
  
  for(i=0;i<=p ;i++) 
    for(j=0;j<S.P.cols();j++)
      bpts(i,j) = S.P(i,j);
  
  while(b<m){ // Big loop thru knot vector
    i=b ;
    while(b<m && S.U[b] >= S.U[b+1]) // for some odd reasons... == doesn't work
      b++ ;
    mul = b-i+1 ; 
    mh += mul+t ;
    ub = S.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){ // Insert knot to get Bezier segment
      numer = ub-ua ;
      for(k=p;k>mul;k--){
        alphas[k-mul-1] = numer/(S.U[a+k]-ua) ;
      }
      for(j=1;j<=r;j++){
        save = r-j ; s = mul+j ;
        for(k=p;k>=s;k--){
          for(colJ=0;colJ<S.P.cols();colJ++){
            bpts(k,colJ) = alphas[k-s]*bpts(k,colJ)+(1.0-alphas[k-s])*bpts(k-1,colJ) ;}
        }
        if(Nextbpts.rows()>0)
          for(colJ=0;colJ<S.P.cols();colJ++){
            Nextbpts(save,colJ) = bpts(p,colJ) ;}
      }
    }
    
    for(i=lbz;i<=ph;i++){ // Degree elevate Bezier,  only the points lbz,...,ph are used
      for(colJ=0;colJ<S.P.cols();colJ++){
        ebpts(i,colJ) = 0.0 ;}
      mpi = minimum(p,i) ;
      for(j=maximum(0,i-t); j<=mpi ; j++)
        for(colJ=0;colJ<S.P.cols();colJ++){
          ebpts(i,colJ) += bezalfs(i,j)*bpts(j,colJ) ;}
    }
    
    if(oldr>1){ // Must remove knot u=c.U[a] oldr times
      // if(oldr>2) // Alphas on the right do not change
      //        alfj = (ua-nc.U[kind-1])/(ub-nc.U[kind-1]) ;
      first = kind-2 ; last = kind ;
      den = ub-ua ;
      bet = (ub-U[kind-1])/den ;
      for(int tr=1; tr<oldr; tr++){ // Knot removal loop
        i = first ; j = last ;
        kj = j-kind+1 ;
        while(j-i>tr){ // Loop and compute the new control points for one removal step
          if(i<cind){
            alf=(ub-U[i])/(ua-U[i]) ;
            for(colJ=0;colJ<S.P.cols();colJ++){
              P(i,colJ) = alf*P(i,colJ) + (1.0-alf)*P(i-1,colJ) ;}
          }
          if( j>= lbz){
            if(j-tr <= kind-ph+oldr){
              gam = (ub-U[j-tr])/den ;
              for(colJ=0;colJ<S.P.cols();colJ++){
                ebpts(kj,colJ) = gam*ebpts(kj,colJ) + (1.0-gam)*ebpts(kj+1,colJ) ;}
            }
            else{
              for(colJ=0;colJ<S.P.cols();colJ++){
                ebpts(kj,colJ) = bet*ebpts(kj,colJ)+(1.0-bet)*ebpts(kj+1,colJ) ;}
            }
          }
          ++i ; --j; --kj ;
        }
        --first ; ++last ;
      }
    }
  
    if(a!=p) // load the knot u=c.U[a]
      for(i=0;i<ph-oldr; i++){
        U[kind++] = ua ; 
      }
    
    for(j=lbz; j<=rbz ; j++) { // load control points onto nc
      for(colJ=0;colJ<S.P.cols();colJ++)
        P(cind,colJ) = ebpts(j,colJ) ; 
      ++cind ;
    }
    
    if(b<m){ // Set up for next pass thru loop
      for(colJ=0;colJ<S.P.cols();colJ++){
        for(j=0;j<r;j++)
          bpts(j,colJ) = Nextbpts(j,colJ) ;
        for(j=r;j<=p;j++)
          bpts(j,colJ) = S.P(b-p+j,colJ) ;
      }
      a=b ; 
      b++ ;
      ua = ub ;
    }
    else{
      for(i=0;i<=ph;i++)
        U[kind+i] = ub ;
    }
  }
  // Resize to the proper number of control points  
  resizeKeep(mh-ph,S.P.cols(),ph,S.degV) ;
}

template <class T, int N>
void NurbsSurface<T,N>::degreeElevateV(int t) {
  if(t<=0){
    return ;
  }

  NurbsSurface<T,N> S(*this) ;

  int i,j,k ;
  int n = S.ctrlPnts().cols()-1;
  int p = S.degV ;
  int m = n+p+1;
  int ph = p+t ;
  int ph2 = ph/2 ;
  Matrix<T> bezalfs(p+t+1,p+1) ; // coefficients for degree elevating the Bezier segment
  Matrix< HPoint_nD<T,N> > bpts(p+1,S.P.rows()) ; // pth-degree Bezier control points of the current segment
  Matrix< HPoint_nD<T,N> > ebpts(p+t+1,S.P.rows()) ; // (p+t)th-degree Bezier control points of the  current segment
  Matrix< HPoint_nD<T,N> > Nextbpts(p-1,S.P.rows()) ; // leftmost control points of the next Bezier segment
  //Vector<T> alphas(p-1) ; // knot instertion alphas.
  T* alphas = (T*) alloca((p-1)*sizeof(T)) ;
  
  // Compute the Binomial coefficients
  Matrix<T> Bin(ph+1,ph2+1) ;
  binomialCoef(Bin) ;
  
  // Compute Bezier degree elevation coefficients
  T inv,mpi ;
  bezalfs(0,0) = bezalfs(ph,p) = 1.0 ;
  for(i=1;i<=ph2;i++){
    inv= 1.0/Bin(ph,i) ;
    mpi = minimum(p,i) ;
    for(j=maximum(0,i-t); j<=mpi; j++){
      bezalfs(i,j) = inv*Bin(p,j)*Bin(t,i-j) ;
    }
  }
  
  for(i=ph2+1;i<ph ; i++){
    mpi = minimum(p,i) ;
    for(j=maximum(0,i-t); j<=mpi ; j++)
      bezalfs(i,j) = bezalfs(ph-i,p-j) ;
  }
  
  resize(S.P.rows(),S.P.cols()+S.P.cols()*t,S.degU,ph) ; // Allocate more control points than necessary
  
  int rowJ ;
  int mh = ph ;
  int kind = ph+1 ;
  T va = S.V[0] ;
  T vb = 0.0 ;
  int r=-1 ; 
  int oldr ;
  int a = p ;
  int b = p+1 ; 
  int cind = 1 ;
  int rbz,lbz = 1 ; 
  int mul,save,s;
  T alf ;
  int first, last, kj ;
  T den,bet,gam,numer ;
  
  for(i=0;i<S.P.rows();i++)
    P(i,0) = S.P(i,0) ;
  for(i=0; i <= ph ; i++){
    V[i] = va ;
  }
  
  // Initialize the first Bezier segment
  
  for(i=0;i<=p ;i++) 
    for(j=0;j<S.P.rows();j++)
      bpts(i,j) = S.P(j,i);
  
  while(b<m){ // Big loop thru knot vector
    i=b ;
    while(b<m && S.V[b] >= S.V[b+1]) // for some odd reasons... == doesn't work
      b++ ;
    mul = b-i+1 ; 
    mh += mul+t ;
    vb = S.V[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){ // Insert knot to get Bezier segment
      numer = vb-va ;
      for(k=p;k>mul;k--){
        alphas[k-mul-1] = numer/(S.V[a+k]-va) ;
      }
      for(j=1;j<=r;j++){
        save = r-j ; s = mul+j ;
        for(k=p;k>=s;k--){
          for(rowJ=0;rowJ<S.P.rows();rowJ++){
            bpts(k,rowJ) = alphas[k-s]*bpts(k,rowJ)+(1.0-alphas[k-s])*bpts(k-1,rowJ) ;}
        }
        if(Nextbpts.rows()>0)
          for(rowJ=0;rowJ<S.P.rows();rowJ++){
            Nextbpts(save,rowJ) = bpts(p,rowJ) ;}
      }
    }
    
    for(i=lbz;i<=ph;i++){ // Degree elevate Bezier,  only the points lbz,...,ph are used
      for(rowJ=0;rowJ<S.P.rows();rowJ++){
        ebpts(i,rowJ) = 0.0 ;}
      mpi = minimum(p,i) ;
      for(j=maximum(0,i-t); j<=mpi ; j++)
        for(rowJ=0;rowJ<S.P.rows();rowJ++){
          ebpts(i,rowJ) += bezalfs(i,j)*bpts(j,rowJ) ;}
    }
    
    if(oldr>1){ // Must remove knot V=S.V[a] oldr times
      // if(oldr>2) // Alphas on the right do not change
      //        alfj = (va-nc.U[kind-1])/(vb-V[kind-1]) ;
      first = kind-2 ; last = kind ;
      den = vb-va ;
      bet = (vb-V[kind-1])/den ;
      for(int tr=1; tr<oldr; tr++){ // Knot removal loop
        i = first ; j = last ;
        kj = j-kind+1 ;
        while(j-i>tr){ // Loop and compute the new control points for one removal step
          if(i<cind){
            alf=(vb-V[i])/(va-V[i]) ;
            for(rowJ=0;rowJ<S.P.rows();rowJ++){
              P(rowJ,i) = alf*P(rowJ,i) + (1.0-alf)*P(rowJ,i-1) ;}
          }
          if( j>= lbz){
            if(j-tr <= kind-ph+oldr){
              gam = (vb-V[j-tr])/den ;
              for(rowJ=0;rowJ<S.P.rows();rowJ++){
                ebpts(kj,rowJ) = gam*ebpts(kj,rowJ) + (1.0-gam)*ebpts(kj+1,rowJ) ;}
            }
            else{
              for(rowJ=0;rowJ<S.P.rows();rowJ++){
                ebpts(kj,rowJ) = bet*ebpts(kj,rowJ)+(1.0-bet)*ebpts(kj+1,rowJ) ;}
            }
          }
          ++i ; --j; --kj ;
        }
        --first ; ++last ;
      }
    }
  
    if(a!=p) // load the knot v=S.V[a]
      for(i=0;i<ph-oldr; i++){
        V[kind++] = va ; 
      }
    
    for(j=lbz; j<=rbz ; j++) { // load control points onto nc
      for(rowJ=0;rowJ<S.P.rows();rowJ++)
        P(rowJ,cind) = ebpts(j,rowJ) ; 
      ++cind ;
    }
    
    if(b<m){ // Set up for next pass thru loop
      for(rowJ=0;rowJ<S.P.rows();rowJ++){
        for(j=0;j<r;j++)
          bpts(j,rowJ) = Nextbpts(j,rowJ) ;
        for(j=r;j<=p;j++)
          bpts(j,rowJ) = S.P(rowJ,b-p+j) ;
      }
      a=b ; 
      b++ ;
      va = vb ;
    }
    else{
      for(i=0;i<=ph;i++)
        V[kind+i] = vb ;
    }
  }
  // Resize to the proper number of control points  
  resizeKeep(S.P.rows(),mh-ph,S.degU,ph) ; 
}
template <class T, int N>
int NurbsSurface<T,N>::findMultU(int r) const {
  int s=1 ;
  for(int i=r;i>degU+1;i--)
    if(U[i]<=U[i-1])
      s++ ;
    else
      return s ;
  return s ;
}

template <class T, int N>
int NurbsSurface<T,N>::findMultV(int r) const {
  int s=1 ;
  for(int i=r;i>degV+1;i--)
    if(V[i]<=V[i-1])
      s++ ;
    else
      return s ;
  return s ;
}




template <class T, int N>
void NurbsSurface<T,N>::findSpan(T u, T v, int& spanU, int& spanV) const{
  spanU = findSpanU(u) ;
  spanV = findSpanV(v) ;
}

template <class T, int N>
int NurbsSurface<T,N>::findSpanU(T u) const{
  if(u>=U[P.rows()]) 
    return P.rows()-1 ;
  if(u<=U[degU])
    return degU ;

  //AF
  int low = 0 ;
  int high = P.rows()+1 ; 
  int 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 <class T, int N>
int NurbsSurface<T,N>::findSpanV(T v) const{
  if(v>=V[P.cols()]) 
    return P.cols()-1 ;
  if(v<=V[degV])
    return degV ;

  //AF
  int low  = 0 ;
  int high = P.cols()+1 ; 
  int mid = (low+high)/2 ;

  while(v<V[mid] || v>= V[mid+1]){
    if(v<V[mid])
      high = mid ;
    else
      low = mid ;
    mid = (low+high)/2 ;
  }
  return mid ;
}

template <class T, int N>
void NurbsSurface<T,N>::basisFuns(T u, T v, int spanU, int spanV, Vector<T>& Nu, Vector<T> &Nv) const{
  basisFunsU(u,spanU,Nu) ;
  basisFunsV(v,spanV,Nv) ;
}

template <class T, int nD>
void NurbsSurface<T,nD>::basisFunsU(T u, int span, Vector<T>& N) const {
  //Vector<T> left(degU+1), right(degU+1) ;
  T* left = (T*) alloca((degU+1)*sizeof(T)) ;
  T* right = (T*) alloca((degU+1)*sizeof(T)) ;
  T temp,saved ;
   

  N.resize(degU+1) ;

  N[0] = 1.0 ;
  for(int j=1; j<= degU ; j++){
    left[j] = u-U[span+1-j] ;
    right[j] = U[span+j]-u ;
    saved = 0.0 ;
    for(int 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 <class T, int nD>
void NurbsSurface<T,nD>::basisFunsV(T v, int span, Vector<T>& N) const {
  //Vector<T> left(degV+1), right(degV+1) ;
  T* left = (T*) alloca((degV+1)*sizeof(T)) ;
  T* right = (T*) alloca((degV+1)*sizeof(T)) ;
  T temp,saved ;
   

  N.resize(degV+1) ;

  N[0] = 1.0 ;
  for(int j=1; j<= degV ; j++){
    left[j] = v-V[span+1-j] ;
    right[j] = V[span+j]-v ;
    saved = 0.0 ;
    for(int 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 <class T, int nD>
void NurbsSurface<T,nD>::deriveAtH(T u, T v, int d, Matrix< HPoint_nD<T,nD> > &skl) const {
  int k,l,du,dv;
  skl.resize(d+1,d+1) ;

  du = minimum(d,degU) ;
  for(k=degU+1;k<=d;++k)
    for(l=0;l<=d-k;++l)
      skl(k,l) = 0.0 ;
  dv=minimum(d,degV) ;
  for(l=degV+1;l<=d;++l)
    for(k=0;k<=d-l;++k)
      skl(k,l) = 0.0 ;
  int uspan = findSpanU(u) ;
  int vspan = findSpanV(v) ;
  Matrix<T> Nu,Nv ;
  nurbsDersBasisFuns(du,u,uspan,degU,U,Nu) ;
  nurbsDersBasisFuns(dv,v,vspan,degV,V,Nv) ;

  Vector< HPoint_nD<T,nD> > temp(degV+1) ;
  int dd,r,s ;
  for(k=0;k<=du;++k){
    for(s=0;s<=degV;++s){
      temp[s] = 0.0 ;
      for(r=0;r<=degU;++r)
        temp[s] += Nu(k,r)*P(uspan-degU+r,vspan-degV+s) ;
    }
    dd = minimum(d-k,dv) ;
    for(l=0;l<=dd;++l){
      skl(k,l) = 0.0 ;
      for(s=0;s<=degV;++s)
        skl(k,l) += Nv(l,s)*temp[s] ;
    }
  }
}

template <class T, int N>
void NurbsSurface<T,N>::deriveAt(T u, T v, int d, Matrix< Point_nD<T,N> > &skl) const {
  int k,l ;
  Matrix< HPoint_nD<T,N> > ders ;
  Point_nD<T,N> pv,pv2 ;
  
  skl.resize(d+1,d+1) ;

  deriveAtH(u,v,d,ders) ;
  
  Matrix<T> Bin(d+1,d+1) ;
  binomialCoef(Bin) ;
  int i,j ; 

  for(k=0;k<=d;++k){
    for(l=0;l<=d-k;++l){
      pv.x() = ders(k,l).x() ;
      pv.y() = ders(k,l).y() ;
      pv.z() = ders(k,l).z() ;
      for(j=1;j<=l;j++)
        pv -= Bin(l,j)*ders(0,j).w()*skl(k,l-j) ;
      for(i=1;i<=k;i++){
        pv -= Bin(k,i)*ders(i,0).w()*skl(k-i,l) ;
        pv2 = 0.0 ;
        for(j=1;j<=l;j++)
          pv2 += Bin(l,j)*ders(i,j).w()*skl(k-i,l-j) ;
        pv -= Bin(k,i)*pv2 ;
      }
      skl(k,l) = pv/ders(0,0).w() ;
    }
  }
}

template <class T, int N>
Point_nD<T,N> NurbsSurface<T,N>::normal(T u, T v) const {
  Matrix< Point_nD<T,N> > ders ;

  deriveAt(u,v,1,ders) ;

  return crossProduct(ders(1,0),ders(0,1)) ;
}

template <class T, int N>
HPoint_nD<T,N> NurbsSurface<T,N>::operator()(T u, T v) const{
  int uspan = findSpanU(u) ;
  int vspan = findSpanV(v) ;
  Vector<T> Nu,Nv ;
  Vector< HPoint_nD<T,N> > temp(degV+1)  ;

  basisFuns(u,v,uspan,vspan,Nu,Nv) ;

  int l;
  for(l=0;l<=degV;l++){
    temp[l] =0.0 ;
    for(int k=0;k<=degU;k++){
      temp[l] += Nu[k]*P(uspan-degU+k,vspan-degV+l) ;
    }
  }
  HPoint_nD<T,N> sp(0,0,0,0) ;
  for(l=0;l<=degV;l++){
    sp += Nv[l]*temp[l] ;
  }
  return sp ; 
}

inline
int max3(int a,int b, int c){
  int m = a ;
  if(m <b)
    m = b ;
  if(m <c)
    m = c ;
  return m ;
}

template <class T, int N>
void gordonSurface(NurbsCurveArray<T,N>& lU, NurbsCurveArray<T,N>& lV, const Matrix< Point_nD<T,N> >& intersections, NurbsSurface<T,N>& gS){
  NurbsSurface<T,N> sU,sV,sI ;
  sU.skinU(lU,3) ;
  sV.skinV(lV,3) ;
  sI.globalInterp(intersections,3,3) ;

  int du = max3(sU.degU,sV.degU,sI.degU) ;
  int dv = max3(sU.degV,sV.degV,sI.degV) ;
  

  NurbsSurface<T,N> SU,SV,SI ;
  degreeElevate(sU,du-sU.degU,dv-sU.degV,SU) ;
  degreeElevate(sV,du-sV.degU,dv-sV.degV,SV) ;
  degreeElevate(sI,du-sI.degU,dv-sI.degV,SI) ;


  Vector<T> U,V ;
  U = knotUnion(SU.knotU(),SV.knotU()) ;
  U = knotUnion(U,SI.knotU()) ;
  V = knotUnion(SU.knotV(),SV.knotV()) ;
  V = knotUnion(V,SI.knotV()) ;

  SU.mergeKnots(U,V) ;
  SV.mergeKnots(U,V) ;
  SI.mergeKnots(U,V) ;
  
  gS = SU ;

  for(int i=0;i<gS.P.rows();i++)
    for(int j=0;j<gS.P.cols();j++){
      gS.P(i,j) += SV.P(i,j) ;
      gS.P(i,j) -= SI.P(i,j) ;
    }
}

template <class T>
inline void insert(T u, Vector<T>& v){
  int i ;
  if(u<v[0] || u>v[v.n()-1])
    return ;
  v.resize(v.n()+1) ;
  i = v.n()-1 ;
  while(v[i-1]>u){
    v[i] = v[i-1] ;
    --i ;
  }
  v[i] = u ;
}


template <class T, int N>
void NurbsSurface<T,N>::sweep(const NurbsCurve<T,N>& Trj, const NurbsCurve<T,N>& C, const NurbsCurve<T,N>& Sv, int K, int useAy, int invAz) {
  int i,j,k,m,q,ktv,nsect ;
  q = Trj.degree() ;
  ktv = Trj.knot().n() ;
  nsect = K  ;
  
  V.resize(Trj.knot().n()) ;
  V = Trj.knot() ;

  if(ktv <= nsect+q){
    m = nsect+q-ktv+1 ;
    // insert m knots into T(v)
    // locations are not critical so inserting in the middle of the 
    // biggest span is used.
    for(i=0;i<m;++i){
      T md,mt ;
      T mu = -1;
      md = 0 ;
      for(j=1;j<V.n();++j){
        mt = V[j]-V[j-1] ;
        if(mt>md){
          md = mt ;
          mu = (V[j]+V[j-1])/2.0 ;
        }
      }
      insert(mu,V) ;
    }
  }
  else{
    if(ktv>nsect){ // must increase the number of instances of C(u)
      nsect = ktv-q-1 ; 
    }
  }

  Vector<T> v; 

  // Compute the parameters by averaging the knots
  v.resize(nsect) ;
  v[0] = Trj.knot(Trj.degree()) ;
  v[nsect-1] = Trj.knot(Trj.knot().n()-Trj.degree()-1) ;
  for(k=1;k<nsect-1;++k){ 
    v[k] = 0.0 ;
    for(i=k+1;i<k+q+1;++i) 
      v[k] += V[i] ;
    v[k] /= q ;
  }

  resize(C.ctrlPnts().n(),nsect,C.degree(),Trj.degree()) ;

  U = C.knot() ;

  // setup Bv ;
  Vector< Point_nD<T,N> > B ;
  NurbsCurve<T,N> Bv ;
  Vector< Point_nD<T,N> > Td ;

  B.resize(v.n()) ;
  Trj.deriveAt(v[0],1,Td) ;

  B[0] = Td[1] ;
  if(Td[1].y() ==0)
    B[0] = crossProduct(Point_nD<T,N>(0,1,0),B[0]) ;
  else
    B[0] = crossProduct(Point_nD<T,N>(1,0,0),B[0]) ;
  B[0] /= norm(B[0]) ;


  for(i=1;i<v.n();++i){
    Trj.deriveAt(v[i],1,Td);
    Point_nD<T,N> Ti(Td[1]) ;
    Ti = Ti/norm(Ti) ;
    Point_nD<T,N> bi ;
    bi = B[i-1]-(B[i-1]*Ti)*Ti ;
    B[i] = bi/norm(bi) ;
  }

  Bv.globalInterp(B,v,minimum(3,B.n()-1)) ;

  Vector< HPoint_nD<T,N> > Q(C.ctrlPnts().n()) ; 
  for(k=0;k<nsect;++k){
    // scale the control points by Sv(v[k])
    for(i=0;i<Q.n();++i){
      HPoint_nD<T,N> sk(Sv(v[k])) ;
      Q[i].x() = sk.x()*C.ctrlPnts(i).x() ;
      Q[i].y() = sk.y()*C.ctrlPnts(i).y() ;
      Q[i].z() = sk.z()*C.ctrlPnts(i).z() ;
      Q[i].w() = sk.w()*C.ctrlPnts(i).w() ;
    }
    // compute o(v[k])
    Point_nD<T,N> o = Trj.pointAt(v[k]) ;
    //T w = T(v[k]).w() ;

    // compute x(v[k])
    Trj.deriveAt(v[k],1,Td) ;

    Point_nD<T,N> x = Td[1]/norm(Td[1]) ;

    // compute z(v[k])
    Point_nD<T,N> z = Bv.pointAt(v[k]) ; 
    z /= norm(z) ;

    // compute y(v[k]) 
    Point_nD<T,N> y = crossProduct(z,x) ;

    /*
    // compute the transform matrix
    double az = M_PI+atan2(y.y(),y.x()) ;
    double ax = -atan2(z.y(),z.z()) ;
    double ay = 0 ;
    if(useAy){
      ay = atan2(x.z(),x.x()) ;
    }
    if(invAz){
      az = -1.0*az ;
    }
    MatrixRT_DOUBLE A(ax,ay,az,o.x(),o.y(),o.z()) ;
    */

    MatrixRT_DOUBLE R ; // M(4,4)
    R(0,0) = y.x();
    R(1,0) = y.y();
    R(2,0) = y.z();
    R(0,1) = x.x();
    R(1,1) = x.y();
    R(2,1) = x.z();
    R(0,2) = z.x();
    R(1,2) = z.y();
    R(2,2) = z.z();
    //R(3,3) = 1.0; 

    MatrixRT_DOUBLE Tx ;
    Tx.translate(o.x(),o.y(),o.z());
    //MatrixRT_DOUBLE R(M) ;
    MatrixRT_DOUBLE A ;
    A = Tx * R ;

    for(i=0;i<Q.n();++i){
      P(i,k) = A*(Q[i]) ;
    }
  }

  // interpolate along V
  Q.resize(P.cols()) ;
  NurbsCurve<T,N> R;
  for(i=0;i<P.rows();++i){
    for(k=0;k<P.cols();++k)
      Q[k] = P(i,k) ;
    R.globalInterpH(Q,v,V,degV) ;
    for(k=0;k<P.cols();++k)
      P(i,k) = R.ctrlPnts(k) ;
  }
  
}

template <class T, int N>
void NurbsSurface<T,N>::sweep(const NurbsCurve<T,N>& Trj, const NurbsCurve<T,N>& C, int K, int useAy,int invAz) {
  // setup Sv 
  Vector< HPoint_nD<T,N> > p(2) ;
  p[0] = HPoint_nD<T,N>(1,1,1,1) ;
  p[1] = HPoint_nD<T,N>(1,1,1,1) ;
  Vector<T> u(4) ;
  u[0] = u[1] = 0.0 ; u[2] = u[3] = 1.0 ;

  NurbsCurve<T,N> Sv(p,u,1) ;
  
  sweep(Trj,C,Sv,K,useAy,invAz) ;
}

template <class T, int N>
void NurbsSurface<T,N>::transform(const MatrixRT<T>& A){
  for(int i=0;i<P.rows();++i)
    for(int j=0;j<P.cols();++j)
      P(i,j) = A*P(i,j) ;
}

template <class T, int N>
void NurbsSurface<T,N>::refineKnots(const Vector<T>& nU, const Vector<T>& nV){
  refineKnotU(nU) ;
  refineKnotV(nV) ;
}

template <class T, int N>
void NurbsSurface<T,N>::refineKnotU(const Vector<T>& X) {
  if(X.n()<=0)
    return ;
  int n = P.rows()-1 ;
  int p = degU;
  int m = n+p+1 ;
  int a,b ;
  int r = X.n()-1 ;
  NurbsSurface<T,N> nS ;
  
  nS = *this ;
  nS.resize(r+1+n+1,P.cols(),degU,degV) ;

  a = findSpanU(X[0]) ;
  b = findSpanU(X[r]) ;
  ++b ;

  int j,col ;
  for(col=0;col<P.cols();col++){
    for(j=0; j<=a-p ; j++)
      nS.P(j,col) = P(j,col);
    for(j = b-1 ; j<=n ; j++)
      nS.P(j+r+1,col) = P(j,col) ;
  }
  for(j=0; j<=a ; j++)
    nS.U[j] = U[j] ;
  for(j=b+p ; j<=m ; j++)
    nS.U[j+r+1] = U[j] ;
  int i = b+p-1 ; 
  int k = b+p+r ;
  for(j=r; j>=0 ; j--){
    while(X[j] <= U[i] && i>a){
      for(col=0;col<P.cols();col++)
        nS.P(k-p-1,col) = P(i-p-1,col) ;
      nS.U[k] = U[i] ;
      --k ;
      --i ;
    }
    for(col=0;col<P.cols();col++)
        nS.P(k-p-1,col) = nS.P(k-p,col) ;
    for(int l=1; l<=p ; l++){
      int ind = k-p+l ;
      T alpha = nS.U[k+l] - X[j] ;
      if(alpha==0.0)
        for(col=0;col<P.cols();col++)
          nS.P(ind-1,col) = nS.P(ind,col) ;
      else
        alpha /= nS.U[k+l]-U[i-p+l] ;
      for(col=0;col<P.cols();col++)
        nS.P(ind-1,col) = alpha*nS.P(ind-1,col) + (1.0-alpha)*nS.P(ind,col) ;
    }
    nS.U[k] = X[j] ;
    --k ;
  }
  *this = nS ;
}

template <class T, int N>
void NurbsSurface<T,N>::refineKnotV(const Vector<T>& X) {
  if(X.n()<=0)
    return ;
  int n = P.cols()-1 ;
  int p = degV;
  int m = n+p+1 ;
  int a,b ;
  int r = X.n()-1 ;
  NurbsSurface<T,N> nS ;
  
  try {
    nS = *this ;
    nS.resize(P.rows(),r+1+n+1,degU,degV) ;
  }
  catch(...){
    cerr << "Out of memory\n" ; 
  }

  a = findSpanV(X[0]) ;
  b = findSpanV(X[r]) ;
  ++b ;

  int j,row ;
  for(row=0;row<P.rows();row++){
    for(j=0; j<=a-p ; j++)
      nS.P(row,j) = P(row,j);
    for(j = b-1 ; j<=n ; j++)
      nS.P(row,j+r+1) = P(row,j) ;
  }
  for(j=0; j<=a ; j++)
    nS.V[j] = V[j] ;
  for(j=b+p ; j<=m ; j++)
    nS.V[j+r+1] = V[j] ;
  int i = b+p-1 ; 
  int k = b+p+r ;
  for(j=r; j>=0 ; j--){
    while(X[j] <= V[i] && i>a){
      for(row=0;row<P.rows();row++)
        nS.P(row,k-p-1) = P(row,i-p-1) ;
      nS.V[k] = V[i] ;
      --k ;
      --i ;
    }
    for(row=0;row<P.rows();row++)
        nS.P(row,k-p-1) = nS.P(row,k-p) ;
    for(int l=1; l<=p ; l++){
      int ind = k-p+l ;
      T alpha = nS.V[k+l] - X[j] ;
      if(alpha==0.0)
        for(row=0;row<P.rows();row++)
          nS.P(row,ind-1) = nS.P(row,ind) ;
      else
        alpha /= nS.V[k+l]-V[i-p+l] ;
      for(row=0;row<P.rows();row++)
        nS.P(row,ind-1) = alpha*nS.P(row,ind-1) + (1.0-alpha)*nS.P(row,ind) ;
    }
    nS.V[k] = X[j] ;
    --k ;
  }
  *this = nS ;
}


template <class T, int N>
void NurbsSurface<T,N>::mergeKnots(const Vector<T>& nU, const Vector<T>& nV) {
  mergeKnotU(nU) ;
  mergeKnotV(nV) ;
}

template <class T, int N>
void NurbsSurface<T,N>::mergeKnotU(const Vector<T>& X){
  int i,ia,ib ;
  // Find the knots to insert
  Vector<T> I(U.n()) ;
  int done = 0 ;
  i = ia = ib = 0 ;
  while(!done) {
    if(X[ib] == U[ia]){
      ++ib ; ++ia ;
    }
    else{
      I[i++] = X[ib] ;
      ib++ ;
    }
    done = (ia>=U.n() || ib >= X.n()) ;
  }
  I.resize(i) ;

  if(I.n()>0)
    refineKnotU(I) ;
}

template <class T, int N>
void NurbsSurface<T,N>::mergeKnotV(const Vector<T>& X){
  int i,ia,ib ;
  // Find the knots to insert
  Vector<T> I(V.n()) ;
  int done = 0 ;
  i = ia = ib = 0 ;
  while(!done) {
    if(X[ib] == V[ia]){
      ++ib ; ++ia ;
    }
    else{
      I[i++] = X[ib] ;
      ib++ ;
    }
    done = (ia>=V.n() || ib >= X.n()) ;
  }
  I.resize(i) ;

  if(I.n()>0)
    refineKnotV(I) ;
}

template <class T, int N>
int NurbsSurface<T,N>::read(ifstream &fin){
  if(!fin) {
    return 0 ;
  }
  int nu,nv,du,dv;
  char *type ;
  type = new char[3] ;
  if(!fin.read(type,sizeof(char)*3)) { delete []type ; return 0 ;}
  int r1 = strncmp(type,"ns3",3) ;
  int r2 = strncmp(type,"ns4",3) ;
  if(!(r1 || r2)) 
    return 0 ;
  int st ;
  char stc ;
  if(!fin.read((char*)&stc,sizeof(char))) { delete []type ; return 0 ;}
  if(!fin.read((char*)&nu,sizeof(int))) { delete []type ; return 0 ;}
  if(!fin.read((char*)&nv,sizeof(int))) { delete []type ; return 0 ;}
  if(!fin.read((char*)&du,sizeof(int))) { delete []type ; return 0 ;}
  if(!fin.read((char*)&dv,sizeof(int))) { delete []type ; return 0 ;}

  st = stc - '0' ; 
  if(st != sizeof(T)){ // not of the same type size
    delete []type ;
    return 0 ; 
  }

  resize(nu,nv,du,dv) ;
  
  if(!fin.read((char*)U.memory(),sizeof(T)*U.n())) { delete []type ; return 0 ;}
  if(!fin.read((char*)V.memory(),sizeof(T)*V.n())) { delete []type ; return 0 ;}
     

  T *p,*p2 ;
  if(!r1){
    p = new T[3*nu*nv] ;
    if(!fin.read((char*)p,sizeof(T)*3*nu*nv)) { delete []type ; return 0 ;}
    p2 = p ;
    for(int i=0;i<nu;i++)
      for(int j=0;j<nv;j++){
        P(i,j).x() = *(p++) ;
        P(i,j).y() = *(p++) ;
        P(i,j).z() = *(p++) ;
        P(i,j).w() = 1.0 ;
      }
    delete []p2 ;
  }
  else{
    p = new T[4*nu*nv] ;
    if(!fin.read((char*)p,sizeof(T)*4*nu*nv)) { delete []type ; return 0 ;}
    p2 = p ;
    for(int i=0;i<nu;i++)
      for(int j=0;j<nv;j++){
        P(i,j).x() = *(p++) ;
        P(i,j).y() = *(p++) ;
        P(i,j).z() = *(p++) ;
        P(i,j).w() = *(p++) ;
      }
    delete []p2 ;
  }

  delete []type ;
  return 1 ;
}

template <class T, int N>
int NurbsSurface<T,N>::read(const char* filename){
  ifstream fin(filename) ;
  if(!fin) {
    return 0 ;
  }

  return read(fin) ;
}

template <class T, int N>
int NurbsSurface<T,N>::write(ofstream &fout) const {
  if(!fout)
    return 0 ;
  int prows = P.rows();
  int pcols = P.cols();
  char st = '0' + sizeof(T) ;
  if(!fout.write((char*)&"ns4",sizeof(char)*3)) return 0 ;
  if(!fout.write((char*)&st,sizeof(char))) return 0 ; 
  if(!fout.write((char*)&prows,sizeof(int))) return 0 ;
  if(!fout.write((char*)&pcols,sizeof(int))) return 0 ;
  if(!fout.write((char*)&degU,sizeof(int))) return 0 ;
  if(!fout.write((char*)&degV,sizeof(int))) return 0 ;
  if(!fout.write((char*)U.memory(),sizeof(T)*U.n())) return 0 ;
  if(!fout.write((char*)V.memory(),sizeof(T)*V.n())) return 0 ;
  
  T *p,*p2 ;
  p = new T[P.rows()*P.cols()*4] ;
  p2 = p ;
  for(int i=0;i<P.rows();i++) 
    for(int j=0;j<P.cols();j++){
      *p = P(i,j).x() ; p++ ;
      *p = P(i,j).y() ; p++ ;
      *p = P(i,j).z() ; p++ ;
      *p = P(i,j).w() ; p++ ;
    }
  if(!fout.write((char*)p2,sizeof(T)*P.rows()*P.cols()*4)) return 0 ;
  delete []p2 ;
  return 1 ;
}

template <class T, int N>
int NurbsSurface<T,N>::write(const char* filename) const {
  ofstream fout(filename) ;    
  if(!fout)
    return 0 ;
  return write(fout);
}

template <class T, int N>
NurbsSurface<T,N>& NurbsSurface<T,N>::transpose(void){
  Vector<T> t(U) ;
  int d ;
  U = V ;
  V = t ;
  d = degU ;
  degU = degV ;
  degV = d ;
  P = ::transpose(P) ;
  return *this ;
}

template <class T, int N>
int NurbsSurface<T,N>::movePoint(T u, T v, const Point_nD<T,N>& delta){
  // setup B
  Matrix_DOUBLE B(1,(degU+1)*(degV+1)) ;
  int i,j,k ;
  
  int spanU,spanV ;

  Vector<T> Ru,Rv ;

  B.reset(0.0) ;
  
  findSpan(u,v,spanU,spanV) ;
  nurbsBasisFuns(u,spanU,degU,U,Ru) ;
  nurbsBasisFuns(v,spanV,degV,V,Rv) ;
  for(j=0;j<=degU;++j){
    for(k=0;k<=degV;++k){
      B(0,j*(degV+1)+k) 
        = (double)Ru[j]*(double)Rv[k] ;
    }
  }
  
  Matrix_DOUBLE A  ;
  Matrix_DOUBLE Bt(::transpose(B)) ;
  Matrix_DOUBLE BBt ;

  BBt = inverse(B*Bt) ;
  A = Bt*BBt ;

  Matrix_DOUBLE dD(1,3) ;

  for(j=0;j<3;++j)
    dD(0,j) = (double)delta.data[j] ;

  Matrix_DOUBLE dP ;

  dP = A*dD ;

  i= 0 ;

  for(j=0;j<=degU;++j){
    for(k=0;k<=degV;++k){
      T w = P(spanU-degU+j,spanV-degV+k).w() ;
      P(spanU-degU+j,spanV-degV+k).x() += dP(i,0)*w ;
      P(spanU-degU+j,spanV-degV+k).y() += dP(i,1)*w ;
      P(spanU-degU+j,spanV-degV+k).z() += dP(i,2)*w ;
      ++i ;
    }
  }

  return 1 ;
}

template <class T, int N>
int NurbsSurface<T,N>::movePoint(const Vector<T>& ur, const Vector<T>& vr, const Vector< Point_nD<T,N> >& D, const Vector_INT& Du, const Vector_INT& Dv) {
  Vector_INT Dk(Du.n()),Dl(Dv.n()) ;
  BasicArray<Coordinate> fixCP(0) ;
  Dk.reset(0) ;
  Dl.reset(0) ;
  return movePoint(ur,vr,D,Du,Dv,Dk,Dl,fixCP) ;
}

template <class T, int N>
int NurbsSurface<T,N>::movePoint(const Vector<T>& ur, const Vector<T>& vr, const Vector< Point_nD<T,N> >& D, const Vector_INT& Du, const Vector_INT& Dv, const Vector_INT& Dk, const Vector_INT& Dl) {
  BasicArray<Coordinate> fixCP(0) ;
  return movePoint(ur,vr,D,Du,Dv,Dk,Dl,fixCP) ;
}



template <class T, int N>
int NurbsSurface<T,N>::movePoint(const Vector<T>& ur, const Vector<T>& vr, const Vector< Point_nD<T,N> >& D, const Vector_INT& Du, const Vector_INT& Dv, const Vector_INT& Dk, const Vector_INT& Dl, const BasicArray<Coordinate>& fixCP) {
  int i,j,k,n ;

  if(D.n() != Du.n() || D.n() !=Du.n() || D.n() != Dv.n() || D.n() != Dv.n()){
#ifdef USE_EXCEPTION
    throw NurbsInputError();
#else
    Error err("movePoint(ur,D,Dr,Dk,fixCP)");
    err << "The D,Dr,Dk vectors are not of the same size\n" ;
    err << "D.n()= " << D.n() << ", Du.n() = " << Du.n() 
        << ", Dk.n() = " << Dk.n() << ", Dv.n() = " << Dv.n() 
        << ", Dl.n() = " << Dl.n() << endl ;
    err.fatal() ;
#endif
  }  
  
  // setup B
  Matrix_DOUBLE B ;
  
  B.resize(D.n(),P.rows()*P.cols()) ;
  
  int spanU,spanV ;

  Matrix<T> Ru,Rv ;

  B.reset(0.0) ;
  
  for(i=0;i<D.rows();++i){
    findSpan(ur[Du[i]],vr[Dv[i]],spanU,spanV) ;
    nurbsDersBasisFuns(Dk[i],ur[Du[i]],spanU,degU,U,Ru) ;
    nurbsDersBasisFuns(Dl[i],vr[Dv[i]],spanV,degV,V,Rv) ;
    for(j=0;j<=degU;++j){
      for(k=0;k<=degV;++k){
        B(i,(spanU-degU+j)*P.cols()+spanV-degV+k) 
          = (double)Ru(Dk[i],j)*(double)Rv(Dl[i],k) ;
      }
    }
  }
  
  // optimize B
  Vector_INT remove(B.cols()) ;
  BasicArray<Coordinate> map(B.cols()) ;
  remove.reset((int)1.0) ;

  for(j=0;j<B.cols();++j){
    for(i=0;i<B.rows();++i)
      if((B(i,j)*B(i,j))>1e-10){
        remove[j] = 0 ;
        break ;
      }
  }

  for(i=0;i<fixCP.n();++i){
    remove[fixCP[i].i*P.cols()+fixCP[i].j] = 1 ;
  }
  
  n = 0 ;
  for(i=0;i<B.cols();++i){
    if(!remove[i]){
      map[n].i = i/P.cols() ;
      map[n].j = i%P.cols() ;
      ++n ;
    }
  }

  map.resize(n) ;

  Matrix_DOUBLE Bopt(B.rows(),n) ;
  for(j=0;j<n;++j){
    for(i=0;i<B.rows();++i)
      Bopt(i,j) = B(i,map[j].i*P.cols()+map[j].j) ;
  }

  Matrix_DOUBLE A  ;
  Matrix_DOUBLE Bt(::transpose(Bopt)) ;
  Matrix_DOUBLE BBt ;

  BBt = inverse(Bopt*Bt) ;

  A = Bt*BBt ;

  Matrix_DOUBLE dD(D.rows(),3) ;

  for(i=0;i<D.rows();++i){
    const Point_nD<T,N>& d = D[i] ; // this makes the SGI compiler happy
    for(j=0;j<3;++j)
      dD(i,j) = (double)d.data[j] ;
  }

  Matrix_DOUBLE dP ;

  dP = A*dD ;

  for(i=0;i<map.n();++i){
    P(map[i].i,map[i].j).x() += dP(i,0)*P(map[i].i,map[i].j).w() ;
    P(map[i].i,map[i].j).y() += dP(i,1)*P(map[i].i,map[i].j).w() ;
    P(map[i].i,map[i].j).z() += dP(i,2)*P(map[i].i,map[i].j).w() ;
  }

  return 1 ;
}


template <class T, int N>
void projectToLine(const Point_nD<T,N>& S, const Point_nD<T,N>& Trj, const Point_nD<T,N>& pnt, Point_nD<T,N>& p) {

  Point_nD<T,N> a = pnt-S ;  
  //p = S+ norm(a)*cos(angle(a,Trj))*Trj.unitLength() ;
  T fraction, denom ;
  denom = norm2(Trj) ; 
  fraction = (denom == 0.0) ? 0.0 : (Trj*a) / denom ; 
  p =  fraction * Trj ; 
  p += S ; 
}

template <class T, int N>
void NurbsSurface<T,N>::makeFromRevolution(const NurbsCurve<T,N>& profile, const Point_nD<T,N>& S, const Point_nD<T,N>& Tvec, double theta){
  double angle,dtheta ;
  int narcs ;
  int i,j ;

  //while(theta>2.0*M_PI)
  //  theta -= 2.0*M_PI ;

  if(theta>2.0*M_PI)
    theta = 2.0*M_PI ;

  if(theta <= 0)
    theta = 0 ;

  if(theta<M_PI/2.0)    
    narcs = 1 ;
  else{
    if(theta<M_PI)
      narcs = 2 ;
    else{
      if(theta<1.5*M_PI)
        narcs = 3 ;
      else
        narcs = 4 ;
    }
  }
  dtheta = theta/(double)narcs ;

  int n = 2*narcs+1 ; // n control points ;
  resize(n,profile.ctrlPnts().n(),2,profile.degree()) ;

  switch(narcs){
  case 1: break ;
  case 2: 
    U[3] = U[4] = 0.5 ;
    break ;
  case 3:
    U[3] = U[4] = 1.0/3.0 ;
    U[5] = U[6] = 2.0/3.0 ;
    break ;
  case 4:
    U[3] = U[4] = 0.25 ;
    U[5] = U[6] = 0.50 ;  
    U[7] = U[8] = 0.75 ;
    break ;    
  }

  j = 3 + 2*(narcs-1) ;  // loading the end knots
  for(i=0;i<3;++j,++i){
    U[i] = 0.0 ;
    U[j] = 1.0 ;
  }

  V = profile.knot() ;

  double wm = cos(dtheta/2.0) ; // dtheta/2.0 is base angle

  Vector_DOUBLE cosines(narcs+1),sines(narcs+1) ;
  angle = 0 ;
  for(i=1;i<=narcs;++i){
    angle = dtheta*(double)i ;
    cosines[i] = cos(angle) ;
    sines[i] = sin(angle) ;
  }

  Point_nD<T,N> P0,T0,P2,T2,P1 ;

  for(j=0;j<P.cols();++j){
    Point_nD<T,N> O ;
    T wj = profile.ctrlPnts(j).w() ;

    projectToLine(S,Tvec,project(profile.ctrlPnts(j)),O) ;
    //projectToLine(S,Tvec,profile.ctrlPnts(j).projectW(),O) ;
    Point_nD<T,N> X,Y ;
    X = project(profile.ctrlPnts(j))-O ;
    //X = profile.ctrlPnts(j).projectW() - O ; 

    double r = norm(X) ;

    if(r < 1e-7){
      for(i=0;i<P.rows();++i){
        P(i,j) = O ;
        P(i,j) *= wj ;
      }
      continue ;
    }

    T b1 = X.norm2() ;
    T b2 = X.norm() ;
    T b3 = sqrt(b1) ;
    T b4 = norm(X) ;
   
    X = X.unitLength() ;
    Y = crossProduct(Tvec,X) ;
    Y = Y.unitLength() ;

    P0 = project(profile.ctrlPnts(j)) ;
    //P0 = profile.ctrlPnts(j).projectW() ;
    P(0,j) = profile.ctrlPnts(j) ;

    T0 = Y ;
    int index = 0 ;

    for(i=1;i<=narcs;++i){
      angle = dtheta*(double)i ;
      P2 = O+ r*cosines[i]*X + r*sines[i]*Y ;  
      //P2 = O+ r*cos(angle)*X + r*sin(angle)*Y ;  
      P(index+2,j) = P2 ;
      P(index+2,j) *= wj ;
      T2 = -sines[i]*X + cosines[i]*Y ;
      //T2 = -sin(angle)*X + cos(angle)*Y ;
      intersectLine(P0,T0,P2,T2,P1) ;
      P(index+1,j) = P1 ;
      P(index+1,j) *= wm*wj ;
      index += 2 ;
      P0 = P2 ;
      T0 = T2 ;
    }
  }
}

template <class T, int N>
void NurbsSurface<T,N>::makeFromRevolution(const NurbsCurve<T,N>& profile, const Point_nD<T,N>& S, const Point_nD<T,N>& Tvec){
  double angle,dtheta ;
  int narcs ;
  int i,j ;

  int n = 9 ; // n control points ;
  resize(n,profile.ctrlPnts().n(),2,profile.degree()) ;

  U[0] = U[1] = U[2] = 0 ;
  U[3] = U[4] = 0.25 ;
  U[5] = U[6] = 0.50 ;  
  U[7] = U[8] = 0.75 ;
  U[9] = U[10] = U[11] = 1 ;

  V = profile.knot() ;


  const T wm = T(0.707106781185) ; // sqrt(2)/2

  for(j=0;j<P.cols();++j){
    Point_nD<T,N> O ;
    T wj = profile.ctrlPnts(j).w() ;

    projectToLine(S,Tvec,project(profile.ctrlPnts(j)),O) ;
    //projectToLine(S,Tvec,profile.ctrlPnts(j).projectW(),O) ;
    Point_nD<T,N> X,Y ;
    X = project(profile.ctrlPnts(j))-O ;
    //X = profile.ctrlPnts(j).projectW() - O ; 

    double r = norm(X) ;

    if(r < 1e-7){
      for(i=0;i<P.rows();++i){
        P(i,j) = O ;
        P(i,j) *= wj ;
      }
      continue ;
    }

    X = X.unitLength() ;
    Y = crossProduct(Tvec,X) ;
    Y = Y.unitLength() ;

    P(0,j) = O + r*X ;
    P(1,j) = O + r*X + r*Y ; 
    P(2,j) = O + r*Y ; 
    P(3,j) = O - r*X + r*Y ; 
    P(4,j) = O - r*X ;
    P(5,j) = O - r*X - r*Y ; 
    P(6,j) = O - r*Y ; 
    P(7,j) = O + r*X - r*Y ; 
    P(8,j) = P(0,j) ; 
    P(0,j) *= wj ;
    P(1,j) *= wj*wm ;
    P(2,j) *= wj ;
    P(3,j) *= wj*wm ;
    P(4,j) *= wj ;
    P(5,j) *= wj*wm ;
    P(6,j) *= wj ;
    P(7,j) *= wj*wm ;
    P(8,j) *= wj ;
  }
}

template <class T, int N>
void NurbsSurface<T,N>::makeFromRevolution(const NurbsCurve<T,N>& profile){
  int j ; 
  int n = 9 ; // n control points ;
  resize(n,profile.ctrlPnts().n(),2,profile.degree()) ;

  U[0] = U[1] = U[2] = 0 ;
  U[3] = U[4] = 0.25 ;
  U[5] = U[6] = 0.50 ;  
  U[7] = U[8] = 0.75 ;
  U[9] = U[10] = U[11] = 1 ;

  V = profile.knot() ;


  const T wm = T(0.707106781185) ; // sqrt(2)/2

  for(j=0;j<P.cols();++j){
    T wj = profile.ctrlPnts(j).w() ;

    Point_nD<T,N> p = project(profile.ctrlPnts(j)) ;
    T r = T(sqrt(p.x()*p.x()+p.y()*p.y()));

    T wjm = wj*wm ; 
    T rwjm = r*wj*wm ;
    T rwj = r*wj ; 
    p.z() *= wj ; 
    
    P(0,j) = HPoint_nD<T,N>(rwj,0,p.z(),wj) ;          
    P(1,j) = HPoint_nD<T,N>(rwjm,rwjm,p.z()*wm,wjm) ;   
    P(2,j) = HPoint_nD<T,N>(0,rwj,p.z(),wj) ;          
    P(3,j) = HPoint_nD<T,N>(-rwjm,rwjm,p.z()*wm,wjm) ;  
    P(4,j) = HPoint_nD<T,N>(-rwj,0,p.z(),wj) ;         
    P(5,j) = HPoint_nD<T,N>(-rwjm,-rwjm,p.z()*wm,wjm) ; 
    P(6,j) = HPoint_nD<T,N>(0,-rwj,p.z(),wj) ;         
    P(7,j) = HPoint_nD<T,N>(rwjm,-rwjm,p.z()*wm,wjm) ;  
    P(8,j) = HPoint_nD<T,N>(rwj,0,p.z(),wj) ;          
  }
}

template <class T, int D>
void NurbsSurface<T,D>::isoCurveU(T u, NurbsCurve<T,D>& c) const {
  c.resize(P.cols(),degV) ;

  c.modKnot(V) ;

  if(u>U[U.n()-1])
    u = U[U.n()-1] ;
  if(u<U[0])
    u = U[0] ;



  int span = findSpanU(u) ;
  
  Vector<T> N ;
  basisFunsU(u,span,N) ;

  HPoint_nD<T,D> p ;
  for(int i=0;i<P.cols();++i){
    p = 0 ;
    for(int j=0;j<=degU;++j)
      p += N[j]*P(span-degU+j,i) ;
    c.modCP(i,p) ;
  }
  
}

template <class T, int D>
void NurbsSurface<T,D>::isoCurveV(T v, NurbsCurve<T,D>& c) const {
  c.resize(P.rows(),degU) ;

  c.modKnot(U) ;

  if(v>V[V.n()-1])
    v = V[V.n()-1] ;
  if(v<V[0])
    v = V[0] ;

  int span = findSpanV(v) ;
  
  Vector<T> N ;
  basisFunsV(v,span,N) ;
  
  HPoint_nD<T,D> p ;
  for(int i=0;i<P.rows();++i){
    p = 0  ;
    for(int j=0;j<=degV;++j)
      p += N[j]*P(i,span-degV+j) ;
    c.modCP(i,p) ;
  }
}

template <class T, int N>
int NurbsSurface<T,N>::decompose(NurbsSurfaceArray<T,N>& S) const {
  int i,m,a,b,nb,mult,j,r,save,s,k,row,col ;
  T numer,alpha ;


  //Vector<T> alphas(degU+1) ;
  T* alphas = (T*) alloca((maximum(degU,degV)+1)*sizeof(T)) ;
  // all the surfaces will have the same knot vector in both the U and V
  // direction
  Vector<T> nU ;
  nU.resize(2*(degU+1)) ;
  for(i=0;i<nU.n()/2;++i)
    nU[i] = 0 ;
  for(i=nU.n()/2;i<nU.n();++i)
    nU[i] = 1 ;
  
  Vector<T> nV ;
  nV.resize(2*(degV+1)) ;
  for(i=0;i<nV.n()/2;++i)
    nV[i] = 0 ;
  for(i=nV.n()/2;i<nV.n();++i)
    nV[i] = 1 ;

  NurbsSurfaceArray<T,N> Su ;
  Su.resize(P.rows()-degU) ; // )*(P.cols()-degV)) ;
  for(i=0;i<Su.n();i++){
    Su[i].resize(degU+1,P.cols(),degU,degV) ;
    Su[i].U = nU ;
    Su[i].V = V ;
  }
  
  m = P.rows()+degU ;
  a = degU ;
  b = degU+1 ;
  nb = 0 ;

  for(i=0;i<=degU;i++)
    for(col=0;col<P.cols();col++)
      Su[nb].P(i,col) = P(i,col) ;
  while(b<m){
    i = b ;
    while(b<m && U[b+1] <= U[b]) b++ ;
    mult = b-i+1 ;
    if(mult<degU){
      numer = U[b]-U[a] ; // the enumerator of the alphas
      for(j=degU;j>mult;j--) // compute and store the alphas
        alphas[j-mult-1] = numer/(U[a+j]-U[a]) ;
      r = degU-mult ; // insert knot r times
      for(j=1;j<=r;j++){
        save=r-j;
        s=mult+j; // this many new points
        for(k=degU;k>=s;k--){
          alpha = alphas[k-s] ;
          for(col=0;col<P.cols();++col)
            Su[nb].P(k,col) = alpha*Su[nb].P(k,col)+(1.0-alpha)*Su[nb].P(k-1,col);
        }
        if(b<m) // control point of next patch
          for(col=0;col<P.cols();++col)
            Su[nb+1].P(save,col) = Su[nb].P(degU,col) ;
      }
    }
    ++nb ;
    if(b<m){ // initialize for next segment
      for(i=degU-mult; i<= degU ; ++i)
        for(col=0;col<P.cols();++col)
          Su[nb].P(i,col) = P(b-degU+i,col) ;
      a=b ;
      ++b ;
    }
  }
  Su.resize(nb) ;
  
  S.resize(Su.n()*(P.cols()-degV)) ;

  for(i=0;i<S.n();i++){
    S[i].resize(degU+1,degV+1,degU,degV) ;
    S[i].U = nU ;
    S[i].V = nV ;
  }
  
  nb = 0 ;

  for(int np=0;np<Su.n();++np){
    for(i=0;i<=degU;i++)
      for(j=0;j<=degV;++j)
        S[nb].P(i,j) = Su[np].P(i,j) ;
    m = P.cols()+degV ;
    a = degV ;
    b = degV+1 ;
    while(b<m){
      i = b ;
      while(b<m && V[b+1] <= V[b]) b++ ;
      mult = b-i+1 ;
      if(mult<degV){
        numer = V[b]-V[a] ; // the enumerator of the alphas
        for(j=degV;j>mult;j--) // compute and store the alphas
          alphas[j-mult-1] = numer/(V[a+j]-V[a]) ;
        r = degV-mult ; // insert knot r times
        for(j=1;j<=r;j++){
          save=r-j;
          s=mult+j; // this many new points
          for(k=degV;k>=s;k--){
            alpha = alphas[k-s] ;
            for(row=0;row<=degU;++row)
              S[nb].P(row,k) = alpha*S[nb].P(row,k)+(1.0-alpha)*S[nb].P(row,k-1);
          }
          if(b<m) // control point of next patch
            for(row=0;row<=degU;++row)
              S[nb+1].P(row,save) = S[nb].P(row,degV) ;
        }
      }
      ++nb ;
      if(b<m){ // initialize for next patch
        for(i=degV-mult; i<= degV ; ++i)
          for(row=0;row<=degU;++row)
            S[nb].P(row,i) = Su[np].P(row,b-degV+i) ;
        a=b ;
        ++b ;
      }
    }
  }

  S.resize(nb) ;

  return Su.n() ;
}

template <class T, int N>
int NurbsSurface<T,N>::writePOVRAY(ostream& povray, int patch_type, double flatness, int num_u_steps, int num_v_steps) const {
  if(degU>3 || degV>3){
#ifdef USE_EXCEPTION
    throw NurbsInputError();
#else
    Error err("NurbsSurface<T,N> writePOVRAY") ;
    err << "The degree of the surface is higher than 3!\n"
        << "A povrary file can not be generated.\n" ;
    err.warning() ;
    return 0;
#endif
  }

  NurbsSurface<T,N> S(*this) ;

  S.degreeElevate(3-degU,3-degV) ;

  NurbsSurfaceArray<T,N> Sa ;
  S.decompose(Sa) ;
  int bad = 0 ; 

  povray << "//\n//Generated for POV-Ray(tm) 3.0 by Phil's NURBS library\n" ;
  povray << "//   http://yukon.genie.uottawa.ca/info/soft/nurbs\n//\n" ;

  for(int i=0;i<Sa.n();++i){
    povray << "bicubic_patch {\n\ttype " << patch_type << endl ;
    povray << "\tflatness " << flatness << endl ;
    povray << "\tu_steps " << num_u_steps << endl ;
    povray << "\tv_steps " << num_v_steps << endl ;
    for(int j=0;j<4;++j){
      for(int k=0;k<4;++k){
        Point_nD<T,N> p = project(Sa[i].ctrlPnts(j,k)) ;
        if(Sa[i].ctrlPnts(j,k).w()>1.0001 || Sa[i].ctrlPnts(j,k).w()<0.9999)
          bad = 1 ;
        povray << "\t<" << p.x() << ", " << p.y() << ", " << p.z() << ">" ;
        if(j==3 && k==3)
          povray << "\n}\n\n" ;
        else
          povray << ",\n " ;
      }
      povray << endl ;
    }
  }
  if(bad){
#ifdef USE_EXCEPTION
    throw NurbsWarning();
#else
    Error err("NurbsSurface<T,N> writePOVRAY") ;
    err << "Warning: at least one of the control point was not rational\n" ;
    err << "The resulting surface will NOT be the same as the one which\n" ;
    err << "generated it.\n" ;
    err.warning() ;
#endif
  }
  return bad ;
}


template <class T, int N>
int NurbsSurface<T,N>::writePOVRAY(T tolerance, ostream& povray,const Color& color, int smooth, T ambient, T diffuse) const {
  BasicList<Point_nD<T,N> > points ;
  BasicList<int> connect ;
  BasicList<Point_nD<T,N> > norm ;

  if(smooth)
    tesselate(tolerance,points,connect,&norm) ;
  else
    tesselate(tolerance,points,connect) ;
    
  povray << "mesh {\n" ;
  
  BasicNode<int> *node ;

  BasicNode<Point_nD<T,N> > *nodeP ;
  BasicNode<Point_nD<T,N> > *nodeN ;
  nodeP = points.goToFirst() ;
  nodeN = 0 ;

  if(smooth)
    nodeN = norm.goToFirst() ;

  Vector< Point_nD<T,N> > Pts(points.size()) ;
  Vector< Point_nD<T,N> > Norm(norm.size()) ;
  for(int i=0;i<Pts.n();++i){
    Pts[i] = *nodeP->data ;
    nodeP = points.goToNext() ;
    if(smooth){
      Norm[i] = *nodeN->data ;
      nodeN = norm.goToNext() ;
    }
  }

  node = connect.goToFirst() ;
  while(node){
    if(smooth)
      povray << "\tsmooth_triangle { " ;
    else
      povray << "\ttriangle { " ;
    povray << "< " << Pts[*node->data].x() << ", " << Pts[*node->data].y() << ", " << Pts[*node->data].z() << "> , " ;
    if(smooth)
      povray << "< " << Norm[*node->data].x() << ", " << Norm[*node->data].y() << ", " << Norm[*node->data].z() << "> , " ;
    node = connect.goToNext() ; 
    povray << "< " << Pts[*node->data].x() << ", " << Pts[*node->data].y() << ", " << Pts[*node->data].z() << "> , " ;
    if(smooth)
      povray << "< " << Norm[*node->data].x() << ", " << Norm[*node->data].y() << ", " << Norm[*node->data].z() << "> , " ;
    node = connect.goToNext() ;
    povray << "< " << Pts[*node->data].x() << ", " << Pts[*node->data].y() << ", " << Pts[*node->data].z() << "> " ;
    if(smooth)
      povray << ", < " << Norm[*node->data].x() << ", " << Norm[*node->data].y() << ", " << Norm[*node->data].z() << "> " ;
    node = connect.goToNext() ; // skip the -1 value
    node = connect.goToNext() ;
    povray << "}\n" ;
  }
  
  povray << "\ttexture{ pigment { rgb < " << (double)color.r/255.0 << ", " << (double)color.g/255.0 << ", " << (double)color.b/255.0 << " >}\n" ;
  povray << "\t\tfinish { ambient " << ambient << " diffuse " << diffuse << " }\n\t}\n" ;
  
  povray << "}\n\n" ;
  return povray.good() ;
}

template <class T, int N>
int NurbsSurface<T,N>::writePOVRAY(const char *filename, const Color& col, const Point_nD<T,N>& cView, const Point_nD<T,N>& up, int patch_type, double flatness, int num_u_steps, int num_v_steps) const {
  ofstream fout(filename) ;
  if(!fout)
    return 0 ;
  
  Point_nD<T,N> view(T(-1.0)*cView) ;

  fout << "//\n//Generated for POV-Ray(tm) 3.0 by Phil's NURBS library\n//\n" ;
  fout << "\n#include \"colors.inc\"\n" ;

  // we want the camera to look at the center of the object
  // and be able to view the whole object
  // we use and angle of 36 to view the object
  // and position the rest according to this.
  Point_nD<T,N> minP, maxP ;
  minP.x() = extremum(1,coordX) ;
  minP.y() = extremum(1,coordY) ;
  minP.z() = extremum(1,coordZ) ;
  maxP.x() = extremum(0,coordX) ;
  maxP.y() = extremum(0,coordY) ;
  maxP.z() = extremum(0,coordZ) ;

  Point_nD<T,N> lookAt  ;
  lookAt.x() = (minP.x()+maxP.x())/2.0 ;
  lookAt.y() = (minP.y()+maxP.y())/2.0 ;
  lookAt.z() = (minP.z()+maxP.z())/2.0 ;

  Point_nD<T,N> camera1, camera2 ;

  Point_nD<T,N> q1 = minP-lookAt ;
  Point_nD<T,N> q2 = maxP-lookAt ;
  T D1 = absolute(dot(q1,view))/norm(view) ;
  T D2 = absolute(dot(q2,view))/norm(view) ;

  T a1 = norm(q1)*cos(angle(view,q1)) ;
  T a2 = norm(q2)*cos(angle(view,q2)) ;

  T b1 = D1/tan(18.0*M_PI/180.0) ;
  T b2 = D2/tan(18.0*M_PI/180.0) ; // this gives the 36 degree angle

  camera1 = lookAt+(a1+b1)*view.unitLength() ;
  camera2 = lookAt+(a2+b2)*view.unitLength() ;

  Point_nD<T,N> right ;

  right = crossProduct(view,up) ; // inversed because pov-ray uses a left-handed system

  fout << "camera {\n\tlocation <" ;
  if(norm2(camera1-lookAt)>norm2(camera2-lookAt))
    fout << camera1.x() << ", " << camera1.y() << ", " << camera1.z() << ">\n" ;
  else
    fout << camera2.x() << ", " << camera2.y() << ", " << camera2.z() << ">\n" ;
  fout << "\tup < " << up.x() << ", " << up.y() << ", " << up.z() << ">\n" ;
  fout << "\tright < " << right.x() << ", " << right.y() << ", " << right.z() << ">\n" ;
  fout << "\tlook_at < " << lookAt.x() << ", " << lookAt.y() << ", " << lookAt.z() << ">\n\tangle 36\n}\n\n" ;

  fout << "union {\n" ;
  writePOVRAY(fout,patch_type,flatness,num_u_steps,num_v_steps) ;
  fout << " texture {\n\tpigment {\n\t\tcolor rgb < " << (double)col.r/255.0 << 
    ", " << (double)col.g/255.0 << ", " << (double)col.b/255.0 << "> \n" <<
    "\t}\n\tfinish { \n\t\tambient .2\n\t\tdiffuse .6\n\t}\n }\n" ;
  fout << "\n}\n" ;

  fout << "light_source { < " ;
  if(norm2(camera1-lookAt)>norm2(camera2-lookAt))
    fout << camera1.x()+view.x() << ", " << camera1.y()+view.y() << ", " << camera1.z()+view.z() << "> color White}\n\n" ;
  else
    fout << camera2.x()+view.x() << ", " << camera2.y()+view.y() << ", " << camera2.z()+view.z() << "> color White}\n\n" ;
  
  return fout.good() ;
}


template <class T, int N>
int NurbsSurface<T,N>::writePOVRAY(T tolerance, const char *filename, const Color& col, const Point_nD<T,N>& cView, const Point_nD<T,N>& up, int smooth, T ambient, T diffuse) const {
  ofstream fout(filename) ;
  if(!fout)
    return 0 ;
  
  Point_nD<T,N> view(T(-1)*cView) ;

  fout << "//\n//Generated for POV-Ray(tm) 3.0 by Phil's NURBS library\n//\n" ;
  fout << "\n#include \"colors.inc\"\n" ;

  // we want the camera to look at the center of the object
  // and be able to view the whole object
  // we use and angle of 36 to view the object
  // and position the rest according to this.
  Point_nD<T,N> minP, maxP ;
  minP.x() = extremum(1,coordX) ;
  minP.y() = extremum(1,coordY) ;
  minP.z() = extremum(1,coordZ) ;
  maxP.x() = extremum(0,coordX) ;
  maxP.y() = extremum(0,coordY) ;
  maxP.z() = extremum(0,coordZ) ;

  Point_nD<T,N> lookAt  ;
  lookAt.x() = (minP.x()+maxP.x())/2.0 ;
  lookAt.y() = (minP.y()+maxP.y())/2.0 ;
  lookAt.z() = (minP.z()+maxP.z())/2.0 ;

  Point_nD<T,N> camera1, camera2 ;

  Point_nD<T,N> q1 = minP-lookAt ;
  Point_nD<T,N> q2 = maxP-lookAt ;
  T D1 = absolute(dot(q1,view))/norm(view) ;
  T D2 = absolute(dot(q2,view))/norm(view) ;

  T a1 = norm(q1)*cos(angle(view,q1)) ;
  T a2 = norm(q2)*cos(angle(view,q2)) ;

  T b1 = D1/tan(18.0*M_PI/180.0) ;
  T b2 = D2/tan(18.0*M_PI/180.0) ; // this gives the 36 degree angle

  camera1 = lookAt+(a1+b1)*view.unitLength() ;
  camera2 = lookAt+(a2+b2)*view.unitLength() ;

  Point_nD<T,N> right ;

  right = crossProduct(view,up) ; // inversed because pov-ray uses a left-handed system

  fout << "camera {\n\tlocation <" ;
  if(norm2(camera1-lookAt)>norm2(camera2-lookAt))
    fout << camera1.x() << ", " << camera1.y() << ", " << camera1.z() << ">\n" ;
  else
    fout << camera2.x() << ", " << camera2.y() << ", " << camera2.z() << ">\n" ;
  fout << "\tup < " << up.x() << ", " << up.y() << ", " << up.z() << ">\n" ;
  fout << "\tright < " << right.x() << ", " << right.y() << ", " << right.z() << ">\n" ;
  fout << "\tlook_at < " << lookAt.x() << ", " << lookAt.y() << ", " << lookAt.z() << ">\n\tangle 36\n}\n\n" ;

  writePOVRAY(tolerance,fout,col,smooth,ambient,diffuse) ;
  

  fout << "light_source { < " ;
  if(norm2(camera1-lookAt)>norm2(camera2-lookAt))
    fout << camera1.x()+view.x() << ", " << camera1.y()+view.y() << ", " << camera1.z()+view.z() << "> color White}\n\n" ;
  else
    fout << camera2.x()+view.x() << ", " << camera2.y()+view.y() << ", " << camera2.z()+view.z() << "> color White}\n\n" ;
  
  return fout.good() ;
}




template <class T, int N>
int NurbsSurface<T,N>::writeRIB(ostream& rib) const {
  rib << "NuPatch " << P.rows() << ' ' << degU+1 << " [ " ; 
  int k;
  // I have to loop since RIB is new line sensitive and 
  // by default a vector adds a new line at the end when using cout
  for(k=0;k<U.n();++k) 
    rib << U[k] << ' ' ;
  rib << " ] " << U[0] << ' ' << U[U.n()-1] << ' ' << P.cols() << ' ' << degV+1 << " [ " ;
  for(k=0;k<V.n();++k)
    rib << V[k] << ' ' ;
  rib << " ] " << V[0] << ' ' << V[V.n()-1] << " \"Pw\" [ " ;
  for(int j=0;j<P.cols();++j)
    for(int i=0;i<P.rows();++i)
      rib << P(i,j) ;
  rib << " ]\n" ;
  return rib.good() ;
}

template <class T, int N>
int NurbsSurface<T,N>::writeRIB(const char* filename, const Color& col, const Point_nD<T,N>& view) const {
  ofstream fout(filename) ;
  if(!fout)
    return 0;
  
  // The following code is based on Listing 8.2 from the RenderMan Companion
  // http://pete.cs.caltech.edu/RMR/Pixar/ch8/listing8_2.c
  int i,j ;

  // aimZ
  MatrixRT<double> Rx ;
  double xzlen, yzlen, yrot, xrot;

  //
  // The initial rotation about the y axis is given by the projection of
  // the direction vector onto the x,z plane: the x and z components
  // of the direction. 
  xzlen = sqrt(view.x()*view.x()+view.z()*view.z()) ;
  if(xzlen == 0)    
    yrot = (view.y() < 0) ? M_PI : 0;
  else
    yrot = acos(view.y()/xzlen);
  
  //
  // The second rotation, about the x axis, is given by the projection on
  // the y,z plane of the y-rotated direction vector: the original y
  // component, and the rotated x,z vector from above. 
  //
  yzlen = sqrt(view.y()*view.y()+xzlen*xzlen);
  xrot = acos(xzlen/yzlen);       /* yzlen should never be 0 */

  // A rotation around y first 
  if (view.y() > 0){
    if(view.x()>0)
      Rx.rotate(xrot,yrot,0.0) ;
    else
      Rx.rotate(xrot,-yrot,0.0) ;
  }
  else{
    if(view.x()>0)
      Rx.rotate(-xrot,yrot,0.0) ;
    else
      Rx.rotate(-xrot,-yrot,0.0) ;
  }

  Point_nD<T,N> minP, maxP ;
  minP.x() = extremum(1,coordX) ;
  minP.y() = extremum(1,coordY) ;
  minP.z() = extremum(1,coordZ) ;
  maxP.x() = extremum(0,coordX) ;
  maxP.y() = extremum(0,coordY) ;
  maxP.z() = extremum(0,coordZ) ;

  Point_nD<T,N> lookAt  ;
  lookAt.x() = (minP.x()+maxP.x())/2.0 ;
  lookAt.y() = (minP.y()+maxP.y())/2.0 ;
  lookAt.z() = (minP.z()+maxP.z())/2.0 ;

  Point_nD<T,N> camera1, camera2 ;

  Point_nD<T,N> q1 = minP-lookAt ;
  Point_nD<T,N> q2 = maxP-lookAt ;
  T D1 = absolute(dot(q1,view))/norm(view) ;
  T D2 = absolute(dot(q2,view))/norm(view) ;

  T a1 = norm(q1)*cos(angle(view,q1)) ;
  T a2 = norm(q2)*cos(angle(view,q2)) ;

  T b1 = D1/tan(18.0*M_PI/180.0) ;
  T b2 = D2/tan(18.0*M_PI/180.0) ; // this gives the 36 degree angle

  camera1 = lookAt+(a1+b1)*view.unitLength() ;
  camera2 = lookAt+(a2+b2)*view.unitLength() ;

  Point_nD<T,N> camera ;

  if(norm(camera1-lookAt)>norm(camera2-lookAt))
    camera = camera1 ;
  else
    camera = camera2 ;

  char front[1024] ;

  char *ext ;
  ext = strstr(filename,".rib") ;
  if(ext){
    for(i=0;i<1024;++i){
      if(&filename[i] == ext)
        break ;
      else
        front[i] = filename[i] ;
      if(!front[i])
        break ;
    }
  }
  else{
    strcpy(front,filename) ;
  }
  

  fout << "##RenderMan RIB-Structure 1.0\n" ;
  fout << "#" << filename << endl;
  fout << "Format 400 400 1\n";
  fout << "Display \"" << front << ".tif\" \"file\" \"rgba\"\n" ;
  fout << "Projection \"perspective\" \"fov\" [36]\n" ;
  fout << "Translate 0 0 " << norm(camera-lookAt) << endl ;
  fout << "Option \"render\" \"prmanspecular\" [1]\n" ;

  fout << "\nWorldBegin\n" ;
  fout << "LightSource \"ambientlight\" 0 \"intensity\" [0.3]\n" ;
  fout << "LightSource \"distantlight\" 1 \"to\" [ " << view.x() << ' ' << view.y() << ' ' << view.z()  << "]\n" ;
  //fout << "LightSource \"spotlight\" 1 \"intensity\" [300] \"from\" [ " ;
  //fout << camera.x()+view.x() << ' ' << camera.y()+view.y() << ' ' << camera.z()+view.z() << " ]\n\n" ;
  fout << "AttributeBegin\nSurface \"plastic\"\n" ;
  fout << "Color [ " << (double)col.r/255.0 << ' ' << (double)col.g/255.0 << ' ' << double(col.b)/255.0 << "]\n" ;

  fout << "Transform [ " ;

  for(j=0;j<4;++j)
    for(i=0;i<4;++i)
      fout << Rx(i,j) << ' ' ;
  fout << "]\n" ;
  fout << "Translate " << -lookAt.x() << ' ' << -lookAt.y() << ' ' << -lookAt.z() << endl ;


  writeRIB(fout) ;

  fout << "AttributeEnd\nWorldEnd\n\n" ;

  return fout.good() ;
}


template <class T, int N>
void NurbsSurface<T,N>::tesselate(T tolerance, BasicList<Point_nD<T,N> > &points, BasicList<int> &connect, BasicList<Point_nD<T,N> > *Norm) const {

#ifdef USE_EXCEPTION
  throw NurbsError() ;
#else
  Error err("NurbsSurface<T,N>::tesselate");
  err << "The tesselate member function is deprecated. Please use\n"
    "the NurbsSubSurface class member functions instead.\n" ; 
  err.fatal();
#endif
}


template <class T, int N>
void NurbsSurface<T,N>::makeSphere(const Point_nD<T,N>& O, T r) {
  NurbsCurve<T,N> c ;

  const T wm = T(0.707106781185) ;  // sqrt(2)/2

  c.resize(5,2) ; 

  c.modCP(0,HPoint_nD<T,N>(0,0,r,1)) ; 
  c.modCP(1,HPoint_nD<T,N>(-r*wm,0,r*wm,wm)) ; 
  c.modCP(2,HPoint_nD<T,N>(-r,0,0,1)) ; 
  c.modCP(3,HPoint_nD<T,N>(-r*wm,0,-r*wm,wm)) ; 
  c.modCP(4,HPoint_nD<T,N>(0,0,-r,1)) ; 
  
  Vector<T> k(5+2+1) ;
  k[0] = k[1] = k[2] = 0 ;
  k[3] = k[4] = 0.5 ;
  k[5] = k[6] = k[7] = 1 ; 
  
  c.modKnot(k) ; 

  makeFromRevolution(c) ; 
  MatrixRT<T> Tx ;
  Tx.translate(O.x(),O.y(),O.z()) ;
  transform(Tx) ;
}

template <class T, int N>
int NurbsSurface<T,N>::writePS(const char* filename, int nu, int nv, const Point_nD<T,N>& camera, const Point_nD<T,N>& lookAt, int cp,T magFact,T dash) const {
  NurbsCurveArray<T,N> Ca ;
  if(nu<=0 || nv<=0)
    return 0 ;


  // We need to find the rotation matrix to have z axis along nv
  Point_nD<T,N> np  = lookAt-camera ; 
  np = np.unitLength() ; 
  np *= -1 ;
  
  T rx = atan2(np.z(),np.y()) ;
  T ry = atan2(np.z(),np.x()) ;
  
  MatrixRT<T> Tx(rx,ry,0,camera.x(),camera.y(),camera.z()) ;
  //MatrixRT<T,N> Sc ; Sc.scale(1,1,T(norm(lookAt-camera))/plane) ;
  //MatrixRT<T,N> Tg(Sc*Tx) ; 

  Ca.resize(nu+nv+2) ; 
  int i ; 
  for(i=0;i<=nu;++i){
    T u = U[0]+T(i)*(U[U.n()-1]-U[0])/T(nu) ;
    isoCurveU(u , Ca[i]) ;
    Ca[i].transform(Tx) ; 
  }
  for(;i<=nv+nu+1;++i){
    T v = V[0]+T(i-nu-1)*(V[V.n()-1]-V[0])/T(nv) ;
    isoCurveV(v , Ca[i]) ;
    Ca[i].transform(Tx) ; 
  }


  return Ca.writePS(filename,cp,magFact,dash) ;
}

template <class T, int N>
int NurbsSurface<T,N>::writePSp(const char*, int nu, int nv, const Point_nD<T,N>& camera, const Point_nD<T,N>& lookAt, const Vector< Point_nD<T,N> >&,const Vector< Point_nD<T,N> >&, int cp,T magFact,T dash) const {
  cerr << "Not implemented. Not sure what is different between this and writePS\n";
  return 0;
}

template <class T, int N>
ostream& NurbsSurface<T,N>::print(ostream& o) const {
  o << "Degree: " << degU << ' ' << degV << endl;
  o << "U : " << U << endl;
  o << "V: " << V << endl ;
  o << "matrix size: " << P.rows() << ' ' << P.cols() << endl ;
  o << P << endl;
  return o;
}

template <class T, int N>
int surfMeshParamsClosedU(const Matrix< Point_nD<T,N> >& Q, Vector<T>& uk, Vector<T>& vl, int degU){

  int n,m,k,l,num ;
  double d,total ;
  Vector<T> cds(Q.rows()) ;

  n = Q.rows() ;
  m = Q.cols() ;
  uk.resize(n) ;
  vl.resize(m) ;
  num = m ;
  
  // Compute the uk
  uk.reset(0) ;

  for(l=0;l<m;l++){
    total = 0.0 ;
    for(k=1;k<=n-degU;k++){
      cds[k] = norm(Q(k,l)-Q(k-1,l)) ;
      total += cds[k] ;
    }
    for(k=n-degU+1; k<n ;k++)
      cds[k] = norm(Q(k,l)-Q(k-1,l)) ;

    if(total==0.0) 
      num-- ;
    else {
      d = 0.0 ;
      for(k=1;k<n;k++){
        d += cds[k] ;
        uk[k] += d/total ;
      }
    }
  }
  if(num==0) 
    return 0 ;
  for(k=1;k<n;k++)
    uk[k] /= num ;

  // Compute the vl
  vl.reset(0) ;
  cds.resize(m) ;

  num = n ;

  for(k=0;k<n;k++){
    total = 0.0 ;
    for(l=1;l<m;l++){
      cds[l] = norm(Q(k,l)-Q(k,l-1)) ;
      total += cds[l] ;
    }
    if(total==0.0) 
      num-- ;
    else {
      d = 0.0 ;
      for(l=1;l<m;l++){
        d += cds[l] ;
        vl[l] += d/total ;
      }
    }
  }
  if(num==0) 
    return 0 ;
  for(l=1;l<m-1;l++)
    vl[l] /= num ;
  vl[m-1] = 1.0 ;


  return 1 ;
}

template <class T, int N>
int surfMeshParamsClosedUH(const Matrix< HPoint_nD<T,N> >& Q, Vector<T>& uk, Vector<T>& vl, int degU){

  int n,m,k,l,num ;
  double d,total ;
  Vector<T> cds(Q.rows()) ;

  n = Q.rows() ;
  m = Q.cols() ;
  uk.resize(n) ;
  vl.resize(m) ;
  num = m ;
  
  // Compute the uk
  uk.reset(0) ;

  for(l=0;l<m;l++){
    total = 0.0 ;
    // Normalization factor
    for(k=1;k<=n-degU;k++){
      cds[k] = norm(Q(k,l)-Q(k-1,l)) ;
      total += cds[k] ;
    }
    for(k=n-degU+1; k<n ;k++)
      cds[k] = norm(Q(k,l)-Q(k-1,l)) ;
    if(total==0.0) 
      num-- ;
    else {
      d = 0.0 ;
      for(k=1;k<n;k++){
        d += cds[k] ;
        uk[k] += d/total ;
      }
    }
  }
  if(num==0) 
    return 0 ;
  for(k=1;k<n;k++)
    uk[k] /= num ;

  // Compute the vl
  vl.reset(0) ;
  cds.resize(m) ;

  num = n ;

  for(k=0;k<n;k++){
    total = 0.0 ;
    for(l=1;l<m;l++){
      cds[l] = norm(Q(k,l)-Q(k,l-1)) ;
      total += cds[l] ;
    }
    if(total==0.0) 
      num-- ;
    else {
      d = 0.0 ;
      for(l=1;l<m;l++){
        d += cds[l] ;
        vl[l] += d/total ;
      }
    }
  }
  if(num==0) 
    return 0 ;
  for(l=1;l<m-1;l++)
    vl[l] /= num ;
  vl[m-1] = 1.0 ;

  return 1 ;
}


template <class T, int N>
void NurbsSurface<T,N>::globalInterpClosedU(const Matrix< Point_nD<T,N> >& Q, int pU, int pV){
  Vector<T> vk,uk ;

  resize(Q.rows(),Q.cols(),pU,pV) ;

  surfMeshParamsClosedU(Q,uk,vk,pU) ;
  knotAveragingClosed(uk,pU,U) ;
  knotAveraging(vk,pV,V) ;


  Vector< HPoint_nD<T,N> > Pts(Q.cols()) ;
  NurbsCurve<T,N> R ;
  
  int i,j ;
  for(i=0;i<Q.rows();i++){
    for(j=0;j<Q.cols();j++)
      Pts[j] = Q(i,j) ;
    R.globalInterpH(Pts,vk,V,pV) ;
    for(j=0;j<Q.cols();j++)
      P(i,j) = R.ctrlPnts(j) ;
  }
  
  Pts.resize(Q.rows()) ;
  for(j=0;j<Q.cols();j++){
    for(i=0;i<Q.rows();i++)
      Pts[i] = P(i,j) ;
    
    R.globalInterpClosedH(Pts,uk,U,pU);
    for(i=0;i<Q.rows();i++)
      P(i,j) = R.ctrlPnts(i) ;
  }

}

template <class T, int N>
void NurbsSurface<T,N>::globalInterpClosedUH(const Matrix< HPoint_nD<T,N> >& Q, int pU, int pV){
  Vector<T> vk,uk ;

  resize(Q.rows(),Q.cols(),pU,pV) ;

  surfMeshParamsClosedUH(Q,uk,vk,pU) ;
  knotAveragingClosed(uk,pU,U) ;
  knotAveraging(vk,pV,V) ;
  
  Vector< HPoint_nD<T,N> > Pts(Q.rows()) ;
  NurbsCurve<T,N> R ;
  
  int i,j ;

  for(j=0;j<Q.cols();j++){
    for(i=0;i<Q.rows();i++)
      Pts[i] = Q(i,j) ;
    R.globalInterpH(Pts,uk,U,pU);
    for(i=0;i<Q.rows();i++)
      P(i,j) = R.ctrlPnts(i) ;
  }

  Pts.resize(Q.cols()) ;
  for(i=0;i<Q.rows();i++){
    for(j=0;j<Q.cols();j++)
      Pts[j] = P(i,j) ;
    R.globalInterpClosedH(Pts,vk,V,pV) ;
    for(j=0;j<Q.cols();j++)
      P(i,j) = R.ctrlPnts(j) ;
  }
}


template <class T, int N>
void NurbsSurface<T,N>::leastSquaresClosedU(const Matrix< Point_nD<T,N> >& Q, int pU, int pV, int nU, int nV){
  Vector<T> vk,uk ;

  resize(nU+pU,nV,pU,pV) ;

  surfMeshParamsClosedU(Q,uk,vk,pU) ;


  Vector< HPoint_nD<T,N> > Pts(Q.rows()) ;
  NurbsCurve<T,N> R ;
  
  int i,j ;

  Matrix< HPoint_nD<T,N> > P2 ;

  P2.resize(nU+pU,Q.cols()) ;

  for(j=0;j<Q.cols();j++){
    for(i=0;i<Q.rows();i++)
      Pts[i] = Q(i,j) ;
    R.leastSquaresClosedH(Pts,pU,nU,uk);
    for(i=0;i<P.rows();i++)
      P2(i,j) = R.ctrlPnts(i) ;
    if(j==0)
      U = R.knot() ;
  }

  Pts.resize(Q.cols()) ;
  for(i=0;i<P.rows();i++){
    for(j=0;j<Q.cols();j++)
      Pts[j] = P2(i,j) ;
    R.leastSquaresH(Pts,pV,nV,vk) ;
    for(j=0;j<P.cols();j++)
      P(i,j) = R.ctrlPnts(j) ;
    if(i==0)
      V = R.knot() ;
  }  

}

template <class T, int N>
int NurbsSurface<T,N>::writeOOGL(const char* filename,
                                 T fDu,
                                 T fDv,
                                 T fBu, T fBv,
                                 T fEu, T fEv,
                                 bool  bDrawCP) const {
  ofstream fout(filename) ;
  if(!fout)
    return 0 ;
  
  // Write file header
  fout << "{ LIST \n";
  fout << "\t{ appearance { shading smooth } \n " ;
  fout << "\tNMESH" << endl;
  T Nu = (fEu-fBu)/fDu + 1;
  T Nv = (fEv-fBv)/fDv + 1;
  fout << "\t"<< Nu << ' ' << Nv << endl;

  // Write mesh vertexes
  Point_nD<T,N> Sp, Np;
  T u,v;
  for (u = fBu; u<fEu+fDu/2; u+=fDu)
    for (v = fBv; v<fEv+fDv/2; v+=fDv){
      Sp = pointAt(u,v);
      Np = normal(u,v);
      Np = (norm(Np)!=0)?Np.unitLength():Point_nD<T,N>(0.0);
      fout << "\t" << Sp << "\t " << Np << endl;
    }
  fout << "\t}" << endl << std::flush;

  // Write the control points 
  if (bDrawCP){
    fout << "\t{ " << endl; 
    fout << "\t  appearance {shading smooth linewidth 5 } " << endl;

    fout << "\t" << "SKEL" << endl; 
    fout <<  P.rows()*P.cols() << ' ' << P.rows()*P.cols() << endl;
    for (int i = 0; i<P.rows(); i++)
      for (int j = 0; j<P.cols(); j++) 
        fout << "\t" << project(P(i,j)) << endl;    
    fout << endl;
    for (int i = 0; i<P.rows()*P.cols(); i++)
      fout << "\t" << "1 " << i << " 0 0 1 0.5 " << endl;
      fout << "\t" << " }" << endl;
  }
  fout << "} " << endl;
  fout << std::flush;

  return 1 ;
}

template <class T, int N>
int NurbsSurface<T,N>::writeDisplayQUADMESH(const char* filename, int iNu,int iNv,const Color& color,T fA, T fO) const 
{

  T fDu = 1/T(iNu);
  T fDv = 1/T(iNv-1);

  ofstream fout(filename) ;
  if(!fout)
    return 0 ;

  // Save the object type
  const char QUADMESH='q'+ ('A' - 'a');
  fout << QUADMESH << ' ';       ;

  // Compute surface properties
  T a, d, s, se, t;

  // Ambient reflectance coefficient
  a  = 0.3;
  // Diffusse reflectance coefficient
  d  = 0.3;
  // Specularity reflectance coeficcient
  s  = 0.4;
  // Specularity reflectance exponent
  se = 10;
  // Opacity
  t  = fO;

  // Save surface properties
  fout << a << ' ' 
       << d << ' ' 
       << s << ' ' 
       << se << ' ' 
       << t  << ' ';

  // Save mesh dimensions  
  fout << iNu << ' ' 
       << iNv << ' ' ;
  
  // Save wrapp status in each direction (v,u)  
  fout << "F T ";

  // New line
  fout << endl ;
  
  // Surface color RGBA (one color for the whole surface)
  T fR= T(color.r)/255;
  T fG= T(color.g)/255;
  T fB= T(color.b)/255;

  /* Colour flag = ONE_COLOUR */
  fout << 0 << ' ' ;
  /* The colour */
  fout << fR << ' '
       << fG << ' '
       << fB << ' '
       << fA << endl;
  
  // New line
  fout << endl ;

  // Now the list of 3D points
  T u,v;
  Point_nD<T,N> Sp;
  for (v = 0; v<1+fDv/2; v+=fDv)
    for (u = 0; u<1-fDu/2; u+=fDu){
      // The change in sign and the swap of y and z coordinates is
      // for conversion to MINC format.
      Sp = -(T)1.0 * pointAt(u,v) ;
      fout << Sp.x() << ' ' << Sp.z() << ' ' << Sp.y() << endl;
    }

  // New line
  fout << endl ;

  // Now the normal vectors
  Point_nD<T,N> Np;
  for (v = 0; v<1+fDv/2; v+=fDv)
    for (u = 0; u<1-fDu/2; u+=fDu){
      Np = normal(u,v);
      Np = (norm(Np)!=0)?Np.unitLength():Point_nD<T,N>(0.0);
      fout << Np.x() << ' ' << Np.z() << ' ' << Np.y() << endl;
    }

  // New line
  fout << endl ;


  return 1;
}

template <class T, int N>
int NurbsSurface<T,N>::writeOOGL(const char* filename) const {
  ofstream fout(filename) ;
  if(!fout)
    return 0 ;

  // Write file header
  int iPointDim = 4;
  fout << "BEZ" << degU << degV << iPointDim << endl;

  // Decompose surface in its Bezier Patches
  NurbsSurfaceArray<T,N> S;
  NurbsSurface<T,N>      surface(*this);
  surface.decompose(S);

  // Write patch vertexes
  for (int iPatch = 0; iPatch < S.n(); iPatch++){
    for(int iu = 0; iu < degU + 1; iu++){
      for(int iv = 0; iv < degV + 1; iv++)
        fout << S[iPatch].ctrlPnts(iu,iv).x() << ' ' 
             << S[iPatch].ctrlPnts(iu,iv).y() << ' '  
             << S[iPatch].ctrlPnts(iu,iv).z() << ' '
             << S[iPatch].ctrlPnts(iu,iv).w() << endl;
    }
    fout << endl;    
  }
  fout << std::flush;

  return 1 ;
}


template <class T, int N>
void wrapPointMatrix(const Matrix< Point_nD<T,N> >& Q, int d, int dir,  Matrix< Point_nD<T,N> >& Qw){
  int i, row, col;

  Qw = Q;

  if (dir==0){
    //    cout << " Wrapping in U dir " << endl << std::flush ;
    Qw.resizeKeep(Q.rows()+d,Q.cols());
    for (col=0; col < Q.cols(); col++)
      for (i=0; i<d; i++)
        Qw(i+Q.rows(),col) = Q(i,col);
  }
  else{
    //    cout << " Wrapping in V dir " << endl << std::flush ;
    Qw.resizeKeep(Q.rows(),Q.cols()+d);
    for (row=0; row < Q.rows(); row++)
      for (i=0; i<d; i++)
        Qw(row,i+Q.cols()) = Q(row,i);
  }
  
} 

template <class T, int N>
void NurbsSurface<T,N>::dersBasisFuns(T u, T v, int dU, int dV, int uspan, int vspan, Matrix<T> & Niku, Matrix<T> & Njkv ) const {
  // Get derivatives
  nurbsDersBasisFuns(dU,u,uspan,degU,U,Niku) ;
  nurbsDersBasisFuns(dV,v,vspan,degV,V,Njkv) ;
}

template <class T, int N>
void NurbsSurface<T,N>::makeTorus(const Point_nD<T,N>& O, T R, T r) {
  // These define the shape of a unit torus centered about the origin. 
  T xvalues[] = { 0.0, -1.0, -1.0, -1.0, 0.0, 1.0, 1.0, 1.0, 0.0 };
  T yvalues[] = { 1.0, 1.0, 0.0, -1.0, -1.0, -1.0, 0.0, 1.0, 1.0 };
  T zvalues[] = { 0.0, 1.0, 1.0, 1.0, 0.0, -1.0, -1.0, -1.0, 0.0 };
  T offsets[] = { -1.0, -1.0, 0.0, 1.0, 1.0, 1.0, 0.0, -1.0, -1.0 };

  // Piecewise Bezier knot vector for a quadratic curve with four segments 
  T knotsMem[] = { 0, 0, 0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75 , 1, 1, 1 };
  Vector<T> knots(knotsMem,12) ;

  resize(9,9,2,2);

  int i, j;

  double r2over2 = sqrt( 2.0 ) / 2.0;
  double weight;

  for (i = 0; i < 9; i++){
    for (j = 0; j < 9; j++) {
      HPoint_nD<T,N> hp ;
      weight = ((j & 1) ? r2over2 : 1.0) * ((i & 1) ? r2over2 : 1.0);
      // Notice how the weights are pre-multiplied with the x, y and z values
      P(i,j).x() = xvalues[j]* (R + offsets[i] * r) * weight;
      P(i,j).y() = yvalues[j]* (R + offsets[i] * r) * weight;
      P(i,j).z() = (zvalues[i] * r) * weight;
      P(i,j).w() = weight;
    }
  }

  // The knot vectors define piecewise Bezier segments 
  // (the same in both U and V).
  
  U = knots ;
  V = knots ; 

  MatrixRT<T> Tx ;
  Tx.translate(O.x(),O.y(),O.z()) ;
  transform(Tx) ;  
}


} // end namespace

#endif // #define PLIB_NURBS_NURBSS_SOURCE

---