Attachment #78077 for bug #120303

View | Details | Raw Unified | Return to bug 120303
Collapse All | Expand All




   echo >&2 "$0: script expects -patch|-unpatch as argument"
   exit 1
   -patch) patch -f -p0 < $0;;
   -unpatch) patch -f -R -p0 < $0;;
   *)
       echo >&2 "$0: script expects -patch|-unpatch as argument"
       exit 1
++ image/image.cpp

void MatrixImage<T>::drawLine(int i1, int j1, int i2, int j2, T color){
  int i,j ;
  double mx,b ;
  if(i1<0 || j1<0 || i1>this->rows() || j1>=this->cols()  ){
#ifdef USE_EXCEPTION
    throw OutOfBound2D(i1,j1,0,this->rows()-1,0,this->cols()-1) ;
#else
    Error error("MatrixImage<T>::drawLine") ;
    error << "Error in drawing line\n Invalid index ("<< i1 << ", " << j1 << ") to ( " << i2 << ", " << j2 << ") \n" ;

#endif
    return ;
  }
  if(i2 <0 || j2<0 || i2>this->rows() || j2>=this->cols() ){
#ifdef USE_EXCEPTION
    throw OutOfBound2D(i2,j2,0,this->rows()-1,0,this->cols()-1) ;
#else
    Error error("MatrixImage<T>::drawLine") ;
    error << "Error in drawing line\n Invalid index ("<< i1 << ", " << j1 << ") to ( " << i2 << ", " << j2 << ") \n" ;

  // check if line is vertical
  if(j1==j2){
    for(i=minimum(i1,i2);i<=maximum(i1,i2);i++)
     this->operator()(i,j1) = color ;
    return ;
  }
  mx = (double)(i1-i2)/(double)(j1-j2) ;

    if(i1>i2){
      for(i=i1;i>=i2;i--){
	j = int(((double)i-b)/mx) ;
	this->operator()(i,j) = color ;
      }
    }
    else{
      for(i=i1;i<=i2;i++){
	j = (int)((i-b)/mx) ;
	this->operator()(i,j) = color ;
      }
    }
  }

    if(j1>j2){
      for(j=j1;j>=j2;j--){
	i = (int)(mx*j+b) ;
	this->operator()(i,j) = color ;
      }
    }
    else{
      for(j=j1;j<=j2;j++){
	i = (int)(mx*j+b) ;
	this->operator()(i,j) = color ;
      }
    }
  }

void MatrixImage<T>::drawPoint(int i, int j, double r , T color){
  for(int y=i-int(ceil(r)) ; y<i+int(ceil(r)) ; y++)
    for(int x = j-int(ceil(r)) ; x<j+int(ceil(r)) ; x++){
      if(y>=0 && y<this->rows() && x>=0 && x<this->cols()){
	if( ((y-i)*(y-i)+(x-j)*(x-j))<= r*r)
	  this->operator()(y,x) = color ;
      }
    }
}

*/
template <class T>
void MatrixImage<T>::store(Matrix<T>& a){
  if(a.rows() != this->rows() || a.cols() != this->cols()) {
    a.resize(this->rows(),this->cols()) ;
  }
  T *aptr, *bptr ;
  int size,i ;
  aptr = &a(0,0)-1 ;
  bptr = this->m-1 ;
  size = this->cols()*this->rows() ;
  for(i=0;i<size;i++)
    *(++aptr) = *(++bptr) ;  
}
++ matrix/cvector.h

    CVector(const BasicArray<T>& v) : Vector<T>(v), index(0)  {;}
    virtual ~CVector() {}
    
    T& operator[](const int i) { return this->x[i%this->sze]; }
    T  operator[](const int i) const   { return this->x[i%this->sze]; }
    
    void put(T v) { this->x[index] = v ; index = (index+1)%this->sze; }
    
  protected:
    int index ;
++ matrix/list.h

