1 /* This file is part of Evoral.
2 * Copyright (C) 2008 David Robillard <http://drobilla.net>
3 * Copyright (C) 2000-2008 Paul Davis
5 * Evoral is free software; you can redistribute it and/or modify it under the
6 * terms of the GNU General Public License as published by the Free Software
7 * Foundation; either version 2 of the License, or (at your option) any later
10 * Evoral is distributed in the hope that it will be useful, but WITHOUT ANY
11 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for details.
14 * You should have received a copy of the GNU General Public License along
15 * with this program; if not, write to the Free Software Foundation, Inc.,
16 * 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
27 #include <glibmm/threads.h>
29 #include "pbd/control_math.h"
31 #include "evoral/Curve.hpp"
32 #include "evoral/ControlList.hpp"
40 Curve::Curve (const ControlList& cl)
55 if ((npoints = _list.events().size()) > 2) {
57 /* Compute coefficients needed to efficiently compute a constrained spline
58 curve. See "Constrained Cubic Spline Interpolation" by CJC Kruger
59 (www.korf.co.uk/spline.pdf) for more details.
62 vector<double> x(npoints);
63 vector<double> y(npoints);
65 ControlList::EventList::const_iterator xx;
67 for (i = 0, xx = _list.events().begin(); xx != _list.events().end(); ++xx, ++i) {
68 x[i] = (double) (*xx)->when;
69 y[i] = (double) (*xx)->value;
72 double lp0, lp1, fpone;
74 lp0 = (x[1] - x[0])/(y[1] - y[0]);
75 lp1 = (x[2] - x[1])/(y[2] - y[1]);
80 fpone = 2 / (lp1 + lp0);
85 for (i = 0, xx = _list.events().begin(); xx != _list.events().end(); ++xx, ++i) {
87 double xdelta; /* gcc is wrong about possible uninitialized use */
88 double xdelta2; /* ditto */
89 double ydelta; /* ditto */
94 xdelta = x[i] - x[i-1];
95 xdelta2 = xdelta * xdelta;
96 ydelta = y[i] - y[i-1];
99 /* compute (constrained) first derivatives */
105 fplast = ((3 * (y[1] - y[0]) / (2 * (x[1] - x[0]))) - (fpone * 0.5));
107 /* we don't store coefficients for i = 0 */
111 } else if (i == npoints - 1) {
115 fpi = ((3 * ydelta) / (2 * xdelta)) - (fplast * 0.5);
119 /* all other segments */
121 double slope_before = ((x[i+1] - x[i]) / (y[i+1] - y[i]));
122 double slope_after = (xdelta / ydelta);
124 if (slope_after * slope_before < 0.0) {
125 /* slope changed sign */
128 fpi = 2 / (slope_before + slope_after);
132 /* compute second derivative for either side of control point `i' */
134 fppL = (((-2 * (fpi + (2 * fplast))) / (xdelta))) +
135 ((6 * ydelta) / xdelta2);
137 fppR = (2 * ((2 * fpi) + fplast) / xdelta) -
138 ((6 * ydelta) / xdelta2);
140 /* compute polynomial coefficients */
144 d = (fppR - fppL) / (6 * xdelta);
145 c = ((x[i] * fppL) - (x[i-1] * fppR))/(2 * xdelta);
150 xim12 = x[i-1] * x[i-1]; /* "x[i-1] squared" */
151 xim13 = xim12 * x[i-1]; /* "x[i-1] cubed" */
152 xi2 = x[i] * x[i]; /* "x[i] squared" */
153 xi3 = xi2 * x[i]; /* "x[i] cubed" */
155 b = (ydelta - (c * (xi2 - xim12)) - (d * (xi3 - xim13))) / xdelta;
159 (*xx)->create_coeffs();
160 (*xx)->coeff[0] = y[i-1] - (b * x[i-1]) - (c * xim12) - (d * xim13);
174 Curve::rt_safe_get_vector (double x0, double x1, float *vec, int32_t veclen) const
176 Glib::Threads::RWLock::ReaderLock lm(_list.lock(), Glib::Threads::TRY_LOCK);
181 _get_vector (x0, x1, vec, veclen);
187 Curve::get_vector (double x0, double x1, float *vec, int32_t veclen) const
189 Glib::Threads::RWLock::ReaderLock lm(_list.lock());
190 _get_vector (x0, x1, vec, veclen);
194 Curve::_get_vector (double x0, double x1, float *vec, int32_t veclen) const
196 double rx, lx, hx, max_x, min_x;
198 int32_t original_veclen;
205 if ((npoints = _list.events().size()) == 0) {
206 /* no events in list, so just fill the entire array with the default value */
207 for (int32_t i = 0; i < veclen; ++i) {
208 vec[i] = _list.descriptor().normal;
214 for (int32_t i = 0; i < veclen; ++i) {
215 vec[i] = _list.events().front()->value;
220 /* events is now known not to be empty */
222 max_x = _list.events().back()->when;
223 min_x = _list.events().front()->when;
226 /* totally past the end - just fill the entire array with the final value */
227 for (int32_t i = 0; i < veclen; ++i) {
228 vec[i] = _list.events().back()->value;
234 /* totally before the first event - fill the entire array with
237 for (int32_t i = 0; i < veclen; ++i) {
238 vec[i] = _list.events().front()->value;
243 original_veclen = veclen;
247 /* fill some beginning section of the array with the
248 initial (used to be default) value
251 double frac = (min_x - x0) / (x1 - x0);
252 int64_t fill_len = (int64_t) floor (veclen * frac);
254 fill_len = min (fill_len, (int64_t)veclen);
256 for (i = 0; i < fill_len; ++i) {
257 vec[i] = _list.events().front()->value;
264 if (veclen && x1 > max_x) {
266 /* fill some end section of the array with the default or final value */
268 double frac = (x1 - max_x) / (x1 - x0);
269 int64_t fill_len = (int64_t) floor (original_veclen * frac);
272 fill_len = min (fill_len, (int64_t)veclen);
273 val = _list.events().back()->value;
275 for (i = veclen - fill_len; i < veclen; ++i) {
282 lx = max (min_x, x0);
283 hx = min (max_x, x1);
287 const double lpos = _list.events().front()->when;
288 const double lval = _list.events().front()->value;
289 const double upos = _list.events().back()->when;
290 const double uval = _list.events().back()->value;
292 /* dx that we are using */
294 const double dx_num = hx - lx;
295 const double dx_den = veclen - 1;
296 const double lower = _list.descriptor().lower;
297 const double upper = _list.descriptor().upper;
299 /* gradient of the line */
300 const double m_num = uval - lval;
301 const double m_den = upos - lpos;
302 /* y intercept of the line */
303 const double c = uval - (m_num * upos / m_den);
305 switch (_list.interpolation()) {
306 case ControlList::Logarithmic:
307 for (int i = 0; i < veclen; ++i) {
308 const double fraction = (lx - lpos + i * dx_num / dx_den) / m_den;
309 vec[i] = interpolate_logarithmic (lval, uval, fraction, lower, upper);
312 case ControlList::Exponential:
313 for (int i = 0; i < veclen; ++i) {
314 const double fraction = (lx - lpos + i * dx_num / dx_den) / m_den;
315 vec[i] = interpolate_gain (lval, uval, fraction, upper);
318 case ControlList::Discrete:
319 // any discrete vector curves somewhere?
321 case ControlList::Curved:
322 // fallthrough, no 2 point spline
324 for (int i = 0; i < veclen; ++i) {
325 vec[i] = (lx * (m_num / m_den) + m_num * i * dx_num / (m_den * dx_den)) + c;
330 double fraction = (lx - lpos) / (upos - lpos);
331 switch (_list.interpolation()) {
332 case ControlList::Logarithmic:
333 vec[0] = interpolate_logarithmic (lval, uval, fraction, _list.descriptor().lower, _list.descriptor().upper);
335 case ControlList::Exponential:
336 vec[0] = interpolate_gain (lval, uval, fraction, _list.descriptor().upper);
338 case ControlList::Discrete:
339 // any discrete vector curves somewhere?
341 case ControlList::Curved:
342 // fallthrough, no 2 point spline
344 vec[0] = interpolate_linear (lval, uval, fraction);
360 dx = (hx - lx) / (veclen - 1);
363 for (i = 0; i < veclen; ++i, rx += dx) {
364 vec[i] = multipoint_eval (rx);
369 Curve::multipoint_eval (double x) const
371 pair<ControlList::EventList::const_iterator,ControlList::EventList::const_iterator> range;
373 ControlList::LookupCache& lookup_cache = _list.lookup_cache();
375 if ((lookup_cache.left < 0) ||
376 ((lookup_cache.left > x) ||
377 (lookup_cache.range.first == _list.events().end()) ||
378 ((*lookup_cache.range.second)->when < x))) {
380 ControlEvent cp (x, 0.0);
382 lookup_cache.range = equal_range (_list.events().begin(), _list.events().end(), &cp, ControlList::time_comparator);
385 range = lookup_cache.range;
389 a) x is an existing control point, so first == existing point, second == next point
393 b) x is between control points, so range is empty (first == second, points to where
398 if (range.first == range.second) {
400 /* x does not exist within the list as a control point */
402 lookup_cache.left = x;
404 if (range.first == _list.events().begin()) {
405 /* we're before the first point */
406 // return default_value;
407 return _list.events().front()->value;
410 if (range.second == _list.events().end()) {
411 /* we're after the last point */
412 return _list.events().back()->value;
415 ControlEvent* after = (*range.second);
417 ControlEvent* before = (*range.second);
419 double vdelta = after->value - before->value;
422 return before->value;
425 double tdelta = x - before->when;
426 double trange = after->when - before->when;
428 switch (_list.interpolation()) {
429 case ControlList::Discrete:
430 return before->value;
431 case ControlList::Logarithmic:
432 return interpolate_logarithmic (before->value, after->value, tdelta / trange, _list.descriptor().lower, _list.descriptor().upper);
433 case ControlList::Exponential:
434 return interpolate_gain (before->value, after->value, tdelta / trange, _list.descriptor().upper);
435 case ControlList::Curved:
437 ControlEvent* ev = after;
439 return ev->coeff[0] + (ev->coeff[1] * x) + (ev->coeff[2] * x2) + (ev->coeff[3] * x2 * x);
441 // no break, fallthru
443 return before->value + (vdelta * (tdelta / trange));
447 /* x is a control point in the data */
448 /* invalidate the cached range because its not usable */
449 lookup_cache.left = -1;
450 return (*range.first)->value;
453 } // namespace Evoral