00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef PLIB_NURBS_TRI_SPLINE_SOURCE
00027 #define PLIB_NURBS_TRI_SPLINE_SOURCE
00028
00029 #include "tri_spline.h"
00030
00033 namespace PLib {
00034
00035 int triangularNumber(int deg){
00036 return ((deg+1)*(deg+2))/2 ;
00037 }
00038
00039 static int row_start[5][5] = {
00040 {0,-1,-1,-1,-1},
00041 {0,2,-1,-1,-1},
00042 {0,3,5,-1,-1},
00043 {0,4,7,9,-1},
00044 {0,5,9,12,14}} ;
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081 inline int baryCoord(const int Deg, int i, int j, int k){
00082 return row_start[Deg][j] + i ;
00083 }
00084
00085 inline void reverseBaryCoord(int deg, int b, int &i, int &j ,int &k)
00086 {
00087 j = deg ;
00088 while(row_start[deg][j]>b)
00089 --j ;
00090 i = b - row_start[deg][j] ;
00091 k = deg - i - j ;
00092 }
00093
00094
00095 template < class T>
00096 inline T basis(int deg, T* B, int i, int j, int k)
00097 {
00098 if(i<0 || j<0 || k<0)
00099 return 0 ;
00100 if(deg<=0)
00101 return 1;
00102 return B[baryCoord(deg,i,j,k)] ;
00103 }
00104
00105 template <class T, int D>
00106 TriangularBSpline<T,D>::TriangularBSpline(int degree) : cp(triangularNumber(degree)), deg(degree){
00107 ;
00108 }
00109
00111 template <class T, int D>
00112 Point_nD<T,D>& TriangularBSpline<T,D>::b(int i, int j, int )
00113 {
00114
00115 return cp[baryCoord(deg,i,j,-1)];
00116 }
00117
00119 template <class T, int D>
00120 Point_nD<T,D> TriangularBSpline<T,D>::b(int i, int j, int ) const
00121 {
00122
00123 return cp[baryCoord(deg,i,j,-1)];
00124 }
00125
00127 template <class T, int D>
00128 Point_nD<T,D> TriangularBSpline<T,D>::operator()(T u, T v) const
00129 {
00130 T w = 1 - u - v;
00131 T u2 = u * u;
00132 T v2 = v * v;
00133 T w2 = w * w;
00134 switch(deg){
00135 case 2:
00136 return (w2*b(0,0,2) + (2*u*w)*b(1,0,1) + u2*b(2,0,0) +
00137 (2*v*w)*b(0,1,1) + (2*u*v)*b(1,1,0) + v2*b(0,2,0));
00138 case 4:{
00139 T u3 = u2 * u;
00140 T u4 = u3 * u;
00141 T v3 = v2 * v;
00142 T v4 = v3 * v;
00143 T w3 = w2 * w;
00144 T w4 = w3 * w;
00145 return (w4*b(0,0,4) + (4*u*w3)*b(1,0,3) + (6*u2*w2)*b(2,0,2) +
00146 (4*u3*w)*b(3,0,1) + u4*b(4,0,0) + (4*v*w3)*b(0,1,3) +
00147 (12*u*v*w2)*b(1,1,2) + (12*u2*v*w)*b(2,1,1) + (4*u3*v)*b(3,1,0) +
00148 (6*v2*w2)*b(0,2,2) + (12*u*v2*w)*b(1,2,1) + (6*u2*v2)*b(2,2,0) +
00149 (4*v3*w)*b(0,3,1) + (4*u*v3)*b(1,3,0) + v4*b(0,4,0));
00150 }
00151 default:
00152 return Point_nD<T,D>(0,0,0) ;
00153 }
00154 }
00155
00156 template <class T, int D>
00157 RTriangularBSpline<T,D>::RTriangularBSpline(int degree) : cp(triangularNumber(degree)), deg(degree){
00158 ;
00159 }
00160
00162 template <class T, int D>
00163 HPoint_nD<T,D>& RTriangularBSpline<T,D>::b(int i, int j, int )
00164 {
00165
00166 return cp[baryCoord(deg,i,j,-1)];
00167 }
00168
00170 template <class T, int D>
00171 HPoint_nD<T,D> RTriangularBSpline<T,D>::b(int i, int j, int ) const
00172 {
00173
00174 return cp[baryCoord(deg,i,j,-1)];
00175 }
00176
00177
00179 template <class T, int D>
00180 HPoint_nD<T,D> RTriangularBSpline<T,D>::operator()(T u, T v) const
00181 {
00182 T w = T(1) - u - v;
00183
00184
00185
00186
00187
00188
00189
00190 T *tmp1 = new T[triangularNumber(deg)] ;
00191 T *tmp2 = new T[triangularNumber(deg)] ;
00192
00193 T *B = tmp1 ;
00194 T *Bold = tmp2 ;
00195
00196 int switchBtmp = -1 ;
00197
00198 B[0] = Bold[0] = 1 ;
00199
00200 int i,j,k ;
00201 for(int p = 1 ; p <= deg ; ++p){
00202 if(switchBtmp>0){
00203 B = tmp2 ;
00204 Bold = tmp1 ;
00205 }
00206 else{
00207 B = tmp1 ;
00208 Bold = tmp2 ;
00209 }
00210 switchBtmp *= -1 ;
00211
00212 for(int s = 0; s<triangularNumber(p) ; ++s){
00213 reverseBaryCoord(p,s,i,j,k) ;
00214 B[baryCoord(p,i,j,k)] =
00215 u*basis(p-1,Bold,i-1,j,k) +
00216 v*basis(p-1,Bold,i,j-1,k) +
00217 w*basis(p-1,Bold,i,j,k-1) ;
00218 }
00219 }
00220
00221 HPoint_nD<T,D> result(0,0,0,0) ;
00222
00223 for(int s=0;s<triangularNumber(deg);++s){
00224 result += B[s]*cp[s] ;
00225 }
00226
00227 return result ;
00228 }
00229
00230 template <class T, int D>
00231 int RTriangularBSpline<T,D>::writeVRML(const char* filename, const Color& color, int Nu, int Nv, int Nw) const{
00232
00233 ofstream fout(filename) ;
00234
00235 if(!fout)
00236 return 0 ;
00237 return writeVRML(fout,color,Nu,Nv,Nw) ;
00238
00239 }
00240
00241 template <class T, int D>
00242 int RTriangularBSpline<T,D>::writeVRML(ostream &fout, const Color& color, int Nu, int Nv, int Nw) const{
00243
00244
00245 fout << "#VRML V1.0 ascii\n" ;
00246 fout << "# Generated by Phil's NURBS library\n" ;
00247 fout << "\nSeparator {\n" << "\tMaterialBinding { value PER_VERTEX_INDEXED }\n" ;
00248 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" ;
00249 fout << "\tCoordinate3 {\n" ;
00250 fout << "\t\tpoint [\n" ;
00251
00252 T u,v,du,dv,w,dw ;
00253
00254 const T uS = 0 ;
00255 const T uE = 1 ;
00256 const T vS = 0 ;
00257 const T vE = 1 ;
00258 const T wS = 0 ;
00259 const T wE = 1 ;
00260
00261 if(Nu<=1)
00262 Nu = 2 ;
00263 if(Nv<=1)
00264 Nv = 2 ;
00265 if(Nw<=1)
00266 Nw = 2 ;
00267
00268 u = uS ;
00269 v = vS ;
00270 w = wS ;
00271 du = (uE-uS)/(T)(Nu) ;
00272 dv = (vE-vS)/(T)(Nv) ;
00273 dw = (wE-wS)/(T)(Nw) ;
00274
00275 Point_nD<T,3> p1,p2,p3,p4 ;
00276
00277 int n =0 ;
00278
00279 for(int i=0;i<Nu;++i){
00280 u = uS + T(i)*du ;
00281 for(int j=0;j<Nv;++j){
00282 v = vS + T(j)*dv ;
00283 if(u+v>=T(1)-du)
00284 break ;
00285 p1 = project(operator()(u,v)) ;
00286 p2 = project(operator()(u,v+dv)) ;
00287 p3 = project(operator()(u+du,v+dv)) ;
00288 p4 = project(operator()(u+du,v)) ;
00289 fout << "\t\t\t" << p1.x() << ' ' << p1.y() << ' ' << p1.z() << ",\n" ;
00290 fout << "\t\t\t" << p2.x() << ' ' << p2.y() << ' ' << p2.z() << ",\n" ;
00291 fout << "\t\t\t" << p3.x() << ' ' << p3.y() << ' ' << p3.z() << ",\n" ;
00292 fout << "\t\t\t" << p4.x() << ' ' << p4.y() << ' ' << p4.z() << ",\n" ;
00293 ++n ;
00294 }
00295 }
00296
00297 for(int i=0;i<Nu;++i){
00298 u = uS + T(i)*du ;
00299 v = T(1)-u ;
00300 p1 = project(operator()(u,v)) ;
00301 p2 = project(operator()(u,v-du)) ;
00302 p3 = project(operator()(u+du,v-du)) ;
00303 fout << "\t\t\t" << p1.x() << ' ' << p1.y() << ' ' << p1.z() << ",\n" ;
00304 fout << "\t\t\t" << p3.x() << ' ' << p3.y() << ' ' << p3.z() << ",\n" ;
00305 fout << "\t\t\t" << p2.x() << ' ' << p2.y() << ' ' << p2.z() << ",\n" ;
00306
00307 }
00308
00309
00310 fout << "\t\t]\n" ;
00311 fout << "\t}\n" ;
00312
00313 fout << "\tIndexedFaceSet{\n" ;
00314 fout << "\t\tcoordIndex[\n" ;
00315
00316 for(int k=0;k<n;++k){
00317 fout << "\t\t\t" << 4*k << ", " << 4*k+1 << ", " << 4*k+2 << ", -1,\n" ;
00318 fout << "\t\t\t" << 4*k << ", " << 4*k+2 << ", " << 4*k+3 << ", -1,\n" ;
00319 }
00320
00321 for(int k=0;k<Nu;++k){
00322 fout << "\t\t\t" << 3*k+4*n<< ", " << 3*k+1+4*n << ", " << 3*k+2+4*n << ", -1,\n" ;
00323 }
00324
00325 fout << "\t\t]\n" ;
00326 fout << "\t}\n" ;
00327
00328 fout << "}\n" ;
00329
00330 return 1 ;
00331
00332
00333 }
00334
00335 template < class T, int D>
00336 void convert(const NurbsSurface<T,D>& surf, RTriangularBSpline<T,D> &t1, RTriangularBSpline<T,D> &t2) {
00337
00338 if(surf.degreeU() != surf.degreeV()){
00339 #ifdef USE_EXCEPTION
00340 throw NurbsError();
00341 #else
00342 Error error("convert");
00343 error << "The surface must have have the same degree in U and V.\n" ;
00344 error.fatal();
00345 #endif
00346 return ;
00347 }
00348
00349 const Matrix<HPoint_nD<T,D> > &p = surf.ctrlPnts() ;
00350
00351 switch(surf.degreeU()){
00352 case 1:
00353 t1.setDegree(2) ;
00354 t2.setDegree(2) ;
00355
00356
00357 t1.b(0,0,2) = p(0,0);
00358 t1.b(1,0,1) = 0.5 * (p(0,0) + p(0,1));
00359 t1.b(2,0,0) = p(0,1);
00360
00361 t1.b(0,1,1) = 0.5 * (p(0,0) + p(1,0));
00362 t1.b(1,1,0) = 0.5 * (p(0,0) + p(1,1));
00363
00364 t1.b(0,2,0) = p(1,0);
00365
00366
00367 t2.b(0,0,2) = p(1,1);
00368 t2.b(1,0,1) = 0.5 * (p(1,1) + p(0,1));
00369 t2.b(2,0,0) = p(0,1);
00370
00371 t2.b(0,1,1) = 0.5 * (p(1,1) + p(1,0));
00372 t2.b(1,1,0) = 0.5 * (p(0,0) + p(1,1));
00373
00374 t2.b(0,2,0) = p(1,0);
00375 break ;
00376
00377 case 2:
00378 t1.setDegree(4) ;
00379 t2.setDegree(4) ;
00380
00381
00382 t1.b(0,0,4) = p(0,0);
00383 t1.b(1,0,3) = 0.5 * (p(0,0) + p(0,1));
00384 t1.b(2,0,2) = (p(0,0) + 4.0 * p(0,1) + p(0,2)) / 6.0;
00385 t1.b(3,0,1) = 0.5 * (p(0,1) + p(0,2));
00386 t1.b(4,0,0) = p(0,2);
00387
00388 t1.b(0,1,3) = 0.5 * (p(0,0) + p(1,0));
00389 t1.b(1,1,2) = (p(0,0) + p(1,1)) / 3.0 + (p(0,1) + p(1,0)) / 6.0;
00390 t1.b(2,1,1) = (p(0,0) + p(1,2)) / 6.0 + (p(0,1) + p(1,1)) / 3.0;
00391 t1.b(3,1,0) = 0.5 * (p(0,1) + p(1,2));
00392
00393 t1.b(0,2,2) = (p(0,0) + 4.0 * p(1,0) + p(2,0)) / 6.0;
00394 t1.b(1,2,1) = (p(0,0) + p(2,1)) / 6.0 + (p(1,0) + p(1,1)) / 3.0;
00395 t1.b(2,2,0) = (p(0,0) + 4.0 * p(1,1) + p(2,2)) / 6.0;
00396
00397 t1.b(0,3,1) = 0.5 * (p(1,0) + p(2,0));
00398 t1.b(1,3,0) = 0.5 * (p(1,0) + p(2,1));
00399
00400 t1.b(0,4,0) = p(2,0);
00401
00402
00403 t2.b(0,0,4) = p(2,2);
00404 t2.b(1,0,3) = 0.5 * (p(2,2) + p(1,2));
00405 t2.b(2,0,2) = (p(2,2) + 4.0 * p(1,2) + p(0,2)) / 6.0;
00406 t2.b(3,0,1) = 0.5 * (p(1,2) + p(0,2));
00407 t2.b(4,0,0) = p(0,2);
00408
00409 t2.b(0,1,3) = 0.5 * (p(2,2) + p(2,1));
00410 t2.b(1,1,2) = (p(2,2) + p(1,1)) / 3.0 + (p(1,2) + p(2,1)) / 6.0;
00411 t2.b(2,1,1) = (p(2,2) + p(0,1)) / 6.0 + (p(1,2) + p(1,1)) / 3.0;
00412 t2.b(3,1,0) = 0.5 * (p(0,1) + p(1,2));
00413
00414 t2.b(0,2,2) = (p(2,2) + 4.0 * p(2,1) + p(2,0)) / 6.0;
00415 t2.b(1,2,1) = (p(2,2) + p(1,0)) / 6.0 + (p(2,1) + p(1,1)) / 3.0;
00416 t2.b(2,2,0) = (p(2,2) + 4.0 * p(1,1) + p(0,0)) / 6.0;
00417
00418 t2.b(0,3,1) = 0.5 * (p(2,1) + p(2,0));
00419 t2.b(1,3,0) = 0.5 * (p(1,0) + p(2,1));
00420
00421 t2.b(0,4,0) = p(2,0);
00422 break ;
00423
00424 default:
00425 #ifdef USE_EXCEPTION
00426 throw NurbsError();
00427 #else
00428 Error error("convert");
00429 error << "There is no routine to convert a NURBS surface of degree" << surf.degreeU() << endl ;
00430 error.fatal();
00431 #endif
00432
00433 }
00434 }
00435
00436 template <class T, int D>
00437 void RTriangularBSpline<T,D>::setDegree(int d){
00438 cp.resize(triangularNumber(d));
00439 deg = d ;
00440 }
00441
00442 }
00443
00444 #endif // #define PLIB_NURBS_TRI_SPLINE_SOURCE