surface.cpp Source File

00001 /*=============================================================================
00002         File: surface.cpp
00003      Purpose:       
00004     Revision: $Id: surface.cpp,v 1.3 2003/01/13 19:42:29 philosophil Exp $
00005   Created by: Philippe Lavoie          (3 Oct, 1996)
00006  Modified by: 
00007 
00008  Copyright notice:
00009           Copyright (C) 1996-1999 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 
00026 #ifndef PLIB_NURBS_SURFACE_SOURCE
00027 #define PLIB_NURBS_SURFACE_SOURCE
00028 
00029 
00030 #include <surface.h>
00031 #include "matrixMat.h"
00032 #include <math.h>
00033 
00036 namespace PLib {
00037 
00071 template <class T, int N>
00072 T ParaSurface<T,N>::minDist2(const Point_nD<T,N>& p, T& guessU, T& guessV,T error,T s,int sep,int maxIter, T um, T uM, T vm, T vM) const {
00073   T d,d1,d2 ;
00074   Point_nD<T,N> p2 ;
00075   p2 = pointAt(guessU,guessV) ;
00076   d = norm2(p-p2) ;
00077   d2 = d1 = 0 ;
00078   int niter = 0 ;
00079   T u1,u2 ;
00080   T v1,v2 ;
00081   T step ;
00082   step = 2.0*s/(T)sep ;
00083   u1 = guessU-s ;
00084   u2 = guessU+s ;
00085   v1 = guessV-s ;
00086   v2 = guessV+s ;
00087   while(d>error && niter<maxIter) {
00088     if(u1<um)
00089       u1=um;
00090     if(u2>uM)
00091       u2 = uM ;
00092     if(v1<vm)
00093       v1 = vm ;
00094     if(v2>vM)
00095       v2 = vM ;
00096     T u,v ;
00097     d2 = d1 ;
00098     for(u=u1;u<u2;u+=step)
00099       for(v=v1;v<v2;v+=step){
00100         p2 = pointAt(u,v) ;
00101         d1 = norm2(p-p2) ;
00102         if(d1<d){
00103           d = d1 ;
00104           guessU = u ;
00105           guessV = v ;
00106         }
00107       }
00108     s /= 2.0 ;
00109     u1 = guessU - s ;
00110     u2 = guessU + s ;
00111     v1 = guessV - s ;
00112     v2 = guessV + s ;
00113     step = 2.0*s/(T)sep ;
00114     if(d-d2==0) niter = maxIter ;
00115     if(step<error) niter = maxIter ;
00116     ++niter;
00117   }
00118   return d ;
00119 }
00120 
00154 template <class T, int N>
00155 T ParaSurface<T,N>::minDist2b(const Point_nD<T,N>& p, T& guessU, T& guessV,T error,T s,int sep,int maxIter, T um, T uM, T vm, T vM) const {
00156   T d,d1,d2 ;
00157   Point_nD<T,N> p2 ;
00158   p2 = pointAt(guessU,guessV) ;
00159   d = norm2(p-p2) ;
00160   d2 = d1 = 0 ;
00161   int niter = 0 ;
00162   T u1,u2 ;
00163   T v1,v2 ;
00164   T step ;
00165   step = 2.0*s/(T)sep ;
00166   u1 = guessU-s ;
00167   u2 = guessU+s ;
00168   v1 = guessV-s ;
00169   v2 = guessV+s ;
00170   while(d>error && niter<maxIter) {
00171     if(u1<um)
00172       u1=um;
00173     if(u2>uM)
00174       u2 = uM ;
00175     if(v1<vm)
00176       v1 = vm ;
00177     if(v2>vM)
00178       v2 = vM ;
00179     T u,v ;
00180     d2 = d1 ;
00181     for(u=u1;u<=u2;u+=step)
00182       for(v=v1;v<=v2;v+=step){
00183         p2 = pointAt(u,v) ;
00184         d1 = norm2(p-p2) ;
00185         if(d1<d){
00186           d = d1 ;
00187           guessU = u ;
00188           guessV = v ;
00189         }
00190       }
00191     s = step ;
00192     u1 = guessU - s ;
00193     u2 = guessU + s ;
00194     v1 = guessV - s ;
00195     v2 = guessV + s ;
00196     step = s/2.0 ; // *s/(T)sep ;
00197     if(d-d2==0) niter = maxIter ;
00198     if(step<error) niter = maxIter ;
00199     ++niter;
00200   }
00201   return d ;
00202   
00203 }
00204 
00205 template <class T, int N>
00206 inline T to2power_xy(const Point_nD<T,N> &p){
00207   return (p.x()*p.x())+(p.y()*p.y()) ;
00208 }
00209 
00249 template <class T, int N>
00250 T ParaSurface<T,N>::minDist2xy(const Point_nD<T,N>& p, T& guessU, T& guessV,T error,T dU,T s,int sepU,int sepV,int maxIter, T um, T uM, T vm, T vM) const {
00251   T d,d1,d2,dz ;
00252   Point_nD<T,N> p2 ;
00253   p2 = pointAt(guessU,guessV) ;
00254   d = to2power_xy(p-p2) ;
00255   dz = to2power(p.z()-p2.z()) ;
00256   d2 = d1 = 0 ;
00257   int niter = 0 ;
00258   T u1,u2 ;
00259   T v1,v2 ;
00260   T stepU,stepV ;
00261   if(sepU>0){
00262     stepU = 2.0*s/(T)sepU ;
00263     u1 = guessU-s ;
00264     u2 = guessU+s ;
00265   }
00266   else{
00267     stepU = s ;
00268     u1 = guessU ;
00269     u2 = guessU ;
00270   }
00271   if(sepV>0){
00272     stepV = 2.0*s/(T)sepV ;
00273     v1 = guessV-s ;
00274     v2 = guessV+s ;
00275   }
00276   else{
00277     stepV = s ;
00278     v1 = guessV ;
00279     v2 = guessV ;
00280   }
00281   while(d>error && niter<maxIter) {
00282     if(u1<um)
00283       u1=um;
00284     if(u2>uM)
00285       u2 = uM ;
00286     if(v1<vm)
00287       v1 = vm ;
00288     if(v2>vM)
00289       v2 = vM ;
00290     T u,v ;
00291     d2 = d1 ;
00292     for(u=u1;u<=u2;u+=stepU)
00293       for(v=v1;v<=v2;v+=stepV){
00294         p2 = pointAt(u,v) ;
00295         d1 = to2power_xy(p-p2) ;
00296         if(d1<d){
00297           d = d1 ;
00298           dz = to2power(p.z()-p2.z()) ;
00299           guessU = u ;
00300           guessV = v ;
00301         }
00302       }
00303     if(d-d2==0) niter = maxIter ;
00304     if(stepU<dU) niter = maxIter ;
00305     if(stepV<dU) niter = maxIter ;
00306     s = stepU*0.55 ;
00307     u1 = guessU - s ;
00308     u2 = guessU + s ;
00309     s = stepV*0.55 ;
00310     v1 = guessV - s ;
00311     v2 = guessV + s ;
00312     stepU *= 0.5 ;
00313     stepV *= 0.5 ;
00314     ++niter;
00315   }
00316   return dz ;
00317   
00318 }
00319 
00343 template <class T, int N>
00344 int ParaSurface<T,N>::writeVRML(const char* filename,const Color& color,int Nu,int Nv, T uS,T uE,T vS, T vE) const {
00345   ofstream fout(filename) ;
00346 
00347   if(!fout)
00348     return 0 ;
00349   return writeVRML(fout,color,Nu,Nv,uS,uE,vS,vE) ;
00350 }
00351 
00352 
00376 template <class T, int N>
00377 int ParaSurface<T,N>::writeVRML(ostream &fout,const Color& color,int Nu,int Nv, T uS,T uE,T vS, T vE) const {
00378 
00379   fout << "#VRML V1.0 ascii\n" ;
00380   fout << "#  Generated by Phil's NURBS library\n" ;
00381   fout << "\nSeparator {\n" << "\tMaterialBinding { value PER_VERTEX_INDEXED }\n"  ;
00382     fout << "\tMaterial{\n\t\tambientColor 0.25 0.25 0.25\n\t\tdiffuseColor " << float(color.r/255.0) << ' ' << float(color.g/255.0) << ' ' << float(color.b/255.0) << "\n\t}\n" ;
00383   fout << "\tCoordinate3 {\n" ;
00384   fout << "\t\tpoint [\n" ;
00385 
00386   T u,v,du,dv ;
00387 
00388   if(Nu<=1)
00389     Nu = 2 ;  // Should I put a warning message ?
00390   if(Nv<=1)
00391     Nv = 2 ;  // Should I put a warning message ?
00392 
00393   u = uS ;
00394   v = vS ;
00395   du = (uE-uS)/(T)(Nu-1) ;
00396   dv = (vE-vS)/(T)(Nv-1) ;
00397 
00398   int i,j ;
00399 
00400 
00401   for(i=0;i<Nu;++i){
00402     v = vS ;
00403     for(j=0 ;j<Nv;++j){
00404       Point_nD<T,N> p ;
00405       p = pointAt(u,v) ;
00406       fout << "\t\t\t" << p.x() << ' ' << p.y() << ' ' << p.z() << ",\n" ;
00407       v += dv ;
00408     }
00409     u += du ;
00410   }
00411   
00412   fout << "\t\t]\n" ; // point [
00413   fout << "\t}\n" ; // cordinate3
00414 
00415   fout << "\tIndexedFaceSet{\n" ;
00416   fout << "\t\tcoordIndex[\n" ;
00417   
00418   for(i=0;i<Nu-1;++i){
00419     for(j=0;j<Nv-1;++j) {
00420       fout << "\t\t\t" << i*Nv+j << ", " << i*Nv+j+1 << ", " << (i+1)*Nv+j << ", -1,\n";
00421       fout << "\t\t\t" << i*Nv+j+1 << ", " << (i+1)*Nv+j+1<< ", " << (i+1)*Nv+j  << ", -1,\n";
00422     }
00423   }
00424   fout << "\t\t]\n" ;
00425   fout << "\t}\n" ; // IndexedFaceSet
00426 
00427   fout << "}\n" ;
00428 
00429   return 1 ;
00430 }
00431 
00432 template <class T>
00433 inline T compare(int findMin, T a, T b){
00434   if(findMin)
00435     return minimum(a,b) ;
00436   return maximum(a,b) ;
00437 }
00438 
00439 // VCC includes curve.cpp which defines this function already
00440 #ifndef INCLUDE_TEMPLATE_SOURCE
00441 template <class T,int N>
00442 inline T coordValue(CoordinateType coord, const Point_nD<T,N>& p){
00443   switch(coord){
00444   case coordX: return p.x() ; break ;
00445   case coordY: return p.y() ; break ;
00446   case coordZ: return p.z() ; break ;
00447   }
00448   return 0.0 ; // elliminates warning messages
00449 }
00450 #endif
00451 
00478 template <class T, int N>
00479 T ParaSurface<T,N>::extremum(int findMin, CoordinateType coord, T minDu, int sepU, int sepV, int maxIter, T um, T uM, T vm, T vM) const {
00480   T c,du,dv,d1,d2,result,guessU,guessV ;
00481   T minDv = minDu ;
00482   Point_nD<T,N> p2 ;
00483 
00484   guessU = 0 ;
00485   guessV = 0 ;
00486 
00487   // check for corner values...
00488   // because of the nature of += these value might be missed 
00489   // in the loop
00490   p2 = pointAt(um,vm) ;
00491   c = coordValue(coord,p2) ;
00492   p2 = pointAt(um,vM) ;
00493   c = compare(findMin,c,coordValue(coord,p2)) ;
00494   p2 = pointAt(uM,vm) ;
00495   c = compare(findMin,c,coordValue(coord,p2)) ;
00496   p2 = pointAt(uM,vM) ;
00497   c = compare(findMin,c,coordValue(coord,p2)) ; 
00498   result = c ;
00499   du = minDu*10.0 ;
00500   dv = minDv*10.0 ;
00501   d2 = d1 = 0 ;
00502   int niter = 0 ;
00503   T u1,u2 ;
00504   T v1,v2 ;
00505   T stepU,stepV ;
00506   T s ;
00507   s = uM - um ;
00508   stepU = s/(T)(sepU+1) ;
00509   stepV = s/(T)(sepV+1) ;
00510   u1 = um ;
00511   u2 = uM ;
00512   v1 = vm ;
00513   v2 = vM ;
00514   while((du>minDu || dv>minDv) && niter<maxIter) {
00515     if(u1<um)
00516       u1=um;
00517     if(u2>uM)
00518       u2 = uM ;
00519     if(v1<vm)
00520       v1 = vm ;
00521     if(v2>vM)
00522       v2 = vM ;
00523     T u,v ;
00524     d2 = c ;
00525     du = guessU ;
00526     dv = guessV ;
00527     for(u=u1;u<=u2;u+=stepU)
00528       for(v=v1;v<=v2;v+=stepV){
00529         p2 = pointAt(u,v) ;
00530         if(findMin){
00531           d1 = minimum(c,coordValue(coord,p2)) ;
00532           if(d1<c){
00533             c = d1 ;
00534             guessU = u ;
00535             guessV = v ;
00536             result = d1 ;
00537           }
00538         }
00539         else{
00540           d1 = maximum(c,coordValue(coord,p2)) ;
00541           if(d1>c){
00542             c = d1 ;
00543           guessU = u ;
00544           guessV = v ;
00545           result = d1 ;
00546           }
00547         }
00548       }
00549     s /= 2.0 ;
00550     u1 = guessU - s ;
00551     u2 = guessU + s ;
00552     v1 = guessV - s ;
00553     v2 = guessV + s ;
00554     stepU = 2.0*s/(T)sepU ;
00555     stepV = 2.0*s/(T)sepV ;
00556     if((c-d2)==0.0) niter = maxIter ;
00557     if(stepU<minDu) niter = maxIter ;
00558     if(stepV<minDv) niter = maxIter ;
00559     du = absolute(guessU-du) ;
00560     dv = absolute(guessV-dv) ;
00561     niter++;
00562   }
00563   return result ;
00564 }
00565 
00613 template <class T, int N>
00614 int ParaSurface<T,N>::projectOn(const Point_nD<T,N>& p, T& u, T& v, int maxI, const T um, const T uM, const T vm, const T vM) const {
00615   int i = 0 ;
00616   Point_nD<T,N> xu,xv,x,t ;
00617   Matrix< Point_nD<T,N> >ders ;
00618   Matrix_DOUBLE A(6,2) ;
00619   Matrix_DOUBLE B(6,1) ;
00620   Matrix_DOUBLE X(2,1) ;
00621 
00622 
00623 
00624   while(i < maxI){
00625     // make sure u and v are in bounds
00626     boundTo(u,um,uM);
00627     boundTo(v,vm,vM);
00628 
00629     // compute the over determined solution to the equation
00630     deriveAt(u,v,1,ders) ;
00631     xu = ders(1,0) ;
00632     xv = ders(0,1) ;
00633     x = ders(0,0) ;
00634     
00635     t = xu*xu ;
00636     A(0,0) = t.x() ;
00637     A(1,0) = t.y() ;
00638     A(2,0) = t.z() ;
00639 
00640     t = xu*xv ;
00641     A(0,1) = t.x() ;
00642     A(1,1) = t.y() ;
00643     A(2,1) = t.z() ;
00644 
00645     t = xv*xu ;
00646     A(3,0) = t.x() ;
00647     A(4,0) = t.y() ;
00648     A(5,0) = t.z() ;
00649 
00650     t = xv*xv ;
00651     A(3,1) = t.x() ;
00652     A(4,1) = t.y() ;
00653     A(5,1) = t.z() ;
00654 
00655     t = (p-x)*xu ;
00656     B(0,0) = t.x() ;
00657     B(1,0) = t.y() ;
00658     B(2,0) = t.z() ;
00659 
00660     t = (p-x)*xv ;
00661     B(3,0) = t.x() ;
00662     B(4,0) = t.y() ;
00663     B(5,0) = t.z() ;
00664     
00665 
00666     SVDMatrix<double> svd(A) ;
00667     if(!svd.solve(B,X))  // the matrix is singular
00668       return 0 ;
00669 
00670     // X now contains the du and dv ;
00671 
00672     if(T(X(0,0)) == T(0) && (T)X(1,0) == T(0)){
00673       // we are hopefully done
00674       return 1 ;
00675     }
00676     
00677     u += (T)X(0,0) ;
00678     v += (T)X(1,0) ;
00679     ++i ;
00680   }
00681 
00682   return 0 ;
00683 }
00684 
00708 template <class T, int N>
00709 int ParaSurface<T,N>::intersectWith(const ParaSurface<T,N> &S, Point_nD<T,N>& p, T& u, T& v, T& s, T& t, int maxI, T um, T uM, T vm, T vM) const{
00710   Point_nD<T,N> xu,xv,x1,x2,tmp,n1,n2,n12 ;
00711 
00712   Matrix< Point_nD<T,N> > ders ; 
00713   Matrix_DOUBLE A(6,2) ; 
00714   Matrix_DOUBLE B(6,1) ; 
00715   Matrix_DOUBLE X(2,1) ; 
00716   Matrix_DOUBLE A2(3,3) ;
00717   Matrix_DOUBLE B2(3,1) ; 
00718   Matrix_DOUBLE X2(3,1) ; 
00719 
00720   int done = 0 ; 
00721   int i = 0 ; 
00722   while(i<maxI){
00723     // make sure u and v are in bounds
00724     boundTo(u,um,uM) ; 
00725     boundTo(v,vm,vM) ; 
00726     boundTo(s,um,uM) ; 
00727     boundTo(t,vm,vM) ; 
00728 
00729     deriveAt(u,v,1,ders) ; 
00730     xu = ders(1,0) ; 
00731     xv = ders(0,1) ; 
00732     x1 = ders(0,0) ; 
00733     n1 = crossProduct(xu,xv) ; 
00734     
00735     tmp = xu*xu ; 
00736     A(0,0) = tmp.x() ;
00737     A(1,0) = tmp.y() ;
00738     A(2,0) = tmp.z() ;
00739 
00740     tmp = xu*xv ;
00741     A(0,1) = tmp.x() ;
00742     A(1,1) = tmp.y() ;
00743     A(2,1) = tmp.z() ;
00744 
00745     tmp = xv*xu ;
00746     A(3,0) = tmp.x() ;
00747     A(4,0) = tmp.y() ;
00748     A(5,0) = tmp.z() ;
00749 
00750     tmp = xv*xv ;
00751     A(3,1) = tmp.x() ;
00752     A(4,1) = tmp.y() ;
00753     A(5,1) = tmp.z() ;
00754 
00755     tmp = (p-x1)*xu ;
00756     B(0,0) = tmp.x() ;
00757     B(1,0) = tmp.y() ;
00758     B(2,0) = tmp.z() ;
00759 
00760     tmp = (p-x1)*xv ;
00761     B(3,0) = tmp.x() ;
00762     B(4,0) = tmp.y() ;
00763     B(5,0) = tmp.z() ;
00764     
00765 
00766     SVDMatrix<double> svd(A) ;
00767     if(!svd.solve(B,X))  // the matrix is singular
00768       return 0 ;
00769 
00770     // X now contains the du and dv ;
00771 
00772     if(T(X(0,0)) == T(0) && (T)X(1,0) == T(0)){
00773       done = 1 ; // we are hopefully done
00774       //return 1 ;
00775     }
00776     
00777     u += (T)X(0,0) ;
00778     v += (T)X(1,0) ;
00779 
00780     // now compute the ds and dt for the other surface
00781 
00782     S.deriveAt(s,t,1,ders) ; 
00783     xu = ders(1,0) ; 
00784     xv = ders(0,1) ; 
00785     x2 = ders(0,0) ; 
00786     n2 = crossProduct(xu,xv) ; 
00787     
00788     tmp = xu*xu ; 
00789     A(0,0) = tmp.x() ;
00790     A(1,0) = tmp.y() ;
00791     A(2,0) = tmp.z() ;
00792 
00793     tmp = xu*xv ;
00794     A(0,1) = tmp.x() ;
00795     A(1,1) = tmp.y() ;
00796     A(2,1) = tmp.z() ;
00797 
00798     tmp = xv*xu ;
00799     A(3,0) = tmp.x() ;
00800     A(4,0) = tmp.y() ;
00801     A(5,0) = tmp.z() ;
00802 
00803     tmp = xv*xv ;
00804     A(3,1) = tmp.x() ;
00805     A(4,1) = tmp.y() ;
00806     A(5,1) = tmp.z() ;
00807 
00808     tmp = (p-x2)*xu ;
00809     B(0,0) = tmp.x() ;
00810     B(1,0) = tmp.y() ;
00811     B(2,0) = tmp.z() ;
00812 
00813     tmp = (p-x2)*xv ;
00814     B(3,0) = tmp.x() ;
00815     B(4,0) = tmp.y() ;
00816     B(5,0) = tmp.z() ;
00817     
00818 
00819     svd.decompose(A) ;
00820     if(!svd.solve(B,X))  // the matrix is singular
00821       return 0 ;
00822 
00823     // X now contains the ds and dt ;
00824 
00825     if(T(X(0,0)) == T(0) && (T)X(1,0) == T(0)){
00826       if(done) {// we are hopefully done
00827         cerr << i << endl ; 
00828         return 1 ;
00829       }
00830       done = 0 ; 
00831     }
00832     
00833     s += (T)X(0,0) ;
00834     t += (T)X(1,0) ;
00835 
00836     // Finally we need to refine the point p 
00837     // this is one step behind refinement of F(u,v) and G(s,t)
00838 
00839     if(i>0){
00840       A2(0,0) = n1.x() ; 
00841       A2(0,1) = n1.y() ; 
00842       A2(0,2) = n1.z() ; 
00843       A2(1,0) = n2.x() ; 
00844       A2(1,1) = n2.y() ; 
00845       A2(1,2) = n2.z() ; 
00846       n12 = crossProduct(n1,n2) ; 
00847       A2(2,0) = n12.x() ; 
00848       A2(2,1) = n12.y() ; 
00849       A2(2,2) = n12.z() ; 
00850       
00851       B2(0,0) = x1*n1 ; 
00852       B2(1,0) = x2*n2 ; 
00853       const T alpha = 0.3 ;
00854       const T beta = 0.3 ; 
00855       const T lambda = 0.4 ; 
00856       B2(2,0) = (alpha*x1 + beta *x2 + lambda * p)*n12 ; 
00857 
00858       SVDMatrix<double> svd3(A2) ;
00859       if(!svd3.solve(B2,X2))  // the matrix is singular
00860         return 0 ;
00861       
00862       // X2 now contains the p2 ;
00863       p.x() = X2(0,0) ; 
00864       p.y() = X2(1,0) ; 
00865       p.z() = X2(2,0) ; 
00866     }
00867 
00868     ++i ;    
00869   }
00870   
00871   return 1 ; 
00872 }
00873 
00874 
00891 template <class T, int N>
00892 int ParaSurface<T,N>::intersectWith(const ParaSurface<T,N> &S, InterPoint<T,N> &iter, int maxI, T um, T uM, T vm, T vM) const{
00893   Point_nD<T,N> xu,xv,x1,x2,tmp,n1,n2,n12 ;
00894 
00895   Matrix< Point_nD<T,N> > ders ; 
00896   Matrix_DOUBLE A(6,2) ; 
00897   Matrix_DOUBLE B(6,1) ; 
00898   Matrix_DOUBLE X(2,1) ; 
00899   Matrix_DOUBLE A2(3,3) ;
00900   Matrix_DOUBLE B2(3,1) ; 
00901   Matrix_DOUBLE X2(3,1) ; 
00902 
00903   int done = 0 ; 
00904   int i = 0 ; 
00905   
00906   Point_nD<T,N> &p = iter.point ; 
00907   T &u = iter.paramA.u ; 
00908   T &v = iter.paramA.v ; 
00909   T &s = iter.paramB.u ; 
00910   T &t = iter.paramB.v ; 
00911 
00912   while(i<maxI){
00913     // make sure u and v are in bounds
00914     boundTo(u,um,uM) ; 
00915     boundTo(v,vm,vM) ; 
00916     boundTo(s,um,uM) ; 
00917     boundTo(t,vm,vM) ; 
00918 
00919     deriveAt(u,v,1,ders) ; 
00920     xu = ders(1,0) ; 
00921     xv = ders(0,1) ; 
00922     x1 = ders(0,0) ; 
00923     n1 = crossProduct(xu,xv) ; 
00924     
00925     tmp = xu*xu ; 
00926     A(0,0) = tmp.x() ;
00927     A(1,0) = tmp.y() ;
00928     A(2,0) = tmp.z() ;
00929 
00930     tmp = xu*xv ;
00931     A(0,1) = tmp.x() ;
00932     A(1,1) = tmp.y() ;
00933     A(2,1) = tmp.z() ;
00934 
00935     tmp = xv*xu ;
00936     A(3,0) = tmp.x() ;
00937     A(4,0) = tmp.y() ;
00938     A(5,0) = tmp.z() ;
00939 
00940     tmp = xv*xv ;
00941     A(3,1) = tmp.x() ;
00942     A(4,1) = tmp.y() ;
00943     A(5,1) = tmp.z() ;
00944 
00945     tmp = (p-x1)*xu ;
00946     B(0,0) = tmp.x() ;
00947     B(1,0) = tmp.y() ;
00948     B(2,0) = tmp.z() ;
00949 
00950     tmp = (p-x1)*xv ;
00951     B(3,0) = tmp.x() ;
00952     B(4,0) = tmp.y() ;
00953     B(5,0) = tmp.z() ;
00954     
00955 
00956     SVDMatrix<double> svd(A) ;
00957     if(!svd.solve(B,X)){  // the matrix is singular
00958       iter.tangent = crossProduct(n1,n2).unitLength() ; 
00959       return 0 ;
00960     }
00961 
00962     // X now contains the du and dv ;
00963 
00964     if(T(X(0,0)) == T(0) && (T)X(1,0) == T(0)){
00965       done = 1 ; // we are hopefully done
00966       //return 1 ;
00967     }
00968     
00969     u += (T)X(0,0) ;
00970     v += (T)X(1,0) ;
00971 
00972     // now compute the ds and dt for the other surface
00973 
00974     S.deriveAt(s,t,1,ders) ; 
00975     xu = ders(1,0) ; 
00976     xv = ders(0,1) ; 
00977     x2 = ders(0,0) ; 
00978     n2 = crossProduct(xu,xv) ; 
00979     
00980     tmp = xu*xu ; 
00981     A(0,0) = tmp.x() ;
00982     A(1,0) = tmp.y() ;
00983     A(2,0) = tmp.z() ;
00984 
00985     tmp = xu*xv ;
00986     A(0,1) = tmp.x() ;
00987     A(1,1) = tmp.y() ;
00988     A(2,1) = tmp.z() ;
00989 
00990     tmp = xv*xu ;
00991     A(3,0) = tmp.x() ;
00992     A(4,0) = tmp.y() ;
00993     A(5,0) = tmp.z() ;
00994 
00995     tmp = xv*xv ;
00996     A(3,1) = tmp.x() ;
00997     A(4,1) = tmp.y() ;
00998     A(5,1) = tmp.z() ;
00999 
01000     tmp = (p-x2)*xu ;
01001     B(0,0) = tmp.x() ;
01002     B(1,0) = tmp.y() ;
01003     B(2,0) = tmp.z() ;
01004 
01005     tmp = (p-x2)*xv ;
01006     B(3,0) = tmp.x() ;
01007     B(4,0) = tmp.y() ;
01008     B(5,0) = tmp.z() ;
01009     
01010 
01011     svd.decompose(A) ;
01012     if(!svd.solve(B,X)){  // the matrix is singular
01013       iter.tangent = crossProduct(n1,n2).unitLength() ; 
01014       return 0 ;
01015     }
01016 
01017     // X now contains the ds and dt ;
01018 
01019     if(T(X(0,0)) == T(0) && (T)X(1,0) == T(0)){
01020       if(done) {// we are hopefully done
01021         iter.tangent = crossProduct(n1,n2).unitLength() ; 
01022         return 1 ;
01023       }
01024       done = 0 ; 
01025     }
01026     
01027     s += (T)X(0,0) ;
01028     t += (T)X(1,0) ;
01029 
01030     // Finally we need to refine the point p 
01031     // this is one step behind refinement of F(u,v) and G(s,t)
01032 
01033     if(i>0){
01034       A2(0,0) = n1.x() ; 
01035       A2(0,1) = n1.y() ; 
01036       A2(0,2) = n1.z() ; 
01037       A2(1,0) = n2.x() ; 
01038       A2(1,1) = n2.y() ; 
01039       A2(1,2) = n2.z() ; 
01040       n12 = crossProduct(n1,n2) ; 
01041       A2(2,0) = n12.x() ; 
01042       A2(2,1) = n12.y() ; 
01043       A2(2,2) = n12.z() ; 
01044       
01045       B2(0,0) = x1*n1 ; 
01046       B2(1,0) = x2*n2 ; 
01047       const T alpha = 0.3 ;
01048       const T beta = 0.3 ; 
01049       const T lambda = 0.4 ; 
01050       B2(2,0) = (alpha*x1 + beta *x2 + lambda * p)*n12 ; 
01051 
01052       SVDMatrix<double> svd3(A2) ;
01053       if(!svd3.solve(B2,X2)){  // the matrix is singular
01054         iter.tangent = crossProduct(n1,n2).unitLength() ; 
01055         return 0 ;
01056       }
01057       
01058       // X2 now contains the p2 ;
01059       p.x() = X2(0,0) ; 
01060       p.y() = X2(1,0) ; 
01061       p.z() = X2(2,0) ; 
01062     }
01063 
01064     ++i ;    
01065   }
01066   
01067   iter.tangent = crossProduct(n1,n2).unitLength() ; 
01068   return 1 ; 
01069 }
01070 
01071 template <class T, int N>
01072 inline int isNear(const SurfParam<T,N> &a, const SurfParam<T,N>& b, double tol=1e-5){
01073   double d2 = ((double)a.u-(double)b.u)*((double)a.u-(double)b.u);
01074   d2 += ((double)a.v-(double)b.v)*((double)a.v-(double)b.v) ;
01075   if(d2<tol*tol)
01076     return 1 ;
01077   return 0 ; 
01078 }
01079 
01080 template <class T, int N>
01081 inline int isNear(const InterPoint<T,N>& a, const InterPoint<T,N> &b){
01082   return isNear(a.paramA,b.paramA) || isNear(a.paramB,b.paramB) ; 
01083 }
01084 
01085 template <class T, int N>
01086 inline int onBoundary(const SurfParam<T,N>& a, T m = 0, T M=1){
01087   if(a.u<=m)
01088     return 1 ;
01089   if(a.u>=M)
01090     return 1 ; 
01091   return 0 ; 
01092 }
01093 
01094 template <class T, int N>
01095 void intersectSurfaces(const ParaSurface<T,N> &surfA, const ParaSurface<T,N> &surfB, BasicList<InterPoint<T,N> > &points, int maxI, T um, T uM, T vm, T vM){
01096   points.reset() ;
01097 
01098   Point_nD<T,N> p ;
01099 
01100   p = surfA.pointAt(0.5,0.5) ;
01101 
01102   InterPoint<T,N> I0,I,Ilast ;
01103 
01104   I0.point = p ;
01105 
01106   surfA.intersectWith(surfB,I0,maxI,um,uM,vm,vM) ;
01107 
01108   I = Ilast = I0 ; 
01109 
01110   T d ;
01111 
01112   d = 0.01 ; // should be dependant on the control points locations
01113   T direction = 1 ; 
01114 
01115   int closed = 0 ; 
01116   int reach_bound = 0 ;
01117 
01118   const T error  = 0.1; 
01119   const T up_bound = 1.5 ;
01120 
01121   int n = 0 ; 
01122 
01123   while(1){
01124     points.add(I) ; 
01125     I.point += direction*d*I.tangent ; 
01126     surfA.intersectWith(surfB,I,maxI,um,uM,vm,vM) ; 
01127     if(isNear(I0,I)){
01128       closed = 1 ; 
01129       break ; 
01130     }
01131     if(onBoundary(I.paramA) || onBoundary(I.paramB)){
01132       reach_bound =1 ;
01133       break ; 
01134     }
01135     d = norm(Ilast.point-I.point) ; 
01136     d *= error/acos(I.tangent*Ilast.tangent/up_bound) ;
01137     if(d<0.01) // putting a lower boundary
01138       d = 0.01 ; 
01139     if(d>100)
01140       break ; 
01141     Ilast = I ; 
01142     cout << I.point << "\t" << I.paramA.u << "," << I.paramA.v << ":" << I.paramB.u << "," << I.paramB.v << "\t" << I.tangent << "|" << d << endl ; 
01143     ++n ; 
01144     if(n>100) 
01145       break ; 
01146   }
01147   if(reach_bound){
01148     // repeat but in the other direction
01149     direction *= -1 ;
01150     I = Ilast = I0 ; 
01151     n = 0 ; 
01152     while(1){
01153       points.add(I) ; 
01154       I.point += direction*d*I.tangent ; 
01155       surfA.intersectWith(surfB,I,maxI,um,uM,vm,vM) ; 
01156       if(isNear(I0,I)){
01157         closed = 1 ; 
01158         break ; 
01159       }
01160       if(onBoundary(I.paramA) || onBoundary(I.paramB)){
01161         reach_bound =1 ;
01162         break ; 
01163       }
01164       d = norm(Ilast.point-I.point) ; 
01165       d *= error/acos(I.tangent*Ilast.tangent/up_bound) ;
01166       if(d<0.01) // putting a lower boundary
01167         d = 0.01 ; 
01168       if(d>100)
01169         break ; 
01170       Ilast = I ; 
01171       cout << I.point << "\t" << I.paramA.u << "," << I.paramA.v << ":" << I.paramB.u << "," << I.paramB.v << "\t" << I.tangent << "|" << d << endl ; 
01172       ++n ; 
01173       if(n>100) 
01174         break ; 
01175     }    
01176   }
01177 }
01178 
01202 template <class T, int N>
01203 int ParaSurface<T,N>::writeVRML97(const char* filename,const Color& color,int Nu,int Nv, T uS,T uE,T vS, T vE) const {
01204   ofstream fout(filename) ;
01205 
01206   if(!fout)
01207     return 0 ;
01208   return writeVRML97(fout,color,Nu,Nv,uS,uE,vS,vE) ;
01209 }
01210 
01211 
01235 template <class T, int N>
01236 int ParaSurface<T,N>::writeVRML97(ostream &fout,const Color& color,int Nu,int Nv, T uS,T uE,T vS, T vE) const {
01237 
01238   fout << "#VRML V2.0 utf8\n" ;
01239   fout << "#  Generated by Phil's NURBS library\n" ;
01240   fout << "\nGroup {\n" ;
01241   fout << "\n  children [\n" ; 
01242   //fout << "\tDEF PS SphereSensor {}\n" ;
01243   fout << "\tDEF T Transform {\n"; 
01244   fout << "\t  children [\n" ;
01245   fout << "\t\tShape {\n" ;
01246   fout << "\t\t appearance Appearance {\n" ; 
01247   fout << "\t\t\tmaterial Material{ diffuseColor " << float(color.r/255.0) << ' ' << float(color.g/255.0) << ' ' << float(color.b/255.0) << " }\n" ;
01248   fout << "\t\t }\n" ; 
01249   fout << "\t\t geometry IndexedFaceSet {\n" ;
01250   fout << "\t\t\tsolid FALSE\n" ; 
01251   fout << "\t\t\tcoord Coordinate {\n" ;
01252   fout << "\t\t\t point [\n" ;
01253 
01254   T u,v,du,dv ;
01255 
01256   if(Nu<=1)
01257     Nu = 2 ;  // Should I put a warning message ?
01258   if(Nv<=1)
01259     Nv = 2 ;  // Should I put a warning message ?
01260 
01261   u = uS ;
01262   v = vS ;
01263   du = (uE-uS)/(T)(Nu-1) ;
01264   dv = (vE-vS)/(T)(Nv-1) ;
01265 
01266   int i,j ;
01267 
01268   Point_nD<T,N> p_min = pointAt(u,v) ;
01269   Point_nD<T,N> p_max = pointAt(u,v) ; 
01270 
01271   for(i=0;i<Nu;++i){
01272     v = vS ;
01273     for(j=0 ;j<Nv;++j){
01274       Point_nD<T,N> p ;
01275       p = pointAt(u,v) ;
01276       fout << "\t\t\t\t" << p.x() << ' ' << p.y() << ' ' << p.z() << ",\n" ;
01277       v += dv ;
01278       if(p.x()<p_min.x()) p_min.x() = p.x();
01279       if(p.y()<p_min.y()) p_min.y() = p.y();
01280       if(p.z()<p_min.z()) p_min.z() = p.z();
01281       if(p.x()>p_max.x()) p_max.x() = p.x();
01282       if(p.y()>p_max.y()) p_max.y() = p.y();
01283       if(p.z()>p_max.z()) p_max.z() = p.z();
01284     }
01285     u += du ;
01286   }
01287   
01288   fout << "\t\t\t ]\n" ; // point [
01289   fout << "\t\t\t}\n" ; // coord
01290 
01291   fout << "\t\t\t coordIndex [\n" ;
01292   
01293   for(i=0;i<Nu-1;++i){
01294     for(j=0;j<Nv-1;++j) {
01295       fout << "\t\t\t\t" << i*Nv+j << ", " << i*Nv+j+1 << ", " << (i+1)*Nv+j << ", -1,\n";
01296       fout << "\t\t\t\t" << i*Nv+j+1 << ", " << (i+1)*Nv+j+1<< ", " << (i+1)*Nv+j  << ", -1,\n";
01297     }
01298   }
01299   fout << "\t\t\t ]\n" ;
01300   fout << "\t\t\t}\n" ; // IndexedFaceSet
01301   fout << "\t\t}\n" ;
01302   fout << "\t ]\n" ; 
01303   fout << "\t}\n" ;
01304   fout << "  ]\n" ; 
01305   fout << "}\n" ; 
01306 
01307   Point_nD<T,N> p_mid((p_max.x()+p_min.x())/T(2),
01308                       (p_max.y()+p_min.y())/T(2),
01309                       (p_max.z()+p_min.z())/T(2));
01310   
01311   T x_axis = p_max.x() - p_min.x() ;
01312   T y_axis = p_max.y() - p_min.y() ; 
01313 
01314   T axis = (x_axis< y_axis) ? y_axis : x_axis ;
01315   axis *= T(2) ;
01316 
01317   fout << "Viewpoint {\n\t position " << p_mid.x() << ' ' << p_mid.y() << ' ' << p_max.z()+axis << "\n\t description \"top\"\n}\n" ; 
01318 
01319 
01320   fout << "NavigationInfo { type \"EXAMINE\" }\n" ; 
01321 
01322   return 1 ;
01323 }
01324 
01325 
01326 
01327 
01328 } // end namespace
01329 
01330 #endif // #define PLIB_NURBS_SURFACE_SOURCE