nurbs.h Source File

00001 /*=============================================================================
00002         File: nurbs.h
00003      Purpose:       
00004     Revision: $Id: nurbs.h,v 1.3 2003/01/15 03:03:59 philosophil Exp $
00005   Created by: Philippe Lavoie          (3 Oct, 1996)
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 #ifndef _nurbs_nurbs_h_
00026 #define _nurbs_nurbs_h_
00027 
00028 #include "matrixRT.h"
00029 
00030 #include "curve.h"
00031 #include "matrix.h"
00032 #include "matrixMat.h"
00033 #include "error.h"
00034 #include "image.h"
00035 
00038 namespace PLib {
00039 
00040   template <class T, int N> class NurbsSurface ;
00041   template <class T, int N> class NurbsCurve ;
00042   template <class T, int N> class NurbsCurveArray ;
00043 
00044   template <class T, int N> void generateCompatibleCurves(NurbsCurveArray<T,N> &ca);
00045 
00046 
00057   template <class T, int N>
00058     class NurbsCurve : public ParaCurve<T,N>{
00059     public:
00060       NurbsCurve() ;
00061       NurbsCurve(const NurbsCurve<T,N>& nurb) ;
00062       NurbsCurve(const Vector< HPoint_nD<T,N> >& P1, const Vector<T> &U1, int deg=3) ;
00063       NurbsCurve(const Vector< Point_nD<T,N> >& P1, const Vector<T> &W, const Vector<T> &U1, int deg=3) ;
00064       virtual ~NurbsCurve() { } // empty destructor
00065       
00066       // Reference to internal data
00067       int degree() const 
00068         { return deg_ ; }
00069       const Vector< HPoint_nD<T,N> >& ctrlPnts() const 
00070         { return P; }
00071       const HPoint_nD<T,N> ctrlPnts(int i) const 
00072         { return P[i] ; }
00073       const Vector<T>& knot() const 
00074         { return U ; }
00075       T knot(int i) const 
00076         { return U[i]; }
00077 
00078       // basic functions
00079       void resize(int n, int Deg) ;
00080       virtual void reset(const Vector< HPoint_nD<T,N> >& P1, const Vector<T> &U1, int deg) ;
00081       virtual NurbsCurve& operator=(const NurbsCurve<T,N>&) ;
00082       
00083       // Evaluation functions
00084       virtual HPoint_nD<T,N> operator()(T u) const;
00085       HPoint_nD<T,N> hpointAt(T u) const  
00086         { return operator()(u) ; }
00087       HPoint_nD<T,N> hpointAt(T u, int span) const ; 
00088       friend HPoint_nD<T,N> C(T u, const NurbsCurve<T,N>& nurb) {return nurb(u) ; } 
00089       friend Point_nD<T,N> Cp(T u, const NurbsCurve<T,N>& nurb) {return project(nurb(u)) ;} 
00090       
00091       
00092       // derivative functions
00093       void deriveAtH(T u, int, Vector< HPoint_nD<T,N> >&) const;
00094       void deriveAt(T u, int, Vector< Point_nD<T,N> >&) const;
00095       void deriveAtH(T u, int, int, Vector< HPoint_nD<T,N> >&) const;
00096       void deriveAt(T u, int, int, Vector< Point_nD<T,N> >&) const;
00097       Point_nD<T,N> derive3D(T u, int d) const;
00098       HPoint_nD<T,N> derive(T u, int d) const;
00099       Point_nD<T,N> normal(T u, const Point_nD<T,N>& v) const ;
00100       
00101       HPoint_nD<T,N> firstD(T u) const ; 
00102       HPoint_nD<T,N> firstD(T u, int span) const ; 
00103       Point_nD<T,N> firstDn(T u) const ; 
00104       Point_nD<T,N> firstDn(T u, int span) const ; 
00105       
00106       // Basis functions
00107       T basisFun(T u, int i, int p=-1) const ;
00108       void basisFuns(T u, int span, Vector<T>& N) const ;
00109       void dersBasisFuns(int n,T u, int span, Matrix<T>& N) const;
00110       
00111       // Knot functions
00112       T minKnot() const 
00113         { return U[0] ; } 
00114       T maxKnot() const 
00115         { return U[U.n()-1] ; }
00116       int findSpan(T u) const ;
00117       void findMultSpan(T u, int& r, int& s) const;
00118       int findMult(int r) const;
00119       int findKnot(T u) const;
00120       T getRemovalBnd(int r, int s) const;
00121       
00122       void removeKnot(int r, int s, int num);
00123       void removeKnotsBound(const Vector<T>& ub, Vector<T>& ek, T E) ;
00124       int knotInsertion(T u, int r,NurbsCurve<T,N>& nc);
00125       void refineKnotVector(const Vector<T>& X);
00126       void refineKnotVectorClosed(const Vector<T>& X);
00127       void mergeKnotVector(const Vector<T> &Um);
00128       
00129       void clamp() ; 
00130       void unclamp();
00131 
00132       // Curve fitting functions
00133       int leastSquares(const Vector< Point_nD<T,N> >& Q, int degC, int n ) ;
00134       int leastSquares(const Vector< Point_nD<T,N> >& Q, int degC, int n, const Vector<T>& ub);
00135       int leastSquaresH(const Vector< HPoint_nD<T,N> >& Q, int degC, int n, const Vector<T>& ub);
00136       int leastSquares(const Vector< Point_nD<T,N> >& Q, int degC, int n, const Vector<T>& ub, const Vector<T>& knot);
00137       int leastSquaresH(const Vector< HPoint_nD<T,N> >& Q, int degC, int n, const Vector<T>& ub, const Vector<T>& knot);
00138       
00139       int leastSquaresClosed(const Vector< Point_nD<T,N> >& Q, int degC, int n ) ;
00140       int leastSquaresClosed(const Vector< Point_nD<T,N> >& Q, int degC, int n, const Vector<T>& ub);
00141       int leastSquaresClosedH(const Vector< HPoint_nD<T,N> >& Q, int degC, int n, const Vector<T>& ub);
00142       int leastSquaresClosed(const Vector< Point_nD<T,N> >& Q, int degC, int n, const Vector<T>& ub, const Vector<T>& knot);
00143       int leastSquaresClosedH(const Vector< HPoint_nD<T,N> >& Q, int degC, int n, const Vector<T>& ub, const Vector<T>& knot);
00144       
00145       void globalApproxErrBnd(Vector< Point_nD<T,N> >& Q, int deg, T E);
00146       void globalApproxErrBnd(Vector< Point_nD<T,N> >& Q, Vector<T>& ub, int deg, T E);
00147       void globalApproxErrBnd2(Vector< Point_nD<T,N> >& Q, int degC, T E);
00148       void globalApproxErrBnd3(Vector< Point_nD<T,N> >& Q, int degC, T E);
00149       void globalApproxErrBnd3(Vector< Point_nD<T,N> >& Q, const Vector<T> &ub, int degC, T E);
00150       
00151       void globalInterp(const Vector< Point_nD<T,N> >& Q, int d);
00152       void globalInterp(const Vector< Point_nD<T,N> >& Q, const Vector<T>& ub, int d);
00153       void globalInterpH(const Vector< HPoint_nD<T,N> >& Q, int d);
00154       void globalInterpH(const Vector< HPoint_nD<T,N> >& Q, const Vector<T>& U, int d);
00155       void globalInterpH(const Vector< HPoint_nD<T,N> >& Q, const Vector<T>& ub, const Vector<T>& U, int d);
00156       
00157       void globalInterpClosed(const Vector< Point_nD<T,N> >& Qw, int d);
00158       void globalInterpClosed(const Vector< Point_nD<T,N> >& Qw, const Vector<T>& ub, int d);
00159       void globalInterpClosedH(const Vector< HPoint_nD<T,N> >& Qw, int d);
00160       void globalInterpClosedH(const Vector< HPoint_nD<T,N> >& Qw, const Vector<T>& U, int d);
00161       void globalInterpClosedH(const Vector< HPoint_nD<T,N> >& Qw, const Vector<T>& ub, const Vector<T>& U, int d);
00162       void globalInterpClosed(const Vector< Point_nD<T,N> >& Qw, const Vector<T>& ub, const Vector<T>& Uc, int d);
00163       
00164       void globalInterpD(const Vector< Point_nD<T,N> >& Q, const Vector< Point_nD<T,N> >& D, int d, int unitD, T a=1.0);
00165       
00166       void projectTo(const Point_nD<T,N>& p, T guess, T& u, Point_nD<T,N>& r, T e1=0.001, T e2=0.001,int maxTry=100) const;
00167       
00168       
00169       T length(T eps=0.001,int n=100) const ; 
00170       T lengthIn(T us, T ue, T eps=0.001, int n=100) const ; 
00171       T lengthF(T) const ;
00172       T lengthF(T,int) const ;
00173       
00174       // Generate type of curve
00175       void makeCircle(const Point_nD<T,N>& O, const Point_nD<T,N>& X, const Point_nD<T,N>& Y, T r, double as, double ae);
00176       void makeCircle(const Point_nD<T,N>& O, T r, double as, double ae);
00177       void makeCircle(const Point_nD<T,N>& O, T r);
00178       void makeLine(const Point_nD<T,N>& P0, const Point_nD<T,N>& P1, int d) ; 
00179       virtual void degreeElevate(int t);
00180       
00181 #ifndef HAVE_ISO_FRIEND_DECL
00182       friend void generateCompatibleCurves (NurbsCurveArray<T,N> &ca);
00183 #else
00184       friend void generateCompatibleCurves <>(NurbsCurveArray<T,N> &ca);
00185 #endif
00186 
00187       void decompose(NurbsCurveArray<T,N>& c) const;
00188       void decomposeClosed(NurbsCurveArray<T,N>& c) const ;
00189       
00190       int splitAt(T u, NurbsCurve<T,N>& cl, NurbsCurve<T,N>& cu) const;
00191       int mergeOf(const NurbsCurve<T,N>& cl, const NurbsCurve<T,N> &cu) ;
00192       
00193       // Modifies the NURBS curve
00194       void transform(const MatrixRT<T>& A) ;
00195       void modCP(int i,const HPoint_nD<T,N>& a) { P[i] = a ; } // To manipulate the value of the control point $P[i]$
00196       void modCPby(int i,const HPoint_nD<T,N>& a) { P[i] += a ; } // To manipulate the value of the control point $P[i]$
00197       virtual void modKnot(const Vector<T>& knotU) { if(knotU.n()-deg_-1==P.n()) U = knotU ; }  // to change the values of the knot vector only if the size is compatible with P.n
00198       
00199       int movePoint(T u, const Point_nD<T,N>& delta) ;
00200       int movePoint(T u, const BasicArray< Point_nD<T,N> >& delta) ;
00201       int movePoint(const BasicArray<T>& ur, const BasicArray< Point_nD<T,N> >& D);
00202       int movePoint(const BasicArray<T>& ur, const BasicArray< Point_nD<T,N> >& D, const BasicArray_INT& Dr, const BasicArray_INT& Dk) ;
00203       int movePoint(const BasicArray<T>& ur, const BasicArray< Point_nD<T,N> >& D, const BasicArray_INT& Dr, const BasicArray_INT& Dk, const BasicArray_INT& fixCP);
00204       
00205       void setTangent(T u, const Point_nD<T,N>& T0) ;
00206       void setTangentAtEnd(const Point_nD<T,N>& T0, const Point_nD<T,N>& T1) ;
00207 
00208       // I/O functions
00209       int read(const char*) ;
00210       int write(const char*) const ;
00211       virtual int read(ifstream &fin) ;
00212       int write(ofstream &fout) const ;
00213       int writePS(const char*,int cp=0,T magFact=T(-1),T dash=T(5), bool bOpen=true) const ;
00214       int writePSp(const char*,const Vector< Point_nD<T,N> >&,const Vector< Point_nD<T,N> >&, int cp=0,T magFact=0.0,T dash=5.0, bool bOpen=true) const ;
00215       
00216       int writeVRML(ostream &fout,T radius,int K, const Color& color,int Nu,int Nv, T u_s, T u_e) const ;
00217       int writeVRML(const char* filename,T radius,int K, const Color& color,int Nu,int Nv, T u_s, T u_e) const ;
00218       int writeVRML(const char* filename,T radius=1,int K=5, const Color& color=whiteColor,int Nu=20,int Nv=20) const { return writeVRML(filename,radius,K,color,Nu,Nv,U[0],U[U.n()-1]) ; } // writes the curve to a VRML file
00219       int writeVRML(ostream& fout,T radius=1,int K=5, const Color& color=whiteColor,int Nu=20,int Nv=20) const { return writeVRML(fout,radius,K,color,Nu,Nv,U[0],U[U.n()-1]) ; } // writes the curve to a VRML file
00220 
00221       int writeVRML97(const char* filename,T radius,int K, const Color& color,int Nu,int Nv, T u_s, T u_e) const ;
00222       int writeVRML97(ostream &fout,T radius,int K, const Color& color,int Nu,int Nv, T u_s, T u_e) const ;
00223       int writeVRML97(const char* filename,T radius=1,int K=5, const Color& color=whiteColor,int Nu=20,int Nv=20) const { return writeVRML97(filename,radius,K,color,Nu,Nv,U[0],U[U.n()-1]) ; } // writes the curve to a VRML file
00224       int writeVRML97(ostream& fout,T radius=1,int K=5, const Color& color=whiteColor,int Nu=20,int Nv=20) const { return writeVRML97(fout,radius,K,color,Nu,Nv,U[0],U[U.n()-1]) ; } // writes the curve to a VRML file
00225 
00226       int writeDisplayLINE(const char* filename, int iNu, const Color& color=blueColor,T fA=1) const ;
00227       int writeDisplayLINE(const char* filename,const Color& color, int iNu,T u_s, T u_e) const;
00228       void drawImg(Image_UBYTE& Img,unsigned char color=255,T step=0.01) ;
00229       void drawImg(Image_Color& Img,const Color& color,T step=0.01) ;
00230       void drawAaImg(Image_Color& Img, const Color& color, int precision=3,int alpha=1) ;
00231       void drawAaImg(Image_Color& Img, const Color& color, const NurbsCurve<T,3>& profile, int precision=3,int alpha=1) ;
00232       NurbsSurface<T,3> drawAaImg(Image_Color& Img, const Color& color, const NurbsCurve<T,3>& profile, const NurbsCurve<T,3> &scaling, int precision=3,int alpha=1) ;
00233       
00234       BasicList<Point_nD<T,N> > tesselate(T tolerance, BasicList<T> *uk) const ;
00235       
00236     protected:
00237       Vector< HPoint_nD<T,N> > P; // the vector of control points
00238       Vector<T> U ;  // the knot vector
00239       int deg_ ;  // the degree of the NURBS curve
00240     };
00241   
00242   typedef NurbsCurve<float,3> NurbsCurvef ;
00243   typedef NurbsCurve<double,3> NurbsCurved ;
00244   typedef NurbsCurve<float,2> NurbsCurve_2Df ;
00245   typedef NurbsCurve<double,2> NurbsCurve_2Dd ;
00246   
00247 } // end namespace
00248 
00249 typedef PLib::NurbsCurve<float,3> PlNurbsCurvef ;
00250 typedef PLib::NurbsCurve<double,3> PlNurbsCurved ;
00251 typedef PLib::NurbsCurve<float,2> PlNurbsCurve_2Df ;
00252 typedef PLib::NurbsCurve<double,2> PlNurbsCurve_2Dd ;
00253 
00256 namespace PLib {
00257   
00258   template <class T, int N>
00259     T chordLengthParam(const Vector< Point_nD<T,N> >& Q, Vector<T> &ub);
00260   template <class T, int N>
00261     T chordLengthParamH(const Vector< HPoint_nD<T,N> >& Q, Vector<T> &ub);
00262   template <class T, int N>
00263     T chordLengthParamClosed(const Vector< Point_nD<T,N> >& Q, Vector<T> &ub, int deg);
00264   template <class T, int N>
00265     T chordLengthParamClosedH(const Vector< HPoint_nD<T,N> >& Q, Vector<T> &ub, int deg);
00266   template <class T>
00267     void binomialCoef(Matrix<T>& Bin) ;
00268   template <class T>
00269     Vector<T> knotUnion(const Vector<T>& Ua, const Vector<T>& Ub);
00270   template <class T>
00271     T nurbsBasisFun(T u, int i, int p, const Vector<T>& U) ; 
00272   template <class T>
00273     void nurbsBasisFuns(T u, int span, int deg, const Vector<T>& U, Vector<T>& N);
00274   template <class T>
00275     void nurbsDersBasisFuns(int n, T u, int span, int deg, const Vector<T>& U, Matrix<T>& ders) ;
00276   
00277   template <class T, int N> int intersectLine(const Point_nD<T,N>& p1, const Point_nD<T,N>& t1, const Point_nD<T,N>& p2, const Point_nD<T,N>& t2, Point_nD<T,N>& p);
00278 #ifndef HAVE_TEMPLATE_OF_TEMPLATE
00279   template <> int intersectLine(const Point_nD<float,2>& p1, const Point_nD<float,2>& t1, const Point_nD<float,2>& p2, const Point_nD<float,2>& t2, Point_nD<float,2>& p);
00280   template <> int intersectLine(const Point_nD<double,2>& p1, const Point_nD<double,2>& t1, const Point_nD<double,2>& p2, const Point_nD<double,2>& t2, Point_nD<double,2>& p);
00281 #endif
00282   
00283   template <class T> void knotAveraging(const Vector<T>& uk, int deg, Vector<T>& U) ;
00284   template <class T> void knotAveragingClosed(const Vector<T>& uk, int deg, Vector<T>& U) ;
00285   template <class T> void knotApproximationClosed( Vector<T>& U, const  Vector<T>& ub, int n, int p);
00286   template <class T> void averagingKnots(const Vector<T>& U, int deg, Vector<T>& uk);
00287   template <class T> int findSpan(T u, const Vector<T>& U, int deg);
00288   template <class T> int maxInfluence(int i, const Vector<T>& U, int p, T &u);
00289   
00290   template <class T> void to3D(const NurbsCurve<T,2>&, NurbsCurve<T,3>&);
00291   template <class T> void to3D(const NurbsCurve<T,3>&, NurbsCurve<T,3>&);
00292   template <class T> void to2D(const NurbsCurve<T,3>&, NurbsCurve<T,2>&);
00293   
00294   template <class T, int N> 
00295     void wrapPointVector(const Vector<Point_nD<T,N> >& Q, int d, Vector<Point_nD<T,N> >& Qw);
00296   template <class T, int N> 
00297     void wrapPointVectorH(const Vector<HPoint_nD<T,N> >& Q, int d, Vector<HPoint_nD<T,N> >& Qw);
00298   
00299   
00317   template <class T, int N>
00318     inline int getCoordinates(const Point_nD<T,N>& p, int& i, int& j, int rows, int cols){
00319     i = int(rint(p.y())) ;
00320     j = int(rint(p.x())) ;
00321     if(i>=rows) return 0 ;
00322     if(j>=cols) return 0 ;
00323     if(i<0) return 0 ;
00324     if(j<0) return 0 ;
00325     return 1 ;
00326   }
00327   
00336   template <class T, int N>
00337     class NurbsCurveArray {
00338     public:
00339       int n () const { return sze ; }
00340       NurbsCurveArray(NurbsCurve<T,N>* Ca, int size) ;
00341       NurbsCurveArray() { C = 0 ; sze = 0 ; rsize = 0 ;}
00342       virtual ~NurbsCurveArray(){ if(C){ for(int i=0;i<rsize;i++) delete C[i];  delete []C ; }}
00343       
00344       virtual NurbsCurve<T,N>& operator[](int i) { return *(C[i]) ; }
00345       virtual NurbsCurve<T,N> operator[](int i) const { return *(C[i]) ; }
00346       
00347       virtual void resize(int s) ;
00348       void init(NurbsCurve<T,N>* Ca, int size) ;
00349       int read(const char *filename);
00350       int write(const char *filename);
00351       int writePS(const char*,int cp=0,T magFact=T(-1),T dash=T(5), bool bOpen=true) const ;
00352       int writePSp(const char*,const Vector< Point_nD<T,N> >&,const Vector< Point_nD<T,N> >&, int cp=0,T magFact=0.0,T dash=5.0, bool bOpen=true) const ;
00353       
00354     protected:
00355       NurbsCurve<T,N>& curve(int i) { return *(C[i]) ; }
00356       NurbsCurve<T,N> curve(int i) const { return *(C[i]) ; }
00357       
00358       int sze ; // the number of NURBS curves in the array
00359       int rsize ; // the number of space allocated for the array
00360       NurbsCurve<T,N>** C ; // An array of pointers to NURBS curves
00361     };
00362   
00363   typedef NurbsCurveArray<float,3> NurbsCurveArrayf ;
00364   typedef NurbsCurveArray<double,3> NurbsCurveArrayd ;
00365   typedef NurbsCurveArray<float,2> NurbsCurveArray_2Df ;
00366   typedef NurbsCurveArray<double,2> NurbsCurveArray_2Dd ;
00367   
00368 } // end namespace 
00369 
00370 typedef PLib::NurbsCurveArray<float,3> PlNurbsCurveArrayf ;
00371 typedef PLib::NurbsCurveArray<double,3> PlNurbsCurveArrayd ;
00372 typedef PLib::NurbsCurveArray<float,2> PlNurbsCurveArray_2Df ;
00373 typedef PLib::NurbsCurveArray<double,2> PlNurbsCurveArray_2Dd ;
00374 
00375 
00376 
00377 #ifdef INCLUDE_TEMPLATE_SOURCE
00378 #include "nurbs.cpp"
00379 #include "nurbsArray.cpp"
00380 #endif
00381 
00382 #endif 
00383