{
public:
  BasicList() ;
  BasicList(const BasicList<T>& a) ;
  ~BasicList() { reset() ; }



  Modified by:
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
template <class T>
BasicList<T>::BasicList(const BasicList<T>& a): BasicNode<T>() {
  first_ = last_ = 0 ;
  current = first_ ;
  *this = a ; 
++ matrix/matrix.cpp

  if ( this == &a )
    return *this;
  
  if ( a.rows() != this->rows() || a.cols() != this->cols() ){
    resize(a.rows(),a.cols()) ;
  }
  
  int sze = this->rows()*this->cols() ;
  T *ptr, *aptr ;
  ptr = this->m-1 ;
  aptr = a.m-1 ;
  
  for (i = sze; i > 0; --i)
    *(++ptr) = *(++aptr) ;
  
  this->by_columns = a.by_columns;
  
  return *this;
}

{
  int rwz,coz,i,j;
  
  if ( this->rows() % a.rows() != 0 || this->cols() % a.cols() != 0 || this->rows() < a.rows() || this->cols() < a.cols() )
    {
#ifdef USE_EXCEPTION
      throw WrongSize2D(this->rows(),this->cols(),a.rows(),a.cols()) ;
#else
      Error error("Matrix<T>::submatrix");
      error << "Matrix and submatrix incommensurate" ;

#endif
    }
  
  if ( sr >= this->rows()/a.rows() || sr < 0 || sc >= this->cols()/a.cols() || sc < 0 )
    {
#ifdef USE_EXCEPTION
      throw OutOfBound2D(sr,sc,0,this->rows()/a.rows()-1,0,this->cols()/a.cols()-1) ;
#else
      Error error("Matrix<T>::submatrix");
      error << "Submatrix location out of bounds.\nrowblock " << sr << ", " << rows()/a.rows() << " colblock " << sc << ", " << a.cols() << endl ;

  aptr = a.m - 1;
  for ( i = a.rows()-1; i >= 0; --i )
    {
      ptr = &this->m[(i+rwz)*this->cols()+coz]-1 ;
      for ( j = a.cols(); j > 0; --j)
	*(++ptr) = *(++aptr) ;
    }  

  // Assign matrix a to this matrix at (i,j)
  int i, j;
  
  if ( (rw + a.rows()) > this->rows() || ( cl + a.cols()) > this->cols()) {
#ifdef USE_EXCEPTION
    throw MatrixErr();
#else

  T *pptr,*aptr ;
  aptr = a.m-1 ;
  for ( i = 0; i<a.rows(); ++i) {
    pptr = &this->m[(i+rw)*this->cols()+cl]-1 ;
    for ( j = 0; j < a.cols(); ++j)
      *(++pptr) = *(++aptr);
  }

Matrix<T> Matrix<T>::get(int rw, int cl, int nr, int nc) const
{
  Matrix<T> getmat(nr,nc) ;
  if ( (rw+nr) > this->rows() || (cl+nc) > this->cols()) {
#ifdef USE_EXCEPTION
    throw MatrixErr();
#else

  T *pptr,*aptr ;
  aptr = getmat.m-1;
  for (i = 0; i < nr; ++i) {
    pptr = &this->m[(i+rw)*this->cols()+cl]-1 ;
    for ( j = 0; j < nc; ++j)
      *(++aptr) = *(++pptr) ;
  }

  double sum, maxsum;
  int init=0 ;
  T *pptr ;
  pptr = this->m-1 ;
  maxsum = 0 ; // Silence the warning message
  for(i=0;i<this->rows();++i){
    sum = 0 ;
    for ( j = 0; j < this->cols(); ++j) 
      sum += *(++pptr) ;
    if(init)
      maxsum = (maxsum>sum) ? maxsum : sum;

{
  int i, iend;
  
  iend = this->rows();
  if ( iend > this->cols() )
    iend = this->cols();
  
  for (i = iend-1; i >=0; --i)
    this->elem(i,i) = a;

}


template <class T>
Vector<T> Matrix<T>::getDiag(){
  int i, iend;
  Vector<T> vec(minimum(this->rows(),this->cols())) ;
  iend = minimum(this->rows(),this->cols());
  for (i = iend-1; i >=0; --i)
      vec[i] = this->elem(i,i);
  return vec ;
}


Matrix<T>& Matrix<T>::operator+=(double a)
{
  T *p1 ;
  p1 = this->m-1 ;
  const int size = this->rows()*this->cols() ;
  for(int i=size; i>0; --i)
    *(++p1) += a ;  
  return *this ;

Matrix<T>& Matrix<T>::operator-=(double a)
{
  T *p1 ;
  p1 = this->m-1 ;
  const int size = this->rows()*this->cols() ;
  for(int i=size; i>0; --i)
    *(++p1) -= a ;  
  return *this ;

Matrix<T>& Matrix<T>::operator*=(double a)
{
  T *p1 ;
  p1 = this->m-1 ;
  const int size = this->rows()*this->cols() ;
  for(int i=size; i>0; --i)
    *(++p1) *= a ;  
  return *this ;

Matrix<T>& Matrix<T>::operator/=(double a)
{
  T *p1 ;
  p1 = this->m-1 ;
  const int size = this->rows()*this->cols() ;
  for(int i=size; i>0; --i)
    *(++p1) /= a ;  
  return *this ;

template <class T> 
Matrix<T>& Matrix<T>::operator+=(const Matrix<T> &a)
{
  if ( a.rows() != this->rows() || a.cols() != this->cols() )
    {
#ifdef USE_EXCEPTION
      throw WrongSize2D(this->rows(),this->cols(),a.rows(),a.cols());
#else
      Error error("Matrix<T>::operator+=") ;
      if ( rows() != a.rows() )

  int i, sze ;
  T *aptr,*sptr ;
  aptr = a.m - 1 ;
  sptr = this->m - 1 ;
  sze = this->rows()*this->cols() ;
  for (i = sze; i > 0; --i){
      *(++sptr) += *(++aptr) ;
  }

template <class T> 
Matrix<T>& Matrix<T>::operator-=(const Matrix<T> &a)
{
  if ( a.rows() != this->rows() || a.cols() != this->cols() )
    {
#ifdef USE_EXCEPTION
      throw WrongSize2D(this->rows(),this->cols(),a.rows(),a.cols());
#else
      Error error("Matrix<T>::operator-=") ;
      if ( rows() != a.rows() )

  int i, size;
  T *aptr,*sptr ;
  aptr = a.m - 1 ;
  sptr = this->m - 1 ;
  size = this->rows()*this->cols() ;
  for (i = size; i > 0; --i){
      *(++sptr) -= *(++aptr) ;
  }

template <class T>
T Matrix<T>::trace() const
{
  int size = this->rows();
  T sum = (T)0;
  
  if ( size > this->cols() )
    size = this->cols();
  
  for (int d = 0; d < size; ++d)
    sum += this->elem(d,d) ;
  
  return sum;
}

template <class T>
Matrix<T> Matrix<T>::herm() const
{
  int i, j, r = this->cols(), c = this->rows();
  Matrix<T> adj(r,c);
  
  for (i = 0; i < r; ++i)
    for (j = 0; j < c; ++j)
      adj.elem(i,j) = this->elem(j,i) ;

  return adj;


template <class T>
Matrix<T> Matrix<T>::flop() const
{					
  Matrix<T> f(this->rows(),this->cols()) ;
  for(int i=this->rows()-1;i>=0;--i)
    for(int j=this->cols()-1;j>=0;--j)
      {
	f(i,j) = elem(i,this->cols()-j-1);
      }
  return f; 
}

{					
  // same as hermitian for real Matrix<T>
  int i, j;
  const int& r = this->cols();
  const int& c = this->rows();
  Matrix<T> adj(r,c);
  
  for (i = r-1; i >=0; --i)
    for (j = c-1; j >=0; --j)
      adj.elem(i,j) = this->elem(j,i) ;
  
  
  return adj; 

int Matrix<T>::read(char* filename) {
  ifstream fin(filename) ;
  if(!fin) {
    this->resize(1,1) ;
    return 0 ;
  }
  int r,c ;

  if(r) return 0 ;
  if(!fin.read((char*)&r,sizeof(int))) return 0 ;
  if(!fin.read((char*)&c,sizeof(int))) return 0 ;
  this->resize(r,c) ;
  if(!fin.read((char*)this->m,sizeof(T)*r*c)) return 0 ;

  delete []type ;
  return 1 ;

int Matrix<T>::read(char* filename,int r, int c) {
  ifstream fin(filename) ;
  if(!fin) {
    this->resize(1,1) ;
    return 0 ;
  }
  this->resize(r,c) ;
  if(!fin.read((char*)this->m,sizeof(T)*r*c)) return 0 ;

  return 1 ;
}

  if(!fout)
    return 0 ;
  int r,c ;
  r = this->rows() ; c = this->cols() ;
  if(!fout.write((char*)&"matrix",sizeof(char)*6)) return 0 ;
  if(!fout.write((char*)&r,sizeof(int))) return 0 ;
  if(!fout.write((char*)&c,sizeof(int))) return 0 ;
  if(!fout.write((char*)this->m,sizeof(T)*r*c)) return 0 ;
  return 1;
}


  ofstream fout(filename) ;
  if(!fout)
    return 0 ;
  if(!fout.write((char*)this->m,sizeof(T)*this->rows()*this->cols())) return 0 ;
  return 1;
}

++ matrix/vector.cpp

  if(this==&b)
    return *this ;

  if ( this->n() != b.n())
    {
      resize(b.n()) ;
    }

  this->sze = b.n() ;
  T *pa, *pb ;
  pa = this->x-1 ;
  pb = b.x-1 ;
  for(int i=this->n();i>0;--i){
    *(++pa) = *(++pb) ;
  }
  return *this;

template <class T>
Vector<T>& Vector<T>::operator=(const BasicArray<T> &b)
{
  if ( this->size() != b.size())
    {
      resize(b.size()) ;
    }
  T *ptr ;
  ptr = this->x - 1 ;
  for(int i=this->size()-1;i>=0;--i)
     *(++ptr) = b[i] ;

  return *this;

template <class T>
T Vector<T>::operator=(const T d)
{
  const int sz = this->size(); 
  T *ptr ;
  ptr = this->x-1 ;
  for (int i = sz; i > 0; --i)
    *(++ptr) = d ;


template <class T>
Vector<T>& Vector<T>::operator+=(const Vector<T> &a)
{
  if ( a.size() != this->size())
    {
#ifdef USE_EXCEPTION
      throw WrongSize(this->size(),a.size()) ;
#else
      Error error("Vector<T>::operator+=(Vector<T>&)");
      error << "Vector<T> a += Vector<T> b different sizes, a = " << size() << ", b = " << a.size() ;
      error.fatal() ;
#endif
    }
  const int sz = this->size();
  T *ptr,*aptr ;
  ptr = this->x-1 ;
  aptr = a.x-1 ;
  for (int i = sz; i >0; --i)
    *(++ptr) += *(++aptr) ;

template <class T>
Vector<T>& Vector<T>::operator-=(const Vector<T> &a)
{
  if ( a.size() != this->size())
    {
#ifdef USE_EXCEPTION
      throw WrongSize(this->size(),a.size()) ;
#else
      Error error("Vector<T>::operator-=(Vector<T>&)");
      error << "Vector<T> a -= Vector<T> b different sizes, a = " << size() << ", b = " << a.size() ;

#endif
    }
  
  const int sz = this->size(); 
  T *ptr,*aptr ;
  ptr = this->x-1 ;
  aptr = a.x-1 ;
  for (int i = sz; i > 0; --i)
    *(++ptr) -= *(++aptr) ;

    }

  T *aptr,*bptr ;
  aptr = &this->x[i]-1 ;
  bptr = b.x-1 ;
  for ( int j = b.rows(); j > 0; --j)
      *(++aptr) = *(++bptr) ;


  Vector<T> subvec(l) ;
  T *aptr, *bptr ;
  aptr = &this->x[i] - 1 ;
  bptr = subvec.x -1 ;
  for ( int j = l; j > 0; --j)
    *(++bptr) = *(++aptr) ;

*/
template <class T>
int Vector<T>::minIndex() const {
  T min = this->x[0] ;
  int index = 0 ;

  for(int i=1;i<this->n();i++){
    if(this->x[i]<=min){
      min = this->x[i] ;
      index = i ;
    }
  }

  T a ;
  T *v1,*v2  ;

  ir = this->sze-1 ;
  l = 0 ;
  
  while(1){
    if(ir-l<M){ // perform an insertion sort when the array is small enough
      v1 = &this->x[l] ;
      for(j=l+1;j<=ir;++j){
	a = *(++v1) ;
	v2 = v1 ;

    }
    else{
      k=(l+ir) >> 1 ;
      swap(this->x[k],this->x[l+1]) ;
      if(this->x[l+1] > this->x[ir]){
	swap(this->x[l+1],this->x[ir]) ;
      }
      if(this->x[l]> this->x[ir]){
	swap(this->x[l],this->x[ir]) ;
      }
      if(this->x[l+1] > this->x[l]){
	swap(this->x[l+1],this->x[l]) ;
      }
      i=l+1 ;
      j=ir ;
      a=this->x[l] ;
      v1 = &this->x[i] ;
      v2 = &this->x[j] ;
      while(1){
	while(*v1 < a) { ++i ; ++v1 ; }
	while(*v2 > a) { --j ; --v2 ; }
	if(j<i) break ;
	if(*v1 == *v2)  // both are equal to a...
	  break ;
	swap(this->x[i],this->x[j]) ;
      }
      this->x[l] = this->x[j] ;
      this->x[j] = a ;
      jstack += 2 ;
      if(jstack>=Nstack){
	istack.resize(istack.n()+Nstack) ; // increase the vector size

  int jstack=0;
  T a ;

  ir = this->sze-1 ;
  l = 0 ;
  
  index.resize(this->sze) ;
  for(i=0;i<index.n();++i)
    index[i] = i ;


    if(ir-l<M){ // perform an insertion sort when the array is small enough
      for(j=l+1;j<=ir;++j){
	indext = index[j] ;
	a = this->x[indext] ;
	for(i=j-1;i>=0;--i){
	  if(this->x[index[i]] <= a) break ;
	  index[i+1] = index[i] ;
	}
	index[i+1] = indext ;

    else{
      k=(l+ir) >> 1 ;
      swap(index[k],index[l+1]) ;
      if(this->x[index[l+1]] > this->x[index[ir]]){
	swap(index[l+1],index[ir]) ;
      }
      if(this->x[index[l]]> this->x[index[ir]]){
	swap(index[l],index[ir]) ;
      }
      if(this->x[index[l+1]] > this->x[index[l]]){
	swap(index[l+1],index[l]) ;
      }
      i=l+1 ;
      j=ir ;
      indext = index[l] ;
      a=this->x[indext] ;
      while(1){
	while(this->x[index[i]] < a) { ++i ; }
	while(this->x[index[j]] > a) { --j ; }
	if(j<i) break ;
	if(this->x[index[i]] == this->x[index[j]])
	  break ;
	swap(index[i],index[j]) ;
      }
++ matrix/vector.h

  {
  public:
    int rows() const //!< a reference to the size of the vector
      { return this->sze ;}
    Vector() : BasicArray<T>(1) {} //!< Basic constructor
    Vector(const int r) : BasicArray<T>(r) {}
    Vector(const Vector<T>& v) : BasicArray<T>(v) {}
++ numerical/matrixMat.cpp

template <class T>
LUMatrix<T>& LUMatrix<T>::operator=(const LUMatrix<T>& a){
  resize(a.rows(),a.cols()) ;
  for(int i=0;i<this->rows();++i)
    for(int j=0;j<this->cols();++j)
      this->elem(i,j) = a(i,j) ;
  pivot_ = a.pivot_ ;
  return *this ;
}

  //	lu = a;	 must do it by copying or LUFACT will be recursively called !
  for(i=0;i<n;++i)
    for(j=0;j<n;++j)
      this->elem(i,j) = a(i,j) ;

  errval = 0;
  nm1 = n - 1;

	    }
	  pivot_[k] = l;

	  if ( this->elem(l,k) != 0.0 )
	    {			// nonsingular pivot found 
	      if (l != k ){	// interchange needed 
		for (i = k; i < n; i++)
		  {
		    t = this->elem(l,i) ;
		    this->elem(l,i) = this->elem(k,i) ;
		    this->elem(k,i) = t ; 
		  }
		sign = -sign ;
	      }
	      q =  this->elem(k,k);	/* scale row */
	      for (i = kp1; i < n; i++)
		{
		  t = - this->elem(i,k)/q;
		  this->elem(i,k) = t;
		  for (j = kp1; j < n; j++)
		    this->elem(i,j) += t * this->elem(k,j);
		}
	    }
	  else		/* pivot singular */

    }
  
  pivot_[nm1] = nm1;
  if (this->elem(nm1,nm1) == 0.0)
    errval = nm1;  
  return *this;
}

*/
template <class T>
T LUMatrix<T>::determinant(){
  T det = this->elem(0,0) ;
  for(int i=1;i<this->rows();++i)
    det *= this->elem(i,i) ;
  return det * (T)sign ;
}


  T ten;
  int i, j, k, l, kb, kp1, nm1, n, coln;

  if ( this->rows() != this->cols() )
    {
#ifdef USE_EXCEPTION
    throw WrongSize2D(this->rows(),this->cols(),0,0) ;
#else
      Error error("invm");
      error << "matrix inverse, not square: " << rows() << " by " << cols() << endl;

#endif
    }

  n = coln = this->rows();


  inv = *this ;

template <class T>
Matrix<T> LUMatrix<T>::inverse() 
{
  if ( this->rows() != this->cols() )
    {
#ifdef USE_EXCEPTION
      throw WrongSize2D(this->rows(),this->cols(),0,0) ;
#else
      Error error("invm");
      error << "matrix inverse, not square: " << rows() << " by " << cols() << endl;
++ nurbs/d_nurbsSub.cpp

  template class RenderMeshPoints<double> ;

  
  template <> double NurbSurface<double>::epsilon = 1e-6 ;
  template <> double SurfSample<double>::epsilon = 1e-6 ;

  template void DrawSubdivision( NurbSurface<double> *, double tolerance );
  template void DrawEvaluation( NurbSurface<double> * );
++ nurbs/f_nurbsSub.cpp

  template class RenderMeshPoints<float> ;

  
  template <> float NurbSurface<float>::epsilon = 1e-6 ;
  template <> float SurfSample<float>::epsilon = 1e-6 ;

  template void DrawSubdivision( NurbSurface<float> *, float tolerance );
  template void DrawEvaluation( NurbSurface<float> * );
++ nurbs/hnurbsS.cpp

  initBase() ;
  offset.resize(baseSurf.ctrlPnts()) ;

  this->P = baseSurf.ctrlPnts() ;
  this->U = baseSurf.knotU() ;
  this->V = baseSurf.knotV() ;
  this->degU = baseSurf.degreeU() ;
  this->degV = baseSurf.degreeV() ;

  //updateSurface() ;


  baseUpdateN = baseLevel_->modifiedN()-1 ; // Set it so that initBase will run
  initBase() ;
  offset.resize(baseSurf.ctrlPnts()) ;
  this->P = baseSurf.ctrlPnts() ;
  this->U = baseSurf.knotU() ;
  this->V = baseSurf.knotV() ;
  this->degU = baseSurf.degreeU() ;
  this->degV = baseSurf.degreeV() ;
  //updateSurface() ;

}

  rU.resize(0) ;
  rV.resize(0) ;

  offset = this->P ;
}

/*!

  }
  if(baseLevel_){
    if(initBase()){
      this->P = baseSurf.ctrlPnts() ;
      this->U = baseSurf.knotU() ;
      this->V = baseSurf.knotV() ;
      this->degU = baseSurf.degreeU() ;
      this->degV = baseSurf.degreeV() ;
    }
    if(i0>=0 && j0>=0){
      Point_nD<T,N> vecOffset ;

	  offset(i0,j0).y()*jvec(i0,j0) +
	  offset(i0,j0).z()*kvec(i0,j0) ;
      }
      this->P(i0,j0).x() = baseSurf.ctrlPnts()(i0,j0).x()+vecOffset.x() ;
      this->P(i0,j0).y() = baseSurf.ctrlPnts()(i0,j0).y()+vecOffset.y() ;
      this->P(i0,j0).z() = baseSurf.ctrlPnts()(i0,j0).z()+vecOffset.z() ;
    }
    else{
      for(int i=0;i<this->P.rows();++i)
	for(int j=0;j<this->P.cols();++j){
	  if(offset(i,j).x() != 0 || 
	     offset(i,j).y() != 0 || offset(i,j).z() || 0){
	    Point_nD<T,N> vecOffset ;

		offset(i,j).y()*jvec(i,j) +
		offset(i,j).z()*kvec(i,j) ;
	    }
	    this->P(i,j).x() = baseSurf.ctrlPnts()(i,j).x()+vecOffset.x() ;
	    this->P(i,j).y() = baseSurf.ctrlPnts()(i,j).y()+vecOffset.y() ;
	    this->P(i,j).z() = baseSurf.ctrlPnts()(i,j).z()+vecOffset.z() ;
	  }
	}
    }
  }
  else{
    if(i0>=0 && j0>=0)
      this->P(i0,j0) = offset(i0,j0) ;
    else{
      for(int i=0;i<this->P.rows();++i)
	for(int j=0;j<this->P.cols();++j){
	  this->P(i,j) = offset(i,j) ;
	}
    }
  }

      return mod ;
  }

  if(u<this->knotU()[0] || u>this->knotU()[this->knotU().n()-1])
    return -1 ;
  if(v<this->knotV()[0] || v>this->knotU()[this->knotV().n()-1])
    return -1 ;

  int su = findSpanU(u) ;
  int sv = findSpanV(v) ;

  for(int i=0;i<=this->degU;++i)
    for(int j=0;j<=this->degV;++j){
      if(offset(su-this->degU+i,sv+this->degV+j) != HPoint_nD<T,N>(0,0,0,0))
	return level_ ;
    }


template <class T, int N>
void HNurbsSurface<T,N>::splitUV(int nu, int nv, Vector<T> &nU, Vector<T> &nV){

  nU.resize(this->knotU().n()*nu) ;
  nV.resize(this->knotV().n()*nv) ;
  
  int i,j,n ;

  n = 0 ; 
  for(i=1;i<this->knotU().n();++i){
    if(this->knotU()[i] >this->knotU()[i-1]){
      T a = this->knotU()[i-1] ;
      T b = this->knotU()[i] ;


      for(j=0;j<nu;++j){

  nU.resize(n) ;

  n = 0 ;
  for(i=1;i<this->knotV().n();++i){
    if(this->knotV()[i] > this->knotV()[i-1]){
      T a = this->knotV()[i-1] ;
      T b = this->knotV()[i] ;

      for(j=0;j<nv;++j){
	nV[n] = a + (b-a)*T(j+1)/T(nv+1) ;

  int i,j,n ;

  if(su<=0)
    su = this->degU  ;
  if(sv<=0)
    sv = this->degV  ;
  if(su>this->degU+1)
    su = this->degU+1 ;
  if(sv>this->degV+1)
    sv = this->degV+1 ;

  nU.resize(this->knotU().n()*nu*su) ;
  nV.resize(this->knotV().n()*nv*sv) ;
  
  n = 0 ; 
  for(i=1;i<this->knotU().n();++i){
    if(this->knotU()[i] >this->knotU()[i-1]){
      T a = this->knotU()[i-1] ;
      T b = this->knotU()[i] ;


      for(j=0;j<nu;++j){

  nU.resize(n) ;

  n = 0 ;
  for(i=1;i<this->knotV().n();++i){
    if(this->knotV()[i] > this->knotV()[i-1]){
      T a = this->knotV()[i-1] ;
      T b = this->knotV()[i] ;

      for(j=0;j<nv;++j){
	T v = a + (b-a)*T(j+1)/T(nv+1) ;

    if(!fin.read((char*)&du,sizeof(int))) { delete []type ; return 0 ;}
    if(!fin.read((char*)&dv,sizeof(int))) { delete []type ; return 0 ;}
    
    this->resize(nu,nv,du,dv) ;
    
    if(!fin.read((char*)this->U.memory(),sizeof(T)*this->U.n())) { delete []type ; return 0 ;}
    if(!fin.read((char*)this->V.memory(),sizeof(T)*this->V.n())) { delete []type ; return 0 ;}
    
    if(!r1){
      p = new T[3*nu*nv] ;

      p2 = p ;
      for(int i=0;i<nu;i++)
	for(int j=0;j<nv;j++){
	  this->P(i,j).x() = *(p++) ;
	  this->P(i,j).y() = *(p++) ;
	  this->P(i,j).z() = *(p++) ;
	  this->P(i,j).w() = 1.0 ;
	}
      delete []p2 ;
    }

      p2 = p ;
      for(int i=0;i<nu;i++)
	for(int j=0;j<nv;j++){
	  this->P(i,j).x() = *(p++) ;
	  this->P(i,j).y() = *(p++) ;
	  this->P(i,j).z() = *(p++) ;
	  this->P(i,j).w() = *(p++) ;
	}
      delete []p2 ;
    }
    offset = this->P ;
    this->updateSurface() ;
  }
  else { // reading the offset information

  if(!fout)
    return 0 ;
  if(!baseLevel_){
    int prows = this->P.rows();
    int pcols = this->P.cols();
    char st = '0' + sizeof(T) ; 
    if(!fout.write((char*)&"hns4",sizeof(char)*4)) 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*)&this->degU,sizeof(int))) return 0 ;
    if(!fout.write((char*)&this->degV,sizeof(int))) return 0 ;
    if(!fout.write((char*)this->U.memory(),sizeof(T)*this->U.n())) return 0 ;
    if(!fout.write((char*)this->V.memory(),sizeof(T)*this->V.n())) return 0 ;
    
    T *p,*p2 ;
    p = new T[this->P.rows()*this->P.cols()*4] ;
    p2 = p ;
    for(int i=0;i<this->P.rows();i++) 
      for(int j=0;j<this->P.cols();j++){
	*p = offset(i,j).x() ; p++ ;
	*p = offset(i,j).y() ; p++ ;
	*p = offset(i,j).z() ; p++ ;
	*p = offset(i,j).w() ; p++ ;
      }
    if(!fout.write((char*)p2,sizeof(T)*this->P.rows()*this->P.cols()*4)) return 0 ;
    delete []p2 ;
  }
  else{

  int i,j ;
  j = 0 ;
  for(i=0;i<X.n();++i){
    if(X[i]>=this->U[0] && X[i]<=this->U[this->U.n()-1]){
      Xu[j] = X[i] ;
      ++j ;
    }

      nextLevel_->refineKnotU(Xu) ;
    }
    
    NurbsSurface<T,N> osurf(this->degU,this->degV,this->U,this->V,offset) ;
    
    osurf.refineKnotU(Xu) ;
    

  int i,j ;
  j = 0 ;
  for(i=0;i<X.n();++i){
    if(X[i]>=this->V[0] && X[i]<=this->V[this->V.n()-1]){
      Xv[j] = X[i] ;
      ++j ;
    }

      nextLevel_->refineKnotV(Xv) ;
    }
    
    NurbsSurface<T,N> osurf(this->degU,this->degV,this->U,this->V,offset) ;
    
    osurf.refineKnotV(Xv) ;
    

*/
template <class T, int N>
int HNurbsSurface<T,N>::movePointOffset(T u, T v, const Point_nD<T,N>& delta){
  this->P = offset ; 

  // by definition the offset has w = 0 , but this isn't valid for
  // the control points, increasing the w by 1, will generate a valid surface
  if(baseLevel_)
    for(int i=0;i<this->P.rows();++i)
      for(int j=0;j<this->P.cols();++j){
	this->P(i,j).w() += T(1) ; 
      }

  if(NurbsSurface<T,N>::movePoint(u,v,delta)){
    offset = this->P ;
    // need to reset the offset weights
    if(baseLevel_)
      for(int i=0;i<this->P.rows();++i)
	for(int j=0;j<this->P.cols();++j){
	  this->P(i,j).w() -= T(1) ; 
      }
    
    this->P = baseSurf.ctrlPnts() ; 
    updateSurface() ; 
    return 1 ;
  }
++ nurbs/hnurbsS_sp.cpp

*/
template <class T, int N>
void HNurbsSurfaceSP<T,N>::updateMaxU() {
  if(this->degU>3){
#ifdef USE_EXCEPTION
    throw NurbsError();
#else

#endif
  }
  else{
    maxU.resize(this->P.rows()) ;
    maxAtU_.resize(this->P.rows()) ;
    for(int i=0;i<this->P.rows();++i){
      if(!maxInfluence(i,this->U,this->degU,maxAtU_[i]))
	cerr << "Problem in maxInfluence U!\n" ;
      maxU[i] = nurbsBasisFun(maxAtU_[i],i,this->degU,this->U) ;
    }
    
  }

*/
template <class T, int N>
void HNurbsSurfaceSP<T,N>::updateMaxV() {
  if(this->degV>3){
#ifdef USE_EXCEPTION
    throw NurbsError();
#else

#endif
  }
  else{
    maxV.resize(this->P.cols()) ;
    maxAtV_.resize(this->P.cols()) ;
    for(int i=0;i<this->P.cols();++i){
      if(!maxInfluence(i,this->V,this->degV,maxAtV_[i]))
	cerr << "Problem in maxInfluence V!\n" ;
      maxV[i] = nurbsBasisFun(maxAtV_[i],i,this->degV,this->V) ;
    }
    
  }

*/
template <class T, int N>
void HNurbsSurfaceSP<T,N>::modSurfCPby(int i, int j, const HPoint_nD<T,N>& a) {
  this->offset(i,j) +=  a / (maxU[i]*maxV[j]) ; 
  if(this->baseLevel_){
    Point_nD<T,N> vecOffset ; 
    vecOffset = this->offset(i,j).x()*this->ivec(i,j) +
      this->offset(i,j).y()*this->jvec(i,j) +
      this->offset(i,j).z()*this->kvec(i,j) ;
    this->P(i,j).x() = this->baseSurf.ctrlPnts()(i,j).x()+vecOffset.x() ;
    this->P(i,j).y() = this->baseSurf.ctrlPnts()(i,j).y()+vecOffset.y() ;
    this->P(i,j).z() = this->baseSurf.ctrlPnts()(i,j).z()+vecOffset.z() ;
  }
  else
    this->P(i,j) = this->offset(i,j) ; 
}

/*!

void HNurbsSurfaceSP<T,N>::modOnlySurfCPby(int i, int j, const HPoint_nD<T,N>& a){
  int k ; 

  this->P = this->offset ; 

  // by definition the offset has w = 0 , but this isn't valid for
  // the control points, increasing the w by 1, will generate a valid surface
  if(this->baseLevel_)
    for(k=0;k<this->P.rows();++k)
      for(int l=0;l<this->P.cols();++l)
	this->P(k,l).w() += T(1) ; 

  int sizeU, sizeV ;

  sizeU = 2*this->degU+3 ; 
  if(i-this->degU-1<0) sizeU += i-this->degU-1 ; 
  if(i+this->degU+1>=this->P.rows()) sizeU -= i+this->degU+1-this->P.rows() ;

  sizeV = 2*this->degV+3 ;
  if(j-this->degV-1<0) sizeV += j-this->degV-1 ; 
  if(j+this->degV+1>=this->P.cols()) sizeV -= j+this->degV+1-this->P.cols() ;
  
  Vector<T> u(sizeU) ;
  Vector<T> v(sizeV) ;

  int n=0;
  int nu = 0 ;
  int nv = 0 ; 
  for(k=i-this->degU-1;k<=i+this->degU+1;++k){
    if(k<0)
      continue ;
    if(k>=this->P.rows())
      break ; 
    nv = 0 ;
    for(int l=j-this->degV-1;l<=j+this->degV+1;++l){
      if(l<0)
	continue ;
      if(l>=this->P.cols())
	break ; 
      if( k == i && j==l){
	pts[n].x() = a.x() ; 

  pv.resize(n) ; 

  if(NurbsSurface<T,N>::movePoint(u,v,pts,pu,pv)){
    this->offset = this->P ; 
    // an offset shouldn't have a weight value.
    if(this->baseLevel_)
      for(k=0;k<this->P.rows();++k)
	for(int l=0;l<this->P.cols();++l)
	  this->offset(k,l).w() -= T(1) ; 
  }
  updateSurface(); 
}

HNurbsSurfaceSP<T,N>* HNurbsSurfaceSP<T,N>::addLevel(int n, int s) {
  HNurbsSurfaceSP<T,N> *newLevel ;

  if(this->nextLevel_)
    return 0 ;

  Vector<T> newU,newV ;

HNurbsSurfaceSP<T,N>* HNurbsSurfaceSP<T,N>::addLevel() {
  HNurbsSurfaceSP<T,N> *newLevel ;

  if(this->nextLevel_)
    return 0 ;

  newLevel = new HNurbsSurfaceSP<T,N>(this) ;

  levelP = nS.nextLevel() ;

  NurbsSurface<T,N>::operator=(nS) ;
  this->rU = nS.rU ;
  this->rV = nS.rV ;
  this->offset = nS.offset ;

  updateMaxUV() ; 

  this->firstLevel_ = this ;

  if(levelP){
    HNurbsSurfaceSP<T,N> *newLevel ;
    newLevel =  new HNurbsSurfaceSP<T,N>(this) ; 
    newLevel->copy(*levelP) ;
    this->nextLevel_ = newLevel ;
    this->lastLevel_ = this->nextLevel_->lastLevel() ;
  }
  else{
    this->lastLevel_ = this ;
  }

}

template <class T, int N>
void HNurbsSurfaceSP<T,N>::updateSurface(int i0, int j0){
  if(i0>=0 && j0>=0){
    if(this->offset(i0,j0).x()==0.0 && this->offset(i0,j0).y()==0.0 && this->offset(i0,j0).z()==0.0)
      return ;
  }
  if(this->baseLevel_){
    if(this->initBase()){
      this->P = this->baseSurf.ctrlPnts() ;
      this->U = this->baseSurf.knotU() ;
      this->V = this->baseSurf.knotV() ;
      this->degU = this->baseSurf.degreeU() ;
      this->degV = this->baseSurf.degreeV() ;
      updateMaxUV() ; 
    }
    if(i0>=0 && j0>=0){
      Point_nD<T,N> vecOffset ;
      vecOffset = this->offset(i0,j0).x()*this->ivec(i0,j0) +
	this->offset(i0,j0).y()*this->jvec(i0,j0) +
	this->offset(i0,j0).z()*this->kvec(i0,j0) ;
      this->P(i0,j0).x() = this->baseSurf.ctrlPnts()(i0,j0).x()+vecOffset.x() ;
      this->P(i0,j0).y() = this->baseSurf.ctrlPnts()(i0,j0).y()+vecOffset.y() ;
      this->P(i0,j0).z() = this->baseSurf.ctrlPnts()(i0,j0).z()+vecOffset.z() ;
    }
    else{
      for(int i=0;i<this->P.rows();++i)
	for(int j=0;j<this->P.cols();++j){
	  if(this->offset(i,j).x() != 0 || 
	     this->offset(i,j).y() != 0 || this->offset(i,j).z() || 0){
	    Point_nD<T,N> vecOffset ;
	    vecOffset = this->offset(i,j).x()*this->ivec(i,j) +
	      this->offset(i,j).y()*this->jvec(i,j) +
	      this->offset(i,j).z()*this->kvec(i,j) ;
	    this->P(i,j).x() = this->baseSurf.ctrlPnts()(i,j).x()+vecOffset.x() ;
	    this->P(i,j).y() = this->baseSurf.ctrlPnts()(i,j).y()+vecOffset.y() ;
	    this->P(i,j).z() = this->baseSurf.ctrlPnts()(i,j).z()+vecOffset.z() ;
	  }
	}
    }
  }
  else{
    if(i0>=0 && j0>=0)
      this->P(i0,j0) = this->offset(i0,j0) ;
    else{
      for(int i=0;i<this->P.rows();++i)
	for(int j=0;j<this->P.cols();++j){
	  this->P(i,j) = this->offset(i,j) ;
	}
    }
  }

  ++this->updateN ;
}

/*!

  if(!okMax())
    updateMaxUV() ; 
  if(upLevel>=0){
    if(this->level()<=upLevel){
      this->updateSurface() ;
    }
  }

    this->updateSurface() ;
  }

  if(upLevel>this->level() || upLevel<0){
    if(this->nextLevel_)
      ((HNurbsSurfaceSP<T,N>*)this->nextLevel_)->updateLevels(upLevel) ;
  }
}


    if(!fin.read((char*)&du,sizeof(int))) { delete []type ; return 0 ;}
    if(!fin.read((char*)&dv,sizeof(int))) { delete []type ; return 0 ;}
    
    this->resize(nu,nv,du,dv) ;
    
    if(!fin.read((char*)this->U.memory(),sizeof(T)*this->U.n())) { delete []type ; return 0 ;}
    if(!fin.read((char*)this->V.memory(),sizeof(T)*this->V.n())) { delete []type ; return 0 ;}
    
    if(!r1){
      p = new T[3*nu*nv] ;

      p2 = p ;
      for(int i=0;i<nu;i++)
	for(int j=0;j<nv;j++){
	  this->P(i,j).x() = *(p++) ;
	  this->P(i,j).y() = *(p++) ;
	  this->P(i,j).z() = *(p++) ;
	  this->P(i,j).w() = 1.0 ;
	}
      delete []p2 ;
    }

      p2 = p ;
      for(int i=0;i<nu;i++)
	for(int j=0;j<nv;j++){
	  this->P(i,j).x() = *(p++) ;
	  this->P(i,j).y() = *(p++) ;
	  this->P(i,j).z() = *(p++) ;
	  this->P(i,j).w() = *(p++) ;
	}
      delete []p2 ;
    }
    this->offset = this->P ;
    this->updateSurface() ;
  }
  else { // reading the offset information
    int ru,rv ;
    if(!fin.read((char*)&ru,sizeof(int))) { delete []type ; return 0 ;}
    if(!fin.read((char*)&rv,sizeof(int))) { delete []type ; return 0 ;}
    this->rU.resize(ru) ;
    this->rV.resize(rv) ;
    if(this->rU.n()>0)
      if(!fin.read((char*)this->rU.memory(),sizeof(T)*this->rU.n())) { delete []type ; return 0 ;}
    if(this->rV.n()>0)
      if(!fin.read((char*)this->rV.memory(),sizeof(T)*this->rV.n())) { 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 ;}

    p = new T[4*nu*nv] ;
    if(!fin.read((char*)p,sizeof(T)*4*nu*nv)) { delete []type ; return 0 ;}
    p2 = p ;
    this->offset.resize(nu,nv) ;
    for(int i=0;i<nu;i++)
      for(int j=0;j<nv;j++){
	this->offset(i,j).x() = *(p++) ;
	this->offset(i,j).y() = *(p++) ;
	this->offset(i,j).z() = *(p++) ;
	this->offset(i,j).w() = *(p++) ;
      }
    delete []p2 ;    
    --this->baseUpdateN ;
    this->updateSurface() ;
  }

++ nurbs/matrixRT.cpp

  // *this = C.translate(x,y,z)*B.rotate(ax,ay,az) ;
  rotate(ax,ay,az) ;
#ifdef COLUMN_ORDER
  this->m[12] = x ;
  this->m[13] = y ;
  this->m[14] = z ;  
#else
  this->m[3] = x ;
  this->m[7] = y ;
  this->m[11] = z ;
#endif
}


 */
template <class T>
MatrixRT<T>::MatrixRT() : Matrix<T>(4,4) {
  this->reset(0) ;
  this->diag(1.0) ;
}

/*!

  t10 = cos(ax);
  t13 = t4*t6;
#ifdef COLUMN_ORDER
  this->m[0] = t1*t2;
  this->m[4] = -t4*t2;
  this->m[8] = t6;
  this->m[12] = 0 ;
  this->m[1] = t7*t8+t4*t10;
  this->m[5] = -t13*t8+t1*t10;
  this->m[9] = -t2*t8;
  this->m[13] = 0 ;
  this->m[2] = -t7*t10+t4*t8;
  this->m[6] = t13*t10+t1*t8;
  this->m[10] = t2*t10;
  this->m[14] = this->m[3] = this->m[7] = this->m[11] = 0.0 ;
  this->m[15] = 1.0 ;
#else
  this->m[0] = t1*t2;
  this->m[1] = -t4*t2;
  this->m[2] = t6;
  this->m[3] = 0 ;
  this->m[4] = t7*t8+t4*t10;
  this->m[5] = -t13*t8+t1*t10;
  this->m[6] = -t2*t8;
  this->m[7] = 0 ;
  this->m[8] = -t7*t10+t4*t8;
  this->m[9] = t13*t10+t1*t8;
  this->m[10] = t2*t10;
  this->m[11] = this->m[12] = this->m[13] = this->m[14] = 0 ;
  this->m[15] = 1.0 ;
#endif
  return *this ;
}

  t9 = (T)sin((double)ax);
  t17 = t4*t7;
#ifdef COLUMN_ORDER
  this->m[0] = t1*t2;
  this->m[4] = -t4*t5+t8*t9;
  this->m[8] = t4*t9+t8*t5;
  this->m[12] = 0.0 ;
  this->m[1] = t4*t2;
  this->m[5] = t1*t5+t17*t9;
  this->m[9] = -t1*t9+t17*t5;
  this->m[13] = 0.0 ;
  this->m[2] = -t7;
  this->m[6] = t2*t9;
  this->m[10] = t2*t5;
  this->m[14] = this->m[3] = this->m[7] = this->m[11] = 0 ;
  this->m[15] = 1.0 ;
#else
  this->m[0] = t1*t2;
  this->m[1] = -t4*t5+t8*t9;
  this->m[2] = t4*t9+t8*t5;
  this->m[3] = 0.0 ;
  this->m[4] = t4*t2;
  this->m[5] = t1*t5+t17*t9;
  this->m[6] = -t1*t9+t17*t5;
  this->m[7] = 0.0 ;
  this->m[8] = -t7;
  this->m[9] = t2*t9;
  this->m[10] = t2*t5;
  this->m[11] = this->m[12] = this->m[13] = this->m[14] = 0 ;
  this->m[15] = 1.0 ;
#endif
  return *this ;
}

 */
template <class T>
MatrixRT<T>& MatrixRT<T>::translate(T x, T y, T z){
  this->reset(0) ;
  this->diag(1.0) ;
#ifdef COLUMN_ORDER
  this->m[12] = x ;
  this->m[13] = y ;
  this->m[14] = z ;
#else
  this->m[3] = x ;
  this->m[7] = y ;
  this->m[11] = z ;
#endif
  return *this ;
}

 */
template <class T>
MatrixRT<T>& MatrixRT<T>::scale(T x, T y, T z){
  this->reset(0) ;
  this->m[0] = x ;
  this->m[5] = y ;
  this->m[10] = z ;
  this->m[15] = 1.0 ;
  return *this ;
}


    error.fatal() ;
  }
  T *a,*b ;
  a = this->m-1 ;
  b = M[0] - 1 ;
  for(int i=0;i<16;++i){
    *(++a) = *(++b) ;

template <class T>
MatrixRT<T>& MatrixRT<T>::operator=(const MatrixRT<T>& M) {
  T *a,*b ;
  a = this->m-1 ;
  b = M.m - 1 ;
  for(int i=0;i<16;++i){
    *(++a) = *(++b) ;
++ nurbs/nurbs.cpp

  T du,dv ;
  // compute a coarse distance for the curve
  Point_nD<T,N> a,b,c ;
  a = this->pointAt(0.0) ;
  b = this->pointAt(0.5) ;
  c = this->pointAt(1.0) ;

  T distance = norm(b-a) + norm(c-b) ;

++ nurbs/nurbsS.cpp

  // 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() = this->extremum(1,coordX) ;
  minP.y() = this->extremum(1,coordY) ;
  minP.z() = this->extremum(1,coordZ) ;
  maxP.x() = this->extremum(0,coordX) ;
  maxP.y() = this->extremum(0,coordY) ;
  maxP.z() = this->extremum(0,coordZ) ;

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

  // 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() = this->extremum(1,coordX) ;
  minP.y() = this->extremum(1,coordY) ;
  minP.z() = this->extremum(1,coordZ) ;
  maxP.x() = this->extremum(0,coordX) ;
  maxP.y() = this->extremum(0,coordY) ;
  maxP.z() = this->extremum(0,coordZ) ;

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

  }

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

  Point_nD<T,N> lookAt  ;
  lookAt.x() = (minP.x()+maxP.x())/2.0 ;
++ nurbs/nurbsS_sp.cpp

*/
template <class T, int N>
void NurbsSurfaceSP<T,N>::updateMaxU() {
  if(this->degU>3){
#ifdef USE_EXCEPTION
    throw NurbsInputError();
#else

#endif
  }
  else{
    maxU.resize(this->P.rows()) ;
    maxAtU_.resize(this->P.rows()) ;
    for(int i=0;i<this->P.rows();++i){
      if(!maxInfluence(i,this->U,this->degU,maxAtU_[i]))
	cerr << "Problem in maxInfluence U!\n" ;
      maxU[i] = nurbsBasisFun(maxAtU_[i],i,this->degU,this->U) ;
    }
    
  }

*/
template <class T, int N>
void NurbsSurfaceSP<T,N>::updateMaxV() {
  if(this->degV>3){
#ifdef USE_EXCEPTION
    throw NurbsInputError();
#else

#endif
  }
  else{
    maxV.resize(this->P.cols()) ;
    maxAtV_.resize(this->P.cols()) ;
    for(int i=0;i<this->P.cols();++i){
      if(!maxInfluence(i,this->V,this->degV,maxAtV_[i]))
	cerr << "Problem in maxInfluence V!\n" ;
      maxV[i] = nurbsBasisFun(maxAtV_[i],i,this->degV,this->V) ;
    }
    
  }

NurbsSurfaceSP<T,N> NurbsSurfaceSP<T,N>::generateParallel(T d) const {
  NurbsSurfaceSP<T,N> p(*this) ;

  Vector< Point_nD<T,N> > offset(this->P.rows()*this->P.cols()) ;
  Vector<T> ur(this->P.rows()*this->P.cols()) ;
  Vector<T> vr(this->P.rows()*this->P.cols()) ;
  Vector_INT Du(this->P.rows()*this->P.cols()) ;
  Vector_INT Dv(this->P.rows()*this->P.cols()) ;

  Du.reset(0) ;
  Dv.reset(0) ;


  no = 0 ;

  for(i=0;i<this->P.rows();++i)
    for(j=0;j<this->P.cols();++j){
      Point_nD<T,N> norm ;
      norm = normal(maxAtU_[i],maxAtV_[j]) ;
      if(norm.x() == T(0) && 

	  norm /= T(2) ;
	  ok = 1 ;
	}
	if(i==this->P.rows()-1 && j==this->P.cols()-1){
	  norm = normal(maxAtU_[i]-delta,maxAtV_[j]) ;
	  norm += normal(maxAtU_[i],maxAtV_[j]-delta) ;
	  norm /= T(2) ;
	  ok = 1 ;
	}
	if(i==0 && j==this->P.cols()-1){
	  norm = normal(maxAtU_[i]-delta,maxAtV_[j]) ;
	  norm += normal(maxAtU_[i],maxAtV_[j]+delta) ;
	  norm /= T(2) ;
	  ok = 1 ;
	}
	if(i==this->P.rows()-1 && j==0){
	  norm = normal(maxAtU_[i]-delta,maxAtV_[j]) ;
	  norm += normal(maxAtU_[i],maxAtV_[j]+delta) ;
	  norm /= T(2) ;

	  while(norm.x() == T(0) && 
	     norm.y() == T(0) &&
	     norm.z() == T(0)){
	    if( nt*d >(this->U[this->U.n()-1]-this->U[0])){
#ifdef USE_EXCEPTION
	      throw NurbsComputationError();
#else

#endif
	    }
	    T u1,u2,v1,v2 ;
	    if(i==0 || i==this->P.rows()-1){
	      u1 = u2 = maxAtU_[i] ;
	      v1 = maxAtV_[j]+ nt*delta ;
	      v2 = maxAtV_[j]- nt*delta ;
	      if(v1>this->V[this->V.n()-1]) v1 = this->V[this->V.n()-1] ;
	      if(v2<this->V[0]) v2 = this->V[0] ;
	      norm = normal(u1,v1);
	      norm += normal(u2,v2) ;
	      norm /= 2 ; 

	      u1 = maxAtU_[i]- nt*delta ;
	      u2 = maxAtU_[i]+ nt*delta ;
	      v1 = v2 = maxAtV_[j] ;
	      if(u1 < this->U[0]) u1 = this->U[0] ;
	      if(u2 > this->U[this->U.n()-1]) u2 = this->U[this->U.n()-1] ;

	      T u3,v3 ;
	      u3 = maxAtU_[i] ;

	      else
		v3 = maxAtV_[j]- nt*delta ;

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

	      norm = normal(u1,v1);
	      norm += normal(u2,v2) ;


  int sizeU, sizeV ;

  sizeU = 2*this->degU+3 ; 
  if(i-this->degU-1<0) sizeU += i-this->degU-1 ; 
  if(i+this->degU+1>=this->P.rows()) sizeU -= i+this->degU+1-this->P.rows() ;

  sizeV = 2*this->degV+3 ;
  if(j-this->degV-1<0) sizeV += j-this->degV-1 ; 
  if(j+this->degV+1>=this->P.cols()) sizeV -= j+this->degV+1-this->P.cols() ;
  
  Vector<T> u(sizeU) ;
  Vector<T> v(sizeV) ;

  int n=0;
  int nu = 0 ;
  int nv = 0 ; 
  for(int k=i-this->degU-1;k<=i+this->degU+1;++k){
    if(k<0)
      continue ;
    if(k>=this->P.rows())
      break ; 
    nv = 0 ;
    for(int l=j-this->degV-1;l<=j+this->degV+1;++l){
      if(l<0)
	continue ;
      if(l>=this->P.cols())
	break ; 
      if( k == i && j==l){
	pts[n].x() = a.x() ; 
++ nurbs/nurbsS_sp.h



  void modSurfCPby(int i, int j, const HPoint_nD<T,N>& a) //!< Moves a surface point by a value
    { this->P(i,j) +=  a / (maxU[i]*maxV[j]) ;  }
  void modSurfCP(int i, int j, const HPoint_nD<T,N>& a) //!< Moves a surface point to a value
    { modSurfCPby(i,j,a-surfP(i,j)) ;  }

++ nurbs/nurbsSub.cpp

  
  /* Allocate storage for the grid of points generated */
  
  CHECK( pts = new SurfSample<T>*[Granularity+1]);
  CHECK( pts[0] = new SurfSample<T>[(Granularity+1)*(Granularity+1)]);
  
  for (i = 1; i <= Granularity; i++)

  
  if (! *alpha)	/* Must allocate alpha */
    {
      CHECK( *alpha = new T*[k+1]);
      for (i = 0; i <= k; i++)
	CHECK( (*alpha)[i] = new T[(m + n + 1)]);
    }
++ nurbs/nurbs_sp.cpp

*/
template <class T, int N>
void NurbsCurveSP<T,N>::updateMaxU() {
  if(this->deg_>3){
#ifdef USE_EXCEPTION
    throw NurbsInputError();
#else

#endif
  }
  else{
    maxU.resize(this->P.n()) ;
    maxAt_.resize(this->P.n()) ;
    for(int i=0;i<this->P.n();++i){
      if(!maxInfluence(i,this->U,this->deg_,maxAt_[i]))
	cerr << "Problem in maxInfluence U!\n" ;
      if(i>0)
	if(maxAt_[i]<maxAt_[i-1]){

	  error.fatal() ; 
#endif
	}
      maxU[i] = basisFun(maxAt_[i],i,this->deg_) ;
    }
    
  }

*/
template <class T, int N>
void NurbsCurveSP<T,N>::modOnlySurfCPby(int i, const HPoint_nD<T,N>& a){
  Vector<T> u(2*this->deg_+3) ;
  Vector< Point_nD<T,N> > pts(2*this->deg_+3) ; 

  int n=0;
  for(int j=i-this->deg_-1;j<=i+this->deg_+1;++j){
    if(j<0)
      continue ;
    if(j>=this->P.n())
      break ; 
    u[n] = maxAt_[j] ;
    if( j == i){
++ nurbs/nurbs_sp.h

  int read(ifstream &fin) ;

  void modSurfCPby(int i, const HPoint_nD<T,N>& a) 
    { this->P[i] +=  a / maxU[i] ;  }
  void modSurfCP(int i, const HPoint_nD<T,N>& a) 
    { modSurfCPby(i,a-surfP(i)) ;  }


Return to bug 120303