hnurbsS.cpp Source File

00001 /*=============================================================================
00002         File: hnurbsS.cpp
00003      Purpose:       
00004     Revision: $Id: hnurbsS.cpp,v 1.3 2002/05/17 18:24:21 philosophil Exp $
00005   Created by: Philippe Lavoie    (7 October 1997)
00006  Modified by: 
00007 
00008  Copyright notice:
00009           Copyright (C) 1996-1998 Philippe Lavoie
00010  
00011           This library is free software; you can redistribute it and/or
00012           modify it under the terms of the GNU Library General Public
00013           License as published by the Free Software Foundation; either
00014           version 2 of the License, or (at your option) any later version.
00015  
00016           This library is distributed in the hope that it will be useful,
00017           but WITHOUT ANY WARRANTY; without even the implied warranty of
00018           MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019           Library General Public License for more details.
00020  
00021           You should have received a copy of the GNU Library General Public
00022           License along with this library; if not, write to the Free
00023           Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00024 =============================================================================*/
00025 #include <hnurbsS.h>
00026 #include <string.h>
00027 
00030 namespace PLib {
00031 
00038 template <class T, int N>
00039 HNurbsSurface<T,N>::HNurbsSurface():NurbsSurface<T,N>(), fixedOffset(0){
00040 
00041   baseLevel_ = 0 ;
00042   nextLevel_ = 0 ;
00043   firstLevel_ = lastLevel_ = this ;
00044 
00045   rU.resize(0) ;
00046   rV.resize(0) ;
00047   updateN = 0 ;
00048   baseUpdateN = 0 ;
00049 
00050   level_ = 0 ;
00051 }
00052 
00065 template <class T, int N>
00066 HNurbsSurface<T,N>::HNurbsSurface(HNurbsSurface<T,N>* base, const Vector<T>& xU, const Vector<T>& xV):NurbsSurface<T,N>(),fixedOffset(0){
00067 
00068   if(!base){
00069     Error error("HNurbsSurface<T,N> constructor") ;
00070     error << "Initializing a HNurbsSurface<T,N> with a null base pointer!" ;
00071     error.fatal() ;
00072   }
00073   if(base->nextLevel_){
00074     Error error("HNurbsSurface<T,N> constructor") ;
00075     error << "You're trying to replace an existing level, this is not allowed." ;
00076     error.fatal() ;
00077   }
00078 
00079 
00080   nextLevel_ = 0 ;
00081   firstLevel_ = base->firstLevel_ ;
00082   lastLevel_ = this ;
00083   baseLevel_ = (HNurbsSurface<T,N>*)base ;
00084 
00085   // update the information in the previous levels
00086   baseLevel_->nextLevel_ = this ;
00087   HNurbsSurface<T,N>* l ;
00088   l = baseLevel_ ;
00089   while(l){
00090     l->lastLevel_ = this ;
00091     l = l->baseLevel_ ;
00092   }
00093 
00094 
00095   level_ = base->level_ + 1 ;
00096 
00097   rU = xU ;
00098   rV = xV ;
00099 
00100   updateN = 0 ;
00101   baseUpdateN = baseLevel_->modifiedN()-1 ; // Set it so that initBaseLevel will run
00102 
00103   initBase() ;
00104   offset.resize(baseSurf.ctrlPnts()) ;
00105 
00106   P = baseSurf.ctrlPnts() ;
00107   U = baseSurf.knotU() ;
00108   V = baseSurf.knotV() ;
00109   degU = baseSurf.degreeU() ;
00110   degV = baseSurf.degreeV() ;
00111 
00112   //updateSurface() ;
00113 
00114 }
00115 
00128 template <class T, int N>
00129 HNurbsSurface<T,N>::HNurbsSurface(HNurbsSurface<T,N>* base):NurbsSurface<T,N>(),fixedOffset(0){
00130 
00131   if(!base){
00132     Error error("HNurbsSurface<T,N> constructor") ;
00133     error << "Initializing a HNurbsSurface<T,N> with a null base pointer!" ;
00134     error.fatal() ;
00135   }
00136   if(base->nextLevel_){
00137     Error error("HNurbsSurface<T,N> constructor") ;
00138     error << "You're trying to replace an existing level, this is not allowed." ;
00139     error.fatal() ;
00140   }
00141 
00142   nextLevel_ = 0 ;
00143   firstLevel_ = base->firstLevel_ ;
00144   lastLevel_ = this ;
00145   baseLevel_ = (HNurbsSurface<T,N>*)base ;
00146 
00147   // update the information in the previous levels
00148   baseLevel_->nextLevel_ = this ;
00149   HNurbsSurface<T,N>* l ;
00150   l = baseLevel_ ;
00151   while(l){
00152     l->lastLevel_ = this ;
00153     l = l->baseLevel_ ;
00154   }
00155 
00156   level_ = base->level_ + 1 ;
00157 
00158   updateN = 0 ;
00159   rU.resize(0) ;
00160   rV.resize(0) ;
00161 
00162   baseUpdateN = baseLevel_->modifiedN()-1 ; // Set it so that initBase will run
00163   initBase() ;
00164   offset.resize(baseSurf.ctrlPnts()) ;
00165   P = baseSurf.ctrlPnts() ;
00166   U = baseSurf.knotU() ;
00167   V = baseSurf.knotV() ;
00168   degU = baseSurf.degreeU() ;
00169   degV = baseSurf.degreeV() ;
00170   //updateSurface() ;
00171 
00172 }
00173 
00189 template <class T, int N>
00190 HNurbsSurface<T,N>::HNurbsSurface(const NurbsSurface<T,N>& S):NurbsSurface<T,N>(S),fixedOffset(0){
00191 
00192   baseLevel_ = 0 ;
00193   nextLevel_ = 0 ; 
00194   firstLevel_ = lastLevel_ = this ;
00195 
00196   level_ = 0 ;
00197   updateN = 0 ;
00198   baseUpdateN = 0 ;
00199 
00200   rU.resize(0) ;
00201   rV.resize(0) ;
00202 
00203   offset = P ;
00204 }
00205 
00214 template <class T, int N>
00215 HNurbsSurface<T,N>::HNurbsSurface(const HNurbsSurface<T,N>& S):NurbsSurface<T,N>(S),fixedOffset(0) {
00216 
00217   baseLevel_ = nextLevel_ = 0 ;
00218   firstLevel_ = lastLevel_ = this ;
00219 
00220   level_ = 0 ;
00221   updateN = 0 ;
00222   baseUpdateN = 0 ;
00223 
00224   rU = S.rU ; 
00225   rV = S.rV ; 
00226 
00227   offset = S.offset ;
00228 
00229   this->copy(S) ;
00230 
00231 }
00232 
00244 template <class T, int N>
00245 HNurbsSurface<T,N>::HNurbsSurface(HNurbsSurface<T,N> *base, const HNurbsSurface<T,N> &surf):NurbsSurface<T,N>(surf),fixedOffset(0) {
00246 
00247   if(!base){
00248     Error error("HNurbsSurface<T,N> constructor") ;
00249     error << "Initializing a HNurbsSurface<T,N> with a null base pointer!" ;
00250     error.fatal() ;
00251   }
00252   if(base->nextLevel_){
00253     Error error("HNurbsSurface<T,N> constructor") ;
00254     error << "You're trying to replace an existing level, this is not allowed." ;
00255     error.fatal() ;
00256   }
00257 
00258   baseLevel_ = (HNurbsSurface<T,N>*)base ;
00259   nextLevel_ = 0 ;
00260   lastLevel_ = this ;
00261 
00262   // update the information in the previous levels
00263   baseLevel_->nextLevel_ = this ;
00264   HNurbsSurface<T,N>* l ;
00265   l = baseLevel_ ;
00266   while(l){
00267     l->lastLevel_ = this ;
00268     l = l->baseLevel_ ;
00269   }
00270 
00271   firstLevel_ = base->firstLevel_ ;
00272   level_ = base->level_+1 ;
00273   baseUpdateN = base->modifiedN()-1 ;
00274 
00275   initBase() ; 
00276 
00277   updateN = 0 ;
00278 
00279   this->copy(surf) ;
00280 
00281 }
00282 
00291 template <class T, int N>
00292 void HNurbsSurface<T,N>::copy(const HNurbsSurface<T,N>& nS){
00293   HNurbsSurface<T,N> *levelP ;
00294   levelP = nS.nextLevel() ;
00295 
00296   NurbsSurface<T,N>::operator=(nS) ;
00297   rU = nS.rU ;
00298   rV = nS.rV ;
00299   offset = nS.offset ;
00300   fixedOffset= nS.fixedOffset ;
00301 
00302   firstLevel_ = this ;
00303 
00304   if(levelP){
00305     HNurbsSurface<T,N> *newLevel = new HNurbsSurface(this,*levelP) ;
00306     nextLevel_ = newLevel ;
00307     lastLevel_ = nextLevel_->lastLevel() ;
00308   }
00309   else{
00310     lastLevel_ = this ;
00311   }
00312 
00313 }
00314 
00329 template <class T, int N>
00330 void HNurbsSurface<T,N>::updateSurface(int i0, int j0){
00331   if(i0>=0 && j0>=0){
00332     if(offset(i0,j0).x()==0.0 && offset(i0,j0).y()==0.0 && offset(i0,j0).z()==0.0)
00333       return ;
00334   }
00335   if(baseLevel_){
00336     if(initBase()){
00337       P = baseSurf.ctrlPnts() ;
00338       U = baseSurf.knotU() ;
00339       V = baseSurf.knotV() ;
00340       degU = baseSurf.degreeU() ;
00341       degV = baseSurf.degreeV() ;
00342     }
00343     if(i0>=0 && j0>=0){
00344       Point_nD<T,N> vecOffset ;
00345       if(fixedOffset){
00346         vecOffset = offset(i0,j0).x()*ivec(0,0) +
00347           offset(i0,j0).y()*jvec(0,0) +
00348           offset(i0,j0).z()*kvec(0,0) ;
00349       }
00350       else{
00351         vecOffset = offset(i0,j0).x()*ivec(i0,j0) +
00352           offset(i0,j0).y()*jvec(i0,j0) +
00353           offset(i0,j0).z()*kvec(i0,j0) ;
00354       }
00355       P(i0,j0).x() = baseSurf.ctrlPnts()(i0,j0).x()+vecOffset.x() ;
00356       P(i0,j0).y() = baseSurf.ctrlPnts()(i0,j0).y()+vecOffset.y() ;
00357       P(i0,j0).z() = baseSurf.ctrlPnts()(i0,j0).z()+vecOffset.z() ;
00358     }
00359     else{
00360       for(int i=0;i<P.rows();++i)
00361         for(int j=0;j<P.cols();++j){
00362           if(offset(i,j).x() != 0 || 
00363              offset(i,j).y() != 0 || offset(i,j).z() || 0){
00364             Point_nD<T,N> vecOffset ;
00365             if(fixedOffset){
00366               vecOffset = offset(i,j).x()*ivec(0,0) +
00367                 offset(i,j).y()*jvec(0,0) +
00368                 offset(i,j).z()*kvec(0,0) ;
00369             }
00370             else{
00371               vecOffset = offset(i,j).x()*ivec(i,j) +
00372                 offset(i,j).y()*jvec(i,j) +
00373                 offset(i,j).z()*kvec(i,j) ;
00374             }
00375             P(i,j).x() = baseSurf.ctrlPnts()(i,j).x()+vecOffset.x() ;
00376             P(i,j).y() = baseSurf.ctrlPnts()(i,j).y()+vecOffset.y() ;
00377             P(i,j).z() = baseSurf.ctrlPnts()(i,j).z()+vecOffset.z() ;
00378           }
00379         }
00380     }
00381   }
00382   else{
00383     if(i0>=0 && j0>=0)
00384       P(i0,j0) = offset(i0,j0) ;
00385     else{
00386       for(int i=0;i<P.rows();++i)
00387         for(int j=0;j<P.cols();++j){
00388           P(i,j) = offset(i,j) ;
00389         }
00390     }
00391   }
00392 
00393   ++updateN ;
00394 }
00395 
00396 
00410 template <class T, int N>
00411 int HNurbsSurface<T,N>::initBase(int force){
00412   if(!baseLevel_){
00413     return 0; // this means that the surface is the base level
00414   }
00415 
00416   // make sure none of the base level need updating
00417   if(baseLevel_->initBase())
00418     force = 1 ;
00419   
00420 
00421   // Only initialize the base level if it isn't up to date
00422   if(baseLevel_->modifiedN() == baseUpdateN && !force){
00423     return 0;
00424   }
00425 
00426   baseUpdateN = baseLevel_->modifiedN() ;
00427 
00428   baseSurf = *baseLevel_ ;
00429 
00430   if(rU.n()>0)
00431     baseSurf.refineKnotU(rU) ;
00432   if(rV.n()>0)
00433     baseSurf.refineKnotV(rV) ;
00434 
00435   Vector<T> maxU,maxV ;
00436   int i,j ;
00437   
00438   if(baseSurf.degreeU()>3)
00439     averagingKnots(baseSurf.knotU(),baseSurf.degreeU(),maxU) ;
00440   else{
00441     maxU.resize(baseSurf.ctrlPnts().rows()) ;
00442     for(i=0;i<baseSurf.ctrlPnts().rows();++i)
00443       if(!maxInfluence(i,baseSurf.knotU(),baseSurf.degreeU(),maxU[i]))
00444         cerr << "Problem in maxInfluence U!\n" ;
00445   }
00446 
00447   if(baseSurf.degreeV()>3)
00448     averagingKnots(baseSurf.knotV(),baseSurf.degreeV(),maxV) ;
00449   else{
00450     maxV.resize(baseSurf.ctrlPnts().cols()) ;
00451     for(i=0;i<baseSurf.ctrlPnts().cols();++i)
00452       if(!maxInfluence(i,baseSurf.knotV(),baseSurf.degreeV(),maxV[i]))
00453         cerr << "Problem in maxInfluence V!\n" ;
00454   }
00455 
00456   if(fixedOffset){
00457     if(ivec.rows() != 1 || ivec.cols()!=1){
00458       ivec.resize(1,1);
00459       jvec.resize(1,1);
00460       kvec.resize(1,1);
00461     }
00462   }
00463   else{
00464     ivec.resize(maxU.n(),maxV.n()) ;
00465     jvec.resize(maxU.n(),maxV.n()) ;
00466     kvec.resize(maxU.n(),maxV.n()) ;
00467 
00468     Matrix< Point_nD<T,N> > ders ;
00469     
00470     for(i=0;i<maxU.n();++i)
00471       for(j=0;j<maxV.n();++j){
00472         baseSurf.deriveAt(maxU[i],maxV[j],1,ders) ;
00473         // It is possible that there are no valid derivative at that point
00474         // we then move to a nearby point to get a valid value
00475         Point_nD<T,N> norm = crossProduct(ders(0,1),ders(1,0)) ;
00476         if(norm.x() == T(0) &&
00477            norm.y() == T(0) &&
00478            norm.z() == T(0)){
00479           const T delta = 0.00001 ;
00480           T nt = 1.0 ; 
00481           Matrix< Point_nD<T,N> > dersT(ders) ;
00482           while(norm.x() == T(0) && 
00483                 norm.y() == T(0) &&
00484                 norm.z() == T(0)){
00485             if( nt*delta >(baseSurf.knotU()[baseSurf.knotU().n()-1]-baseSurf.knotU()[0])){
00486               Error error("initBase");
00487               error << "Can't compute the derivative.\n" ;
00488               error.fatal() ;
00489             }
00490             T nd ;
00491             
00492             ders.reset(0) ;
00493             
00494             nd = 0 ;
00495             if(i != 0){
00496               baseSurf.deriveAt(maxU[i] - nt*delta, maxV[j],1, dersT) ;
00497               ders += dersT ;
00498               ++nd ;
00499             }
00500             if(i != maxU.n()-1){
00501               baseSurf.deriveAt(maxU[i] + nt*delta, maxV[j],1, dersT) ;
00502               ders += dersT ;
00503               ++nd ;
00504             }
00505             if(j != 0){
00506               baseSurf.deriveAt(maxU[i], maxV[j] - nt*delta,1, dersT) ;
00507               ders += dersT ;
00508               ++nd ;
00509             }
00510             if(j != maxV.n()-1){
00511               baseSurf.deriveAt(maxU[i], maxV[j] + nt*delta,1, dersT) ;
00512               ders += dersT ;
00513               ++nd ;
00514             }
00515             
00516             if(nd==0){
00517               Error error("initBase");
00518               error << "Can't compute the derivative.\n" ;
00519               error.fatal() ;
00520             }
00521             
00522             ders /= nd ;
00523             
00524             norm = crossProduct(ders(0,1),ders(1,0)) ;
00525             nt *= 10.0 ; 
00526           }
00527         }       
00528         
00529         ivec(i,j) = ders(0,1).unitLength() ;
00530         kvec(i,j) = crossProduct(ders(0,1),ders(1,0)).unitLength() ;
00531         jvec(i,j) = crossProduct(kvec(i,j),ivec(i,j)).unitLength() ;
00532       }
00533   }
00534     
00535   return 1 ;
00536 }
00537 
00549 template <class T, int N>
00550 int HNurbsSurface<T,N>::modifies(T u, T v){
00551   if(nextLevel_){
00552     int mod = nextLevel_->modifies(u,v) ;
00553     if(mod>=0)
00554       return mod ;
00555   }
00556 
00557   if(u<knotU()[0] || u>knotU()[knotU().n()-1])
00558     return -1 ;
00559   if(v<knotV()[0] || v>knotU()[knotV().n()-1])
00560     return -1 ;
00561 
00562   int su = findSpanU(u) ;
00563   int sv = findSpanV(v) ;
00564 
00565   for(int i=0;i<=degU;++i)
00566     for(int j=0;j<=degV;++j){
00567       if(offset(su-degU+i,sv+degV+j) != HPoint_nD<T,N>(0,0,0,0))
00568         return level_ ;
00569     }
00570 
00571   return -1 ;
00572 }
00573 
00587 template <class T, int N>
00588 HPoint_nD<T,N> HNurbsSurface<T,N>::hpointAt(T u, T v, int lod) const {
00589   if(level_==lod)
00590     return NurbsSurface<T,N>::operator()(u,v) ;
00591 
00592   if(nextLevel_)
00593     return nextLevel_->hpointAt(u,v,lod) ; 
00594   
00595   return NurbsSurface<T,N>::operator()(u,v) ;
00596 }
00597 
00609 template <class T, int N>
00610 int HNurbsSurface<T,N>::maxLevel() const{
00611   if(lastLevel_)
00612     return lastLevel_->level_ ;
00613   return level_ ; 
00614 }
00615 
00631 template <class T, int N>
00632 void HNurbsSurface<T,N>::deriveAtH(T u, T v, int d, Matrix< HPoint_nD<T,N> >& ders,int lod) const {
00633   if(level_==lod){
00634     NurbsSurface<T,N>::deriveAtH(u,v,d,ders) ;
00635     return ;
00636   }
00637 
00638   if(nextLevel_){
00639     nextLevel_->deriveAtH(u,v,d,ders,lod) ; 
00640     return ;
00641   }
00642   
00643   NurbsSurface<T,N>::deriveAtH(u,v,d,ders) ;
00644 }
00645 
00661 template <class T, int N>
00662 void HNurbsSurface<T,N>::deriveAt(T u, T v, int d, Matrix< Point_nD<T,N> >& ders, int lod) const {
00663   if(level_==lod){
00664     NurbsSurface<T,N>::deriveAt(u,v,d,ders) ;
00665     return ;
00666   }
00667 
00668   if(nextLevel_){
00669     nextLevel_->deriveAt(u,v,d,ders,lod) ; 
00670     return ;
00671 
00672   NurbsSurface<T,N>::deriveAt(u,v,d,ders) ;
00673   }
00674 }
00675 
00684 template <class T, int N>
00685 HNurbsSurface<T,N>::~HNurbsSurface(){
00686   if(nextLevel_)
00687     delete nextLevel_ ;
00688   lastLevel_ = 0 ;
00689   if(baseLevel_){
00690     baseLevel_->nextLevel_ = 0 ;
00691     baseLevel_->lastLevel_ = baseLevel_ ;
00692   }
00693   baseLevel_ = 0 ;
00694   nextLevel_ = 0 ;
00695   firstLevel_ = 0 ;
00696 }
00697 
00706 template <class T, int N>
00707 void HNurbsSurface<T,N>::updateLevels(int upLevel){
00708   if(upLevel>=0){
00709     if(level_ <= upLevel){
00710       this->updateSurface() ;
00711     }
00712   }
00713   else{
00714     this->updateSurface() ;
00715   }
00716 
00717   if(upLevel>level() || upLevel<0){
00718     if(nextLevel_)
00719       nextLevel_->updateLevels(upLevel) ;
00720   }
00721 }
00722 
00742 template <class T, int N>
00743 void HNurbsSurface<T,N>::splitUV(int nu, int nv, Vector<T> &nU, Vector<T> &nV){
00744 
00745   nU.resize(knotU().n()*nu) ;
00746   nV.resize(knotV().n()*nv) ;
00747   
00748   int i,j,n ;
00749 
00750   n = 0 ; 
00751   for(i=1;i<knotU().n();++i){
00752     if(knotU()[i] >knotU()[i-1]){
00753       T a = knotU()[i-1] ;
00754       T b = knotU()[i] ;
00755 
00756 
00757       for(j=0;j<nu;++j){
00758         nU[n] = a + (b-a)*T(j+1)/T(nu+1) ;
00759         ++n ;
00760       }
00761     }
00762   }  
00763   nU.resize(n) ;
00764 
00765   n = 0 ;
00766   for(i=1;i<knotV().n();++i){
00767     if(knotV()[i] > knotV()[i-1]){
00768       T a = knotV()[i-1] ;
00769       T b = knotV()[i] ;
00770 
00771       for(j=0;j<nv;++j){
00772         nV[n] = a + (b-a)*T(j+1)/T(nv+1) ;
00773         ++n ;
00774       }
00775     }
00776   }  
00777   nV.resize(n) ;
00778 
00779 }
00780  
00802 template <class T, int N>
00803 void HNurbsSurface<T,N>::splitUV(int nu, int su, int nv, int sv, Vector<T> &nU, Vector<T> &nV){
00804 
00805   int i,j,n ;
00806 
00807   if(su<=0)
00808     su = degU  ;
00809   if(sv<=0)
00810     sv = degV  ;
00811   if(su>degU+1)
00812     su = degU+1 ;
00813   if(sv>degV+1)
00814     sv = degV+1 ;
00815 
00816   nU.resize(knotU().n()*nu*su) ;
00817   nV.resize(knotV().n()*nv*sv) ;
00818   
00819   n = 0 ; 
00820   for(i=1;i<knotU().n();++i){
00821     if(knotU()[i] >knotU()[i-1]){
00822       T a = knotU()[i-1] ;
00823       T b = knotU()[i] ;
00824 
00825 
00826       for(j=0;j<nu;++j){
00827         T u = a + (b-a)*T(j+1)/T(nu+1) ;
00828         for(int l=0;l<su;++l){
00829           nU[n] = u ;
00830           ++n ;
00831         }
00832       }
00833     }
00834   }  
00835   nU.resize(n) ;
00836 
00837   n = 0 ;
00838   for(i=1;i<knotV().n();++i){
00839     if(knotV()[i] > knotV()[i-1]){
00840       T a = knotV()[i-1] ;
00841       T b = knotV()[i] ;
00842 
00843       for(j=0;j<nv;++j){
00844         T v = a + (b-a)*T(j+1)/T(nv+1) ;
00845         for(int l=0;l<sv;++l){
00846           nV[n] = v ;
00847           ++n ;
00848         }
00849       }
00850     }
00851   }  
00852   nV.resize(n) ;
00853 
00854 }
00855  
00868 template <class T, int N>
00869 HNurbsSurface<T,N>* HNurbsSurface<T,N>::addLevel(int n) {
00870   HNurbsSurface<T,N> *newLevel ;
00871 
00872   if(nextLevel_)
00873     return 0 ;
00874 
00875   Vector<T> newU,newV ;
00876   
00877   splitUV(n,n,newU,newV) ;
00878 
00879   newLevel = new HNurbsSurface(this,newU,newV) ;
00880 
00881   return newLevel ;
00882 }
00883 
00896 template <class T, int N>
00897 HNurbsSurface<T,N>* HNurbsSurface<T,N>::addLevel() {
00898   HNurbsSurface<T,N> *newLevel ;
00899 
00900   if(nextLevel_)
00901     return 0 ;
00902 
00903   newLevel = new HNurbsSurface(this) ;
00904 
00905   return newLevel ;
00906 }
00907 
00925 template <class T, int N>
00926 int HNurbsSurface<T,N>::isoCurveU(T u, NurbsCurve<T,N>& c,int lod) const {
00927   if(lod>=0 && level_>lod)
00928     return 0;
00929 
00930   if(lod==level_ || lod<0){
00931     NurbsSurface<T,N>::isoCurveU(u,c) ; 
00932     return 1 ;
00933   }
00934   
00935   if(nextLevel_){
00936     return nextLevel_->isoCurveU(u,c,lod) ;
00937   }
00938 
00939   return 0 ;
00940 }
00941 
00958 template <class T, int N>
00959 int HNurbsSurface<T,N>::isoCurveV(T v, NurbsCurve<T,N>& c,int lod) const {
00960   if(lod>=0 && level_>lod)
00961     return 0;
00962 
00963   if(lod==level_ || lod<0){
00964     NurbsSurface<T,N>::isoCurveV(v,c) ; 
00965     return 1 ;
00966   }
00967   
00968   if(nextLevel_){
00969     return nextLevel_->isoCurveV(v,c,lod) ;
00970   }
00971 
00972   return 0 ;
00973 }
00974 
00975 
00976 
00987 template <class T, int N>
00988 int HNurbsSurface<T,N>::read(ifstream &fin){
00989   if(!fin) {
00990     return 0 ;
00991   }
00992   int nu,nv,du,dv;
00993   char *type ;
00994   type = new char[4] ;
00995   if(!fin.read(type,sizeof(char)*4)) { delete []type ; return 0 ;}
00996   int r1 = strncmp(type,"hns3",4) ;
00997   int r2 = strncmp(type,"hns4",4) ;
00998   int r3 = strncmp(type,"hnso",4) ;
00999   if(!(r1 || r2 || r3)) 
01000     return 0 ;
01001   T *p,*p2 ;
01002   char stc ; 
01003   if(!r1 || !r2){
01004     if(!fin.read((char*)&stc,sizeof(char))) { delete []type ; return 0 ;}
01005 
01006     int st = stc - '0' ; 
01007     if(st != sizeof(T)){ // does not have the same data size
01008       delete []type ;
01009       return 0 ; 
01010     }
01011 
01012     if(!fin.read((char*)&nu,sizeof(int))) { delete []type ; return 0 ;}
01013     if(!fin.read((char*)&nv,sizeof(int))) { delete []type ; return 0 ;}
01014     if(!fin.read((char*)&du,sizeof(int))) { delete []type ; return 0 ;}
01015     if(!fin.read((char*)&dv,sizeof(int))) { delete []type ; return 0 ;}
01016     
01017     resize(nu,nv,du,dv) ;
01018     
01019     if(!fin.read((char*)U.memory(),sizeof(T)*U.n())) { delete []type ; return 0 ;}
01020     if(!fin.read((char*)V.memory(),sizeof(T)*V.n())) { delete []type ; return 0 ;}
01021     
01022     if(!r1){
01023       p = new T[3*nu*nv] ;
01024       if(!fin.read((char*)p,sizeof(T)*3*nu*nv)) { delete []type ; return 0 ;}
01025       p2 = p ;
01026       for(int i=0;i<nu;i++)
01027         for(int j=0;j<nv;j++){
01028           P(i,j).x() = *(p++) ;
01029           P(i,j).y() = *(p++) ;
01030           P(i,j).z() = *(p++) ;
01031           P(i,j).w() = 1.0 ;
01032         }
01033       delete []p2 ;
01034     }
01035     else{
01036       p = new T[4*nu*nv] ;
01037       if(!fin.read((char*)p,sizeof(T)*4*nu*nv)) { delete []type ; return 0 ;}
01038       p2 = p ;
01039       for(int i=0;i<nu;i++)
01040         for(int j=0;j<nv;j++){
01041           P(i,j).x() = *(p++) ;
01042           P(i,j).y() = *(p++) ;
01043           P(i,j).z() = *(p++) ;
01044           P(i,j).w() = *(p++) ;
01045         }
01046       delete []p2 ;
01047     }
01048     offset = P ;
01049     this->updateSurface() ;
01050   }
01051   else { // reading the offset information
01052     int ru,rv ;
01053     if(!fin.read((char*)&ru,sizeof(int))) { delete []type ; return 0 ;}
01054     if(!fin.read((char*)&rv,sizeof(int))) { delete []type ; return 0 ;}
01055     rU.resize(ru) ;
01056     rV.resize(rv) ;
01057     if(rU.n()>0)
01058       if(!fin.read((char*)rU.memory(),sizeof(T)*rU.n())) { delete []type ; return 0 ;}
01059     if(rV.n()>0)
01060       if(!fin.read((char*)rV.memory(),sizeof(T)*rV.n())) { delete []type ; return 0 ;}
01061     
01062     if(!fin.read((char*)&nu,sizeof(int))) { delete []type ; return 0 ;}
01063     if(!fin.read((char*)&nv,sizeof(int))) { delete []type ; return 0 ;}
01064 
01065     p = new T[4*nu*nv] ;
01066     if(!fin.read((char*)p,sizeof(T)*4*nu*nv)) { delete []type ; return 0 ;}
01067     p2 = p ;
01068     offset.resize(nu,nv) ;
01069     for(int i=0;i<nu;i++)
01070       for(int j=0;j<nv;j++){
01071         offset(i,j).x() = *(p++) ;
01072         offset(i,j).y() = *(p++) ;
01073         offset(i,j).z() = *(p++) ;
01074         offset(i,j).w() = *(p++) ;
01075       }
01076     delete []p2 ;    
01077     --baseUpdateN ;
01078     this->updateSurface() ;
01079   }
01080 
01081   // now we must read all the level information (if any)
01082   char *ptxt ;
01083   ptxt = new char[7] ;
01084   // the following line is used so that purify doesn't give a warning
01085   ptxt[0] = ptxt[1] = ptxt[2] = ptxt[3] = ptxt[4] = ptxt[5] = ptxt[6] = 0 ;
01086   int mark  = fin.tellg() ;
01087   if(!fin.read(ptxt,sizeof(char)*5))
01088     { delete []ptxt ; delete []type ; return 1 ; }
01089   char *pat ;
01090   pat = strstr(ptxt,"level") ;
01091   if(pat){
01092     // insert a new level
01093     HNurbsSurface<T,N> *newPatch = this->addLevel() ;
01094     if(newPatch){
01095       if(!newPatch->read(fin)) 
01096         return 0 ;      
01097     }
01098     else
01099       return 0 ;
01100   }
01101   else{ 
01102     // something else is after the HNurbs
01103     // restore were the file was before the read
01104     fin.seekg(mark);
01105   }
01106 
01107   delete []ptxt ;
01108   delete []type ;
01109   return 1 ;
01110 }
01111 
01122 template <class T, int N>
01123 int HNurbsSurface<T,N>::read(const char* filename){
01124   ifstream fin(filename) ;
01125   if(!fin) {
01126     return 0 ;
01127   }
01128 
01129   return this->read(fin) ;
01130 }
01131 
01142 template <class T, int N>
01143 int HNurbsSurface<T,N>::write(ofstream &fout) const {
01144   if(!fout)
01145     return 0 ;
01146   if(!baseLevel_){
01147     int prows = P.rows();
01148     int pcols = P.cols();
01149     char st = '0' + sizeof(T) ; 
01150     if(!fout.write((char*)&"hns4",sizeof(char)*4)) return 0 ;
01151     if(!fout.write((char*)&st,sizeof(char))) return 0 ; 
01152     if(!fout.write((char*)&prows,sizeof(int))) return 0 ;
01153     if(!fout.write((char*)&pcols,sizeof(int))) return 0 ;
01154     if(!fout.write((char*)&degU,sizeof(int))) return 0 ;
01155     if(!fout.write((char*)&degV,sizeof(int))) return 0 ;
01156     if(!fout.write((char*)U.memory(),sizeof(T)*U.n())) return 0 ;
01157     if(!fout.write((char*)V.memory(),sizeof(T)*V.n())) return 0 ;
01158     
01159     T *p,*p2 ;
01160     p = new T[P.rows()*P.cols()*4] ;
01161     p2 = p ;
01162     for(int i=0;i<P.rows();i++) 
01163       for(int j=0;j<P.cols();j++){
01164         *p = offset(i,j).x() ; p++ ;
01165         *p = offset(i,j).y() ; p++ ;
01166         *p = offset(i,j).z() ; p++ ;
01167         *p = offset(i,j).w() ; p++ ;
01168       }
01169     if(!fout.write((char*)p2,sizeof(T)*P.rows()*P.cols()*4)) return 0 ;
01170     delete []p2 ;
01171   }
01172   else{
01173     if(!fout.write((char*)&"hnso",sizeof(char)*4)) return 0 ;
01174     int run = rU.n();
01175     int rvn = rV.n();
01176     if(!fout.write((char*)&run,sizeof(int))) return 0 ;
01177     if(!fout.write((char*)&rvn,sizeof(int))) return 0 ;
01178     if(rU.n()>0)
01179       if(!fout.write((char*)rU.memory(),sizeof(T)*rU.n())) return 0 ;
01180     if(rV.n()>0)
01181       if(!fout.write((char*)rV.memory(),sizeof(T)*rV.n())) return 0 ;
01182     int orows = offset.rows();
01183     int ocols = offset.cols() ;
01184     if(!fout.write((char*)&orows,sizeof(int))) return 0 ;
01185     if(!fout.write((char*)&ocols,sizeof(int))) return 0 ;
01186     T *p,*p2 ;
01187 
01188     p = new T[offset.rows()*offset.cols()*4] ;
01189     p2 = p ;
01190     for(int i=0;i<offset.rows();i++) 
01191       for(int j=0;j<offset.cols();j++){
01192         *p = offset(i,j).x() ; p++ ;
01193         *p = offset(i,j).y() ; p++ ;
01194         *p = offset(i,j).z() ; p++ ;
01195         *p = offset(i,j).w() ; p++ ;
01196       }
01197     if(!fout.write((char*)p2,sizeof(T)*offset.rows()*offset.cols()*4)) 
01198       return 0 ;
01199     delete []p2 ;    
01200   }
01201   if(nextLevel_){
01202     if(!fout.write((char*)&"level",sizeof(char)*5)) return 0 ;
01203     if(!nextLevel_->write(fout)) return 0 ;
01204   }
01205   return 1 ;
01206 }
01207 
01218 template <class T, int N>
01219 int HNurbsSurface<T,N>::write(const char* filename) const {
01220   ofstream fout(filename) ;    
01221   if(!fout)
01222     return 0 ;
01223   return write(fout);
01224 }
01225 
01235 //template <class T, int N>
01236 //void HNurbsSurface<T,N>::toNurbs(NurbsSurface<T,N>& surf) const {
01237   
01238 //}
01239 
01240 template <class T>
01241 inline void mergeInto(const Vector<T> &a, Vector<T> &b) {
01242   // merge a in b ;
01243   int n = b.n() ;
01244   b.resize(b.n()+a.n()) ;  
01245   for(int i=0;i<a.n();++i)
01246     b[n+i] = a[i] ;
01247   b.qSort() ;
01248 }
01249 
01259 template <class T, int N>
01260 void HNurbsSurface<T,N>::refineKnots(const Vector<T>& nU, const Vector<T>& nV){
01261   refineKnotU(nU) ;
01262   mergeInto(nU,rU) ;
01263 
01264   initBase(1) ;
01265 
01266   refineKnotV(nV) ;  
01267   mergeInto(nV,rV) ;
01268 }
01269 
01278 template <class T, int N>
01279 void HNurbsSurface<T,N>::refineKnotU(const Vector<T>& X) {
01280   updateSurface() ;
01281   Vector<T> Xu(X) ;
01282   int i,j ;
01283   j = 0 ;
01284   for(i=0;i<X.n();++i){
01285     if(X[i]>=U[0] && X[i]<=U[U.n()-1]){
01286       Xu[j] = X[i] ;
01287       ++j ;
01288     }
01289   }
01290   Xu.resize(j) ;
01291 
01292   if(Xu.n()>0){
01293     if(nextLevel_){
01294       nextLevel_->refineKnotU(Xu) ;
01295     }
01296     
01297     NurbsSurface<T,N> osurf(degU,degV,U,V,offset) ;
01298     
01299     osurf.refineKnotU(Xu) ;
01300     
01301     offset.resize(osurf.ctrlPnts().rows(),osurf.ctrlPnts().cols()) ;
01302     for(i=0;i<offset.rows();++i)
01303       for(j=0;j<offset.cols();++j)
01304         offset(i,j) = osurf.ctrlPnts()(i,j) ;
01305 
01306     if(!baseLevel_){
01307       NurbsSurface<T,N>::refineKnotU(Xu) ;
01308     }
01309   }
01310 }
01311 
01320 template <class T, int N>
01321 void HNurbsSurface<T,N>::refineKnotV(const Vector<T>& X) {
01322   updateSurface() ;
01323   Vector<T> Xv(X) ;
01324   int i,j ;
01325   j = 0 ;
01326   for(i=0;i<X.n();++i){
01327     if(X[i]>=V[0] && X[i]<=V[V.n()-1]){
01328       Xv[j] = X[i] ;
01329       ++j ;
01330     }
01331   }
01332   Xv.resize(j) ;
01333 
01334   if(Xv.n()>0){
01335     if(nextLevel_){
01336       nextLevel_->refineKnotV(Xv) ;
01337     }
01338     
01339     NurbsSurface<T,N> osurf(degU,degV,U,V,offset) ;
01340     
01341     osurf.refineKnotV(Xv) ;
01342     
01343     offset.resize(osurf.ctrlPnts().rows(),osurf.ctrlPnts().cols()) ;
01344     for(i=0;i<offset.rows();++i)
01345       for(j=0;j<offset.cols();++j)
01346         offset(i,j) = osurf.ctrlPnts()(i,j) ;
01347     if(!baseLevel_){
01348       NurbsSurface<T,N>::refineKnotV(Xv) ;
01349     }
01350   }
01351 }
01352 
01371 template <class T, int N>
01372 int HNurbsSurface<T,N>::movePointOffset(T u, T v, const Point_nD<T,N>& delta){
01373   P = offset ; 
01374 
01375   // by definition the offset has w = 0 , but this isn't valid for
01376   // the control points, increasing the w by 1, will generate a valid surface
01377   if(baseLevel_)
01378     for(int i=0;i<P.rows();++i)
01379       for(int j=0;j<P.cols();++j){
01380         P(i,j).w() += T(1) ; 
01381       }
01382 
01383   if(NurbsSurface<T,N>::movePoint(u,v,delta)){
01384     offset = P ;
01385     // need to reset the offset weights
01386     if(baseLevel_)
01387       for(int i=0;i<P.rows();++i)
01388         for(int j=0;j<P.cols();++j){
01389           P(i,j).w() -= T(1) ; 
01390       }
01391     
01392     P = baseSurf.ctrlPnts() ; 
01393     updateSurface() ; 
01394     return 1 ;
01395   }
01396   updateSurface() ; 
01397   return 0 ;
01398 }
01399 
01408 template <class T, int N>
01409 void HNurbsSurface<T,N>::scale(const Point_nD<T,N>& s){
01410   for(int i=0;i<offset.rows();++i)
01411     for(int j=0;j<offset.cols();++j){
01412       offset(i,j).x() *= s.x() ; 
01413       offset(i,j).y() *= s.y() ; 
01414       offset(i,j).z() *= s.z() ; 
01415     }
01416 
01417   if(nextLevel_){
01418     nextLevel_->scale(s) ; 
01419   }
01420 }
01421 
01422 
01436 template <class T, int N>
01437 void HNurbsSurface<T,N>::setFixedOffsetVector(const Point_nD<T,N> &I, const Point_nD<T,N> &J, const Point_nD<T,N>& K){
01438   fixedOffset = 1;
01439   initBase();
01440   ivec(0,0) = I;
01441   jvec(0,0) = J;
01442   kvec(0,0) = K;
01443   updateSurface();
01444 }
01445 
01456 template <class T, int N>
01457 void HNurbsSurface<T,N>::setVariableOffsetVector(){
01458   fixedOffset = 0 ; 
01459   initBase();
01460   updateSurface();
01461 }
01462 
01463 } //end namespace