1 /*=============================================================================
2         File: hnurbsS.cpp
3      Purpose:
4     Revision: $Id: hnurbsS.cpp,v 1.3 2002/05/17 18:24:21 philosophil Exp $
5   Created by: Philippe Lavoie    (7 October 1997)
6  Modified by:
7 
8  Copyright notice:
9           Copyright (C) 1996-1998 Philippe Lavoie
10 
11           This library is free software; you can redistribute it and/or
12           modify it under the terms of the GNU Library General Public
13           License as published by the Free Software Foundation; either
14           version 2 of the License, or (at your option) any later version.
15 
16           This library is distributed in the hope that it will be useful,
17           but WITHOUT ANY WARRANTY; without even the implied warranty of
18           MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
19           Library General Public License for more details.
20 
21           You should have received a copy of the GNU Library General Public
22           License along with this library; if not, write to the Free
23           Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
24 =============================================================================*/
25 #include <hnurbsS.h>
26 #include <string.h>
27 
28 /*!
29  */
30 namespace PLib {
31 
32 /*!
33   \brief The basic constructor
34 
35   \author Philippe Lavoie
36   \date 7 October 1997
37 */
38 template <class T, int N>
HNurbsSurface()39 HNurbsSurface<T,N>::HNurbsSurface():NurbsSurface<T,N>(), fixedOffset(0){
40 
41   baseLevel_ = 0 ;
42   nextLevel_ = 0 ;
43   firstLevel_ = lastLevel_ = this ;
44 
45   rU.resize(0) ;
46   rV.resize(0) ;
47   updateN = 0 ;
48   baseUpdateN = 0 ;
49 
50   level_ = 0 ;
51 }
52 
53 /*!
54   \brief Constructor with a base level
55 
56   \param base  the base level
57   \param xU  the U knots to insert in this level
58   \param xV  the V knots to insert in this level
59 
60   \warning The base pointer must be pointing to a valid object
61 
62   \author Philippe Lavoie
63   \date 7 October 1997
64 */
65 template <class T, int N>
HNurbsSurface(HNurbsSurface<T,N> * base,const Vector<T> & xU,const Vector<T> & xV)66 HNurbsSurface<T,N>::HNurbsSurface(HNurbsSurface<T,N>* base, const Vector<T>& xU, const Vector<T>& xV):NurbsSurface<T,N>(),fixedOffset(0){
67 
68   if(!base){
69     Error error("HNurbsSurface<T,N> constructor") ;
70     error << "Initializing a HNurbsSurface<T,N> with a null base pointer!" ;
71     error.fatal() ;
72   }
73   if(base->nextLevel_){
74     Error error("HNurbsSurface<T,N> constructor") ;
75     error << "You're trying to replace an existing level, this is not allowed." ;
76     error.fatal() ;
77   }
78 
79 
80   nextLevel_ = 0 ;
81   firstLevel_ = base->firstLevel_ ;
82   lastLevel_ = this ;
83   baseLevel_ = (HNurbsSurface<T,N>*)base ;
84 
85   // update the information in the previous levels
86   baseLevel_->nextLevel_ = this ;
87   HNurbsSurface<T,N>* l ;
88   l = baseLevel_ ;
89   while(l){
90     l->lastLevel_ = this ;
91     l = l->baseLevel_ ;
92   }
93 
94 
95   level_ = base->level_ + 1 ;
96 
97   rU = xU ;
98   rV = xV ;
99 
100   updateN = 0 ;
101   baseUpdateN = baseLevel_->modifiedN()-1 ; // Set it so that initBaseLevel will run
102 
103   initBase() ;
104   offset.resize(baseSurf.ctrlPnts()) ;
105 
106   this->P = baseSurf.ctrlPnts() ;
107   this->U = baseSurf.knotU() ;
108   this->V = baseSurf.knotV() ;
109   this->degU = baseSurf.degreeU() ;
110   this->degV = baseSurf.degreeV() ;
111 
112   //updateSurface() ;
113 
114 }
115 
116 /*!
117   \brief Constructor with a base level
118 
119   \param base  the base level
120   \param xU  the U knots to insert in this level
121   \param xV  the V knots to insert in this level
122 
123   \warning The base pointer must be pointing to a valid object
124 
125   \author Philippe Lavoie
126   \date 7 October 1997
127 */
128 template <class T, int N>
HNurbsSurface(HNurbsSurface<T,N> * base)129 HNurbsSurface<T,N>::HNurbsSurface(HNurbsSurface<T,N>* base):NurbsSurface<T,N>(),fixedOffset(0){
130 
131   if(!base){
132     Error error("HNurbsSurface<T,N> constructor") ;
133     error << "Initializing a HNurbsSurface<T,N> with a null base pointer!" ;
134     error.fatal() ;
135   }
136   if(base->nextLevel_){
137     Error error("HNurbsSurface<T,N> constructor") ;
138     error << "You're trying to replace an existing level, this is not allowed." ;
139     error.fatal() ;
140   }
141 
142   nextLevel_ = 0 ;
143   firstLevel_ = base->firstLevel_ ;
144   lastLevel_ = this ;
145   baseLevel_ = (HNurbsSurface<T,N>*)base ;
146 
147   // update the information in the previous levels
148   baseLevel_->nextLevel_ = this ;
149   HNurbsSurface<T,N>* l ;
150   l = baseLevel_ ;
151   while(l){
152     l->lastLevel_ = this ;
153     l = l->baseLevel_ ;
154   }
155 
156   level_ = base->level_ + 1 ;
157 
158   updateN = 0 ;
159   rU.resize(0) ;
160   rV.resize(0) ;
161 
162   baseUpdateN = baseLevel_->modifiedN()-1 ; // Set it so that initBase will run
163   initBase() ;
164   offset.resize(baseSurf.ctrlPnts()) ;
165   this->P = baseSurf.ctrlPnts() ;
166   this->U = baseSurf.knotU() ;
167   this->V = baseSurf.knotV() ;
168   this->degU = baseSurf.degreeU() ;
169   this->degV = baseSurf.degreeV() ;
170   //updateSurface() ;
171 
172 }
173 
174 /*!
175   \brief Constructs a base HNURBS
176 
177   Constructs a base HNURBS. This HNURBS surface is set to level 0.
178   And it corresponds to the NURBS surface.
179 
180   This constructor does not transform the NURBS surface into
181   a HNURBS surface. It only copies the values from the NURBS
182   surface as it's base offset values.
183 
184   \param S  the NURBS surface at level 0
185 
186   \author Philippe Lavoie
187   \date 7 October 1997
188 */
189 template <class T, int N>
HNurbsSurface(const NurbsSurface<T,N> & S)190 HNurbsSurface<T,N>::HNurbsSurface(const NurbsSurface<T,N>& S):NurbsSurface<T,N>(S),fixedOffset(0){
191 
192   baseLevel_ = 0 ;
193   nextLevel_ = 0 ;
194   firstLevel_ = lastLevel_ = this ;
195 
196   level_ = 0 ;
197   updateN = 0 ;
198   baseUpdateN = 0 ;
199 
200   rU.resize(0) ;
201   rV.resize(0) ;
202 
203   offset = this->P ;
204 }
205 
206 /*!
207   \brief The copy constructor
208 
209   \param S  the HNURBS surface to copy
210 
211   \author Philippe Lavoie
212   \date 7 October 1997
213 */
214 template <class T, int N>
HNurbsSurface(const HNurbsSurface<T,N> & S)215 HNurbsSurface<T,N>::HNurbsSurface(const HNurbsSurface<T,N>& S):NurbsSurface<T,N>(S),fixedOffset(0) {
216 
217   baseLevel_ = nextLevel_ = 0 ;
218   firstLevel_ = lastLevel_ = this ;
219 
220   level_ = 0 ;
221   updateN = 0 ;
222   baseUpdateN = 0 ;
223 
224   rU = S.rU ;
225   rV = S.rV ;
226 
227   offset = S.offset ;
228 
229   this->copy(S) ;
230 
231 }
232 
233 /*!
234   \brief A level constructor
235 
236   \param base the base of this level
237   \param the values for this new level
238 
239   \warning The base pointer must be pointing to a valid object
240 
241   \author Philippe Lavoie
242   \date 7 October 1997
243 */
244 template <class T, int N>
HNurbsSurface(HNurbsSurface<T,N> * base,const HNurbsSurface<T,N> & surf)245 HNurbsSurface<T,N>::HNurbsSurface(HNurbsSurface<T,N> *base, const HNurbsSurface<T,N> &surf):NurbsSurface<T,N>(surf),fixedOffset(0) {
246 
247   if(!base){
248     Error error("HNurbsSurface<T,N> constructor") ;
249     error << "Initializing a HNurbsSurface<T,N> with a null base pointer!" ;
250     error.fatal() ;
251   }
252   if(base->nextLevel_){
253     Error error("HNurbsSurface<T,N> constructor") ;
254     error << "You're trying to replace an existing level, this is not allowed." ;
255     error.fatal() ;
256   }
257 
258   baseLevel_ = (HNurbsSurface<T,N>*)base ;
259   nextLevel_ = 0 ;
260   lastLevel_ = this ;
261 
262   // update the information in the previous levels
263   baseLevel_->nextLevel_ = this ;
264   HNurbsSurface<T,N>* l ;
265   l = baseLevel_ ;
266   while(l){
267     l->lastLevel_ = this ;
268     l = l->baseLevel_ ;
269   }
270 
271   firstLevel_ = base->firstLevel_ ;
272   level_ = base->level_+1 ;
273   baseUpdateN = base->modifiedN()-1 ;
274 
275   initBase() ;
276 
277   updateN = 0 ;
278 
279   this->copy(surf) ;
280 
281 }
282 
283 /*!
284   \brief Copies a HNurbs Surface and all it children
285 
286   \param ns  the HNurbs surface to copy
287 
288   \author Philippe Lavoie
289   \date 7 October 1997
290 */
291 template <class T, int N>
copy(const HNurbsSurface<T,N> & nS)292 void HNurbsSurface<T,N>::copy(const HNurbsSurface<T,N>& nS){
293   HNurbsSurface<T,N> *levelP ;
294   levelP = nS.nextLevel() ;
295 
296   NurbsSurface<T,N>::operator=(nS) ;
297   rU = nS.rU ;
298   rV = nS.rV ;
299   offset = nS.offset ;
300   fixedOffset= nS.fixedOffset ;
301 
302   firstLevel_ = this ;
303 
304   if(levelP){
305     HNurbsSurface<T,N> *newLevel = new HNurbsSurface(this,*levelP) ;
306     nextLevel_ = newLevel ;
307     lastLevel_ = nextLevel_->lastLevel() ;
308   }
309   else{
310     lastLevel_ = this ;
311   }
312 
313 }
314 
315 /*!
316   \brief updates the NURBS surface
317 
318   Updates the NURBS surface according to the offset values
319   and its base level. You can update only one control point
320   from the surface if you specify a value for i and j or you
321   can update all the points if i0 or j0 is below 0.
322 
323   \param i0  the row of the control point to update
324   \param j0  the column of the control point to update
325 
326   \author Philippe Lavoie
327   \date 7 October 1997
328 */
329 template <class T, int N>
updateSurface(int i0,int j0)330 void HNurbsSurface<T,N>::updateSurface(int i0, int j0){
331   if(i0>=0 && j0>=0){
332     if(offset(i0,j0).x()==0.0 && offset(i0,j0).y()==0.0 && offset(i0,j0).z()==0.0)
333       return ;
334   }
335   if(baseLevel_){
336     if(initBase()){
337       this->P = baseSurf.ctrlPnts() ;
338       this->U = baseSurf.knotU() ;
339       this->V = baseSurf.knotV() ;
340       this->degU = baseSurf.degreeU() ;
341       this->degV = baseSurf.degreeV() ;
342     }
343     if(i0>=0 && j0>=0){
344       Point_nD<T,N> vecOffset ;
345       if(fixedOffset){
346 	vecOffset = offset(i0,j0).x()*ivec(0,0) +
347 	  offset(i0,j0).y()*jvec(0,0) +
348 	  offset(i0,j0).z()*kvec(0,0) ;
349       }
350       else{
351 	vecOffset = offset(i0,j0).x()*ivec(i0,j0) +
352 	  offset(i0,j0).y()*jvec(i0,j0) +
353 	  offset(i0,j0).z()*kvec(i0,j0) ;
354       }
355       this->P(i0,j0).x() = baseSurf.ctrlPnts()(i0,j0).x()+vecOffset.x() ;
356       this->P(i0,j0).y() = baseSurf.ctrlPnts()(i0,j0).y()+vecOffset.y() ;
357       this->P(i0,j0).z() = baseSurf.ctrlPnts()(i0,j0).z()+vecOffset.z() ;
358     }
359     else{
360       for(int i=0;i<this->P.rows();++i)
361 	for(int j=0;j<this->P.cols();++j){
362 	  if(offset(i,j).x() != 0 ||
363 	     offset(i,j).y() != 0 || offset(i,j).z() || 0){
364 	    Point_nD<T,N> vecOffset ;
365 	    if(fixedOffset){
366 	      vecOffset = offset(i,j).x()*ivec(0,0) +
367 		offset(i,j).y()*jvec(0,0) +
368 		offset(i,j).z()*kvec(0,0) ;
369 	    }
370 	    else{
371 	      vecOffset = offset(i,j).x()*ivec(i,j) +
372 		offset(i,j).y()*jvec(i,j) +
373 		offset(i,j).z()*kvec(i,j) ;
374 	    }
375 	    this->P(i,j).x() = baseSurf.ctrlPnts()(i,j).x()+vecOffset.x() ;
376 	    this->P(i,j).y() = baseSurf.ctrlPnts()(i,j).y()+vecOffset.y() ;
377 	    this->P(i,j).z() = baseSurf.ctrlPnts()(i,j).z()+vecOffset.z() ;
378 	  }
379 	}
380     }
381   }
382   else{
383     if(i0>=0 && j0>=0)
384       this->P(i0,j0) = offset(i0,j0) ;
385     else{
386       for(int i=0;i<this->P.rows();++i)
387 	for(int j=0;j<this->P.cols();++j){
388 	  this->P(i,j) = offset(i,j) ;
389 	}
390     }
391   }
392 
393   ++updateN ;
394 }
395 
396 
397 /*!
398   \brief Initialize the base surface
399 
400   Initialize the base surface from the previous level if it has
401   been modified.
402 
403   \param force  if set, this forces an update of the base surface
404 
405   \return 1 if the base surface is modified, 0 otherwise.
406 
407   \author Philippe Lavoie
408   \date 15 April 1998
409 */
410 template <class T, int N>
initBase(int force)411 int HNurbsSurface<T,N>::initBase(int force){
412   if(!baseLevel_){
413     return 0; // this means that the surface is the base level
414   }
415 
416   // make sure none of the base level need updating
417   if(baseLevel_->initBase())
418     force = 1 ;
419 
420 
421   // Only initialize the base level if it isn't up to date
422   if(baseLevel_->modifiedN() == baseUpdateN && !force){
423     return 0;
424   }
425 
426   baseUpdateN = baseLevel_->modifiedN() ;
427 
428   baseSurf = *baseLevel_ ;
429 
430   if(rU.n()>0)
431     baseSurf.refineKnotU(rU) ;
432   if(rV.n()>0)
433     baseSurf.refineKnotV(rV) ;
434 
435   Vector<T> maxU,maxV ;
436   int i,j ;
437 
438   if(baseSurf.degreeU()>3)
439     averagingKnots(baseSurf.knotU(),baseSurf.degreeU(),maxU) ;
440   else{
441     maxU.resize(baseSurf.ctrlPnts().rows()) ;
442     for(i=0;i<baseSurf.ctrlPnts().rows();++i)
443       if(!maxInfluence(i,baseSurf.knotU(),baseSurf.degreeU(),maxU[i]))
444 	cerr << "Problem in maxInfluence U!\n" ;
445   }
446 
447   if(baseSurf.degreeV()>3)
448     averagingKnots(baseSurf.knotV(),baseSurf.degreeV(),maxV) ;
449   else{
450     maxV.resize(baseSurf.ctrlPnts().cols()) ;
451     for(i=0;i<baseSurf.ctrlPnts().cols();++i)
452       if(!maxInfluence(i,baseSurf.knotV(),baseSurf.degreeV(),maxV[i]))
453 	cerr << "Problem in maxInfluence V!\n" ;
454   }
455 
456   if(fixedOffset){
457     if(ivec.rows() != 1 || ivec.cols()!=1){
458       ivec.resize(1,1);
459       jvec.resize(1,1);
460       kvec.resize(1,1);
461     }
462   }
463   else{
464     ivec.resize(maxU.n(),maxV.n()) ;
465     jvec.resize(maxU.n(),maxV.n()) ;
466     kvec.resize(maxU.n(),maxV.n()) ;
467 
468     Matrix< Point_nD<T,N> > ders ;
469 
470     for(i=0;i<maxU.n();++i)
471       for(j=0;j<maxV.n();++j){
472 	baseSurf.deriveAt(maxU[i],maxV[j],1,ders) ;
473 	// It is possible that there are no valid derivative at that point
474 	// we then move to a nearby point to get a valid value
475 	Point_nD<T,N> norm = crossProduct(ders(0,1),ders(1,0)) ;
476 	if(norm.x() == T(0) &&
477 	   norm.y() == T(0) &&
478 	   norm.z() == T(0)){
479 	  const T delta = 0.00001 ;
480 	  T nt = 1.0 ;
481 	  Matrix< Point_nD<T,N> > dersT(ders) ;
482 	  while(norm.x() == T(0) &&
483 		norm.y() == T(0) &&
484 		norm.z() == T(0)){
485 	    if( nt*delta >(baseSurf.knotU()[baseSurf.knotU().n()-1]-baseSurf.knotU()[0])){
486 	      Error error("initBase");
487 	      error << "Can't compute the derivative.\n" ;
488 	      error.fatal() ;
489 	    }
490 	    T nd ;
491 
492 	    ders.reset(0) ;
493 
494 	    nd = 0 ;
495 	    if(i != 0){
496 	      baseSurf.deriveAt(maxU[i] - nt*delta, maxV[j],1, dersT) ;
497 	      ders += dersT ;
498 	      ++nd ;
499 	    }
500 	    if(i != maxU.n()-1){
501 	      baseSurf.deriveAt(maxU[i] + nt*delta, maxV[j],1, dersT) ;
502 	      ders += dersT ;
503 	      ++nd ;
504 	    }
505 	    if(j != 0){
506 	      baseSurf.deriveAt(maxU[i], maxV[j] - nt*delta,1, dersT) ;
507 	      ders += dersT ;
508 	      ++nd ;
509 	    }
510 	    if(j != maxV.n()-1){
511 	      baseSurf.deriveAt(maxU[i], maxV[j] + nt*delta,1, dersT) ;
512 	      ders += dersT ;
513 	      ++nd ;
514 	    }
515 
516 	    if(nd==0){
517 	      Error error("initBase");
518 	      error << "Can't compute the derivative.\n" ;
519 	      error.fatal() ;
520 	    }
521 
522 	    ders /= nd ;
523 
524 	    norm = crossProduct(ders(0,1),ders(1,0)) ;
525 	    nt *= 10.0 ;
526 	  }
527 	}
528 
529 	ivec(i,j) = ders(0,1).unitLength() ;
530 	kvec(i,j) = crossProduct(ders(0,1),ders(1,0)).unitLength() ;
531 	jvec(i,j) = crossProduct(kvec(i,j),ivec(i,j)).unitLength() ;
532       }
533   }
534 
535   return 1 ;
536 }
537 
538 /*!
539   \brief Specifies the level that modifies the point
540 
541   Specifies what level modifies the point \a (u,v)
542 
543   \param u  the \a u parametric value
544   \param v  the \a v parametric value
545 
546   \author Philippe Lavoie
547   \date 7 October 1997
548 */
549 template <class T, int N>
modifies(T u,T v)550 int HNurbsSurface<T,N>::modifies(T u, T v){
551   if(nextLevel_){
552     int mod = nextLevel_->modifies(u,v) ;
553     if(mod>=0)
554       return mod ;
555   }
556 
557   if(u<this->knotU()[0] || u>this->knotU()[this->knotU().n()-1])
558     return -1 ;
559   if(v<this->knotV()[0] || v>this->knotU()[this->knotV().n()-1])
560     return -1 ;
561 
562   int su = this->findSpanU(u) ;
563   int sv = this->findSpanV(v) ;
564 
565   for(int i=0;i<=this->degU;++i)
566     for(int j=0;j<=this->degV;++j){
567       if(offset(su-this->degU+i,sv+this->degV+j) != HPoint_nD<T,N>(0,0,0,0))
568 	return level_ ;
569     }
570 
571   return -1 ;
572 }
573 
574 /*!
575   \brief finds the homogenous point at (u,v) for a certain level of detail.
576 
577   \param u  the \a u parametric value
578   \param  v  the \a v parametric value
579   \param  lod  the level of detail3
580 
581   \return the homogenous point at (u,v) and at the level of detail specivied by
582                lod.
583 
584   \author Philippe Lavoie
585   \date 7 October 1997
586 */
587 template <class T, int N>
hpointAt(T u,T v,int lod) const588 HPoint_nD<T,N> HNurbsSurface<T,N>::hpointAt(T u, T v, int lod) const {
589   if(level_==lod)
590     return NurbsSurface<T,N>::operator()(u,v) ;
591 
592   if(nextLevel_)
593     return nextLevel_->hpointAt(u,v,lod) ;
594 
595   return NurbsSurface<T,N>::operator()(u,v) ;
596 }
597 
598 /*!
599   \brief the maximum level of detail
600 
601   Finds the maximum level of detail available from this HNURBS
602   surface
603 
604   \return the maximum level of detail.
605 
606   \author Philippe Lavoie
607   \date 7 October 1997
608 */
609 template <class T, int N>
maxLevel() const610 int HNurbsSurface<T,N>::maxLevel() const{
611   if(lastLevel_)
612     return lastLevel_->level_ ;
613   return level_ ;
614 }
615 
616 /*!
617   \brief Finds the derivative of the point \a (u,v)
618 
619   Computes the matrix of derivatives at \a (u,v) .
620   The value of skl(k,l) represents the
621   derivative of the surface $S(u,v)$ with respect to
622   \a u, \a k times and to \a v, \a l times.
623 
624   \param u  the \a u parametric value
625   \param v  the \a v parametric value
626   \param ders  the Matrix of derivatives.
627 
628   \author Philippe Lavoie
629   \date 7 October 1997
630 */
631 template <class T, int N>
deriveAtH(T u,T v,int d,Matrix<HPoint_nD<T,N>> & ders,int lod) const632 void HNurbsSurface<T,N>::deriveAtH(T u, T v, int d, Matrix< HPoint_nD<T,N> >& ders,int lod) const {
633   if(level_==lod){
634     NurbsSurface<T,N>::deriveAtH(u,v,d,ders) ;
635     return ;
636   }
637 
638   if(nextLevel_){
639     nextLevel_->deriveAtH(u,v,d,ders,lod) ;
640     return ;
641   }
642 
643   NurbsSurface<T,N>::deriveAtH(u,v,d,ders) ;
644 }
645 
646 /*!
647   \brief Finds the derivative of the point \a (u,v)
648 
649   Computes the matrix of derivatives at \a u,v.
650   The value of skl(k,l) represents the
651   derivative of the surface \a S(u,v) with respect to
652   \a u, \a k times and to \a v, \a l times.
653 
654   \param u  the \a u parametric value
655   \param v  the \a v parametric value
656   \param ders  the Matrix of derivatives.
657 
658   \author Philippe Lavoie
659   \date 7 October 1997
660 */
661 template <class T, int N>
deriveAt(T u,T v,int d,Matrix<Point_nD<T,N>> & ders,int lod) const662 void HNurbsSurface<T,N>::deriveAt(T u, T v, int d, Matrix< Point_nD<T,N> >& ders, int lod) const {
663   if(level_==lod){
664     NurbsSurface<T,N>::deriveAt(u,v,d,ders) ;
665     return ;
666   }
667 
668   if(nextLevel_){
669     nextLevel_->deriveAt(u,v,d,ders,lod) ;
670     return ;
671 
672   NurbsSurface<T,N>::deriveAt(u,v,d,ders) ;
673   }
674 }
675 
676 /*!
677   \brief Destructor
678 
679   Deletes all the levels.
680 
681   \author Philippe Lavoie
682   \date 7 October 1997
683 */
684 template <class T, int N>
~HNurbsSurface()685 HNurbsSurface<T,N>::~HNurbsSurface(){
686   if(nextLevel_)
687     delete nextLevel_ ;
688   lastLevel_ = 0 ;
689   if(baseLevel_){
690     baseLevel_->nextLevel_ = 0 ;
691     baseLevel_->lastLevel_ = baseLevel_ ;
692   }
693   baseLevel_ = 0 ;
694   nextLevel_ = 0 ;
695   firstLevel_ = 0 ;
696 }
697 
698 /*!
699   \brief Update the surface for all the levels
700 
701   \param upLevel  updates the levels up to this level of detail
702 
703   \author Philippe Lavoie
704   \date 7 October 1997
705 */
706 template <class T, int N>
updateLevels(int upLevel)707 void HNurbsSurface<T,N>::updateLevels(int upLevel){
708   if(upLevel>=0){
709     if(level_ <= upLevel){
710       this->updateSurface() ;
711     }
712   }
713   else{
714     this->updateSurface() ;
715   }
716 
717   if(upLevel>level() || upLevel<0){
718     if(nextLevel_)
719       nextLevel_->updateLevels(upLevel) ;
720   }
721 }
722 
723 /*!
724   \brief Insert n knots betwen each knots
725 
726   Insert nu knots betwen each knots in the U vector and
727   nv knots between each knots in the V vector.
728 
729   This does not perform a split. It just generates
730   a suitable rU and rV vector. It is suggested that
731   splitting should be done for the level above, not
732   the local level.
733 
734   \param nu  the number of new knots between each knots in U.
735   \param nv  the number of new knots between each knots in V.
736   \param nU  the new refinement knot vector in U
737   \param nV  the new refinement knot vector in V
738 
739   \author Philippe Lavoie
740   \date 7 October 1997
741 */
742 template <class T, int N>
splitUV(int nu,int nv,Vector<T> & nU,Vector<T> & nV)743 void HNurbsSurface<T,N>::splitUV(int nu, int nv, Vector<T> &nU, Vector<T> &nV){
744 
745   nU.resize(this->knotU().n()*nu) ;
746   nV.resize(this->knotV().n()*nv) ;
747 
748   int i,j,n ;
749 
750   n = 0 ;
751   for(i=1;i<this->knotU().n();++i){
752     if(this->knotU()[i] >this->knotU()[i-1]){
753       T a = this->knotU()[i-1] ;
754       T b = this->knotU()[i] ;
755 
756 
757       for(j=0;j<nu;++j){
758 	nU[n] = a + (b-a)*T(j+1)/T(nu+1) ;
759 	++n ;
760       }
761     }
762   }
763   nU.resize(n) ;
764 
765   n = 0 ;
766   for(i=1;i<this->knotV().n();++i){
767     if(this->knotV()[i] > this->knotV()[i-1]){
768       T a = this->knotV()[i-1] ;
769       T b = this->knotV()[i] ;
770 
771       for(j=0;j<nv;++j){
772 	nV[n] = a + (b-a)*T(j+1)/T(nv+1) ;
773 	++n ;
774       }
775     }
776   }
777   nV.resize(n) ;
778 
779 }
780 
781 /*!
782   \brief Insert n knots betwen each knots
783 
784   Insert nu knots betwen each knots in the U vector and
785   nv knots between each knots in the V vector.
786 
787   This doesn't not perform a split. It just generates
788   a suitable rU and rV vector. It is suggested that
789   splitting should be done for the level above, not
790   the local level.
791 
792   \param nu  the number of new knots between each knots in U.
793   \param su  the multiplicity of the each new knot in U
794   \param nv  the number of new knots between each knots in V.
795   \param  sv  the multiplicity of the each new knot in V
796   \param nU  the new refinement knot vector in U
797   \param nV  the new refinement knot vector in V
798 
799   \author Philippe Lavoie
800   \date 7 October 1997
801 */
802 template <class T, int N>
splitUV(int nu,int su,int nv,int sv,Vector<T> & nU,Vector<T> & nV)803 void HNurbsSurface<T,N>::splitUV(int nu, int su, int nv, int sv, Vector<T> &nU, Vector<T> &nV){
804 
805   int i,j,n ;
806 
807   if(su<=0)
808     su = this->degU  ;
809   if(sv<=0)
810     sv = this->degV  ;
811   if(su>this->degU+1)
812     su = this->degU+1 ;
813   if(sv>this->degV+1)
814     sv = this->degV+1 ;
815 
816   nU.resize(this->knotU().n()*nu*su) ;
817   nV.resize(this->knotV().n()*nv*sv) ;
818 
819   n = 0 ;
820   for(i=1;i<this->knotU().n();++i){
821     if(this->knotU()[i] >this->knotU()[i-1]){
822       T a = this->knotU()[i-1] ;
823       T b = this->knotU()[i] ;
824 
825 
826       for(j=0;j<nu;++j){
827 	T u = a + (b-a)*T(j+1)/T(nu+1) ;
828 	for(int l=0;l<su;++l){
829 	  nU[n] = u ;
830 	  ++n ;
831 	}
832       }
833     }
834   }
835   nU.resize(n) ;
836 
837   n = 0 ;
838   for(i=1;i<this->knotV().n();++i){
839     if(this->knotV()[i] > this->knotV()[i-1]){
840       T a = this->knotV()[i-1] ;
841       T b = this->knotV()[i] ;
842 
843       for(j=0;j<nv;++j){
844 	T v = a + (b-a)*T(j+1)/T(nv+1) ;
845 	for(int l=0;l<sv;++l){
846 	  nV[n] = v ;
847 	  ++n ;
848 	}
849       }
850     }
851   }
852   nV.resize(n) ;
853 
854 }
855 
856 /*!
857   \brief Adds a level to this HNURBS surface
858 
859   \param n  the number of new knots between each knots.
860 
861   \return a pointer to the new level or 0 if there was an error.
862 
863   \warning returns 0, if there is already a nextlevel.
864 
865   \author Philippe Lavoie
866   \date 7 October 1997
867 */
868 template <class T, int N>
addLevel(int n)869 HNurbsSurface<T,N>* HNurbsSurface<T,N>::addLevel(int n) {
870   HNurbsSurface<T,N> *newLevel ;
871 
872   if(nextLevel_)
873     return 0 ;
874 
875   Vector<T> newU,newV ;
876 
877   splitUV(n,n,newU,newV) ;
878 
879   newLevel = new HNurbsSurface(this,newU,newV) ;
880 
881   return newLevel ;
882 }
883 
884 /*!
885   \brief Adds a level to this HNURBS surface
886 
887   \param n  the number of new knots between each knots.
888 
889   \return a pointer to the new level or 0 if there was an error.
890 
891   \warning returns 0, if there is already a nextlevel.
892 
893   \author Philippe Lavoie
894   \date 7 October 1997
895 */
896 template <class T, int N>
addLevel()897 HNurbsSurface<T,N>* HNurbsSurface<T,N>::addLevel() {
898   HNurbsSurface<T,N> *newLevel ;
899 
900   if(nextLevel_)
901     return 0 ;
902 
903   newLevel = new HNurbsSurface(this) ;
904 
905   return newLevel ;
906 }
907 
908 /*!
909   \brief generates an iso curve in the U direction
910 
911   Generates an iso-parametric curve which goes through the
912   parametric value u along the U direction.
913 
914   \param u  the U parametric value
915   \param c  the iso-parametric curve
916   \param lod  the level of detail to draw the curve with
917 
918   \return 0 if an error occured, 1 otherwise.
919 
920   \warning the parametric value \a u must be in a valid range
921 
922   \author Philippe Lavoie
923   \date 7 October, 1997
924 */
925 template <class T, int N>
isoCurveU(T u,NurbsCurve<T,N> & c,int lod) const926 int HNurbsSurface<T,N>::isoCurveU(T u, NurbsCurve<T,N>& c,int lod) const {
927   if(lod>=0 && level_>lod)
928     return 0;
929 
930   if(lod==level_ || lod<0){
931     NurbsSurface<T,N>::isoCurveU(u,c) ;
932     return 1 ;
933   }
934 
935   if(nextLevel_){
936     return nextLevel_->isoCurveU(u,c,lod) ;
937   }
938 
939   return 0 ;
940 }
941 
942 /*!
943   \brief generates an iso curve in the V direction
944 
945   Generates an iso-parametric curve which goes through the
946   parametric value v along the V direction.
947 
948   \param v  the V parametric value
949   \param c  the iso-parametric curve
950   \param lod  the level of detail to draw the curve with
951 
952   \return 0 if an error occured, 1 otherwise
953   \warning the parametric value \a v must be in a valid range
954 
955   \author Philippe Lavoie
956   \date 7 October, 1997
957 */
958 template <class T, int N>
isoCurveV(T v,NurbsCurve<T,N> & c,int lod) const959 int HNurbsSurface<T,N>::isoCurveV(T v, NurbsCurve<T,N>& c,int lod) const {
960   if(lod>=0 && level_>lod)
961     return 0;
962 
963   if(lod==level_ || lod<0){
964     NurbsSurface<T,N>::isoCurveV(v,c) ;
965     return 1 ;
966   }
967 
968   if(nextLevel_){
969     return nextLevel_->isoCurveV(v,c,lod) ;
970   }
971 
972   return 0 ;
973 }
974 
975 
976 
977 /*!
978   \brief Read a HNURBS surface from a file stream
979 
980   \param fin  the input file stream
981 
982   \return 0 if an error occurs, 1 otherwise
983 
984   \author Philippe Lavoie
985   \date 7 October 1997
986 */
987 template <class T, int N>
read(ifstream & fin)988 int HNurbsSurface<T,N>::read(ifstream &fin){
989   if(!fin) {
990     return 0 ;
991   }
992   int nu,nv,du,dv;
993   char *type ;
994   type = new char[4] ;
995   if(!fin.read(type,sizeof(char)*4)) { delete []type ; return 0 ;}
996   int r1 = strncmp(type,"hns3",4) ;
997   int r2 = strncmp(type,"hns4",4) ;
998   int r3 = strncmp(type,"hnso",4) ;
999   if(!(r1 || r2 || r3))
1000     return 0 ;
1001   T *p,*p2 ;
1002   char stc ;
1003   if(!r1 || !r2){
1004     if(!fin.read((char*)&stc,sizeof(char))) { delete []type ; return 0 ;}
1005 
1006     int st = stc - '0' ;
1007     if(st != sizeof(T)){ // does not have the same data size
1008       delete []type ;
1009       return 0 ;
1010     }
1011 
1012     if(!fin.read((char*)&nu,sizeof(int))) { delete []type ; return 0 ;}
1013     if(!fin.read((char*)&nv,sizeof(int))) { delete []type ; return 0 ;}
1014     if(!fin.read((char*)&du,sizeof(int))) { delete []type ; return 0 ;}
1015     if(!fin.read((char*)&dv,sizeof(int))) { delete []type ; return 0 ;}
1016 
1017     this->resize(nu,nv,du,dv) ;
1018 
1019     if(!fin.read((char*)this->U.memory(),sizeof(T)*this->U.n())) { delete []type ; return 0 ;}
1020     if(!fin.read((char*)this->V.memory(),sizeof(T)*this->V.n())) { delete []type ; return 0 ;}
1021 
1022     if(!r1){
1023       p = new T[3*nu*nv] ;
1024       if(!fin.read((char*)p,sizeof(T)*3*nu*nv)) { delete []type ; return 0 ;}
1025       p2 = p ;
1026       for(int i=0;i<nu;i++)
1027 	for(int j=0;j<nv;j++){
1028 	  this->P(i,j).x() = *(p++) ;
1029 	  this->P(i,j).y() = *(p++) ;
1030 	  this->P(i,j).z() = *(p++) ;
1031 	  this->P(i,j).w() = 1.0 ;
1032 	}
1033       delete []p2 ;
1034     }
1035     else{
1036       p = new T[4*nu*nv] ;
1037       if(!fin.read((char*)p,sizeof(T)*4*nu*nv)) { delete []type ; return 0 ;}
1038       p2 = p ;
1039       for(int i=0;i<nu;i++)
1040 	for(int j=0;j<nv;j++){
1041 	  this->P(i,j).x() = *(p++) ;
1042 	  this->P(i,j).y() = *(p++) ;
1043 	  this->P(i,j).z() = *(p++) ;
1044 	  this->P(i,j).w() = *(p++) ;
1045 	}
1046       delete []p2 ;
1047     }
1048     offset = this->P ;
1049     this->updateSurface() ;
1050   }
1051   else { // reading the offset information
1052     int ru,rv ;
1053     if(!fin.read((char*)&ru,sizeof(int))) { delete []type ; return 0 ;}
1054     if(!fin.read((char*)&rv,sizeof(int))) { delete []type ; return 0 ;}
1055     rU.resize(ru) ;
1056     rV.resize(rv) ;
1057     if(rU.n()>0)
1058       if(!fin.read((char*)rU.memory(),sizeof(T)*rU.n())) { delete []type ; return 0 ;}
1059     if(rV.n()>0)
1060       if(!fin.read((char*)rV.memory(),sizeof(T)*rV.n())) { delete []type ; return 0 ;}
1061 
1062     if(!fin.read((char*)&nu,sizeof(int))) { delete []type ; return 0 ;}
1063     if(!fin.read((char*)&nv,sizeof(int))) { delete []type ; return 0 ;}
1064 
1065     p = new T[4*nu*nv] ;
1066     if(!fin.read((char*)p,sizeof(T)*4*nu*nv)) { delete []type ; return 0 ;}
1067     p2 = p ;
1068     offset.resize(nu,nv) ;
1069     for(int i=0;i<nu;i++)
1070       for(int j=0;j<nv;j++){
1071 	offset(i,j).x() = *(p++) ;
1072 	offset(i,j).y() = *(p++) ;
1073 	offset(i,j).z() = *(p++) ;
1074 	offset(i,j).w() = *(p++) ;
1075       }
1076     delete []p2 ;
1077     --baseUpdateN ;
1078     this->updateSurface() ;
1079   }
1080 
1081   // now we must read all the level information (if any)
1082   char *ptxt ;
1083   ptxt = new char[7] ;
1084   // the following line is used so that purify doesn't give a warning
1085   ptxt[0] = ptxt[1] = ptxt[2] = ptxt[3] = ptxt[4] = ptxt[5] = ptxt[6] = 0 ;
1086   int mark  = fin.tellg() ;
1087   if(!fin.read(ptxt,sizeof(char)*5))
1088     { delete []ptxt ; delete []type ; return 1 ; }
1089   char *pat ;
1090   pat = strstr(ptxt,"level") ;
1091   if(pat){
1092     // insert a new level
1093     HNurbsSurface<T,N> *newPatch = this->addLevel() ;
1094     if(newPatch){
1095       if(!newPatch->read(fin))
1096 	return 0 ;
1097     }
1098     else
1099       return 0 ;
1100   }
1101   else{
1102     // something else is after the HNurbs
1103     // restore were the file was before the read
1104     fin.seekg(mark);
1105   }
1106 
1107   delete []ptxt ;
1108   delete []type ;
1109   return 1 ;
1110 }
1111 
1112 /*!
1113   \brief Reads a HNURBS surface from a file
1114 
1115   \param filename  the filename to read from
1116 
1117   \return 0 if an error occurs, 1 otherwise
1118 
1119   \author Philippe Lavoie
1120   \date 7 October 1997
1121 */
1122 template <class T, int N>
read(const char * filename)1123 int HNurbsSurface<T,N>::read(const char* filename){
1124   ifstream fin(filename) ;
1125   if(!fin) {
1126     return 0 ;
1127   }
1128 
1129   return this->read(fin) ;
1130 }
1131 
1132 /*!
1133   \brief Write a HNURBS surface to a file stream
1134 
1135   \param fout  the output filestream to write to.
1136 
1137   \return 1 on success, 0 on failure
1138 
1139   \author Philippe Lavoie
1140   \date 7 October 1997
1141 */
1142 template <class T, int N>
write(ofstream & fout) const1143 int HNurbsSurface<T,N>::write(ofstream &fout) const {
1144   if(!fout)
1145     return 0 ;
1146   if(!baseLevel_){
1147     int prows = this->P.rows();
1148     int pcols = this->P.cols();
1149     char st = '0' + sizeof(T) ;
1150     if(!fout.write((char*)&"hns4",sizeof(char)*4)) return 0 ;
1151     if(!fout.write((char*)&st,sizeof(char))) return 0 ;
1152     if(!fout.write((char*)&prows,sizeof(int))) return 0 ;
1153     if(!fout.write((char*)&pcols,sizeof(int))) return 0 ;
1154     if(!fout.write((char*)&this->degU,sizeof(int))) return 0 ;
1155     if(!fout.write((char*)&this->degV,sizeof(int))) return 0 ;
1156     if(!fout.write((char*)this->U.memory(),sizeof(T)*this->U.n())) return 0 ;
1157     if(!fout.write((char*)this->V.memory(),sizeof(T)*this->V.n())) return 0 ;
1158 
1159     T *p,*p2 ;
1160     p = new T[this->P.rows()*this->P.cols()*4] ;
1161     p2 = p ;
1162     for(int i=0;i<this->P.rows();i++)
1163       for(int j=0;j<this->P.cols();j++){
1164 	*p = offset(i,j).x() ; p++ ;
1165 	*p = offset(i,j).y() ; p++ ;
1166 	*p = offset(i,j).z() ; p++ ;
1167 	*p = offset(i,j).w() ; p++ ;
1168       }
1169     if(!fout.write((char*)p2,sizeof(T)*this->P.rows()*this->P.cols()*4)) return 0 ;
1170     delete []p2 ;
1171   }
1172   else{
1173     if(!fout.write((char*)&"hnso",sizeof(char)*4)) return 0 ;
1174     int run = rU.n();
1175     int rvn = rV.n();
1176     if(!fout.write((char*)&run,sizeof(int))) return 0 ;
1177     if(!fout.write((char*)&rvn,sizeof(int))) return 0 ;
1178     if(rU.n()>0)
1179       if(!fout.write((char*)rU.memory(),sizeof(T)*rU.n())) return 0 ;
1180     if(rV.n()>0)
1181       if(!fout.write((char*)rV.memory(),sizeof(T)*rV.n())) return 0 ;
1182     int orows = offset.rows();
1183     int ocols = offset.cols() ;
1184     if(!fout.write((char*)&orows,sizeof(int))) return 0 ;
1185     if(!fout.write((char*)&ocols,sizeof(int))) return 0 ;
1186     T *p,*p2 ;
1187 
1188     p = new T[offset.rows()*offset.cols()*4] ;
1189     p2 = p ;
1190     for(int i=0;i<offset.rows();i++)
1191       for(int j=0;j<offset.cols();j++){
1192 	*p = offset(i,j).x() ; p++ ;
1193 	*p = offset(i,j).y() ; p++ ;
1194 	*p = offset(i,j).z() ; p++ ;
1195 	*p = offset(i,j).w() ; p++ ;
1196       }
1197     if(!fout.write((char*)p2,sizeof(T)*offset.rows()*offset.cols()*4))
1198       return 0 ;
1199     delete []p2 ;
1200   }
1201   if(nextLevel_){
1202     if(!fout.write((char*)&"level",sizeof(char)*5)) return 0 ;
1203     if(!nextLevel_->write(fout)) return 0 ;
1204   }
1205   return 1 ;
1206 }
1207 
1208 /*!
1209   \brief write a HNURBS surface to a file
1210 
1211   \param filename  the filename to write to
1212 
1213   \return 1 on success, 0 on failure
1214 
1215   \author Philippe Lavoie
1216   \date 7 October 1997
1217 */
1218 template <class T, int N>
write(const char * filename) const1219 int HNurbsSurface<T,N>::write(const char* filename) const {
1220   ofstream fout(filename) ;
1221   if(!fout)
1222     return 0 ;
1223   return write(fout);
1224 }
1225 
1226 /*!
1227   Generates a NURBS surface.
1228   \param surf <-- the NURBS surface corresponding to the HNURBS surface
1229   \return
1230   \warning
1231 
1232   \author Philippe Lavoie
1233   \date 7 October 1997
1234 */
1235 //template <class T, int N>
1236 //void HNurbsSurface<T,N>::toNurbs(NurbsSurface<T,N>& surf) const {
1237 
1238 //}
1239 
1240 template <class T>
mergeInto(const Vector<T> & a,Vector<T> & b)1241 inline void mergeInto(const Vector<T> &a, Vector<T> &b) {
1242   // merge a in b ;
1243   int n = b.n() ;
1244   b.resize(b.n()+a.n()) ;
1245   for(int i=0;i<a.n();++i)
1246     b[n+i] = a[i] ;
1247   b.qSort() ;
1248 }
1249 
1250 /*!
1251   \brief Refine both knot vectors
1252 
1253   \param nU  the U knot vector to refine from
1254   \param nV  the V knot vector to refine from
1255 
1256   \author Philippe Lavoie
1257   \date 24 January, 1997
1258 */
1259 template <class T, int N>
refineKnots(const Vector<T> & nU,const Vector<T> & nV)1260 void HNurbsSurface<T,N>::refineKnots(const Vector<T>& nU, const Vector<T>& nV){
1261   refineKnotU(nU) ;
1262   mergeInto(nU,rU) ;
1263 
1264   initBase(1) ;
1265 
1266   refineKnotV(nV) ;
1267   mergeInto(nV,rV) ;
1268 }
1269 
1270 /*!
1271   \brief  Refines the U knot vector
1272 
1273   \param X  the knot vector to refine from
1274 
1275   \author Philippe Lavoie
1276   \date 24 January, 1997
1277 */
1278 template <class T, int N>
refineKnotU(const Vector<T> & X)1279 void HNurbsSurface<T,N>::refineKnotU(const Vector<T>& X) {
1280   updateSurface() ;
1281   Vector<T> Xu(X) ;
1282   int i,j ;
1283   j = 0 ;
1284   for(i=0;i<X.n();++i){
1285     if(X[i]>=this->U[0] && X[i]<=this->U[this->U.n()-1]){
1286       Xu[j] = X[i] ;
1287       ++j ;
1288     }
1289   }
1290   Xu.resize(j) ;
1291 
1292   if(Xu.n()>0){
1293     if(nextLevel_){
1294       nextLevel_->refineKnotU(Xu) ;
1295     }
1296 
1297     NurbsSurface<T,N> osurf(this->degU,this->degV,this->U,this->V,offset) ;
1298 
1299     osurf.refineKnotU(Xu) ;
1300 
1301     offset.resize(osurf.ctrlPnts().rows(),osurf.ctrlPnts().cols()) ;
1302     for(i=0;i<offset.rows();++i)
1303       for(j=0;j<offset.cols();++j)
1304 	offset(i,j) = osurf.ctrlPnts()(i,j) ;
1305 
1306     if(!baseLevel_){
1307       NurbsSurface<T,N>::refineKnotU(Xu) ;
1308     }
1309   }
1310 }
1311 
1312 /*!
1313   \brief  Refines the V knot vector
1314 
1315   \param X  the knot vector to refine from
1316 
1317   \author Philippe Lavoie
1318   \date 24 January, 1997
1319 */
1320 template <class T, int N>
refineKnotV(const Vector<T> & X)1321 void HNurbsSurface<T,N>::refineKnotV(const Vector<T>& X) {
1322   updateSurface() ;
1323   Vector<T> Xv(X) ;
1324   int i,j ;
1325   j = 0 ;
1326   for(i=0;i<X.n();++i){
1327     if(X[i]>=this->V[0] && X[i]<=this->V[this->V.n()-1]){
1328       Xv[j] = X[i] ;
1329       ++j ;
1330     }
1331   }
1332   Xv.resize(j) ;
1333 
1334   if(Xv.n()>0){
1335     if(nextLevel_){
1336       nextLevel_->refineKnotV(Xv) ;
1337     }
1338 
1339     NurbsSurface<T,N> osurf(this->degU,this->degV,this->U,this->V,offset) ;
1340 
1341     osurf.refineKnotV(Xv) ;
1342 
1343     offset.resize(osurf.ctrlPnts().rows(),osurf.ctrlPnts().cols()) ;
1344     for(i=0;i<offset.rows();++i)
1345       for(j=0;j<offset.cols();++j)
1346 	offset(i,j) = osurf.ctrlPnts()(i,j) ;
1347     if(!baseLevel_){
1348       NurbsSurface<T,N>::refineKnotV(Xv) ;
1349     }
1350   }
1351 }
1352 
1353 /*!
1354   \brief  Move a point on the surface
1355 
1356   This moves the point \a s(u,v) by delta. As this is a HNURBS
1357   surface. It moves the offset surface by delta, it doesn't move
1358   the surface point per say.
1359 
1360   \param u  the parameter in the u direction
1361   \param v  the parameter in the v direction
1362   \param delta  the displacement of the point \a s(u,v)
1363 
1364   \return 1 if the operation was succesfull, 0 otherwise
1365 
1366   \warning u and v must be in a valid range.
1367 
1368   \author Philippe Lavoie
1369   \date 3 June 1998
1370 */
1371 template <class T, int N>
movePointOffset(T u,T v,const Point_nD<T,N> & delta)1372 int HNurbsSurface<T,N>::movePointOffset(T u, T v, const Point_nD<T,N>& delta){
1373   this->P = offset ;
1374 
1375   // by definition the offset has w = 0 , but this isn't valid for
1376   // the control points, increasing the w by 1, will generate a valid surface
1377   if(baseLevel_)
1378     for(int i=0;i<this->P.rows();++i)
1379       for(int j=0;j<this->P.cols();++j){
1380 	this->P(i,j).w() += T(1) ;
1381       }
1382 
1383   if(NurbsSurface<T,N>::movePoint(u,v,delta)){
1384     offset = this->P ;
1385     // need to reset the offset weights
1386     if(baseLevel_)
1387       for(int i=0;i<this->P.rows();++i)
1388 	for(int j=0;j<this->P.cols();++j){
1389 	  this->P(i,j).w() -= T(1) ;
1390       }
1391 
1392     this->P = baseSurf.ctrlPnts() ;
1393     updateSurface() ;
1394     return 1 ;
1395   }
1396   updateSurface() ;
1397   return 0 ;
1398 }
1399 
1400 /*!
1401   \brief  Scales the object
1402 
1403   \param s  the scaling factor
1404 
1405   \author Philippe Lavoie
1406   \date 11 June 1998
1407 */
1408 template <class T, int N>
scale(const Point_nD<T,N> & s)1409 void HNurbsSurface<T,N>::scale(const Point_nD<T,N>& s){
1410   for(int i=0;i<offset.rows();++i)
1411     for(int j=0;j<offset.cols();++j){
1412       offset(i,j).x() *= s.x() ;
1413       offset(i,j).y() *= s.y() ;
1414       offset(i,j).z() *= s.z() ;
1415     }
1416 
1417   if(nextLevel_){
1418     nextLevel_->scale(s) ;
1419   }
1420 }
1421 
1422 
1423 /*!
1424   \brief  The offset vector are fixed
1425 
1426   Fixes the offset vector direction to a unique value. The offset
1427   vector's direction won't depend on its base layer.
1428 
1429   \param I  the i vector
1430   \param J  the j vector
1431   \param  K  the k vector
1432 
1433   \author Philippe Lavoie
1434   \date 11 June 1998
1435 */
1436 template <class T, int N>
setFixedOffsetVector(const Point_nD<T,N> & I,const Point_nD<T,N> & J,const Point_nD<T,N> & K)1437 void HNurbsSurface<T,N>::setFixedOffsetVector(const Point_nD<T,N> &I, const Point_nD<T,N> &J, const Point_nD<T,N>& K){
1438   fixedOffset = 1;
1439   initBase();
1440   ivec(0,0) = I;
1441   jvec(0,0) = J;
1442   kvec(0,0) = K;
1443   updateSurface();
1444 }
1445 
1446 /*!
1447   \brief  The offset vector are variable
1448 
1449   Fixes the offsset vector direction to a variable value. The value
1450   depends on its base layer.
1451 
1452 
1453   \author Philippe Lavoie
1454   \date 11 June 1998
1455 */
1456 template <class T, int N>
setVariableOffsetVector()1457 void HNurbsSurface<T,N>::setVariableOffsetVector(){
1458   fixedOffset = 0 ;
1459   initBase();
1460   updateSurface();
1461 }
1462 
1463 } //end namespace
1464