nurbsS_sp.cpp Source File

00001 /*=====================================================================
00002         File: nurbsS_sp.cpp
00003      Purpose:       
00004     Revision: $Id: nurbsS_sp.cpp,v 1.2 2002/05/13 21:07:46 philosophil Exp $
00005   Created by: Philippe Lavoie          (8 May, 1998)
00006  Modified by: 
00007 
00008  Copyright notice:
00009           Copyright (C) 1996-1997 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 <nurbsS_sp.h>
00026 //#include "nan.h"
00027 
00030 namespace PLib {
00031 
00044 template <class T, int N>
00045 void NurbsSurfaceSP<T,N>::updateMaxU() {
00046   if(degU>3){
00047 #ifdef USE_EXCEPTION
00048     throw NurbsInputError();
00049 #else
00050     Error error("NurbsSurfaceSP<T,N>::updateMaxU()") ;
00051     error << "This class doesn't support surfaces having a degree of 4 and higher.\n" ;
00052     error.fatal() ;
00053 #endif
00054   }
00055   else{
00056     maxU.resize(P.rows()) ;
00057     maxAtU_.resize(P.rows()) ;
00058     for(int i=0;i<P.rows();++i){
00059       if(!maxInfluence(i,U,degU,maxAtU_[i]))
00060         cerr << "Problem in maxInfluence U!\n" ;
00061       maxU[i] = nurbsBasisFun(maxAtU_[i],i,degU,U) ;
00062     }
00063     
00064   }
00065 }
00066 
00079 template <class T, int N>
00080 void NurbsSurfaceSP<T,N>::updateMaxV() {
00081   if(degV>3){
00082 #ifdef USE_EXCEPTION
00083     throw NurbsInputError();
00084 #else
00085     Error error("NurbsSurfaceSP<T,N>::updateMaxV()") ;
00086     error << "This class doesn't support surfaces having a degree of 4 and higher.\n" ;
00087     error.fatal() ;
00088 #endif
00089   }
00090   else{
00091     maxV.resize(P.cols()) ;
00092     maxAtV_.resize(P.cols()) ;
00093     for(int i=0;i<P.cols();++i){
00094       if(!maxInfluence(i,V,degV,maxAtV_[i]))
00095         cerr << "Problem in maxInfluence V!\n" ;
00096       maxV[i] = nurbsBasisFun(maxAtV_[i],i,degV,V) ;
00097     }
00098     
00099   }
00100 }
00101 
00123 template <class T, int N>
00124 NurbsSurfaceSP<T,N> NurbsSurfaceSP<T,N>::generateParallel(T d) const {
00125   NurbsSurfaceSP<T,N> p(*this) ;
00126 
00127   Vector< Point_nD<T,N> > offset(P.rows()*P.cols()) ;
00128   Vector<T> ur(P.rows()*P.cols()) ;
00129   Vector<T> vr(P.rows()*P.cols()) ;
00130   Vector_INT Du(P.rows()*P.cols()) ;
00131   Vector_INT Dv(P.rows()*P.cols()) ;
00132 
00133   Du.reset(0) ;
00134   Dv.reset(0) ;
00135 
00136   int i,j,no ;
00137 
00138   no = 0 ;
00139 
00140   for(i=0;i<P.rows();++i)
00141     for(j=0;j<P.cols();++j){
00142       Point_nD<T,N> norm ;
00143       norm = normal(maxAtU_[i],maxAtV_[j]) ;
00144       if(norm.x() == T(0) && 
00145          norm.y() == T(0) &&
00146          norm.z() == T(0)){
00147         // normal is undefined there...
00148         // compute an average and find a suitable normal
00149         const T delta = 0.00001 ;
00150         // must handle the corner cases
00151         int ok = 0 ; 
00152         if(i==0 && j==0){
00153           norm = normal(maxAtU_[i]+delta,maxAtV_[j]) ;
00154           norm += normal(maxAtU_[i],maxAtV_[j]+delta) ;
00155           norm /= T(2) ;
00156           ok = 1 ;
00157         }
00158         if(i==P.rows()-1 && j==P.cols()-1){
00159           norm = normal(maxAtU_[i]-delta,maxAtV_[j]) ;
00160           norm += normal(maxAtU_[i],maxAtV_[j]-delta) ;
00161           norm /= T(2) ;
00162           ok = 1 ;
00163         }
00164         if(i==0 && j==P.cols()-1){
00165           norm = normal(maxAtU_[i]-delta,maxAtV_[j]) ;
00166           norm += normal(maxAtU_[i],maxAtV_[j]+delta) ;
00167           norm /= T(2) ;
00168           ok = 1 ;
00169         }
00170         if(i==P.rows()-1 && j==0){
00171           norm = normal(maxAtU_[i]-delta,maxAtV_[j]) ;
00172           norm += normal(maxAtU_[i],maxAtV_[j]+delta) ;
00173           norm /= T(2) ;
00174           ok = 1 ;
00175         }
00176         if(!ok){
00177           T nt = 1.0 ; 
00178           while(norm.x() == T(0) && 
00179              norm.y() == T(0) &&
00180              norm.z() == T(0)){
00181             if( nt*d >(U[U.n()-1]-U[0])){
00182 #ifdef USE_EXCEPTION
00183               throw NurbsComputationError();
00184 #else
00185               Error error("generateParallel");
00186               error << "Can't compute a normal point.\n" ;
00187               error.fatal() ;
00188 #endif
00189             }
00190             T u1,u2,v1,v2 ;
00191             if(i==0 || i==P.rows()-1){
00192               u1 = u2 = maxAtU_[i] ;
00193               v1 = maxAtV_[j]+ nt*delta ;
00194               v2 = maxAtV_[j]- nt*delta ;
00195               if(v1>V[V.n()-1]) v1 = V[V.n()-1] ;
00196               if(v2<V[0]) v2 = V[0] ;
00197               norm = normal(u1,v1);
00198               norm += normal(u2,v2) ;
00199               norm /= 2 ; 
00200             }
00201             else{
00202               u1 = maxAtU_[i]- nt*delta ;
00203               u2 = maxAtU_[i]+ nt*delta ;
00204               v1 = v2 = maxAtV_[j] ;
00205               if(u1 < U[0]) u1 = U[0] ;
00206               if(u2 > U[U.n()-1]) u2 = U[U.n()-1] ;
00207 
00208               T u3,v3 ;
00209               u3 = maxAtU_[i] ;
00210               if(j==0)
00211                 v3 = maxAtV_[j]+ nt*delta ;
00212               else
00213                 v3 = maxAtV_[j]- nt*delta ;
00214 
00215               if(v3<V[0]) v3 = V[0] ;
00216               if(v3>V[V.n()-1]) v3 = V[V.n()-1] ;
00217 
00218               norm = normal(u1,v1);
00219               norm += normal(u2,v2) ;
00220               norm += normal(u3,v3) ;
00221               norm /= 3 ; 
00222             }
00223             nt *= 10.0 ; 
00224           }
00225         }
00226       }
00227       Point_nD<T,N> unit = norm.unitLength();
00228       unit *= d ;
00229       //HPoint_nD<T,N> offset(unit ) ;
00230       //offset.w() = 0.0 ; 
00231       //p.modSurfCPby(i,j,offset) ;
00232       Du[no] = i ;
00233       Dv[no] = j ;
00234       offset[no] = unit ;
00235       ++no ;
00236     }
00237 
00238   p.movePoint(maxAtU_,maxAtV_,offset,Du,Dv) ;
00239 
00240   return p ;
00241 }
00242 
00261 template <class T, int N>
00262 void NurbsSurfaceSP<T,N>::modOnlySurfCPby(int i, int j, const HPoint_nD<T,N>& a){
00263 
00264   int sizeU, sizeV ;
00265 
00266   sizeU = 2*degU+3 ; 
00267   if(i-degU-1<0) sizeU += i-degU-1 ; 
00268   if(i+degU+1>=P.rows()) sizeU -= i+degU+1-P.rows() ;
00269 
00270   sizeV = 2*degV+3 ;
00271   if(j-degV-1<0) sizeV += j-degV-1 ; 
00272   if(j+degV+1>=P.cols()) sizeV -= j+degV+1-P.cols() ;
00273   
00274   Vector<T> u(sizeU) ;
00275   Vector<T> v(sizeV) ;
00276   Vector<Point_nD<T,N> > pts(sizeU*sizeV) ; 
00277   Vector<int> pu(sizeU*sizeV) ;
00278   Vector<int> pv(sizeU*sizeV) ;
00279 
00280   int n=0;
00281   int nu = 0 ;
00282   int nv = 0 ; 
00283   for(int k=i-degU-1;k<=i+degU+1;++k){
00284     if(k<0)
00285       continue ;
00286     if(k>=P.rows())
00287       break ; 
00288     nv = 0 ;
00289     for(int l=j-degV-1;l<=j+degV+1;++l){
00290       if(l<0)
00291         continue ;
00292       if(l>=P.cols())
00293         break ; 
00294       if( k == i && j==l){
00295         pts[n].x() = a.x() ; 
00296         pts[n].y() = a.y() ; 
00297         pts[n].z() = a.z() ; 
00298       }
00299       //else
00300       //pts[n] = Point3D(0,0,0) ;
00301       pu[n] = nu ; 
00302       pv[n] = nv ; 
00303       if(k==i){
00304         v[nv] = maxAtV_[l] ; // only need to initialise this once
00305       }
00306       ++n ;
00307       ++nv ; 
00308     }  
00309     u[nu] = maxAtU_[k] ;
00310     ++nu ; 
00311   }
00312 
00313   u.resize(nu) ;
00314   v.resize(nv) ; 
00315   pts.resize(n) ;
00316   pu.resize(n) ; 
00317   pv.resize(n) ; 
00318 
00319   movePoint(u,v,pts,pu,pv) ;
00320 }
00321 
00322 } // end namespace