00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "stdafx.h"
00021
00022 #include <iostream>
00023
00024 #include "gul_types.h"
00025 #include "gul_vector.h"
00026 #include "gunu_basics.h"
00027 #include "gunu_refine.h"
00028
00029 namespace gunu {
00030
00031 using gul::Ptr;
00032 using gul::point;
00033 using gul::hpoint;
00034 using gul::point2;
00035 using gul::hpoint2;
00036 using gul::point1;
00037 using gul::hpoint1;
00038
00039
00040
00041
00042
00043
00044
00045
00046 template< class T, class HP >
00047 void RefineCurve( int n, int p, const Ptr<T>& U, const Ptr<HP> P,
00048 const Ptr<T> X, int r, Ptr<T>& Ubar, Ptr<HP> Q )
00049 {
00050 int i,j,k,l,m,a,b,ind;
00051 T alfa;
00052
00053 m = n + p + 1;
00054
00055 a = FindSpan( X[0], n, p, U );
00056 b = FindSpan( X[r], n, p, U );
00057 b = b + 1;
00058
00059 for( j = 0; j <= a-p; j++ )
00060 Q[j] = P[j];
00061 for( j = b-1; j <= n; j++ )
00062 Q[j+r+1] = P[j];
00063
00064 for( j = 0; j <= a; j++ )
00065 Ubar[j] = U[j];
00066 for( j = b+p; j <= m; j++ )
00067 Ubar[j+r+1] = U[j];
00068
00069 i = b+p-1;
00070 k = b+p+r;
00071 for( j = r; j >= 0; j-- )
00072 {
00073 while( (X[j] <= U[i]) && (i > a) )
00074 {
00075 Q[k-p-1] = P[i-p-1];
00076 Ubar[k] = U[i];
00077 k--;
00078 i--;
00079 }
00080 Q[k-p-1] = Q[k-p];
00081
00082 for( l = 1; l <= p; l++ )
00083 {
00084 ind = k-p+l;
00085 alfa = Ubar[k+l] - X[j];
00086 if( fabs(alfa) == 0.0 )
00087 Q[ind-1] = Q[ind];
00088 else
00089 {
00090 alfa = alfa / ( Ubar[k+l] - U[i-p+l] );
00091 Q[ind-1] = alfa * Q[ind-1] + ((T)1 - alfa) * Q[ind];
00092 }
00093 }
00094 Ubar[k] = X[j];
00095 k--;
00096 }
00097 }
00098 template void RefineCurve(
00099 int n, int p, const Ptr<float>& U, const Ptr< hpoint<float> > P,
00100 const Ptr<float> X, int r, Ptr<float>& Ubar, Ptr< hpoint<float> > Q );
00101 template void RefineCurve(
00102 int n, int p, const Ptr<float>& U, const Ptr< point<float> > P,
00103 const Ptr<float> X, int r, Ptr<float>& Ubar, Ptr< point<float> > Q );
00104 template void RefineCurve(
00105 int n, int p, const Ptr<float>& U, const Ptr< hpoint2<float> > P,
00106 const Ptr<float> X, int r, Ptr<float>& Ubar, Ptr< hpoint2<float> > Q );
00107 template void RefineCurve(
00108 int n, int p, const Ptr<float>& U, const Ptr< point2<float> > P,
00109 const Ptr<float> X, int r, Ptr<float>& Ubar, Ptr< point2<float> > Q );
00110
00111 template void RefineCurve(
00112 int n, int p, const Ptr<double>& U, const Ptr< hpoint<double> > P,
00113 const Ptr<double> X, int r, Ptr<double>& Ubar, Ptr< hpoint<double> > Q );
00114 template void RefineCurve(
00115 int n, int p, const Ptr<double>& U, const Ptr< point<double> > P,
00116 const Ptr<double> X, int r, Ptr<double>& Ubar, Ptr< point<double> > Q );
00117 template void RefineCurve(
00118 int n, int p, const Ptr<double>& U, const Ptr< hpoint2<double> > P,
00119 const Ptr<double> X, int r, Ptr<double>& Ubar, Ptr< hpoint2<double> > Q);
00120 template void RefineCurve(
00121 int n, int p, const Ptr<double>& U, const Ptr< point2<double> > P,
00122 const Ptr<double> X, int r, Ptr<double>& Ubar, Ptr< point2<double> > Q );
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133 template< class T, class HP >
00134 void RefineSurfaceU(
00135 int nu, int pu, const Ptr<T>& U, int nv, int pv, const Ptr<T>& V,
00136 const Ptr< Ptr< HP > >& P, const Ptr<T>& X, int r,
00137 Ptr<T>& Ubar, Ptr<T>& Vbar, Ptr< Ptr< HP > >& Q )
00138 {
00139 int i,j,k,l,m,a,b,ind,row;
00140 T alfa;
00141
00142 m = nu + pu + 1;
00143
00144
00145 a = FindSpan( X[0], nu, pu, U );
00146 b = FindSpan( X[r], nu, pu, U );
00147 b = b + 1;
00148
00149
00150 for( j = 0; j <= a; j++ )
00151 Ubar[j] = U[j];
00152 for( j = b+pu; j <= m; j++ )
00153 Ubar[j+r+1] = U[j];
00154
00155
00156 for( j = 0; j <= nv + pv + 1; j++ )
00157 Vbar[j] = V[j];
00158
00159
00160 for( row = 0; row <= nv; row++ )
00161 {
00162 for( k = 0; k <= a-pu; k++ )
00163 Q[row][k] = P[row][k];
00164 for( k = b-1; k <= nu; k++ )
00165 Q[row][k+r+1] = P[row][k];
00166 }
00167
00168 i = b+pu-1;
00169 k = b+pu+r;
00170 for( j = r; j >= 0; j-- )
00171 {
00172 while( (X[j] <= U[i]) && (i > a) )
00173 {
00174
00175 Ubar[k] = U[i];
00176
00177 for( row = 0; row <= nv; row++ )
00178 Q[row][k-pu-1] = P[row][i-pu-1];
00179 k -= 1;
00180 i -= 1;
00181 }
00182 for( row = 0; row <= nv; row++ )
00183 Q[row][k-pu-1] = Q[row][k-pu];
00184
00185 for( l = 1; l <= pu; l++ )
00186 {
00187 ind = k-pu+l;
00188 alfa = Ubar[k+l] - X[j];
00189 if( fabs(alfa) == 0.0 )
00190 for( row = 0; row <= nv; row++ )
00191 Q[row][ind-1] = Q[row][ind];
00192 else
00193 {
00194 alfa = alfa / ( Ubar[k+l] - U[i-pu+l] );
00195 for( row = 0; row <= nv; row++ )
00196 {
00197 Q[row][ind-1] = alfa * Q[row][ind-1] + ((T)1 - alfa) * Q[row][ind];
00198 }
00199 }
00200 }
00201 Ubar[k] = X[j];
00202 k--;
00203 }
00204 }
00205
00206 template void RefineSurfaceU(
00207 int nu, int pu, const Ptr<float>& U, int nv, int pv, const Ptr<float>& V,
00208 const Ptr< Ptr< point<float> > >& P, const Ptr<float>& X, int r,
00209 Ptr<float>& Ubar, Ptr<float>& Vbar, Ptr< Ptr< point<float> > >& Q );
00210 template void RefineSurfaceU(
00211 int nu, int pu, const Ptr<float>& U, int nv, int pv, const Ptr<float>& V,
00212 const Ptr< Ptr< hpoint<float> > >& P, const Ptr<float>& X, int r,
00213 Ptr<float>& Ubar, Ptr<float>& Vbar, Ptr< Ptr< hpoint<float> > >& Q );
00214 template void RefineSurfaceU(
00215 int nu, int pu, const Ptr<float>& U, int nv, int pv, const Ptr<float>& V,
00216 const Ptr< Ptr< point1<float> > >& P, const Ptr<float>& X, int r,
00217 Ptr<float>& Ubar, Ptr<float>& Vbar, Ptr< Ptr< point1<float> > >& Q );
00218 template void RefineSurfaceU(
00219 int nu, int pu, const Ptr<float>& U, int nv, int pv, const Ptr<float>& V,
00220 const Ptr< Ptr< hpoint1<float> > >& P, const Ptr<float>& X, int r,
00221 Ptr<float>& Ubar, Ptr<float>& Vbar, Ptr< Ptr< hpoint1<float> > >& Q );
00222
00223 template void RefineSurfaceU(
00224 int nu, int pu, const Ptr<double>& U, int nv, int pv, const Ptr<double>& V,
00225 const Ptr< Ptr< point<double> > >& P, const Ptr<double>& X, int r,
00226 Ptr<double>& Ubar, Ptr<double>& Vbar, Ptr< Ptr< point<double> > >& Q );
00227 template void RefineSurfaceU(
00228 int nu, int pu, const Ptr<double>& U, int nv, int pv, const Ptr<double>& V,
00229 const Ptr< Ptr< hpoint<double> > >& P, const Ptr<double>& X, int r,
00230 Ptr<double>& Ubar, Ptr<double>& Vbar, Ptr< Ptr< hpoint<double> > >& Q );
00231 template void RefineSurfaceU(
00232 int nu, int pu, const Ptr<double>& U, int nv, int pv, const Ptr<double>& V,
00233 const Ptr< Ptr< point1<double> > >& P, const Ptr<double>& X, int r,
00234 Ptr<double>& Ubar, Ptr<double>& Vbar, Ptr< Ptr< point1<double> > >& Q );
00235 template void RefineSurfaceU(
00236 int nu, int pu, const Ptr<double>& U, int nv, int pv, const Ptr<double>& V,
00237 const Ptr< Ptr< hpoint1<double> > >& P, const Ptr<double>& X, int r,
00238 Ptr<double>& Ubar, Ptr<double>& Vbar, Ptr< Ptr< hpoint1<double> > >& Q );
00239
00240
00241
00242
00243
00244
00245
00246 template< class T >
00247 int BezierPositions( int n, int p, const Ptr<T>& U,
00248 Ptr<int>& K, Ptr<int>& M )
00249 {
00250 int m,b,nb,i;
00251
00252
00253 K.reserve_pool(n-p+1);
00254 M.reserve_pool(n-p+1);
00255
00256 m = n+p+1;
00257 b = p+1;
00258 nb = 0;
00259
00260 while( b < m )
00261 {
00262 i = b;
00263 while( (b < m) && (U[b+1] == U[b]) )
00264 b++;
00265
00266 K[nb] = b;
00267 M[nb] = b - i + 1;
00268 nb++;
00269 b++;
00270 }
00271
00272 return nb;
00273 }
00274
00275 template int BezierPositions(
00276 int n, int p, const Ptr<float>& U,
00277 Ptr<int>& K, Ptr<int>& M );
00278 template int BezierPositions(
00279 int n, int p, const Ptr<double>& U,
00280 Ptr<int>& K, Ptr<int>& M );
00281
00282
00283
00284
00285
00286 template< class T, class HP >
00287 void BezierDecomposeCurve(
00288 int nK, const Ptr<int>& K, const Ptr<int>& M,
00289 int n, int p, const Ptr<T>& U, const Ptr<HP>& Pw,
00290 Ptr< Ptr<HP> > *retQw )
00291 {
00292 int m,a,b,nb,i,j,mult,save,r,s,k;
00293 T numer,alpha;
00294 Ptr<T> alphas;
00295 Ptr< Ptr<HP> > Qw;
00296
00297 Qw.reserve_pool(nK+1);
00298 for( nb = 0; nb <= nK; nb++ )
00299 Qw[nb].reserve_pool(p+1);
00300
00301 alphas.reserve_pool(p+1);
00302
00303 m = n + p + 1;
00304 a = p;
00305 b = p+1;
00306 nb = 0;
00307
00308
00309 for( i = 0; i <= p; i++ )
00310 Qw[nb][i] = Pw[i];
00311
00312 while( nb < nK )
00313 {
00314 b = K[nb];
00315 mult = M[nb];
00316
00317 if( mult < p )
00318 {
00319 numer = U[b] - U[a];
00320
00321 for( j = p; j > mult; j-- )
00322 alphas[j-mult-1] = numer / (U[a+j] - U[a]);
00323
00324 r = p - mult;
00325
00326 for( j = 1; j <= r; j++ )
00327 {
00328 save = r - j;
00329
00330 s = mult + j;
00331
00332 for( k = p; k >= s; k-- )
00333 {
00334 alpha = alphas[k-s];
00335
00336 Qw[nb][k] = alpha * Qw[nb][k] + ((T)1 - alpha) * Qw[nb][k-1];
00337 }
00338 if( b < m )
00339 {
00340 Qw[nb+1][save] = Qw[nb][p];
00341 }
00342 }
00343 }
00344 nb = nb + 1;
00345
00346 if( nb < nK )
00347 {
00348 for( i = p - mult; i <= p; i++ )
00349 Qw[nb][i] = Pw[b-p+i];
00350 a = b;
00351 b = b+1;
00352 }
00353 }
00354 }
00355
00356 template void BezierDecomposeCurve(
00357 int nK, const Ptr<int>& K, const Ptr<int>& M,
00358 int n, int p, const Ptr<float>& U, const Ptr< point<float> >& Pw,
00359 Ptr< Ptr< point<float> > > *retQw );
00360 template void BezierDecomposeCurve(
00361 int nK, const Ptr<int>& K, const Ptr<int>& M,
00362 int n, int p, const Ptr<float>& U, const Ptr< hpoint<float> >& Pw,
00363 Ptr< Ptr< hpoint<float> > > *retQw );
00364 template void BezierDecomposeCurve(
00365 int nK, const Ptr<int>& K, const Ptr<int>& M,
00366 int n, int p, const Ptr<float>& U, const Ptr< point2<float> >& Pw,
00367 Ptr< Ptr< point2<float> > > *retQw );
00368 template void BezierDecomposeCurve(
00369 int nK, const Ptr<int>& K, const Ptr<int>& M,
00370 int n, int p, const Ptr<float>& U, const Ptr< hpoint2<float> >& Pw,
00371 Ptr< Ptr< hpoint2<float> > > *retQw );
00372
00373 template void BezierDecomposeCurve(
00374 int nK, const Ptr<int>& K, const Ptr<int>& M,
00375 int n, int p, const Ptr<double>& U, const Ptr< point<double> >& Pw,
00376 Ptr< Ptr< point<double> > > *retQw );
00377 template void BezierDecomposeCurve(
00378 int nK, const Ptr<int>& K, const Ptr<int>& M,
00379 int n, int p, const Ptr<double>& U, const Ptr< hpoint<double> >& Pw,
00380 Ptr< Ptr< hpoint<double> > > *retQw );
00381 template void BezierDecomposeCurve(
00382 int nK, const Ptr<int>& K, const Ptr<int>& M,
00383 int n, int p, const Ptr<double>& U, const Ptr< point2<double> >& Pw,
00384 Ptr< Ptr< point2<double> > > *retQw );
00385 template void BezierDecomposeCurve(
00386 int nK, const Ptr<int>& K, const Ptr<int>& M,
00387 int n, int p, const Ptr<double>& U, const Ptr< hpoint2<double> >& Pw,
00388 Ptr< Ptr< hpoint2<double> > > *retQw );
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414 template< class T, class HP >
00415 void BezierDecomposeSurfaceU(
00416 int nK, const Ptr<int>& K, const Ptr<int>& M,
00417 int nu, int pu, const Ptr<T>& U,
00418 int nv, int pv, const Ptr<T>& V,
00419 const Ptr< Ptr<HP> >& Pw,
00420 Ptr< Ptr< Ptr<HP> > > *retQw )
00421 {
00422 int m,a,b,nb,i,j,mult,save,r,s,k,row;
00423 T numer,alpha;
00424 Ptr<T> alphas;
00425 Ptr< Ptr< Ptr<HP> > > Qw;
00426
00427 Qw.reserve_pool(nK);
00428 for( nb = 0; nb < nK; nb++ )
00429 {
00430 Qw[nb].reserve_pool(nv+1);
00431 for( i = 0; i <= nv; i++ )
00432 Qw[nb][i].reserve_pool(pu+1);
00433 }
00434
00435 alphas.reserve_pool(pu+1);
00436
00437 m = nu + pu + 1;
00438 a = pu;
00439 b = pu+1;
00440 nb = 0;
00441
00442
00443 for( row = 0; row <= nv; row++ )
00444 for( i = 0; i <= pu; i++ )
00445 Qw[nb][row][i] = Pw[row][i];
00446
00447 while( nb < nK )
00448 {
00449 b = K[nb];
00450 mult = M[nb];
00451
00452 if( mult < pu )
00453 {
00454 numer = U[b] - U[a];
00455
00456 for( j = pu; j > mult; j-- )
00457 alphas[j-mult-1] = numer / (U[a+j] - U[a]);
00458
00459 r = pu - mult;
00460
00461 for( j = 1; j <= r; j++ )
00462 {
00463 save = r - j;
00464
00465 s = mult + j;
00466
00467 for( k = pu; k >= s; k-- )
00468 {
00469 alpha = alphas[k-s];
00470
00471 for( row = 0; row <= nv; row ++ )
00472 {
00473 Qw[nb][row][k] =
00474 alpha * Qw[nb][row][k] + ((T)1 - alpha) * Qw[nb][row][k-1];
00475 }
00476 }
00477
00478 if( b < m )
00479 {
00480 for( row = 0; row <= nv; row ++ )
00481 Qw[nb+1][row][save] = Qw[nb][row][pu];
00482 }
00483 }
00484 }
00485
00486 nb = nb + 1;
00487
00488 if( nb < nK )
00489 {
00490 for( row = 0; row <= nv; row++ )
00491 for( i = pu - mult; i <= pu; i++ )
00492 Qw[nb][row][i] = Pw[row][b-pu+i];
00493 a = b;
00494 }
00495 }
00496
00497 *retQw = Qw;
00498 }
00499
00500 template void BezierDecomposeSurfaceU(
00501 int nK, const Ptr<int>& K, const Ptr<int>& M,
00502 int nu, int pu, const Ptr<float>& U, int nv, int pv, const Ptr<float>& V,
00503 const Ptr< Ptr< point<float> > >& Pw,
00504 Ptr< Ptr< Ptr< point<float> > > > *retQw );
00505 template void BezierDecomposeSurfaceU(
00506 int nK, const Ptr<int>& K, const Ptr<int>& M,
00507 int nu, int pu, const Ptr<float>& U, int nv, int pv, const Ptr<float>& V,
00508 const Ptr< Ptr< hpoint<float> > >& Pw,
00509 Ptr< Ptr< Ptr< hpoint<float> > > > *retQw );
00510 template void BezierDecomposeSurfaceU(
00511 int nK, const Ptr<int>& K, const Ptr<int>& M,
00512 int nu, int pu, const Ptr<float>& U, int nv, int pv, const Ptr<float>& V,
00513 const Ptr< Ptr< point1<float> > >& Pw,
00514 Ptr< Ptr< Ptr< point1<float> > > > *retQw );
00515 template void BezierDecomposeSurfaceU(
00516 int nK, const Ptr<int>& K, const Ptr<int>& M,
00517 int nu, int pu, const Ptr<float>& U, int nv, int pv, const Ptr<float>& V,
00518 const Ptr< Ptr< hpoint1<float> > >& Pw,
00519 Ptr< Ptr< Ptr< hpoint1<float> > > > *retQw );
00520
00521 template void BezierDecomposeSurfaceU(
00522 int nK, const Ptr<int>& K, const Ptr<int>& M,
00523 int nu, int pu, const Ptr<double>& U, int nv, int pv, const Ptr<double>& V,
00524 const Ptr< Ptr< point<double> > >& Pw,
00525 Ptr< Ptr< Ptr< point<double> > > > *retQw );
00526 template void BezierDecomposeSurfaceU(
00527 int nK, const Ptr<int>& K, const Ptr<int>& M,
00528 int nu, int pu, const Ptr<double>& U, int nv, int pv, const Ptr<double>& V,
00529 const Ptr< Ptr< hpoint<double> > >& Pw,
00530 Ptr< Ptr< Ptr< hpoint<double> > > > *retQw );
00531 template void BezierDecomposeSurfaceU(
00532 int nK, const Ptr<int>& K, const Ptr<int>& M,
00533 int nu, int pu, const Ptr<double>& U, int nv, int pv, const Ptr<double>& V,
00534 const Ptr< Ptr< point1<double> > >& Pw,
00535 Ptr< Ptr< Ptr< point1<double> > > > *retQw );
00536 template void BezierDecomposeSurfaceU(
00537 int nK, const Ptr<int>& K, const Ptr<int>& M,
00538 int nu, int pu, const Ptr<double>& U, int nv, int pv, const Ptr<double>& V,
00539 const Ptr< Ptr< hpoint1<double> > >& Pw,
00540 Ptr< Ptr< Ptr< hpoint1<double> > > > *retQw );
00541 }