|
MathTL
|
00001 // -*- c++ -*- 00002 00003 // +--------------------------------------------------------------------+ 00004 // | This file is part of MathTL - the Mathematical Template Library | 00005 // | | 00006 // | Copyright (c) 2002-2009 | 00007 // | Thorsten Raasch, Manuel Werner | 00008 // +--------------------------------------------------------------------+ 00009 00010 #ifndef _MATHTL_CARDINAL_SPLINES_H 00011 #define _MATHTL_CARDINAL_SPLINES_H 00012 00013 #include <cmath> 00014 #include <utils/function.h> 00015 #include <numerics/splines.h> 00016 #include <utils/tiny_tools.h> 00017 #include <algebra/piecewise.h> 00018 #include <algebra/polynomial.h> 00019 #include <algebra/sparse_matrix.h> 00020 00021 namespace MathTL 00022 { 00026 template <int d> 00027 class CardinalKnotSequence 00028 : public KnotSequence 00029 { 00030 public: 00031 /* 00032 index of the first knot 00033 */ 00034 int k0() const { return -d+1; } 00035 00039 double knot(const int k) const { return k; } 00040 }; 00041 00046 template <int d> 00047 inline 00048 double EvaluateCardinalBSpline(const double x) 00049 { 00050 return (x < 0 || x >= d ? 0 00051 : (x * EvaluateCardinalBSpline<d-1>(x) 00052 + (d-x) * EvaluateCardinalBSpline<d-1>(x-1)) / (d-1)); 00053 } 00054 00059 template <> 00060 double EvaluateCardinalBSpline<3>(const double x) 00061 { 00062 if (x < 0 || x >= 3) 00063 return 0; 00064 if (x < 1) 00065 return x*x/2; 00066 else 00067 return (x < 2 ? x*(3-x)-1.5 : (x*(x-6)+9)/2); 00068 } 00069 00074 template <> 00075 inline 00076 double EvaluateCardinalBSpline<2>(const double x) 00077 { 00078 return (x < 0 || x >= 2 ? 0 : 1-abs(x-1)); 00079 } 00080 00085 template <> 00086 inline 00087 double EvaluateCardinalBSpline<1>(const double x) 00088 { 00089 return (x < 0 || x >= 1 ? 0 : 1); 00090 } 00091 00095 template <int d> 00096 inline 00097 double EvaluateCardinalBSpline(const int k, const double x) 00098 { 00099 if (x < k || x >= k+d) 00100 return 0.; 00101 return ((x-k) * EvaluateCardinalBSpline<d-1>(k, x) 00102 + (k+d-x) * EvaluateCardinalBSpline<d-1>(k+1, x)) / (d-1); 00103 } 00104 00108 template <> 00109 inline 00110 double EvaluateCardinalBSpline<1>(const int k, const double x) 00111 { 00112 if (x < k || x >= k+1) 00113 return 0.0; 00114 return 1.0; 00115 } 00116 00121 template <int d> 00122 inline 00123 double EvaluateCardinalBSpline_td(const int j, const int k, const double x) 00124 { 00125 #if 0 00126 const double factor(1<<j); 00127 return sqrt(factor) * EvaluateCardinalBSpline<d>(k, factor * x + d/2); 00128 #else 00129 return twotothejhalf(j) * EvaluateCardinalBSpline<d>((1<<j) * x - k + d/2); 00130 #endif 00131 } 00132 00137 template <int d> 00138 inline 00139 double EvaluateCardinalBSpline_x(const double x) 00140 { 00141 return (x < 0 || x >= d ? 0 00142 : EvaluateCardinalBSpline<d-1>(x) - EvaluateCardinalBSpline<d-1>(x-1)); 00143 } 00144 00149 template <> 00150 inline 00151 double EvaluateCardinalBSpline_x<2>(const double x) 00152 { 00153 if (x < 0 || x >= 2) 00154 return 0; 00155 return (x < 1 ? 1: -1); 00156 } 00157 00161 template <> 00162 inline 00163 double EvaluateCardinalBSpline_x<1>(const double x) 00164 { 00165 return 0; 00166 } 00167 00171 template <int d> 00172 inline 00173 double EvaluateCardinalBSpline_x(const int k, const double x) 00174 { 00175 if (x < k || x >= k+d) 00176 return 0; 00177 return EvaluateCardinalBSpline<d-1>(k, x) - EvaluateCardinalBSpline<d-1>(k+1, x); 00178 } 00179 00183 template <int d> 00184 inline 00185 double EvaluateCardinalBSpline_x_sq(const int k, const double x) 00186 { 00187 return EvaluateCardinalBSpline_x<d-1>(k, x) - EvaluateCardinalBSpline_x<d-1>(k+1, x); 00188 } 00189 00193 template <> 00194 inline 00195 double EvaluateCardinalBSpline_x<1>(const int k, const double x) 00196 { 00197 return 0; 00198 } 00199 00204 template <int d> 00205 inline 00206 double EvaluateCardinalBSpline_td_x(const int j, const int k, const double x) 00207 { 00208 #if 0 00209 const double factor(1<<j); 00210 return factor * sqrt(factor) * EvaluateCardinalBSpline_x<d>(k, factor * x + d/2); 00211 #else 00212 return twotothejhalf(3*j) * EvaluateCardinalBSpline_x<d>((1<<j) * x - k + d/2); 00213 #endif 00214 } 00215 00220 template <int d> 00221 inline 00222 double EvaluateCardinalBSpline_td_x_sq(const int j, const int k, const double x) 00223 { 00224 #if 0 00225 const double factor(1<<j); 00226 return factor * factor * sqrt(factor) * EvaluateCardinalBSpline_x_sq<d>(k, factor * x + d/2); 00227 #else 00228 return twotothejhalf(5*j) * EvaluateCardinalBSpline_x_sq<d>(k, (1<<j) * x + d/2); 00229 #endif 00230 } 00231 00232 00238 double EvaluateCardinalBSpline(const int d, const int k, double x) 00239 { 00240 if (x < k) 00241 return 0.; 00242 else 00243 { 00244 if (x >= k+d) 00245 return 0.; 00246 else 00247 { 00248 // we know that x\in\supp N_d(.-k) 00249 if (d == 1) 00250 return 1.; 00251 else 00252 { 00253 // hard-encode case d=2 00254 if (d == 2) 00255 { 00256 if (x < k+1) 00257 return x-k; 00258 else 00259 return 2.-(x-k); 00260 } 00261 else 00262 return ((x-k) * EvaluateCardinalBSpline(d-1, k, x) 00263 + (k+d-x) * EvaluateCardinalBSpline(d-1, k+1, x)) / (d-1); 00264 } 00265 } 00266 } 00267 } 00268 00269 00274 inline double EvaluateCardinalBSpline_td(const int d, const int j, const int k, const double x) 00275 { 00276 const double factor(1<<j); 00277 return sqrt(factor) * EvaluateCardinalBSpline(d, k, factor * x + d/2); 00278 } 00279 00283 inline double EvaluateCardinalBSpline_x(const int d, const int k, const double x) 00284 { 00285 if (d == 1) 00286 return 0.; 00287 else 00288 return EvaluateCardinalBSpline(d-1, k, x) - EvaluateCardinalBSpline(d-1, k+1, x); 00289 } 00290 00295 inline double EvaluateCardinalBSpline_td_x(const int d, const int j, const int k, const double x) 00296 { 00297 const double factor(1<<j); 00298 return factor * sqrt(factor) * EvaluateCardinalBSpline_x(d, k, factor * x + d/2); 00299 } 00300 00301 00302 /*------------------------------------------------------------------* 00303 * n n n - k + 1 n n * 00304 * Binomialkoeffizient ( ) = ( ) --------- , ( ) = ( ) = 1 * 00305 * k k-1 k 0 n * 00306 * * 00307 * Rekursive Berechnung, damit keine zu grossen Zahlen entstehen * 00308 *------------------------------------------------------------------*/ 00309 00310 double binomd (long int n, long int k) 00311 { 00312 long int i; 00313 long int binomialkoeffizient = 1; 00314 00315 if (n < 0 || k < 0 || k > n) return (-1); 00316 00317 if (k == 0 || k == n) return (binomialkoeffizient); 00318 00319 for (i = 1; i <= k; i++) 00320 binomialkoeffizient = binomialkoeffizient * (n - i + 1) / i; 00321 return (double)(binomialkoeffizient); 00322 } 00323 00324 00325 // helper function for cBspline 00326 SparseMatrix<double> localBSplineCoeffs(const int d) 00327 { 00328 int k, mu, nu; 00329 00330 SparseMatrix<double> A(d,d); 00331 if(d==1) 00332 { 00333 A.set_entry(0, 0, 1); 00334 return A; 00335 } 00336 else 00337 { 00338 SparseMatrix<double> B = localBSplineCoeffs(d-1); 00339 // first row (nu=0) 00340 // A.set_entry(0, 0, 0); 00341 for(k=1; k<d-1; k++) 00342 { 00343 // A.set_entry(0, k, 0); 00344 for(mu=0; mu<=d-2; mu++) 00345 A.set_entry(0, k, A.get_entry(0, k) + (pow(k, mu+1)*(B.get_entry(mu, k-1)-B.get_entry(mu, k)) + pow(-1, mu)*B.get_entry(mu, k-1))/(mu+1)); 00346 } 00347 // A.set_entry(0, d-1, 0); 00348 for(mu=0; mu<=d-2; mu++) 00349 A.set_entry(0, d-1, (A.get_entry(0, d-1) + (pow(k, mu+1)*B.get_entry(mu, d-2) + pow(-1, mu)*B.get_entry(mu, d-2))/(mu+1))); 00350 // second to last row (nu>0) 00351 for(nu=1; nu<d; nu++) 00352 { 00353 A.set_entry(nu, 0, B.get_entry(nu-1, 0)/nu); 00354 for(k=1; k<d-1; k++) 00355 { 00356 A.set_entry(nu, k, B.get_entry(nu-1, k)/nu); 00357 for(mu=nu; mu<=d-1; mu++) 00358 A.set_entry(nu, k, A.get_entry(nu, k) - (pow(-1, mu-nu)*binomd(mu, nu)*B.get_entry(mu-1, k-1))/mu); 00359 } 00360 // A.put(nu, d-1) = 0; 00361 for(mu=nu; mu<=d-1; mu++) 00362 A.set_entry(nu, d-1, A.get_entry(nu, d-1) - (pow(-1, mu-nu)*binomd(mu, nu)*B.get_entry(mu-1, d-2))/mu); 00363 } 00364 00365 return A; 00366 } 00367 } 00368 00369 00370 00375 template <int d> 00376 Piecewise<double> ExpandBspline(const int j, const int k) 00377 { 00378 Piecewise<double> r(j); 00379 Polynomial<double> p, q; 00380 p.set_coefficient(0, floor((double) d/2.0) -k); //floor((double) d/2.0) 00381 p.set_coefficient(1, ldexp(1.0, j)); // p(x)=2^jx-k+floor(d/2) 00382 double factor = sqrt(ldexp(1.0, j)); 00383 int m, n; 00384 SparseMatrix<double> coeffs(d, d); 00385 switch(d) // precalculated entries for 1<=d<=4 00386 { 00387 case 1: 00388 { 00389 coeffs.set_entry(0, 0, 1); 00390 } 00391 break; 00392 case 2: 00393 { 00394 coeffs.set_entry(1, 0, 1); 00395 00396 coeffs.set_entry(0, 1, 2); 00397 coeffs.set_entry(1, 1, -1); 00398 } 00399 break; 00400 case 3: 00401 { 00402 coeffs.set_entry(2, 0, 0.5); 00403 00404 coeffs.set_entry(0, 1, -1.5); 00405 coeffs.set_entry(1, 1, 3); 00406 coeffs.set_entry(2, 1, -1); 00407 00408 coeffs.set_entry(0, 2, 4.5); 00409 coeffs.set_entry(1, 2, -3); 00410 coeffs.set_entry(2, 2, 0.5); 00411 } 00412 break; 00413 case 4: 00414 { 00415 coeffs.set_entry(3, 0, 1.0/6.0); 00416 00417 coeffs.set_entry(0, 1, 2.0/3.0); 00418 coeffs.set_entry(1, 1, -2); 00419 coeffs.set_entry(2, 1, 2); 00420 coeffs.set_entry(3, 1, -0.5); 00421 00422 coeffs.set_entry(0, 2, -22.0/3.0); 00423 coeffs.set_entry(1, 2, 10); 00424 coeffs.set_entry(2, 2, -4); 00425 coeffs.set_entry(3, 2, 0.5); 00426 00427 coeffs.set_entry(0, 3, 32.0/3.0); 00428 coeffs.set_entry(1, 3, -8); 00429 coeffs.set_entry(2, 3, 2); 00430 coeffs.set_entry(3, 3, -1.0/6.0); 00431 } 00432 break; 00433 default: 00434 { 00435 //return r; 00436 coeffs = localBSplineCoeffs(d); 00437 } 00438 } 00439 for (n=0; n<=d-1; n++) 00440 { 00441 // coeffs.findc(n, ind); 00442 // CoeffsType col; 00443 for (m=0; m<=d-1; m++) 00444 q.set_coefficient(m, coeffs.get_entry(m, n)); 00445 00446 r.set_local_expansion(n + k - floor((double) d/2.0), q.substitute_into(p)); //(int)floor((double) d/2.0) 00447 } 00448 r *= factor; 00449 00450 return r; 00451 } 00452 00453 00454 00455 00459 template <int d> 00460 class CardinalBSpline : public Function<1> 00461 { 00462 public: 00466 CardinalBSpline() : Function<1>(1) {} 00467 00471 virtual ~CardinalBSpline() {} 00472 00476 inline double value(const Point<1>& p, 00477 const unsigned int component = 0) const 00478 { 00479 return EvaluateCardinalBSpline<d>(0, p(0)); 00480 } 00481 00485 void vector_value(const Point<1> &p, 00486 Vector<double>& values) const 00487 { 00488 values.resize(1, false); 00489 values[0] = value(p); 00490 } 00491 }; 00492 } 00493 00494 #endif
1.7.6.1
