/*============================================================================= File: nurbsS.cpp Purpose: Revision: $Id: nurbsS.cpp,v 1.2 2002/05/13 21:07:46 philosophil Exp $ Created by: Philippe Lavoie (3 Oct, 1996) Modified by: Copyright notice: Copyright (C) 1996-1997 Philippe Lavoie This library is free software; you can redistribute it and/or modify it under the terms of the GNU Library General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This library is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public License for more details. You should have received a copy of the GNU Library General Public License along with this library; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. =============================================================================*/ #include #include #include #include #include "integrate.h" #ifdef USING_VCC #include #endif /*! */ namespace PLib { /*! \brief Default constructor \warning The surface is initialized to invalid values. Use a reset or a read function to set them to correct values. \author Philippe Lavoie \date 24 January, 1997 */ template NurbsSurface::NurbsSurface(): ParaSurface(), U(1),V(1),P(1,1),degU(0),degV(0) { } /*! \brief the copy constructor \param s the NurbsSurface to copy \author Philippe Lavoie \date 24 January, 1997 */ template NurbsSurface::NurbsSurface(const NurbsSurface& s): ParaSurface(), U(s.U), V(s.V), P(s.P), degU(s.degU), degV(s.degV) { } /*! \brief constructor with points in homogenous space \param DegU the degree in the $u$ direction \param DegV the degree in the $v$ direction \param Uk the $u$ knot vector \param Vk the $v$ knot vector \param Cp the matrix of control points in 4D \author Philippe Lavoie \date 24 January, 1997 */ template NurbsSurface::NurbsSurface(int DegU, int DegV, const Vector& Uk, const Vector& Vk, const Matrix< HPoint_nD >& Cp) : ParaSurface(), U(Uk),V(Vk),P(Cp),degU(DegU),degV(DegV) { int bad = 0 ; if(U.n() != P.rows()+degU+1){ #ifdef USE_EXCEPTION throw NurbsSizeError(P.rows(),U.n(),degU) ; #else Error err("NurbsSurface constructor") ; err << "The U knot vector and the number of rows of the control points are incompatible\n" ; err << "P.rows() = " << P.rows() << ", U.n() = " << U.n() << endl ; bad = 1 ; err.fatal() ; #endif } if(V.n() != P.cols()+degV+1){ #ifdef USE_EXCEPTION throw NurbsSizeError(P.cols(),V.n(),degV); #else Error err("NurbsSurface constructor") ; err << "The V knot vector and the number of columns of the control points are incompatible\n" ; err << "P.cols() = " << P.cols() << ", V.n() = " << V.n() << endl ; bad = 1 ; err.fatal() ; #endif } if(bad) exit(-1) ; } /*! \brief constructor with points in 3D \param DegU the degree of the surface in the U direction \param DegV the degree of the surface in the V direction \param Uk the U knot vector \param Vk the V knot vector \param Cp the matrix of 3D control points \param W the weight value for each control points \author Philippe Lavoie \date 24 January, 1997 */ template NurbsSurface::NurbsSurface(int DegU, int DegV, Vector& Uk, Vector& Vk, Matrix< Point_nD >& Cp, Matrix& W) : ParaSurface(), U(Uk),V(Vk),P(Cp.rows(),Cp.cols()),degU(DegU),degV(DegV) { int bad = 0 ; if(U.n() != Cp.rows()+degU+1){ #ifdef USE_EXCEPTION throw NurbsSizeError(P.rows(),U.n(),degU); #else Error err("NurbsSurface constructor") ; err << "The U knot vector and the number of rows of the control points are incompatible\n" ; err << "P.rows() = " << P.rows() << ", U.n() = " << U.n() << endl ; bad = 1 ; err.fatal() ; #endif } if(V.n() != Cp.cols()+degV+1){ #ifdef USE_EXCEPTION throw NurbsSizeError(P.cols(),V.n(),degV); #else Error err("NurbsSurface constructor") ; err << "The V knot vector and the number of columns of the control points are incompatible\n" ; err << "P.cols() = " << P.cols() << ", V.n() = " << V.n() << endl ; bad = 1 ; err.fatal() ; #endif } if(W.rows() != Cp.rows()){ #ifdef USE_EXCEPTION throw NurbsInputError(W.rows(),Cp.rows()) ; #else Error err("NurbsSurface constructor") ; err << "The dimension of the weights are incompatible with the dimension of the control poitns\n" ; err << "W( " << W.rows() << ", " << W.cols() << ") and P( " << P.rows() << ", " << P.cols() << ") \n" ; bad = 1 ; err.fatal() ; #endif } if(W.cols() != Cp.cols()){ #ifdef USE_EXCEPTION throw NurbsInputError(W.cols(),Cp.cols()) ; #else Error err("NurbsSurface constructor") ; err << "The dimension of the weights are incompatible with the dimension of the control poitns\n" ; err << "W( " << W.rows() << ", " << W.cols() << ") and P( " << P.rows() << ", " << P.cols() << ") \n" ; bad = 1 ; err.fatal() ; #endif } for(int i=0;i constructor") ; err << "Error initializing a NurbsSurface.\n" ; err << "The weight W(" << i << ", " << j << ") is equal to 0.\n" ; bad = 1 ; err.fatal() ; break ; #endif } else { P(i,j) = Cp(i,j) ; P(i,j) *= W(i,j) ; } } if(bad) exit(-1) ; } /*! \brief NurbsSurface assignment \param nS the NURBS surface to copy \author Philippe Lavoie \date 24 January, 1997 */ template NurbsSurface& NurbsSurface::operator=(const NurbsSurface& nS){ P = nS.P ; U = nS.U ; V = nS.V ; degU = nS.degU ; degV = nS.degV ; return *this ; } template struct AreaData { T v; T eps; T knotUi; T knotUii; T knotVj; T knotVjj; const NurbsSurface& s ; const BasicArray w ; AreaData(const NurbsSurface& surf,T e, const BasicArray& ww): s(surf),v(T(0)),eps(e),w(ww), knotUi(T(0)), knotUii(T(1)) {} }; template struct OpAreaAuxFcn : public ClassPOvoid { T operator()(T u, void* data){ AreaData* p = (AreaData*)data ; return (p->s).areaF(u,p->v) ; } }; template struct OpAreaFcn : public ClassPOvoid { T operator()(T v, void* data){ static Vector w; T err; OpAreaAuxFcn f; AreaData* Data = (AreaData*)data ; Data->v = v ; return intcc2((ClassPOvoid*)&f,Data, Data->knotUi,Data->knotUii, Data->eps,Data->w,err); } }; /*! \brief Computes the area of the surface Computes an approximation of the area of the surface using a numerical automatic integrator. That integrator uses a Chebyshev Series Expansion to perform its approximation. This is why you can change the value \a n which sets the number of elements in the series. The method is simple, integrate between each span. This is necessary in case the tangant of a point at u_i is undefined. Add the result and return this as the approximation. \param eps the accepted relative error \param n the number of element in the Chebyshev series \return the area of the NURBS surface. \author Alejandro Frangi \date 20 January 1999 */ template T NurbsSurface::area(T eps,int n) const { T a = T(0.0) ; T err ; static Vector bufFcn(0) ; if(bufFcn.n() != n){ bufFcn.resize(n) ; intccini(bufFcn) ; } AreaData data(*this,eps,bufFcn) ; OpAreaFcn op; for(int i=degU;i= U[i+1] || U[i]>=T(1.0)) continue ; data.knotUi = U[i] ; data.knotUii = U[i+1] ; for(int j=degV;j= V[j+1] || V[j]>=T(1.0)) continue ; data.knotVj = V[j] ; data.knotVjj = V[j+1] ; a += intcc2((ClassPOvoid*)&op,(void*)&data, data.knotVj,data.knotVjj, eps,bufFcn,err) ; } } return a ; } /*! \brief Computes the area of the surface inside [u_s,u_e] Computes an approximation of the area of the surface using a numerical automatic integrator. The area is computed for the range [u_s,u_e] That integrator uses a Chebyshev Series Expansion to perform its approximation. This is why you can change the value \a n which sets the number of elements in the series. The method is similar to the one used by area excepted that it needs to check for the range. \param us the starting range \param ue the end of the range \param vs the starting range \param ve the end of the range \param n the number of element in the Chebyshev series \param eps the accepted relative error \return the area of the NURBS surface. \warning ue (ve) must be greater than us (vs) and both must be in a valid range. \author Alejandro Frangi \date 20 January 1999 */ template T NurbsSurface::areaIn(T us, T ue, T vs, T ve, T eps, int n) const { T l = T() ; T err ; T a ; bool bLastU = false; bool bLastV = false; static Vector bufFcn ; if(bufFcn.n() != n){ bufFcn.resize(n) ; intccini(bufFcn) ; } AreaData data(*this,eps,bufFcn) ; OpAreaFcn op; for(int i=degU;i= U[i+1] || U[i] >= T(1)) continue; if(i=U[i] && ue<=U[i+findMultU(i)]){ data.knotUi = us; data.knotUii = ue; bLastU = true; goto Integrate_I; } if(us>=U[i]){ data.knotUi = us; data.knotUii = U[i+findMultU(i)]; bLastU = false; goto Integrate_I; } if(ue<=U[i+1]){ data.knotUi = U[i]; data.knotUii = ue; bLastU = true; goto Integrate_I; } data.knotUi = U[i] ; data.knotUii = U[i+findMultU(i)] ; Integrate_I: for(int j=degV;j= V[j+1] || V[j]>=T(1)) continue ; if(j=V[j] && ve<=V[j+findMultV(j)]){ data.knotVj = vs; data.knotVjj = ve; bLastV = true; goto Integrate_II; } if(vs>=V[j]){ data.knotVj = vs; data.knotVjj = V[j+findMultV(j)]; bLastV = false; goto Integrate_II; } if(ve<=V[j+1]){ data.knotVj = V[j]; data.knotVjj = ve; bLastV = true; goto Integrate_II; } data.knotVj = V[j] ; data.knotVjj = V[j+findMultV(j)] ; Integrate_II: a += intcc2((ClassPOvoid*)&op,(void*)&data,data.knotVj,data.knotVjj,eps,bufFcn,err) ; if (bLastU && bLastV) return a; if (bLastV) break; } } return a ; } // the definitions are in f_nurbs.cpp and d_nurbs.cpp /*! \a area needs to integrate a function over an interval to determine the area of the NURBS surface. Well, this is the function. \param u the parameter \param v the parameter \return the elemental area at (u,v) \author Alejandro Frangi \date 20 January 1999 */ template T NurbsSurface::areaF(T u, T v) const { Matrix > Skl(2,2) ; deriveAt(u,v,1,Skl); T tmp = norm(crossProduct(Skl(1,0),Skl(0,1))); return tmp ; } /*! \brief Determines if the surface is valid Determines if the surface is valid. The routine only verifies if the number of control points in the U and V direction matches the area of the U and V knot vectors. \return 1 if the surface is valid, 0 otherwise \author Philippe Lavoie \date 24 January, 1997 */ template int NurbsSurface::ok() { if(P.rows() <= degU) return 0 ; if(P.cols() <= degV) return 0 ; if(P.rows() != U.n()+degU+1) return 0 ; if(P.cols() != V.n()+degV+1) return 0 ; return 1 ; } /*! \brief Resize the surface Resize the surface. Proper values must be assigned once this function has been called since the resize operator is destructive. \param Pu the number of control points in the U direction \param Pv the number of control points in the V direction \param DegU the degree of the surface in the U direction \param DegV the degree of the surface in the V direction \author Philippe Lavoie \date 24 January, 1997 */ template void NurbsSurface::resize(int Pu, int Pv, int DegU, int DegV){ P.resize(Pu,Pv) ; degU = DegU ; degV = DegV ; U.resize(Pu+DegU+1) ; V.resize(Pv+DegV+1) ; } /*! \brief Resize the surface while keeping the old values. \param Pu the number of control points in the U direction \param Pv the number of control points in the V direction \param DegU the degree of the surface in the U direction \param DegV the degree of the surface in the V direction \author Philippe Lavoie \date 24 January, 1997 */ template void NurbsSurface::resizeKeep(int Pu, int Pv, int DegU, int DegV){ P.resizeKeep(Pu,Pv) ; degU = DegU ; degV = DegV ; U.resize(Pu+DegU+1) ; V.resize(Pv+DegV+1) ; } /*! \brief Generates a NURBS surface from skinning The NURBS surface is generated from skinning. The skinning is performed in the V direction. \param ca an array of NURBS curves \param degV the degree to skin in the V direction \param surf the skinned surface \return 0 if an error occurs, 1 otherwise \warning The number of curves to skin from must be greater than degV \author Philippe Lavoie \date 24 January, 1997 */ template int NurbsSurface::skinV(NurbsCurveArray& ca, int dV) { Vector vk(ca.n()) ; //Vector d(ca.n()) ; T* d ; int i,k,N,K ; if(ca.n() skinV") ; err << "The number of curves are insufficient for the degree of interpolation specified.\n" ; err << "n of curves = " << ca.n() << ", degree of interpolation = " << dV << endl ; err.warning() ; return 0 ; #endif } generateCompatibleCurves(ca) ; K = ca.n() ; N = ca[0].ctrlPnts().n() ; resize(N,K,ca[0].degree(),dV) ; //d.resize(ca[0].ctrlPnts().n()) ; d = new T[ca[0].ctrlPnts().n()] ; for(i=0;i > Q(K) ; NurbsCurve i_curve ; for(i=0;i int NurbsSurface::skinU(NurbsCurveArray& ca, int dU) { Vector uk(ca.n()); //Vector d(ca.n()) ; T* d; int i,k,N,K ; if(ca.n() skinU") ; err << "The number of curves are insufficient for the degree of interpolation specified.\n" ; err << "n of curves = " << ca.n() << ", degree of interpolation = " << dU << endl ; err.warning() ; return 0 ; #endif } generateCompatibleCurves(ca) ; K = ca.n() ; N = ca[0].ctrlPnts().n() ; resize(K,N,dU,ca[0].degree()) ; //d.resize(ca[0].ctrlPnts().n()) ; d = new T[ca[0].ctrlPnts().n()] ; for(i=0;i > Q(K) ; NurbsCurve i_curve ; for(i=0;i int surfMeshParams(const Matrix< Point_nD >& Q, Vector& uk, Vector& vl){ int n,m,k,l,num ; double d,total ; //Vector cds(Q.rows()) ; T* cds = new T[maximum(Q.rows(),Q.cols())] ; // alloca might not have enough space n = Q.rows() ; m = Q.cols() ; uk.resize(n) ; vl.resize(m) ; num = m ; // Compute the uk uk.reset(0) ; for(l=0;l int surfMeshParamsH(const Matrix< HPoint_nD >& Q, Vector& uk, Vector& vl){ int n,m,k,l,num ; double d,total ; //Vector cds(Q.rows()) ; T* cds = new T[maximum(Q.rows(),Q.cols())] ; n = Q.rows() ; m = Q.cols() ; uk.resize(n) ; vl.resize(m) ; num = m ; // Compute the uk uk.reset(0) ; for(l=0;l void NurbsSurface::globalInterp(const Matrix< Point_nD >& Q, int pU, int pV){ Vector vk,uk ; resize(Q.rows(),Q.cols(),pU,pV) ; surfMeshParams(Q,uk,vk) ; knotAveraging(uk,pU,U) ; knotAveraging(vk,pV,V) ; Vector< HPoint_nD > Pts(Q.rows()) ; NurbsCurve R ; int i,j ; for(j=0;j void NurbsSurface::globalInterpH(const Matrix< HPoint_nD >& Q, int pU, int pV){ Vector vk,uk ; resize(Q.rows(),Q.cols(),pU,pV) ; surfMeshParamsH(Q,uk,vk) ; knotAveraging(uk,pU,U) ; knotAveraging(vk,pV,V) ; Vector< HPoint_nD > Pts(Q.rows()) ; NurbsCurve R ; int i,j ; for(j=0;j void NurbsSurface::leastSquares(const Matrix< Point_nD >& Q, int pU, int pV, int nU, int nV){ Vector vk,uk ; resize(nU,nV,pU,pV) ; surfMeshParams(Q,uk,vk) ; Vector< HPoint_nD > Pts(Q.rows()) ; NurbsCurve R ; int i,j ; Matrix< HPoint_nD > P2 ; P2.resize(nU,Q.cols()) ; for(j=0;j void globalSurfInterpXY(const Matrix< Point_nD >& Q, int pU, int pV, NurbsSurface& S) { Vector uk,vk ; T um,uM ; T vm,vM ; um = Q(0,0).y() ; vm = Q(0,0).x() ; uM = Q(Q.rows()-1,0).y() ; vM = Q(0,Q.cols()-1).x() ; uk.resize(Q.rows()) ; vk.resize(Q.cols()) ; uk[0] = 0.0 ; vk[0] = 0.0 ; uk[uk.n()-1] = 1.0 ; vk[vk.n()-1] = 1.0 ; T dU = uM-um ; T dV = vM-vm ; int i ; for(i=1;i void globalSurfInterpXY(const Matrix< Point_nD >& Q, int pU, int pV, NurbsSurface& S, const Vector& uk, const Vector& vk){ Vector V,U ; int i,j ; knotAveraging(uk,pU,U) ; knotAveraging(vk,pV,V) ; Vector< HPoint_nD > P(Q.rows()) ; NurbsCurve R ; S.resize(Q.rows(),Q.cols(),pU,pV) ; S.U = U ; S.V = V ; for(j=0;j void globalSurfApprox(const Matrix< Point_nD >& Q, int pU, int pV, NurbsSurface& S, double error){ NurbsCurveArray R ; Vector< Point_nD > P ; Matrix< HPoint_nD > St ; Vector Ut,Vt ; Vector vk,uk ; surfMeshParams(Q,uk,vk) ; R.resize(Q.cols()) ; P.resize(Q.rows()) ; int i,j ; for(j=0;j void NurbsSurface::degreeElevate(int tU, int tV) { degreeElevateU(tU) ; degreeElevateV(tV) ; } /*! \brief Degree elevate the surface in the U direction \param t elevate the degree in the u direction by this amount. \author Philippe Lavoie \date 24 January, 1997 */ template void NurbsSurface::degreeElevateU(int t) { if(t<=0){ return ; } NurbsSurface S(*this) ; int i,j,k ; int n = S.ctrlPnts().rows()-1; int p = S.degU ; int m = n+p+1; int ph = p+t ; int ph2 = ph/2 ; Matrix bezalfs(p+t+1,p+1) ; // coefficients for degree elevating the Bezier segment Matrix< HPoint_nD > bpts(p+1,S.P.cols()) ; // pth-degree Bezier control points of the current segment Matrix< HPoint_nD > ebpts(p+t+1,S.P.cols()) ; // (p+t)th-degree Bezier control points of the current segment Matrix< HPoint_nD > Nextbpts(p-1,S.P.cols()) ; // leftmost control points of the next Bezier segment Vector alphas(p-1) ; // knot instertion alphas. // Compute the Binomial coefficients Matrix Bin(ph+1,ph2+1) ; binomialCoef(Bin) ; // Compute Bezier degree elevation coefficients T inv,mpi ; bezalfs(0,0) = bezalfs(ph,p) = 1.0 ; for(i=1;i<=ph2;i++){ inv= 1.0/Bin(ph,i) ; mpi = minimum(p,i) ; for(j=maximum(0,i-t); j<=mpi; j++){ bezalfs(i,j) = inv*Bin(p,j)*Bin(t,i-j) ; } } for(i=ph2+1;i= S.U[b+1]) // for some odd reasons... == doesn't work b++ ; mul = b-i+1 ; mh += mul+t ; ub = S.U[b] ; oldr = r ; r = p-mul ; if(oldr>0) lbz = (oldr+2)/2 ; else lbz = 1 ; if(r>0) rbz = ph-(r+1)/2 ; else rbz = ph ; if(r>0){ // Insert knot to get Bezier segment numer = ub-ua ; for(k=p;k>mul;k--){ alphas[k-mul-1] = numer/(S.U[a+k]-ua) ; } for(j=1;j<=r;j++){ save = r-j ; s = mul+j ; for(k=p;k>=s;k--){ for(colJ=0;colJ0) for(colJ=0;colJ1){ // Must remove knot u=c.U[a] oldr times // if(oldr>2) // Alphas on the right do not change // alfj = (ua-nc.U[kind-1])/(ub-nc.U[kind-1]) ; first = kind-2 ; last = kind ; den = ub-ua ; bet = (ub-U[kind-1])/den ; for(int tr=1; trtr){ // Loop and compute the new control points for one removal step if(i= lbz){ if(j-tr <= kind-ph+oldr){ gam = (ub-U[j-tr])/den ; for(colJ=0;colJ void NurbsSurface::degreeElevateV(int t) { if(t<=0){ return ; } NurbsSurface S(*this) ; int i,j,k ; int n = S.ctrlPnts().cols()-1; int p = S.degV ; int m = n+p+1; int ph = p+t ; int ph2 = ph/2 ; Matrix bezalfs(p+t+1,p+1) ; // coefficients for degree elevating the Bezier segment Matrix< HPoint_nD > bpts(p+1,S.P.rows()) ; // pth-degree Bezier control points of the current segment Matrix< HPoint_nD > ebpts(p+t+1,S.P.rows()) ; // (p+t)th-degree Bezier control points of the current segment Matrix< HPoint_nD > Nextbpts(p-1,S.P.rows()) ; // leftmost control points of the next Bezier segment //Vector alphas(p-1) ; // knot instertion alphas. T* alphas = (T*) alloca((p-1)*sizeof(T)) ; // Compute the Binomial coefficients Matrix Bin(ph+1,ph2+1) ; binomialCoef(Bin) ; // Compute Bezier degree elevation coefficients T inv,mpi ; bezalfs(0,0) = bezalfs(ph,p) = 1.0 ; for(i=1;i<=ph2;i++){ inv= 1.0/Bin(ph,i) ; mpi = minimum(p,i) ; for(j=maximum(0,i-t); j<=mpi; j++){ bezalfs(i,j) = inv*Bin(p,j)*Bin(t,i-j) ; } } for(i=ph2+1;i= S.V[b+1]) // for some odd reasons... == doesn't work b++ ; mul = b-i+1 ; mh += mul+t ; vb = S.V[b] ; oldr = r ; r = p-mul ; if(oldr>0) lbz = (oldr+2)/2 ; else lbz = 1 ; if(r>0) rbz = ph-(r+1)/2 ; else rbz = ph ; if(r>0){ // Insert knot to get Bezier segment numer = vb-va ; for(k=p;k>mul;k--){ alphas[k-mul-1] = numer/(S.V[a+k]-va) ; } for(j=1;j<=r;j++){ save = r-j ; s = mul+j ; for(k=p;k>=s;k--){ for(rowJ=0;rowJ0) for(rowJ=0;rowJ1){ // Must remove knot V=S.V[a] oldr times // if(oldr>2) // Alphas on the right do not change // alfj = (va-nc.U[kind-1])/(vb-V[kind-1]) ; first = kind-2 ; last = kind ; den = vb-va ; bet = (vb-V[kind-1])/den ; for(int tr=1; trtr){ // Loop and compute the new control points for one removal step if(i= lbz){ if(j-tr <= kind-ph+oldr){ gam = (vb-V[j-tr])/den ; for(rowJ=0;rowJ int NurbsSurface::findMultU(int r) const { int s=1 ; for(int i=r;i>degU+1;i--) if(U[i]<=U[i-1]) s++ ; else return s ; return s ; } /*! \brief finds the multiplicity of a knot in the V knot \param r the knot to observe \return the multiplicity of the V knot \warning \a r must be a valid knot index \author Philippe Lavoie \date 24 January, 1997 */ template int NurbsSurface::findMultV(int r) const { int s=1 ; for(int i=r;i>degV+1;i--) if(V[i]<=V[i-1]) s++ ; else return s ; return s ; } /*! \brief finds the span in the U and V direction Finds the span in the U and V direction. The spanU is the index \a k for which the parameter \a u is valid in the \a [u_k,u_{k+1}] range. The spanV is the index \a k for which the parameter \a v is valid in the \a [v_k,v_{k+1}] range. \param u find the U span for this parametric value \param v find the V span for this parametric value \param spanU the U span \param spanV the V span \author Philippe Lavoie \date 24 January, 1997 */ template void NurbsSurface::findSpan(T u, T v, int& spanU, int& spanV) const{ spanU = findSpanU(u) ; spanV = findSpanV(v) ; } /*! \brief finds the span in the U direction Finds the span in the U direction. The span is the index \a k for which the parameter \a u is valid in the \a [u_k,u_{k+1}] range. \param u --> find the span for this parametric value \return the span for \a u \author Philippe Lavoie \date 24 January, 1997 \modified 20 January, 1999 (Alejandro Frangi) */ template int NurbsSurface::findSpanU(T u) const{ if(u>=U[P.rows()]) return P.rows()-1 ; if(u<=U[degU]) return degU ; //AF int low = 0 ; int high = P.rows()+1 ; int mid = (low+high)/2 ; while(u= U[mid+1]){ if(u int NurbsSurface::findSpanV(T v) const{ if(v>=V[P.cols()]) return P.cols()-1 ; if(v<=V[degV]) return degV ; //AF int low = 0 ; int high = P.cols()+1 ; int mid = (low+high)/2 ; while(v= V[mid+1]){ if(v void NurbsSurface::basisFuns(T u, T v, int spanU, int spanV, Vector& Nu, Vector &Nv) const{ basisFunsU(u,spanU,Nu) ; basisFunsV(v,spanV,Nv) ; } /*! \brief Finds the non-zero basis function in the U direction \param u the u parametric value \param span the span of u \param N the vector containing the basis functions \author Philippe Lavoie \date 24 January, 1997 */ template void NurbsSurface::basisFunsU(T u, int span, Vector& N) const { //Vector left(degU+1), right(degU+1) ; T* left = (T*) alloca((degU+1)*sizeof(T)) ; T* right = (T*) alloca((degU+1)*sizeof(T)) ; T temp,saved ; N.resize(degU+1) ; N[0] = 1.0 ; for(int j=1; j<= degU ; j++){ left[j] = u-U[span+1-j] ; right[j] = U[span+j]-u ; saved = 0.0 ; for(int r=0 ; r void NurbsSurface::basisFunsV(T v, int span, Vector& N) const { //Vector left(degV+1), right(degV+1) ; T* left = (T*) alloca((degV+1)*sizeof(T)) ; T* right = (T*) alloca((degV+1)*sizeof(T)) ; T temp,saved ; N.resize(degV+1) ; N[0] = 1.0 ; for(int j=1; j<= degV ; j++){ left[j] = v-V[span+1-j] ; right[j] = V[span+j]-v ; saved = 0.0 ; for(int r=0 ; r void NurbsSurface::deriveAtH(T u, T v, int d, Matrix< HPoint_nD > &skl) const { int k,l,du,dv; skl.resize(d+1,d+1) ; du = minimum(d,degU) ; for(k=degU+1;k<=d;++k) for(l=0;l<=d-k;++l) skl(k,l) = 0.0 ; dv=minimum(d,degV) ; for(l=degV+1;l<=d;++l) for(k=0;k<=d-l;++k) skl(k,l) = 0.0 ; int uspan = findSpanU(u) ; int vspan = findSpanV(v) ; Matrix Nu,Nv ; nurbsDersBasisFuns(du,u,uspan,degU,U,Nu) ; nurbsDersBasisFuns(dv,v,vspan,degV,V,Nv) ; Vector< HPoint_nD > temp(degV+1) ; int dd,r,s ; for(k=0;k<=du;++k){ for(s=0;s<=degV;++s){ temp[s] = 0.0 ; for(r=0;r<=degU;++r) temp[s] += Nu(k,r)*P(uspan-degU+r,vspan-degV+s) ; } dd = minimum(d-k,dv) ; for(l=0;l<=dd;++l){ skl(k,l) = 0.0 ; for(s=0;s<=degV;++s) skl(k,l) += Nv(l,s)*temp[s] ; } } } /*! \brief Computes the point and the derivatives of degree \a d and below at \a (u,v) Computes the matrix of derivatives at \a u,v . The value of skl(k,l) represents the derivative of the surface \a S(u,v) with respect to \a u, \a k times and to \a v, \a l times. \param u the u parametric value \param v the v parametric value \param d the derivative is computed up to and including to this value \param skl the matrix containing the derivatives \author Philippe Lavoie \date 24 January, 1997 */ template void NurbsSurface::deriveAt(T u, T v, int d, Matrix< Point_nD > &skl) const { int k,l ; Matrix< HPoint_nD > ders ; Point_nD pv,pv2 ; skl.resize(d+1,d+1) ; deriveAtH(u,v,d,ders) ; Matrix Bin(d+1,d+1) ; binomialCoef(Bin) ; int i,j ; for(k=0;k<=d;++k){ for(l=0;l<=d-k;++l){ pv.x() = ders(k,l).x() ; pv.y() = ders(k,l).y() ; pv.z() = ders(k,l).z() ; for(j=1;j<=l;j++) pv -= Bin(l,j)*ders(0,j).w()*skl(k,l-j) ; for(i=1;i<=k;i++){ pv -= Bin(k,i)*ders(i,0).w()*skl(k-i,l) ; pv2 = 0.0 ; for(j=1;j<=l;j++) pv2 += Bin(l,j)*ders(i,j).w()*skl(k-i,l-j) ; pv -= Bin(k,i)*pv2 ; } skl(k,l) = pv/ders(0,0).w() ; } } } /*! \brief Computes the normal of the surface at \a (u,v) \param u the u parametric value \param v the v parametric value \return the normal at \a (u,v) . \author Philippe Lavoie \date 24 January, 1997 */ template Point_nD NurbsSurface::normal(T u, T v) const { Matrix< Point_nD > ders ; deriveAt(u,v,1,ders) ; return crossProduct(ders(1,0),ders(0,1)) ; } /*! \brief Returns the point on the surface at \a u,v Returns the point on the surface at \a u,v \param u the u parametric value \param v the v parametric value \return The homogenous point at \a u,v \author Philippe Lavoie \date 24 January, 1997 */ template HPoint_nD NurbsSurface::operator()(T u, T v) const{ int uspan = findSpanU(u) ; int vspan = findSpanV(v) ; Vector Nu,Nv ; Vector< HPoint_nD > temp(degV+1) ; basisFuns(u,v,uspan,vspan,Nu,Nv) ; int l; for(l=0;l<=degV;l++){ temp[l] =0.0 ; for(int k=0;k<=degU;k++){ temp[l] += Nu[k]*P(uspan-degU+k,vspan-degV+l) ; } } HPoint_nD sp(0,0,0,0) ; for(l=0;l<=degV;l++){ sp += Nv[l]*temp[l] ; } return sp ; } inline int max3(int a,int b, int c){ int m = a ; if(m void gordonSurface(NurbsCurveArray& lU, NurbsCurveArray& lV, const Matrix< Point_nD >& intersections, NurbsSurface& gS){ NurbsSurface sU,sV,sI ; sU.skinU(lU,3) ; sV.skinV(lV,3) ; sI.globalInterp(intersections,3,3) ; int du = max3(sU.degU,sV.degU,sI.degU) ; int dv = max3(sU.degV,sV.degV,sI.degV) ; NurbsSurface SU,SV,SI ; degreeElevate(sU,du-sU.degU,dv-sU.degV,SU) ; degreeElevate(sV,du-sV.degU,dv-sV.degV,SV) ; degreeElevate(sI,du-sI.degU,dv-sI.degV,SI) ; Vector U,V ; U = knotUnion(SU.knotU(),SV.knotU()) ; U = knotUnion(U,SI.knotU()) ; V = knotUnion(SU.knotV(),SV.knotV()) ; V = knotUnion(V,SI.knotV()) ; SU.mergeKnots(U,V) ; SV.mergeKnots(U,V) ; SI.mergeKnots(U,V) ; gS = SU ; for(int i=0;i inline void insert(T u, Vector& v){ int i ; if(uv[v.n()-1]) return ; v.resize(v.n()+1) ; i = v.n()-1 ; while(v[i-1]>u){ v[i] = v[i-1] ; --i ; } v[i] = u ; } /*! \brief Generates a surface by sweeping a curve along a trajectory Sweeping consists of creating a surface by moving a curve profile through a trajectory. The method uses here consists of using skinning of \a K instances of the curve \a C(u) along \a T(u). The \a K value should be viewed as the minimum number of sections required. The profile curve \a C(u) should lie on the xz-plane. It follows the trajectory curve \a T(u) along its y-axis. The scaling function is used to modify the shape of the curve profile while it's being swept. It can scale in any of the 4 dimensions to obtain the desired effects on the profile. See A10.1 on page 476 of the NURBS book for more details about the implementation. You might have to play with the useAy and invAz variables to obtain a satisfactory result for the sweep operation. This is either because there is an error in the code or because it is the way it is supposed to work. \param T the trajectory of the curve \param C the curve profile to sweep \param Sv a scaling function \param K the minimum number of insertion \param useAy a 0 indicates that rotation angle around the $y$-axis is not computed, otherwise it is computed. Results are usually better with a value of 0. \param invAz a 1 indicates that the computed rotation around the $z$-axis should be inversed, a 0 indicates it stays as it is. \return the NURBS surface representing the sweep of $C$ along $T$ \warning This will not yield a correct value for a closed trajectory curve. \author Philippe Lavoie \date 25 July, 1997 */ template void NurbsSurface::sweep(const NurbsCurve& Trj, const NurbsCurve& C, const NurbsCurve& Sv, int K, int useAy, int invAz) { int i,j,k,m,q,ktv,nsect ; q = Trj.degree() ; ktv = Trj.knot().n() ; nsect = K ; V.resize(Trj.knot().n()) ; V = Trj.knot() ; if(ktv <= nsect+q){ m = nsect+q-ktv+1 ; // insert m knots into T(v) // locations are not critical so inserting in the middle of the // biggest span is used. for(i=0;imd){ md = mt ; mu = (V[j]+V[j-1])/2.0 ; } } insert(mu,V) ; } } else{ if(ktv>nsect){ // must increase the number of instances of C(u) nsect = ktv-q-1 ; } } Vector v; // Compute the parameters by averaging the knots v.resize(nsect) ; v[0] = Trj.knot(Trj.degree()) ; v[nsect-1] = Trj.knot(Trj.knot().n()-Trj.degree()-1) ; for(k=1;k > B ; NurbsCurve Bv ; Vector< Point_nD > Td ; B.resize(v.n()) ; Trj.deriveAt(v[0],1,Td) ; B[0] = Td[1] ; if(Td[1].y() ==0) B[0] = crossProduct(Point_nD(0,1,0),B[0]) ; else B[0] = crossProduct(Point_nD(1,0,0),B[0]) ; B[0] /= norm(B[0]) ; for(i=1;i Ti(Td[1]) ; Ti = Ti/norm(Ti) ; Point_nD bi ; bi = B[i-1]-(B[i-1]*Ti)*Ti ; B[i] = bi/norm(bi) ; } Bv.globalInterp(B,v,minimum(3,B.n()-1)) ; Vector< HPoint_nD > Q(C.ctrlPnts().n()) ; for(k=0;k sk(Sv(v[k])) ; Q[i].x() = sk.x()*C.ctrlPnts(i).x() ; Q[i].y() = sk.y()*C.ctrlPnts(i).y() ; Q[i].z() = sk.z()*C.ctrlPnts(i).z() ; Q[i].w() = sk.w()*C.ctrlPnts(i).w() ; } // compute o(v[k]) Point_nD o = Trj.pointAt(v[k]) ; //T w = T(v[k]).w() ; // compute x(v[k]) Trj.deriveAt(v[k],1,Td) ; Point_nD x = Td[1]/norm(Td[1]) ; // compute z(v[k]) Point_nD z = Bv.pointAt(v[k]) ; z /= norm(z) ; // compute y(v[k]) Point_nD y = crossProduct(z,x) ; /* // compute the transform matrix double az = M_PI+atan2(y.y(),y.x()) ; double ax = -atan2(z.y(),z.z()) ; double ay = 0 ; if(useAy){ ay = atan2(x.z(),x.x()) ; } if(invAz){ az = -1.0*az ; } MatrixRT_DOUBLE A(ax,ay,az,o.x(),o.y(),o.z()) ; */ MatrixRT_DOUBLE R ; // M(4,4) R(0,0) = y.x(); R(1,0) = y.y(); R(2,0) = y.z(); R(0,1) = x.x(); R(1,1) = x.y(); R(2,1) = x.z(); R(0,2) = z.x(); R(1,2) = z.y(); R(2,2) = z.z(); //R(3,3) = 1.0; MatrixRT_DOUBLE Tx ; Tx.translate(o.x(),o.y(),o.z()); //MatrixRT_DOUBLE R(M) ; MatrixRT_DOUBLE A ; A = Tx * R ; for(i=0;i R; for(i=0;i void NurbsSurface::sweep(const NurbsCurve& Trj, const NurbsCurve& C, int K, int useAy,int invAz) { // setup Sv Vector< HPoint_nD > p(2) ; p[0] = HPoint_nD(1,1,1,1) ; p[1] = HPoint_nD(1,1,1,1) ; Vector u(4) ; u[0] = u[1] = 0.0 ; u[2] = u[3] = 1.0 ; NurbsCurve Sv(p,u,1) ; sweep(Trj,C,Sv,K,useAy,invAz) ; } /*! \brief Performs geometrical modifications Each control points will be modified by a rotation-translation matrix. \param A the rotation-translation matrix \author Philippe Lavoie \date 22 August 1997 */ template void NurbsSurface::transform(const MatrixRT& A){ for(int i=0;i void NurbsSurface::refineKnots(const Vector& nU, const Vector& nV){ refineKnotU(nU) ; refineKnotV(nV) ; } /*! \brief Refines the U knot vector \param X the knot vector to refine from \author Philippe Lavoie \date 24 January, 1997 */ template void NurbsSurface::refineKnotU(const Vector& X) { if(X.n()<=0) return ; int n = P.rows()-1 ; int p = degU; int m = n+p+1 ; int a,b ; int r = X.n()-1 ; NurbsSurface nS ; nS = *this ; nS.resize(r+1+n+1,P.cols(),degU,degV) ; a = findSpanU(X[0]) ; b = findSpanU(X[r]) ; ++b ; int j,col ; for(col=0;col=0 ; j--){ while(X[j] <= U[i] && i>a){ for(col=0;col void NurbsSurface::refineKnotV(const Vector& X) { if(X.n()<=0) return ; int n = P.cols()-1 ; int p = degV; int m = n+p+1 ; int a,b ; int r = X.n()-1 ; NurbsSurface nS ; try { nS = *this ; nS.resize(P.rows(),r+1+n+1,degU,degV) ; } catch(...){ cerr << "Out of memory\n" ; } a = findSpanV(X[0]) ; b = findSpanV(X[r]) ; ++b ; int j,row ; for(row=0;row=0 ; j--){ while(X[j] <= V[i] && i>a){ for(row=0;row void NurbsSurface::mergeKnots(const Vector& nU, const Vector& nV) { mergeKnotU(nU) ; mergeKnotV(nV) ; } /*! \brief merges the U knot vector with another one \param X a knot vector \warning The knot vector must be compatible with the U knot vector \author Philippe Lavoie \date 24 January, 1997 */ template void NurbsSurface::mergeKnotU(const Vector& X){ int i,ia,ib ; // Find the knots to insert Vector I(U.n()) ; int done = 0 ; i = ia = ib = 0 ; while(!done) { if(X[ib] == U[ia]){ ++ib ; ++ia ; } else{ I[i++] = X[ib] ; ib++ ; } done = (ia>=U.n() || ib >= X.n()) ; } I.resize(i) ; if(I.n()>0) refineKnotU(I) ; } /*! \brief merges the V knot vector with another one \param X a knot vector \warning The knot vector must be compatible with the V knot vector \author Philippe Lavoie \date 24 January, 1997 */ template void NurbsSurface::mergeKnotV(const Vector& X){ int i,ia,ib ; // Find the knots to insert Vector I(V.n()) ; int done = 0 ; i = ia = ib = 0 ; while(!done) { if(X[ib] == V[ia]){ ++ib ; ++ia ; } else{ I[i++] = X[ib] ; ib++ ; } done = (ia>=V.n() || ib >= X.n()) ; } I.resize(i) ; if(I.n()>0) refineKnotV(I) ; } /*! \brief Read a surface from an input stream \param fin the input file stream \return 0 if an error occurs, 1 otherwise \author Philippe Lavoie \date 24 January, 1997 */ template int NurbsSurface::read(ifstream &fin){ if(!fin) { return 0 ; } int nu,nv,du,dv; char *type ; type = new char[3] ; if(!fin.read(type,sizeof(char)*3)) { delete []type ; return 0 ;} int r1 = strncmp(type,"ns3",3) ; int r2 = strncmp(type,"ns4",3) ; if(!(r1 || r2)) return 0 ; int st ; char stc ; if(!fin.read((char*)&stc,sizeof(char))) { delete []type ; return 0 ;} if(!fin.read((char*)&nu,sizeof(int))) { delete []type ; return 0 ;} if(!fin.read((char*)&nv,sizeof(int))) { delete []type ; return 0 ;} if(!fin.read((char*)&du,sizeof(int))) { delete []type ; return 0 ;} if(!fin.read((char*)&dv,sizeof(int))) { delete []type ; return 0 ;} st = stc - '0' ; if(st != sizeof(T)){ // not of the same type size delete []type ; return 0 ; } resize(nu,nv,du,dv) ; if(!fin.read((char*)U.memory(),sizeof(T)*U.n())) { delete []type ; return 0 ;} if(!fin.read((char*)V.memory(),sizeof(T)*V.n())) { delete []type ; return 0 ;} T *p,*p2 ; if(!r1){ p = new T[3*nu*nv] ; if(!fin.read((char*)p,sizeof(T)*3*nu*nv)) { delete []type ; return 0 ;} p2 = p ; for(int i=0;i int NurbsSurface::read(const char* filename){ ifstream fin(filename) ; if(!fin) { return 0 ; } return read(fin) ; } /*! \brief Write a surface to a file stream \param fout the output filestream to write to. \return 1 on success, 0 on failure \author Philippe Lavoie \date 24 January, 1997 */ template int NurbsSurface::write(ofstream &fout) const { if(!fout) return 0 ; int prows = P.rows(); int pcols = P.cols(); char st = '0' + sizeof(T) ; if(!fout.write((char*)&"ns4",sizeof(char)*3)) return 0 ; if(!fout.write((char*)&st,sizeof(char))) return 0 ; if(!fout.write((char*)&prows,sizeof(int))) return 0 ; if(!fout.write((char*)&pcols,sizeof(int))) return 0 ; if(!fout.write((char*)°U,sizeof(int))) return 0 ; if(!fout.write((char*)°V,sizeof(int))) return 0 ; if(!fout.write((char*)U.memory(),sizeof(T)*U.n())) return 0 ; if(!fout.write((char*)V.memory(),sizeof(T)*V.n())) return 0 ; T *p,*p2 ; p = new T[P.rows()*P.cols()*4] ; p2 = p ; for(int i=0;i int NurbsSurface::write(const char* filename) const { ofstream fout(filename) ; if(!fout) return 0 ; return write(fout); } /*! \brief Transpose the U and V coordinates of a surface Transpose the U and V coordinates of a surface. After this operation the \a (u,v) points correspond to \a (v,u). \return A reference to itself \warning This is not completely tested \author Philippe Lavoie \date 24 January, 1997 */ template NurbsSurface& NurbsSurface::transpose(void){ Vector t(U) ; int d ; U = V ; V = t ; d = degU ; degU = degV ; degV = d ; P = ::transpose(P) ; return *this ; } /*! \brief Moves a point on the surface This moves the point \a s(u,v) by delta. \param u the parameter in the u direction \param v the parameter in the v direction \param delta the displacement of the point at S(u,v) \return 1 if the operation is possible, 0 if the problem is ill defined \e i.e. there isn't enough information to find a unique solution (the system is overdetermined) or that the system has non-independant components. \author Philippe Lavoie \date 24 January 1997 */ template int NurbsSurface::movePoint(T u, T v, const Point_nD& delta){ // setup B Matrix_DOUBLE B(1,(degU+1)*(degV+1)) ; int i,j,k ; int spanU,spanV ; Vector Ru,Rv ; B.reset(0.0) ; findSpan(u,v,spanU,spanV) ; nurbsBasisFuns(u,spanU,degU,U,Ru) ; nurbsBasisFuns(v,spanV,degV,V,Rv) ; for(j=0;j<=degU;++j){ for(k=0;k<=degV;++k){ B(0,j*(degV+1)+k) = (double)Ru[j]*(double)Rv[k] ; } } Matrix_DOUBLE A ; Matrix_DOUBLE Bt(::transpose(B)) ; Matrix_DOUBLE BBt ; BBt = inverse(B*Bt) ; A = Bt*BBt ; Matrix_DOUBLE dD(1,3) ; for(j=0;j<3;++j) dD(0,j) = (double)delta.data[j] ; Matrix_DOUBLE dP ; dP = A*dD ; i= 0 ; for(j=0;j<=degU;++j){ for(k=0;k<=degV;++k){ T w = P(spanU-degU+j,spanV-degV+k).w() ; P(spanU-degU+j,spanV-degV+k).x() += dP(i,0)*w ; P(spanU-degU+j,spanV-degV+k).y() += dP(i,1)*w ; P(spanU-degU+j,spanV-degV+k).z() += dP(i,2)*w ; ++i ; } } return 1 ; } /*! \brief Moves a point with some constraint This will modify the NURBS surface by respecting a certain number of constraints. \a u_r and \a v_r specifies the parameters on which the constraints should be applied. The constraint are defined by \a D_i(u,v) which requires 3 vectors to fully qualify. $D$ specifies the value of the constraint and \a Du and \a Dv are used to specify on which parameter the constraint is applied. ur and vr should be in an increasing order. \param ur the vector of parameters in the u direction \param vr the vector of parameters in the v direction \param D a vector of the value of \a D_i^(k,l)(u,v) \param Du a vector specifying the index of the value of u for \a D_i \param Dv a vector specifying the index of the value of v for \a D_i \return 1 if the operation is possible, 0 if the problem is ill defined \e i.e. there isn't enough information to find a unique solution (the system is overdetermined) or that the system has non-independant components. \warning The vectors defining \a D_i(u,v) should all be of the same size. \author Philippe Lavoie \date 24 January 1997 */ template int NurbsSurface::movePoint(const Vector& ur, const Vector& vr, const Vector< Point_nD >& D, const Vector_INT& Du, const Vector_INT& Dv) { Vector_INT Dk(Du.n()),Dl(Dv.n()) ; BasicArray fixCP(0) ; Dk.reset(0) ; Dl.reset(0) ; return movePoint(ur,vr,D,Du,Dv,Dk,Dl,fixCP) ; } /*! \brief Moves a point with some constraint This will modify the NURBS surface by respecting a certain number of constraints. \a u_r and \a v_r specifies the parameters on which the constraints should be applied. The constraint are defined by \a D_i^{(k,l)}(u,v) which requires 5 vectors to fully qualify. \a D specifies the value of the constraint and \a Du and \a Dv are used to specify on which parameter the constraint is applied and \a Dk and \a Dl specify the partial degree of the constraint. The values in \a D should be ordered in respect with i,k and l. ur and vr should be in an increasing order. \param ur the vector of parameters in the u direction \param vr the vector of parameters in the v direction \param D a vector of the value of \a D_i^(k,l)(u,v) \param Du a vector specifying the index of the value of u for \a D_i \param Dv a vector specifying the index of the value of v for \a D_i \param Dk a vector specifying the value of \a k for \a D_i \param Dl a vector specifying the value of \a l for \a D_i \return 1 if the operation is possible, 0 if the problem is ill defined \e i.e. there isn't enough information to find a unique solution (the system is overdetermined) or that the system has non-independant components. \warning The vectors defining \a D_i^{(k,l)}(u,v) should all be of the same size. \author Philippe Lavoie \date 24 January 1997 */ template int NurbsSurface::movePoint(const Vector& ur, const Vector& vr, const Vector< Point_nD >& D, const Vector_INT& Du, const Vector_INT& Dv, const Vector_INT& Dk, const Vector_INT& Dl) { BasicArray fixCP(0) ; return movePoint(ur,vr,D,Du,Dv,Dk,Dl,fixCP) ; } /*! \brief Moves a point with some constraint This will modify the NURBS surface by respecting a certain number of constraints. \a u_r and \a v_r specifies the parameters on which the constraints should be applied. The constraint are defined by $D_i^{(k,l)}(u,v)$ which requires 5 vectors to fully qualify. \a D specifies the value of the constraint and \a Du and \a Dv are used to specify on which parameter the constraint is applied and \a Dk and \a Dl specify the partial degree of the constraint. A second constraint \a fixCP consists of specifying which control points can not be moved by the routine. The values in D should be ordered in respect with i,k and l. ur and vr should be in an increasing order. \param ur the vector of parameters in the u direction \param vr the vector of parameters in the v direction \param D a vector of the value of \a D_i^(k,l)(u,v) \param Du a vector specifying the index of the value of u for \a D_i \param Dv a vector specifying the index of the value of v for \a D_i \param Dk a vector specifying the value of \a k for \a D_i \param Dl a vector specifying the value of \a l for \a D_i \param fixCP a vector specifying which control points can \e not be modified. \return 1 if the operation is possible, 0 if the problem is ill defined \e i.e. there isn't enough information to find a unique solution (the system is overdetermined) or that the system has non-independant components. \warning The vectors defining $D_i^{(k,l)}(u,v)$ should all be of the same size. \author Philippe Lavoie \date 24 January 1997 */ template int NurbsSurface::movePoint(const Vector& ur, const Vector& vr, const Vector< Point_nD >& D, const Vector_INT& Du, const Vector_INT& Dv, const Vector_INT& Dk, const Vector_INT& Dl, const BasicArray& fixCP) { int i,j,k,n ; if(D.n() != Du.n() || D.n() !=Du.n() || D.n() != Dv.n() || D.n() != Dv.n()){ #ifdef USE_EXCEPTION throw NurbsInputError(); #else Error err("movePoint(ur,D,Dr,Dk,fixCP)"); err << "The D,Dr,Dk vectors are not of the same size\n" ; err << "D.n()= " << D.n() << ", Du.n() = " << Du.n() << ", Dk.n() = " << Dk.n() << ", Dv.n() = " << Dv.n() << ", Dl.n() = " << Dl.n() << endl ; err.fatal() ; #endif } // setup B Matrix_DOUBLE B ; B.resize(D.n(),P.rows()*P.cols()) ; int spanU,spanV ; Matrix Ru,Rv ; B.reset(0.0) ; for(i=0;i map(B.cols()) ; remove.reset((int)1.0) ; for(j=0;j1e-10){ remove[j] = 0 ; break ; } } for(i=0;i& d = D[i] ; // this makes the SGI compiler happy for(j=0;j<3;++j) dD(i,j) = (double)d.data[j] ; } Matrix_DOUBLE dP ; dP = A*dD ; for(i=0;i void projectToLine(const Point_nD& S, const Point_nD& Trj, const Point_nD& pnt, Point_nD& p) { Point_nD a = pnt-S ; //p = S+ norm(a)*cos(angle(a,Trj))*Trj.unitLength() ; T fraction, denom ; denom = norm2(Trj) ; fraction = (denom == 0.0) ? 0.0 : (Trj*a) / denom ; p = fraction * Trj ; p += S ; } /*! \brief Generates a surface of revolution Generates a surface of revolution of a profile around an arbitrary axis (specified by a starting point S and a tangent T) with a certain angle. The angle is specified in radians. \param profile the curves to rotate around the axis \param S a point on the axis \param T the tangent vector of the rotation axis \param theta the angle of the revolution (in radians) \warning If a point is within a distance of 1e-7 from the axis, it will be assumed to be on the axis. The angle theta is only valid in the range [0,2 PI]. \author Philippe Lavoie \date 7 October, 1997 */ template void NurbsSurface::makeFromRevolution(const NurbsCurve& profile, const Point_nD& S, const Point_nD& Tvec, double theta){ double angle,dtheta ; int narcs ; int i,j ; //while(theta>2.0*M_PI) // theta -= 2.0*M_PI ; if(theta>2.0*M_PI) theta = 2.0*M_PI ; if(theta <= 0) theta = 0 ; if(theta P0,T0,P2,T2,P1 ; for(j=0;j O ; T wj = profile.ctrlPnts(j).w() ; projectToLine(S,Tvec,project(profile.ctrlPnts(j)),O) ; //projectToLine(S,Tvec,profile.ctrlPnts(j).projectW(),O) ; Point_nD X,Y ; X = project(profile.ctrlPnts(j))-O ; //X = profile.ctrlPnts(j).projectW() - O ; double r = norm(X) ; if(r < 1e-7){ for(i=0;i void NurbsSurface::makeFromRevolution(const NurbsCurve& profile, const Point_nD& S, const Point_nD& Tvec){ double angle,dtheta ; int narcs ; int i,j ; int n = 9 ; // n control points ; resize(n,profile.ctrlPnts().n(),2,profile.degree()) ; U[0] = U[1] = U[2] = 0 ; U[3] = U[4] = 0.25 ; U[5] = U[6] = 0.50 ; U[7] = U[8] = 0.75 ; U[9] = U[10] = U[11] = 1 ; V = profile.knot() ; const T wm = T(0.707106781185) ; // sqrt(2)/2 for(j=0;j O ; T wj = profile.ctrlPnts(j).w() ; projectToLine(S,Tvec,project(profile.ctrlPnts(j)),O) ; //projectToLine(S,Tvec,profile.ctrlPnts(j).projectW(),O) ; Point_nD X,Y ; X = project(profile.ctrlPnts(j))-O ; //X = profile.ctrlPnts(j).projectW() - O ; double r = norm(X) ; if(r < 1e-7){ for(i=0;i void NurbsSurface::makeFromRevolution(const NurbsCurve& profile){ int j ; int n = 9 ; // n control points ; resize(n,profile.ctrlPnts().n(),2,profile.degree()) ; U[0] = U[1] = U[2] = 0 ; U[3] = U[4] = 0.25 ; U[5] = U[6] = 0.50 ; U[7] = U[8] = 0.75 ; U[9] = U[10] = U[11] = 1 ; V = profile.knot() ; const T wm = T(0.707106781185) ; // sqrt(2)/2 for(j=0;j p = project(profile.ctrlPnts(j)) ; T r = T(sqrt(p.x()*p.x()+p.y()*p.y())); T wjm = wj*wm ; T rwjm = r*wj*wm ; T rwj = r*wj ; p.z() *= wj ; P(0,j) = HPoint_nD(rwj,0,p.z(),wj) ; P(1,j) = HPoint_nD(rwjm,rwjm,p.z()*wm,wjm) ; P(2,j) = HPoint_nD(0,rwj,p.z(),wj) ; P(3,j) = HPoint_nD(-rwjm,rwjm,p.z()*wm,wjm) ; P(4,j) = HPoint_nD(-rwj,0,p.z(),wj) ; P(5,j) = HPoint_nD(-rwjm,-rwjm,p.z()*wm,wjm) ; P(6,j) = HPoint_nD(0,-rwj,p.z(),wj) ; P(7,j) = HPoint_nD(rwjm,-rwjm,p.z()*wm,wjm) ; P(8,j) = HPoint_nD(rwj,0,p.z(),wj) ; } } /*! \brief Generates an iso curve in the U direction Generates an iso-parametric curve which goes through the parametric value u along the U direction. \param u the U parametric value \param c the iso-parametric curve \warning the parametric value $u$ must be in a valid range \author Philippe Lavoie \date 7 October, 1997 */ template void NurbsSurface::isoCurveU(T u, NurbsCurve& c) const { c.resize(P.cols(),degV) ; c.modKnot(V) ; if(u>U[U.n()-1]) u = U[U.n()-1] ; if(u N ; basisFunsU(u,span,N) ; HPoint_nD p ; for(int i=0;i void NurbsSurface::isoCurveV(T v, NurbsCurve& c) const { c.resize(P.rows(),degU) ; c.modKnot(U) ; if(v>V[V.n()-1]) v = V[V.n()-1] ; if(v N ; basisFunsV(v,span,N) ; HPoint_nD p ; for(int i=0;i int NurbsSurface::decompose(NurbsSurfaceArray& S) const { int i,m,a,b,nb,mult,j,r,save,s,k,row,col ; T numer,alpha ; //Vector alphas(degU+1) ; T* alphas = (T*) alloca((maximum(degU,degV)+1)*sizeof(T)) ; // all the surfaces will have the same knot vector in both the U and V // direction Vector nU ; nU.resize(2*(degU+1)) ; for(i=0;i nV ; nV.resize(2*(degV+1)) ; for(i=0;i Su ; Su.resize(P.rows()-degU) ; // )*(P.cols()-degV)) ; for(i=0;imult;j--) // compute and store the alphas alphas[j-mult-1] = numer/(U[a+j]-U[a]) ; r = degU-mult ; // insert knot r times for(j=1;j<=r;j++){ save=r-j; s=mult+j; // this many new points for(k=degU;k>=s;k--){ alpha = alphas[k-s] ; for(col=0;colmult;j--) // compute and store the alphas alphas[j-mult-1] = numer/(V[a+j]-V[a]) ; r = degV-mult ; // insert knot r times for(j=1;j<=r;j++){ save=r-j; s=mult+j; // this many new points for(k=degV;k>=s;k--){ alpha = alphas[k-s] ; for(row=0;row<=degU;++row) S[nb].P(row,k) = alpha*S[nb].P(row,k)+(1.0-alpha)*S[nb].P(row,k-1); } if(b int NurbsSurface::writePOVRAY(ostream& povray, int patch_type, double flatness, int num_u_steps, int num_v_steps) const { if(degU>3 || degV>3){ #ifdef USE_EXCEPTION throw NurbsInputError(); #else Error err("NurbsSurface writePOVRAY") ; err << "The degree of the surface is higher than 3!\n" << "A povrary file can not be generated.\n" ; err.warning() ; return 0; #endif } NurbsSurface S(*this) ; S.degreeElevate(3-degU,3-degV) ; NurbsSurfaceArray Sa ; S.decompose(Sa) ; int bad = 0 ; povray << "//\n//Generated for POV-Ray(tm) 3.0 by Phil's NURBS library\n" ; povray << "// http://yukon.genie.uottawa.ca/info/soft/nurbs\n//\n" ; for(int i=0;i p = project(Sa[i].ctrlPnts(j,k)) ; if(Sa[i].ctrlPnts(j,k).w()>1.0001 || Sa[i].ctrlPnts(j,k).w()<0.9999) bad = 1 ; povray << "\t<" << p.x() << ", " << p.y() << ", " << p.z() << ">" ; if(j==3 && k==3) povray << "\n}\n\n" ; else povray << ",\n " ; } povray << endl ; } } if(bad){ #ifdef USE_EXCEPTION throw NurbsWarning(); #else Error err("NurbsSurface writePOVRAY") ; err << "Warning: at least one of the control point was not rational\n" ; err << "The resulting surface will NOT be the same as the one which\n" ; err << "generated it.\n" ; err.warning() ; #endif } return bad ; } /*! \brief Writes the surface as a mesh of triangles Writes the surface as a mesh of triangles. You might have to change the values for the tolerance to get exactly what you're looking for. \param povray the output stream \param tolerance the tolerance when performing the tesselation \param color the color of the object \param diffuse the diffuse factor \param ambient the ambient factor \param smooth flags which indicates if we generate smooth triangles \return 1 on success, 0 on failure \warning It doesn't work very well. \author Philippe Lavoie \date 8 October, 1997 */ template int NurbsSurface::writePOVRAY(T tolerance, ostream& povray,const Color& color, int smooth, T ambient, T diffuse) const { BasicList > points ; BasicList connect ; BasicList > norm ; if(smooth) tesselate(tolerance,points,connect,&norm) ; else tesselate(tolerance,points,connect) ; povray << "mesh {\n" ; BasicNode *node ; BasicNode > *nodeP ; BasicNode > *nodeN ; nodeP = points.goToFirst() ; nodeN = 0 ; if(smooth) nodeN = norm.goToFirst() ; Vector< Point_nD > Pts(points.size()) ; Vector< Point_nD > Norm(norm.size()) ; for(int i=0;idata ; nodeP = points.goToNext() ; if(smooth){ Norm[i] = *nodeN->data ; nodeN = norm.goToNext() ; } } node = connect.goToFirst() ; while(node){ if(smooth) povray << "\tsmooth_triangle { " ; else povray << "\ttriangle { " ; povray << "< " << Pts[*node->data].x() << ", " << Pts[*node->data].y() << ", " << Pts[*node->data].z() << "> , " ; if(smooth) povray << "< " << Norm[*node->data].x() << ", " << Norm[*node->data].y() << ", " << Norm[*node->data].z() << "> , " ; node = connect.goToNext() ; povray << "< " << Pts[*node->data].x() << ", " << Pts[*node->data].y() << ", " << Pts[*node->data].z() << "> , " ; if(smooth) povray << "< " << Norm[*node->data].x() << ", " << Norm[*node->data].y() << ", " << Norm[*node->data].z() << "> , " ; node = connect.goToNext() ; povray << "< " << Pts[*node->data].x() << ", " << Pts[*node->data].y() << ", " << Pts[*node->data].z() << "> " ; if(smooth) povray << ", < " << Norm[*node->data].x() << ", " << Norm[*node->data].y() << ", " << Norm[*node->data].z() << "> " ; node = connect.goToNext() ; // skip the -1 value node = connect.goToNext() ; povray << "}\n" ; } povray << "\ttexture{ pigment { rgb < " << (double)color.r/255.0 << ", " << (double)color.g/255.0 << ", " << (double)color.b/255.0 << " >}\n" ; povray << "\t\tfinish { ambient " << ambient << " diffuse " << diffuse << " }\n\t}\n" ; povray << "}\n\n" ; return povray.good() ; } /*! \brief Writes a set of povray bicubic patches \param patch_type the type of the patch \param flatness the flatness level \param num_u_steps the minimum number of triangles in the U direction \param num_v_steps the minimum number of triangles in the V direction \return an ostream containing the definition of the surface \warning POVRAY only accepts rational spline patches. Thus you can't have a value other then 1.0 for the weights of your surface. A warning message will be generated if this is the case. POVRAY only deals with surfaces of degree 3. If the surface as a degree below 3 either in the U or V direction it will be elevated to be at 3 in both directions. If the surface as a degree higher then 3, then the function aborts. \author Philippe Lavoie \date 8 October, 1997 */ template