NumeRe v1.1.4
NumeRe: Framework für Numerische Rechnungen
spline.h
Go to the documentation of this file.
1/*
2 * spline.h
3 *
4 * simple cubic spline interpolation library without external
5 * dependencies
6 *
7 * ---------------------------------------------------------------------
8 * Copyright (C) 2011, 2014 Tino Kluge (ttk448 at gmail.com)
9 *
10 * This program is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU General Public License
12 * as published by the Free Software Foundation; either version 2
13 * of the License, or (at your option) any later version.
14 *
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with this program. If not, see <http://www.gnu.org/licenses/>.
22 * ---------------------------------------------------------------------
23 *
24 */
25
26
27#ifndef TK_SPLINE_H
28#define TK_SPLINE_H
29
30#include <cstdio>
31#include <cassert>
32#include <vector>
33#include <algorithm>
34
35
36// unnamed namespace only because the implementation is in this
37// header file and we don't want to export symbols to the obj files
38namespace
39{
40
41 namespace tk
42 {
43
44// band matrix solver
45 class band_matrix
46 {
47 private:
48 std::vector< std::vector<double> > m_upper; // upper band
49 std::vector< std::vector<double> > m_lower; // lower band
50 public:
51 band_matrix() {}; // constructor
52 band_matrix(int dim, int n_u, int n_l); // constructor
53 ~band_matrix() {}; // destructor
54 void resize(int dim, int n_u, int n_l); // init with dim,n_u,n_l
55 int dim() const; // matrix dimension
56 int num_upper() const
57 {
58 return m_upper.size() - 1;
59 }
60 int num_lower() const
61 {
62 return m_lower.size() - 1;
63 }
64 // access operator
65 double& operator () (int i, int j); // write
66 double operator () (int i, int j) const; // read
67 // we can store an additional diogonal (in m_lower)
68 double& saved_diag(int i);
69 double saved_diag(int i) const;
70 void lu_decompose();
71 std::vector<double> r_solve(const std::vector<double>& b) const;
72 std::vector<double> l_solve(const std::vector<double>& b) const;
73 std::vector<double> lu_solve(const std::vector<double>& b,
74 bool is_lu_decomposed = false);
75
76 };
77
78
79// spline interpolation
80 class spline
81 {
82 public:
83 enum bd_type
84 {
85 first_deriv = 1,
86 second_deriv = 2
87 };
88
89 private:
90 std::vector<double> m_x, m_y; // x,y coordinates of points
91 // interpolation parameters
92 // f(x) = a*(x-x_i)^3 + b*(x-x_i)^2 + c*(x-x_i) + y_i
93 std::vector<double> m_a, m_b, m_c; // spline coefficients
94 double m_b0, m_c0; // for left extrapol
95 bd_type m_left, m_right;
96 double m_left_value, m_right_value;
97 bool m_force_linear_extrapolation;
98
99 public:
100 // set default boundary condition to be zero curvature at both ends
101 spline(): m_left(second_deriv), m_right(second_deriv),
102 m_left_value(0.0), m_right_value(0.0),
103 m_force_linear_extrapolation(false)
104 {
105 ;
106 }
107
108 // optional, but if called it has to come be before set_points()
109 void set_boundary(bd_type left, double left_value,
110 bd_type right, double right_value,
111 bool force_linear_extrapolation = false);
112 void set_points(const std::vector<double>& x,
113 const std::vector<double>& y, bool cubic_spline = true);
114 double operator() (double x) const;
115 double deriv(int order, double x) const;
116 std::vector<double> operator[] (size_t i) const;
117 };
118
119
120
121// ---------------------------------------------------------------------
122// implementation part, which could be separated into a cpp file
123// ---------------------------------------------------------------------
124
125
126// band_matrix implementation
127// -------------------------
128
129 band_matrix::band_matrix(int dim, int n_u, int n_l)
130 {
131 resize(dim, n_u, n_l);
132 }
133 void band_matrix::resize(int dim, int n_u, int n_l)
134 {
135 assert(dim > 0);
136 assert(n_u >= 0);
137 assert(n_l >= 0);
138 m_upper.resize(n_u + 1);
139 m_lower.resize(n_l + 1);
140 for (size_t i = 0; i < m_upper.size(); i++)
141 {
142 m_upper[i].resize(dim);
143 }
144 for (size_t i = 0; i < m_lower.size(); i++)
145 {
146 m_lower[i].resize(dim);
147 }
148 }
149 int band_matrix::dim() const
150 {
151 if (m_upper.size() > 0)
152 {
153 return m_upper[0].size();
154 }
155 else
156 {
157 return 0;
158 }
159 }
160
161
162// defines the new operator (), so that we can access the elements
163// by A(i,j), index going from i=0,...,dim()-1
164 double& band_matrix::operator () (int i, int j)
165 {
166 int k = j - i; // what band is the entry
167 assert( (i >= 0) && (i < dim()) && (j >= 0) && (j < dim()) );
168 assert( (-num_lower() <= k) && (k <= num_upper()) );
169 // k=0 -> diogonal, k<0 lower left part, k>0 upper right part
170 if (k >= 0) return m_upper[k][i];
171 else return m_lower[-k][i];
172 }
173 double band_matrix::operator () (int i, int j) const
174 {
175 int k = j - i; // what band is the entry
176 assert( (i >= 0) && (i < dim()) && (j >= 0) && (j < dim()) );
177 assert( (-num_lower() <= k) && (k <= num_upper()) );
178 // k=0 -> diogonal, k<0 lower left part, k>0 upper right part
179 if (k >= 0) return m_upper[k][i];
180 else return m_lower[-k][i];
181 }
182// second diag (used in LU decomposition), saved in m_lower
183 double band_matrix::saved_diag(int i) const
184 {
185 assert( (i >= 0) && (i < dim()) );
186 return m_lower[0][i];
187 }
188 double& band_matrix::saved_diag(int i)
189 {
190 assert( (i >= 0) && (i < dim()) );
191 return m_lower[0][i];
192 }
193
194// LR-Decomposition of a band matrix
195 void band_matrix::lu_decompose()
196 {
197 int i_max, j_max;
198 int j_min;
199 double x;
200
201 // preconditioning
202 // normalize column i so that a_ii=1
203 for (int i = 0; i < this->dim(); i++)
204 {
205 assert(this->operator()(i, i) != 0.0);
206 this->saved_diag(i) = 1.0 / this->operator()(i, i);
207 j_min = std::max(0, i - this->num_lower());
208 j_max = std::min(this->dim() - 1, i + this->num_upper());
209 for (int j = j_min; j <= j_max; j++)
210 {
211 this->operator()(i, j) *= this->saved_diag(i);
212 }
213 this->operator()(i, i) = 1.0; // prevents rounding errors
214 }
215
216 // Gauss LR-Decomposition
217 for (int k = 0; k < this->dim(); k++)
218 {
219 i_max = std::min(this->dim() - 1, k + this->num_lower()); // num_lower not a mistake!
220 for (int i = k + 1; i <= i_max; i++)
221 {
222 assert(this->operator()(k, k) != 0.0);
223 x = -this->operator()(i, k) / this->operator()(k, k);
224 this->operator()(i, k) = -x; // assembly part of L
225 j_max = std::min(this->dim() - 1, k + this->num_upper());
226 for (int j = k + 1; j <= j_max; j++)
227 {
228 // assembly part of R
229 this->operator()(i, j) = this->operator()(i, j) + x * this->operator()(k, j);
230 }
231 }
232 }
233 }
234// solves Ly=b
235 std::vector<double> band_matrix::l_solve(const std::vector<double>& b) const
236 {
237 assert( this->dim() == (int)b.size() );
238 std::vector<double> x(this->dim());
239 int j_start;
240 double sum;
241 for (int i = 0; i < this->dim(); i++)
242 {
243 sum = 0;
244 j_start = std::max(0, i - this->num_lower());
245 for (int j = j_start; j < i; j++) sum += this->operator()(i, j) * x[j];
246 x[i] = (b[i] * this->saved_diag(i)) - sum;
247 }
248 return x;
249 }
250// solves Rx=y
251 std::vector<double> band_matrix::r_solve(const std::vector<double>& b) const
252 {
253 assert( this->dim() == (int)b.size() );
254 std::vector<double> x(this->dim());
255 int j_stop;
256 double sum;
257 for (int i = this->dim() - 1; i >= 0; i--)
258 {
259 sum = 0;
260 j_stop = std::min(this->dim() - 1, i + this->num_upper());
261 for (int j = i + 1; j <= j_stop; j++) sum += this->operator()(i, j) * x[j];
262 x[i] = ( b[i] - sum ) / this->operator()(i, i);
263 }
264 return x;
265 }
266
267 std::vector<double> band_matrix::lu_solve(const std::vector<double>& b,
268 bool is_lu_decomposed)
269 {
270 assert( this->dim() == (int)b.size() );
271 std::vector<double> x, y;
272 if (is_lu_decomposed == false)
273 {
274 this->lu_decompose();
275 }
276 y = this->l_solve(b);
277 x = this->r_solve(y);
278 return x;
279 }
280
281
282
283
284// spline implementation
285// -----------------------
286
287 void spline::set_boundary(spline::bd_type left, double left_value,
288 spline::bd_type right, double right_value,
289 bool force_linear_extrapolation)
290 {
291 assert(m_x.size() == 0); // set_points() must not have happened yet
292 m_left = left;
293 m_right = right;
294 m_left_value = left_value;
295 m_right_value = right_value;
296 m_force_linear_extrapolation = force_linear_extrapolation;
297 }
298
299
300 void spline::set_points(const std::vector<double>& x,
301 const std::vector<double>& y, bool cubic_spline)
302 {
303 assert(x.size() == y.size());
304 assert(x.size() > 2);
305 m_x = x;
306 m_y = y;
307 int n = x.size();
308
309 for (int i = 0; i < n - 1; i++)
310 {
311 assert(m_x[i] <= m_x[i + 1]);
312 }
313
314 double dAvg = 0.0;
315 size_t nPos = (size_t)-1;
316
317 // Remove duplicates by averaging
318 for (size_t i = 0; i < m_x.size()-1; i++)
319 {
320 if (m_x[i] == m_x[i+1])
321 {
322 if (nPos == (size_t)-1)
323 nPos = i;
324
325 dAvg += m_y[i];
326 }
327 else if (nPos != (size_t)-1)
328 {
329 dAvg += m_y[i];
330 dAvg /= i - nPos + 1;
331 m_y[i] = dAvg;
332 m_x.erase(m_x.begin()+nPos, m_x.begin()+i);
333 m_y.erase(m_y.begin()+nPos, m_y.begin()+i);
334 i = nPos;
335 nPos = (size_t)-1;
336 dAvg = 0.0;
337 }
338 }
339
340 if (cubic_spline == true) // cubic spline interpolation
341 {
342 // setting up the matrix and right hand side of the equation system
343 // for the parameters b[]
344 band_matrix A(n, 1, 1);
345 std::vector<double> rhs(n);
346 for (int i = 1; i < n - 1; i++)
347 {
348 A(i, i - 1) = 1.0 / 3.0 * (x[i] - x[i - 1]);
349 A(i, i) = 2.0 / 3.0 * (x[i + 1] - x[i - 1]);
350 A(i, i + 1) = 1.0 / 3.0 * (x[i + 1] - x[i]);
351 rhs[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
352 }
353 // boundary conditions
354 if (m_left == spline::second_deriv)
355 {
356 // 2*b[0] = f''
357 A(0, 0) = 2.0;
358 A(0, 1) = 0.0;
359 rhs[0] = m_left_value;
360 }
361 else if (m_left == spline::first_deriv)
362 {
363 // c[0] = f', needs to be re-expressed in terms of b:
364 // (2b[0]+b[1])(x[1]-x[0]) = 3 ((y[1]-y[0])/(x[1]-x[0]) - f')
365 A(0, 0) = 2.0 * (x[1] - x[0]);
366 A(0, 1) = 1.0 * (x[1] - x[0]);
367 rhs[0] = 3.0 * ((y[1] - y[0]) / (x[1] - x[0]) - m_left_value);
368 }
369 else
370 {
371 assert(false);
372 }
373 if (m_right == spline::second_deriv)
374 {
375 // 2*b[n-1] = f''
376 A(n - 1, n - 1) = 2.0;
377 A(n - 1, n - 2) = 0.0;
378 rhs[n - 1] = m_right_value;
379 }
380 else if (m_right == spline::first_deriv)
381 {
382 // c[n-1] = f', needs to be re-expressed in terms of b:
383 // (b[n-2]+2b[n-1])(x[n-1]-x[n-2])
384 // = 3 (f' - (y[n-1]-y[n-2])/(x[n-1]-x[n-2]))
385 A(n - 1, n - 1) = 2.0 * (x[n - 1] - x[n - 2]);
386 A(n - 1, n - 2) = 1.0 * (x[n - 1] - x[n - 2]);
387 rhs[n - 1] = 3.0 * (m_right_value - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
388 }
389 else
390 {
391 assert(false);
392 }
393
394 // solve the equation system to obtain the parameters b[]
395 m_b = A.lu_solve(rhs);
396
397 // calculate parameters a[] and c[] based on b[]
398 m_a.resize(n);
399 m_c.resize(n);
400 for (int i = 0; i < n - 1; i++)
401 {
402 m_a[i] = 1.0 / 3.0 * (m_b[i + 1] - m_b[i]) / (x[i + 1] - x[i]);
403 m_c[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
404 - 1.0 / 3.0 * (2.0 * m_b[i] + m_b[i + 1]) * (x[i + 1] - x[i]);
405 }
406 }
407 else // linear interpolation
408 {
409 m_a.resize(n);
410 m_b.resize(n);
411 m_c.resize(n);
412 for (int i = 0; i < n - 1; i++)
413 {
414 m_a[i] = 0.0;
415 m_b[i] = 0.0;
416 m_c[i] = (m_y[i + 1] - m_y[i]) / (m_x[i + 1] - m_x[i]);
417 }
418 }
419
420 // for left extrapolation coefficients
421 m_b0 = (m_force_linear_extrapolation == false) ? m_b[0] : 0.0;
422 m_c0 = m_c[0];
423
424 // for the right extrapolation coefficients
425 // f_{n-1}(x) = b*(x-x_{n-1})^2 + c*(x-x_{n-1}) + y_{n-1}
426 double h = x[n - 1] - x[n - 2];
427 // m_b[n-1] is determined by the boundary condition
428 m_a[n - 1] = 0.0;
429 m_c[n - 1] = 3.0 * m_a[n - 2] * h * h + 2.0 * m_b[n - 2] * h + m_c[n - 2]; // = f'_{n-2}(x_{n-1})
430 if (m_force_linear_extrapolation == true)
431 m_b[n - 1] = 0.0;
432 }
433
434 double spline::operator() (double x) const
435 {
436 size_t n = m_x.size();
437 // find the closest point m_x[idx] < x, idx=0 even if x<m_x[0]
438 std::vector<double>::const_iterator it;
439 it = std::lower_bound(m_x.begin(), m_x.end(), x);
440 int idx = std::max( int(it - m_x.begin()) - 1, 0);
441
442 double h = x - m_x[idx];
443 double interpol;
444 if (x < m_x[0])
445 {
446 // extrapolation to the left
447 interpol = (m_b0 * h + m_c0) * h + m_y[0];
448 }
449 else if (x > m_x[n - 1])
450 {
451 // extrapolation to the right
452 interpol = (m_b[n - 1] * h + m_c[n - 1]) * h + m_y[n - 1];
453 }
454 else
455 {
456 // interpolation
457 interpol = ((m_a[idx] * h + m_b[idx]) * h + m_c[idx]) * h + m_y[idx];
458 }
459 return interpol;
460 }
461
462 double spline::deriv(int order, double x) const
463 {
464 assert(order > 0);
465
466 size_t n = m_x.size();
467 // find the closest point m_x[idx] < x, idx=0 even if x<m_x[0]
468 std::vector<double>::const_iterator it;
469 it = std::lower_bound(m_x.begin(), m_x.end(), x);
470 int idx = std::max( int(it - m_x.begin()) - 1, 0);
471
472 double h = x - m_x[idx];
473 double interpol;
474 if (x < m_x[0])
475 {
476 // extrapolation to the left
477 switch (order)
478 {
479 case 1:
480 interpol = 2.0 * m_b0 * h + m_c0;
481 break;
482 case 2:
483 interpol = 2.0 * m_b0 * h;
484 break;
485 default:
486 interpol = 0.0;
487 break;
488 }
489 }
490 else if (x > m_x[n - 1])
491 {
492 // extrapolation to the right
493 switch (order)
494 {
495 case 1:
496 interpol = 2.0 * m_b[n - 1] * h + m_c[n - 1];
497 break;
498 case 2:
499 interpol = 2.0 * m_b[n - 1];
500 break;
501 default:
502 interpol = 0.0;
503 break;
504 }
505 }
506 else
507 {
508 // interpolation
509 switch (order)
510 {
511 case 1:
512 interpol = (3.0 * m_a[idx] * h + 2.0 * m_b[idx]) * h + m_c[idx];
513 break;
514 case 2:
515 interpol = 6.0 * m_a[idx] * h + 2.0 * m_b[idx];
516 break;
517 case 3:
518 interpol = 6.0 * m_a[idx];
519 break;
520 default:
521 interpol = 0.0;
522 break;
523 }
524 }
525 return interpol;
526 }
527
528 std::vector<double> spline::operator[] (size_t i) const
529 {
530 std::vector<double> vCoeffs;
531 if (i > m_a.size())
532 {
533 vCoeffs.push_back(NAN);
534 vCoeffs.push_back(NAN);
535 vCoeffs.push_back(NAN);
536 vCoeffs.push_back(NAN);
537 }
538 else // y0 + c*x + b*x^2 + a*x^3
539 {
540 vCoeffs.push_back(m_y[i]);
541 vCoeffs.push_back(m_c[i]);
542 vCoeffs.push_back(m_b[i]);
543 vCoeffs.push_back(m_a[i]);
544 }
545 return vCoeffs;
546 }
547
548
549 } // namespace tk
550
551
552} // namespace
553
554#endif /* TK_SPLINE_H */
tk
Definition: spline.h:42
min
#define min(a, b)
Definition: resampler.cpp:34
max
#define max(a, b)
Definition: resampler.cpp:30