Library of Assembled Shared Sources
spline_cubic.inl
Go to the documentation of this file.
1/** @file
2 * @author Bram de Greve (bram@cocamware.com)
3 * @author Tom De Muer (tom@cocamware.com)
4 *
5 * *** BEGIN LICENSE INFORMATION ***
6 *
7 * The contents of this file are subject to the Common Public Attribution License
8 * Version 1.0 (the "License"); you may not use this file except in compliance with
9 * the License. You may obtain a copy of the License at
10 * http://lass.sourceforge.net/cpal-license. The License is based on the
11 * Mozilla Public License Version 1.1 but Sections 14 and 15 have been added to cover
12 * use of software over a computer network and provide for limited attribution for
13 * the Original Developer. In addition, Exhibit A has been modified to be consistent
14 * with Exhibit B.
15 *
16 * Software distributed under the License is distributed on an "AS IS" basis, WITHOUT
17 * WARRANTY OF ANY KIND, either express or implied. See the License for the specific
18 * language governing rights and limitations under the License.
19 *
20 * The Original Code is LASS - Library of Assembled Shared Sources.
21 *
22 * The Initial Developer of the Original Code is Bram de Greve and Tom De Muer.
23 * The Original Developer is the Initial Developer.
24 *
25 * All portions of the code written by the Initial Developer are:
26 * Copyright (C) 2004-2011 the Initial Developer.
27 * All Rights Reserved.
28 *
29 * Contributor(s):
30 *
31 * Alternatively, the contents of this file may be used under the terms of the
32 * GNU General Public License Version 2 or later (the GPL), in which case the
33 * provisions of GPL are applicable instead of those above. If you wish to allow use
34 * of your version of this file only under the terms of the GPL and not to allow
35 * others to use your version of this file under the CPAL, indicate your decision by
36 * deleting the provisions above and replace them with the notice and other
37 * provisions required by the GPL License. If you do not delete the provisions above,
38 * a recipient may use your version of this file under either the CPAL or the GPL.
39 *
40 * *** END LICENSE INFORMATION ***
41 */
42
43
44
45#ifndef LASS_GUARDIAN_OF_INCLUSION_NUM_SPLINE_CUBIC_INL
46#define LASS_GUARDIAN_OF_INCLUSION_NUM_SPLINE_CUBIC_INL
47
48#include "num_common.h"
49#include "spline_cubic.h"
50#include "num_cast.h"
51#include "impl/matrix_solve.h"
52#include "../stde/extended_iterator.h"
53
54#include <cstddef>
55
56namespace lass
57{
58namespace num
59{
60
61// --- public --------------------------------------------------------------------------------------
62
63/** construct an empty spline
64 */
65template <typename S, typename D, typename T>
66SplineCubic<S, D, T>::SplineCubic()
67{
68}
69
70
71
72/** construct a spline from a range of control/data pairs.
73 *
74 * Each node consists of a control value and a data value. This contstructor accepts a single
75 * range [iFirst, iLast) of control/data pairs. The iterator type should have two fields
76 * @c first and @c second that contain respectively the control and data values. This is
77 * choosen so that a std::pair can be used as a representatin of the control/data pair.
78 *
79 * @pre
80 * @arg [iFirst, iLast) is a valid range.
81 * @arg @c PairInputIterator has a member @c first containing the control value
82 * @arg @c PairInputIterator has a member @c second containing the data value
83 *
84 * @par complexity:
85 * O(D * log(N)) with
86 * @arg D = a figure that indicates the complexity of operations on data values.
87 * Is most probably linear with the dimension of the data value
88 * @arg N = number of nodes
89 */
90template <typename S, typename D, typename T>
91template <typename PairInputIterator>
92SplineCubic<S, D, T>::SplineCubic(PairInputIterator iFirst, PairInputIterator iLast)
93{
94 while (iFirst != iLast)
95 {
96 Node node;
97 node.x = iFirst->first;
98 node.d = iFirst->second;
99 nodes_.push_back(node);
100 ++iFirst;
101 }
102 init();
103}
104
105
106
107/** construct a spline from seperate ranges.
108 *
109 * Each node consists of a control value and a data value. This contstructor accepts seperate
110 * ranges for control and data values. The control values are given by the range
111 * [iFirstControl , iLastControl). Of the range of the data values, only the iterator
112 * iFirstData to the the first element is given.
113 *
114 * @pre
115 * @arg [iFirstControl, iLastControl) is a valid range.
116 * @arg [iFirstData, iFirstData + (iLastControl - iFirstControl)) is a valid range.
117 *
118 * @par complexity:
119 * O(D * log(N)) with
120 * @arg D = a figure that indicates the complexity of operations on data values.
121 * Is most probably linear with the dimension of the data value
122 * @arg N = number of nodes
123 */
124template <typename S, typename D, typename T>
125template <typename ScalarInputIterator, typename DataInputIterator>
126SplineCubic<S, D, T>::SplineCubic(ScalarInputIterator iFirstControl,
127 ScalarInputIterator iLastControl,
128 DataInputIterator iFirstData)
129{
130 while (iFirstControl != iLastControl)
131 {
132 Node node;
133 node.x = *iFirstControl++;
134 node.d = *iFirstData++;
135 nodes_.push_back(node);
136 }
137 init();
138}
139
140
141
142/** Get the linear interpolated data value that corresponds with constrol value @a iX.
143 *
144 * @pre this->isEmpty() == false
145 *
146 * @par Complexity:
147 * O(D * log(N)) with
148 * @arg D = a figure that indicates the complexity of operations on data values.
149 * Is most probably linear with the dimension of the data value
150 * @arg N = number of nodes
151 */
152template <typename S, typename D, typename T>
153const typename SplineCubic<S, D, T>::TData
154SplineCubic<S, D, T>::operator ()(TScalar iX) const
155{
156 LASS_ASSERT(!isEmpty());
157
158 const TNodeConstIterator n = findNode(iX);
159 const TScalar s = iX - n->x;
160 TData result(n->d);
161 TDataTraits::multiplyAccumulate(result, n->c, s);
162 TDataTraits::multiplyAccumulate(result, n->b, num::sqr(s));
163 TDataTraits::multiplyAccumulate(result, n->a, num::cubic(s));
164 return result;
165}
166
167
168
169/** Get the first derivative of data value that corresponds with constrol value @a iX.
170 *
171 * @pre this->isEmpty() == false
172 *
173 * @par Complexity:
174 * O(D * log(N)) with
175 * @arg D = a figure that indicates the complexity of operations on data values.
176 * Is most probably linear with the dimension of the data value
177 * @arg N = number of nodes
178 */
179template <typename S, typename D, typename T>
180const typename SplineCubic<S, D, T>::TData
181SplineCubic<S, D, T>::derivative(TScalar iX) const
182{
183 LASS_ASSERT(!isEmpty());
184
185 const TNodeConstIterator n = findNode(iX);
186 const TScalar s = iX - n->x;
187 TData result(n->c);
188 TDataTraits::multiplyAccumulate(result, n->b, 2 * s);
189 TDataTraits::multiplyAccumulate(result, n->a, 3 * num::sqr(s));
190 return result;
191}
192
193
194
195/** Get the second derivative of data value that corresponds with constrol value @a iX.
196 *
197 * @pre this->isEmpty() == false
198 *
199 * @par Complexity:
200 * O(D * log(N)) with
201 * @arg D = a figure that indicates the complexity of operations on data values.
202 * Is most probably linear with the dimension of the data value
203 * @arg N = number of nodes
204 */
205template <typename S, typename D, typename T>
206const typename SplineCubic<S, D, T>::TData
207SplineCubic<S, D, T>::derivative2(TScalar iX) const
208{
209 LASS_ASSERT(!isEmpty());
210
211 const TNodeConstIterator n = findNode(iX);
212 const TScalar s = iX - n->x;
213 TData result(n->b);
214 TDataTraits::multiplyAccumulate(result, n->a, 6 * s);
215 return result;
216}
217
218
219
220/** Get the integrated data value between control points @a iA and @a iB.
221 *
222 * @pre this->isEmpty() == false
223 *
224 * @par Complexity:
225 * O(D * M * log(N)) with
226 * @arg D = a figure that indicates the complexity of operations on data values.
227 * Is most probably linear with the dimension of the data value
228 * @arg M = number of nodes between @a iA and @a iB.
229 * @arg N = total number of nodes in spline
230 */
231template <typename S, typename D, typename T>
232const typename SplineCubic<S, D, T>::TData
233SplineCubic<S, D, T>::integral(TScalar iBegin, TScalar iEnd) const
234{
235 LASS_ASSERT(!isEmpty());
236
237 TNodeConstIterator first = findNode(iBegin);
238 TNodeConstIterator last = findNode(iEnd);
239 if (first == last)
240 {
241 const TScalar s1 = iBegin - first->x;
242 const TScalar s2 = iEnd - first->x;
243 TData result(first->d);
244 TDataTraits::scale(result, s2 - s1);
245 TDataTraits::multiplyAccumulate(result, first->c, (num::sqr(s2) - num::sqr(s1)) / 2);
246 TDataTraits::multiplyAccumulate(result, first->b, (num::cubic(s2) - num::cubic(s1)) / 3);
247 TDataTraits::multiplyAccumulate(result, first->a, (num::sqr(num::sqr(s2)) - num::sqr(num::sqr(s1))) / 4);
248 return result;
249 }
250 else
251 {
252 TScalar multiplier = 1;
253 if (iEnd < iBegin)
254 {
255 std::swap(iBegin, iEnd);
256 std::swap(first, last);
257 multiplier = -1;
258 }
259
260 TNodeConstIterator next = stde::next(first);
261
262 const TScalar s1 = iBegin - first->x;
263 const TScalar s2 = next->x - first->x;
264 TData result(first->d);
265 TDataTraits::scale(result, s2 - s1);
266 TDataTraits::multiplyAccumulate(result, first->c, (num::sqr(s2) - num::sqr(s1)) / 2);
267 TDataTraits::multiplyAccumulate(result, first->b, (num::cubic(s2) - num::cubic(s1)) / 3);
268 TDataTraits::multiplyAccumulate(result, first->a, (num::sqr(num::sqr(s2)) - num::sqr(num::sqr(s1))) / 4);
269 first = next++;
270
271 while (first != last)
272 {
273 const TScalar s = next->x - first->x;
274 TDataTraits::multiplyAccumulate(result, first->d, s);
275 TDataTraits::multiplyAccumulate(result, first->c, num::sqr(s) / 2);
276 TDataTraits::multiplyAccumulate(result, first->b, num::cubic(s) / 3);
277 TDataTraits::multiplyAccumulate(result, first->a, num::sqr(num::sqr(s)) / 4);
278 first = next++;
279 }
280
281 const TScalar s = iEnd - first->x;
282 TDataTraits::multiplyAccumulate(result, first->d, s);
283 TDataTraits::multiplyAccumulate(result, first->c, num::sqr(s) / 2);
284 TDataTraits::multiplyAccumulate(result, first->b, num::cubic(s) / 3);
285 TDataTraits::multiplyAccumulate(result, first->a, num::sqr(num::sqr(s)) / 4);
286
287 TDataTraits::scale(result, multiplier);
288 return result;
289 }
290}
291
292
293
294/** return true if the spline contains any nodes at all.
295 *
296 * @par Complexity:
297 * O(1)
298 */
299template <typename S, typename D, typename T>
300bool SplineCubic<S, D, T>::isEmpty() const
301{
302 return nodes_.empty();
303}
304
305
306
307/** return the range of control values for which the spline can interpolate.
308 * @par complexity:
309 * O(1)
310 */
311template <typename S, typename D, typename T>
312const typename SplineCubic<S, D, T>::TControlRange
313SplineCubic<S, D, T>::controlRange() const
314{
315 return TControlRange(nodes_.front().x, nodes_.back().x);
316}
317
318
319
320// --- private -------------------------------------------------------------------------------------
321
322template <typename S, typename D, typename T>
323void SplineCubic<S, D, T>::init()
324{
325 typedef std::vector<TScalar> TVector;
326
327 // are there any nodes at all? we need at least two!
328 //
329 const size_t n = nodes_.size();
330 if (n < 2)
331 {
332 LASS_THROW("A cubic spline needs at least two nodes! This one only has '"
333 << static_cast<unsigned>(n) << "'.");
334 }
335
336 // - check for matching data dimensions
337 // - check for increasing control components
338 // - precalculate h_i = x_(i+1) - x_i
339 //
340 dataDimension_ = TDataTraits::dimension(nodes_[0].d);
341 TVector h;
342 for (size_t i = 1; i < n; ++i)
343 {
344 h.push_back(nodes_[i].x - nodes_[i - 1].x);
345
346 if (h.back() <= 0)
347 {
348 LASS_THROW("Nodes in cubic spline must have absolutely increasing control components.");
349 }
350 if (TDataTraits::dimension(nodes_[i].d) != dataDimension_)
351 {
352 LASS_THROW("All data elements in cubic spline must have same dimension.");
353 }
354 }
355
356 if (n == 2)
357 {
358 // linear interpolation between node 0 and 1.
359 TDataTraits::zero(nodes_[0].a, dataDimension_);
360 TDataTraits::zero(nodes_[0].b, dataDimension_);
361 nodes_[0].c = nodes_[1].d;
362 TDataTraits::multiplyAccumulate(nodes_[0].c, nodes_[0].d, -1);
363 TDataTraits::scale(nodes_[0].c, num::inv(h[0]));
364 nodes_[1].a = nodes_[0].a;
365 nodes_[1].b = nodes_[0].b;
366 nodes_[1].c = nodes_[0].c;
367 return;
368 }
369
370 // --- init some elements ---
371
372 for (size_t i = 0; i < n; ++i)
373 {
374 TDataTraits::zero(nodes_[i].b, dataDimension_);
375 }
376
377 // --- precalculate coefficients. ---
378
379 // find coefficients for tridiagonal matrix
380 //
381 const size_t numUnknowns = n - 2;
382
383 //*
384 TVector ma(numUnknowns);
385 TVector mb(numUnknowns);
386 TVector mc(numUnknowns);
387 TVector temp(numUnknowns);
388 TVector unknowns(numUnknowns);
389
390 mb[0] = 2 * (h[0] + h[1]);
391 mc[0] = h[1];
392 for (size_t i = 1; i < numUnknowns - 1; ++i)
393 {
394 ma[i] = h[i];
395 mb[i] = 2 * (h[i] + h[i + 1]);
396 mc[i] = h[i + 1];
397 }
398 ma[numUnknowns - 1] = h[numUnknowns - 1];
399 mb[numUnknowns - 1] = 2 * (h[numUnknowns - 1] + h[numUnknowns]);
400
401 for (size_t k = 0; k < dataDimension_; ++k)
402 {
403 for (size_t i = 0; i < numUnknowns; ++i)
404 {
405 // d_i / h_i - d_(i+1) * (h_i + h_(i+1)) / (h_i * h_(i+1)) + d_(i+2) / h_(i+1)
406 const TScalar d0 = TDataTraits::get(nodes_[i].d, k);
407 const TScalar d1 = TDataTraits::get(nodes_[i + 1].d, k);
408 const TScalar d2 = TDataTraits::get(nodes_[i + 2].d, k);
409 unknowns[i] = 3 * (d0 / h[i] - d1 * (h[i] + h[i + 1]) / (h[i] * h[i + 1]) + d2 / h[i + 1]);
410 }
411 if (!num::impl::solveTridiagonal<TScalar>(\
412 ma.begin(), mb.begin(), mc.begin(), unknowns.begin(), temp.begin(), numCast<std::ptrdiff_t>(numUnknowns)))
413 {
414 LASS_THROW("serious logic error, could not solve equation, contact [Bramz]");
415 }
416 for (size_t i = 0; i < numUnknowns; ++i)
417 {
418 TDataTraits::set(nodes_[i + 1].b, k, unknowns[i]);
419 }
420 }
421
422 // find parameters for splines
423 //
424 for (size_t i = 0; i < n - 1; ++i)
425 {
426 Node& node = nodes_[i];
427 const Node& nextNode = nodes_[i + 1];
428
429 // a_i = (b_(i+1) - b_i) / (3h_i)
430 node.a = nextNode.b;
431 TDataTraits::multiplyAccumulate(node.a, node.b, -1);
432 TDataTraits::scale(node.a, num::inv(3 * h[i]));
433
434 // c_i = (y_(i+1) - y_i) / h - (M_(i+1) + 2M_i) * h/6
435 node.c = nextNode.d;
436 TDataTraits::multiplyAccumulate(node.c, node.d, -1);
437 TDataTraits::scale(node.c, num::inv(h[i]));
438 TDataTraits::multiplyAccumulate(node.c, nextNode.b, -h[i] / 3);
439 TDataTraits::multiplyAccumulate(node.c, node.b, -2 * h[i] / 3);
440 }
441
442 Node& node = nodes_[n - 1];
443 const Node& prevNode = nodes_[n - 2];
444 node.c = prevNode.c;
445 TDataTraits::multiplyAccumulate(node.c, prevNode.b, 2 * h[n - 2]);
446 TDataTraits::multiplyAccumulate(node.c, prevNode.a, 3 * num::sqr(h[n - 2]));
447 node.a = prevNode.a;
448
449}
450
451
452/** binary search to find node that belongs to iX
453 *
454 * @return
455 * @arg the index @a i if @a iX is in the interval [@c nodes_[i].x, @c nodes_[i+1].x)
456 * @arg 0 if @a iX is smaller than @c nodes_[0].x</tt>
457 * @arg @c nodes_.size()-2 if @a iX is greater than @c nodes_[nodes_.size()-1].x
458 *
459 * @par complexity:
460 * O(log N)
461 */
462template <typename S, typename D, typename T>
463const typename SplineCubic<S, D, T>::TNodeConstIterator
464SplineCubic<S, D, T>::findNode(TScalar iX) const
465{
466 const typename TNodes::size_type size = nodes_.size();
467 LASS_ASSERT(nodes_.size() >= 2);
468
469 // early outs
470 //
471 if (iX < nodes_[1].x)
472 {
473 return nodes_.begin();
474 }
475 if (iX >= nodes_[size - 2].x)
476 {
477 return stde::prev(nodes_.end(), 2);
478 }
479
480 // binary search
481 //
482 TNodeConstIterator first = nodes_.begin();
483 TNodeConstIterator last = nodes_.end();
484 while (stde::next(first) != last)
485 {
486 TNodeConstIterator middle = stde::next(first, std::distance(first, last) / 2);
487 LASS_ASSERT(middle != first && middle != last);
488
489 if (middle->x <= iX)
490 {
491 first = middle;
492 }
493 else
494 {
495 last = middle;
496 }
497 LASS_ASSERT(first != last);
498 }
499
500 LASS_ASSERT(first->x <= iX && last->x > iX);
501 return first;
502}
503
504
505
506}
507
508}
509
510#endif
511
512// EOF
extended_iterator.h
lass::num::sqr
T sqr(const T &x)
return x * x
Definition basic_ops.h:162
lass::num::inv
T inv(const T &x)
return x ^ -1
Definition basic_ops.h:178
lass::num::cubic
T cubic(const T &x)
return x * x * x
Definition basic_ops.h:170
lass::python::methods::next
const lass::python::impl::IterNextSlot next("__next__", Py_tp_iternext)
__next__ method (iterator next)
matrix_solve.h
lass::num
numeric types and traits.
Definition basic_ops.h:70
lass
Library for Assembled Shared Sources.
Definition config.h:53
num_cast.h
num_common.h
spline_cubic.h