13
0
livetrax/libs/evoral/test/CurveTest.cpp
Colin Fletcher cb3961d953 Add a test for the constrained cubic interpolation of Evoral::Curve
Add a test, based on the worked example in www.korf.co.uk/spline.pdf, for
the constrained cubic spline interpolation.

The delta values for the float comparisons are rather arbitrary, I'm sorry
to say: they're basically chosen so that everything passes.
2015-02-13 12:25:51 +00:00

410 lines
14 KiB
C++
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#include "CurveTest.hpp"
#include "evoral/ControlList.hpp"
#include "evoral/Curve.hpp"
#include <stdlib.h>
CPPUNIT_TEST_SUITE_REGISTRATION (CurveTest);
using namespace Evoral;
// linear y = Y0 + YS * x ; with x = i * (X1 - X0) + X0; and i = [0..1023]
#define VEC1024LINCMP(X0, X1, Y0, YS) \
cl->curve ().get_vector ((X0), (X1), vec, 1024); \
for (int i = 0; i < 1024; ++i) { \
char msg[64]; \
snprintf (msg, 64, "at i=%d (x0=%.1f, x1=%.1f, y0=%.1f, ys=%.3f)", \
i, X0, X1, Y0, YS); \
CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE ( \
msg, \
(Y0) + i * (YS), vec[i], \
1e-24 \
); \
}
void
CurveTest::twoPointLinear ()
{
float vec[1024];
boost::shared_ptr<Evoral::ControlList> cl = TestCtrlList();
cl->create_curve ();
cl->set_interpolation (ControlList::Linear);
// add two points to curve
cl->fast_simple_add ( 0.0 , 2048.0);
cl->fast_simple_add (8192.0 , 4096.0);
cl->curve ().get_vector (1024.0, 2047.0, vec, 1024);
VEC1024LINCMP (1024.0, 2047.0, 2304.f, .25f);
VEC1024LINCMP (2048.0, 2559.5, 2560.f, .125f);
VEC1024LINCMP ( 0.0, 4092.0, 2048.f, 1.f);
// greetings to tartina
cl->curve ().get_vector (2048.0, 2048.0, vec, 1);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 2048..2048", 2560.f, vec[0]);
/* XXX WHAT DO WE EXPECT WITH veclen=1 AND x1 > x0 ? */
#if 0
/* .. interpolated value at (x1+x0)/2 */
cl->curve ().get_vector (2048.0, 2049.0, vec, 1);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 2048-2049", 2560.125f, vec[0]);
cl->curve ().get_vector (2048.0, 2056.0, vec, 1);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 2048-2049", 2561.f, vec[0]);
#else
/* .. value at x0 */
cl->curve ().get_vector (2048.0, 2049.0, vec, 1);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 , 2048..2049", 2560.f, vec[0]);
cl->curve ().get_vector (2048.0, 2056.0, vec, 1);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 , 2048..2049", 2560.f, vec[0]);
#endif
cl->curve ().get_vector (2048.0, 2048.0, vec, 2);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=2 , 2048..2048 @ 0", 2560.f, vec[0]);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=2 , 2048..2048 @ 1", 2560.f, vec[1]);
cl->curve ().get_vector (2048.0, 2056.0, vec, 2);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=2 , 2048..2056 @ 0", 2560.f, vec[0]);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=2 , 2048..2056 @ 0", 2562.f, vec[1]);
cl->curve ().get_vector (2048.0, 2056.0, vec, 3);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 , 2048..2056 @ 0", 2560.f, vec[0]);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 , 2048..2056 @ 1", 2561.f, vec[1]);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 , 2048..2056 @ 2", 2562.f, vec[2]);
/* check out-of range..
* we expect the first and last value - no interpolation
*/
cl->curve ().get_vector (-1, -1, vec, 1);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ -1", 2048.f, vec[0]);
cl->curve ().get_vector (9999.0, 9999.0, vec, 1);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 9999", 4096.f, vec[0]);
cl->curve ().get_vector (-999.0, 0, vec, 13);
for (int i = 0; i < 13; ++i) {
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=13 @ -999..0", 2048.f, vec[i]);
}
cl->curve ().get_vector (9998.0, 9999.0, vec, 8);
for (int i = 0; i < 8; ++i) {
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=8 @ 9998..9999", 4096.f, vec[i]);
}
}
void
CurveTest::threePointLinear ()
{
float vec[4];
boost::shared_ptr<Evoral::ControlList> cl = TestCtrlList();
cl->create_curve ();
cl->set_interpolation (ControlList::Linear);
// add 3 points to curve
cl->fast_simple_add ( 0.0 , 2.0);
cl->fast_simple_add ( 100.0 , 4.0);
cl->fast_simple_add ( 200.0 , 0.0);
cl->curve ().get_vector (50.0, 60.0, vec, 1);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 50", 3.f, vec[0]);
cl->curve ().get_vector (100.0, 100.0, vec, 1);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 100", 4.f, vec[0]);
cl->curve ().get_vector (150.0, 150.0, vec, 1);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 150", 2.f, vec[0]);
cl->curve ().get_vector (130.0, 150.0, vec, 3);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 130..150 @ 0", 2.8f, vec[0]);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 130..150 @ 2", 2.4f, vec[1]);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 130..150 @ 3", 2.0f, vec[2]);
cl->curve ().get_vector (80.0, 160.0, vec, 3);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 80..160 @ 0", 3.6f, vec[0]);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 80..160 @ 2", 3.2f, vec[1]);
CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 80..160 @ 3", 1.6f, vec[2]);
}
void
CurveTest::threePointDiscete ()
{
boost::shared_ptr<Evoral::ControlList> cl = TestCtrlList();
cl->set_interpolation (ControlList::Discrete);
// add 3 points to curve
cl->fast_simple_add ( 0.0 , 2.0);
cl->fast_simple_add ( 100.0 , 4.0);
cl->fast_simple_add ( 200.0 , 0.0);
CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
cl->set_interpolation (ControlList::Linear);
CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
CPPUNIT_ASSERT_EQUAL(3.2, cl->unlocked_eval(120.));
CPPUNIT_ASSERT_EQUAL(1.6, cl->unlocked_eval(160.));
}
void
CurveTest::ctrlListEval ()
{
boost::shared_ptr<Evoral::ControlList> cl = TestCtrlList();
cl->fast_simple_add ( 0.0 , 2.0);
cl->set_interpolation (ControlList::Discrete);
CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(120.));
CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(160.));
cl->set_interpolation (ControlList::Linear);
CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(120.));
CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(160.));
cl->fast_simple_add ( 100.0 , 4.0);
cl->set_interpolation (ControlList::Discrete);
CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
cl->set_interpolation (ControlList::Linear);
CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
cl->fast_simple_add ( 200.0 , 0.0);
cl->set_interpolation (ControlList::Discrete);
CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
cl->set_interpolation (ControlList::Linear);
CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
CPPUNIT_ASSERT_EQUAL(3.2, cl->unlocked_eval(120.));
CPPUNIT_ASSERT_EQUAL(1.6, cl->unlocked_eval(160.));
cl->fast_simple_add ( 300.0 , 8.0);
cl->set_interpolation (ControlList::Discrete);
CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
CPPUNIT_ASSERT_EQUAL(0.0, cl->unlocked_eval(250.));
CPPUNIT_ASSERT_EQUAL(8.0, cl->unlocked_eval(999.));
cl->set_interpolation (ControlList::Linear);
CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
CPPUNIT_ASSERT_EQUAL(3.2, cl->unlocked_eval(120.));
CPPUNIT_ASSERT_EQUAL(1.6, cl->unlocked_eval(160.));
CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(250.));
CPPUNIT_ASSERT_EQUAL(8.0, cl->unlocked_eval(999.));
cl->fast_simple_add ( 400.0 , 9.0);
cl->set_interpolation (ControlList::Discrete);
CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
CPPUNIT_ASSERT_EQUAL(0.0, cl->unlocked_eval(250.));
CPPUNIT_ASSERT_EQUAL(8.0, cl->unlocked_eval(350.));
CPPUNIT_ASSERT_EQUAL(9.0, cl->unlocked_eval(999.));
cl->set_interpolation (ControlList::Linear);
CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
CPPUNIT_ASSERT_EQUAL(3.2, cl->unlocked_eval(120.));
CPPUNIT_ASSERT_EQUAL(1.6, cl->unlocked_eval(160.));
CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(250.));
CPPUNIT_ASSERT_EQUAL(8.5, cl->unlocked_eval(350.));
CPPUNIT_ASSERT_EQUAL(9.0, cl->unlocked_eval(999.));
}
void
CurveTest::constrainedCubic ()
{
struct point {
int x, y;
};
static const struct point data[] = {
/* values from worked example in www.korf.co.uk/spline.pdf */
{ 0, 30 },
{ 10, 130 },
{ 30, 150 },
{ 50, 150 },
{ 70, 170 },
{ 90, 220 },
{ 100, 320 },
};
int32_t type = 0;
Evoral::Parameter p(type);
Evoral::ParameterDescriptor pd;
Evoral::ControlList l(p,pd);
size_t i;
l.set_interpolation(Evoral::ControlList::Curved);
for (i=0; i<sizeof(data)/sizeof(data[0]); i++) {
l.add (data[i].x, data[i].y);
}
Evoral::Curve curve(l);
float f[121];
curve.get_vector(-10, 110, f, 121);
const float *g = &f[10]; /* so g starts at x==0 */
/* given points - should be exactly equal */
CPPUNIT_ASSERT_EQUAL( 30.0f, g[-10]);
CPPUNIT_ASSERT_EQUAL( 30.0f, g[ 0]);
CPPUNIT_ASSERT_EQUAL(130.0f, g[ 10]);
CPPUNIT_ASSERT_EQUAL(150.0f, g[ 30]);
CPPUNIT_ASSERT_EQUAL(150.0f, g[ 40]);
CPPUNIT_ASSERT_EQUAL(150.0f, g[ 50]);
CPPUNIT_ASSERT_EQUAL(320.0f, g[100]);
CPPUNIT_ASSERT_EQUAL(320.0f, g[110]);
/*
First segment, i=1, for 0 <= x <= 10
f'1(x1) = 2/((x2 x1)/(y2 y1) + (x1 x0)/(y1 y0))
= 2/((30 10)/(150 130) + (10 0)/(130 30))
= 1.8181
f'1(x0) = 3/2*(y1 y0)/(x1 x0) - f'1(x1)/2
= 3/2*(130 30)/(10 0) 1.818/2
= 14.0909
f"1(x0) = -2*(f'1(x1) + 2* f'1(x0))/(x1 x0) + 6*(y1 y0)/ (x1 x0)^2
= -2*(1.8181 + 2*14.0909)/(10 0) + 6*(130 30)/(10 0)^2
= 0
f"1(x1) = 2*(2*f'1(x1) + f'1(x0))/(x1 x0) - 6*(y1 y0)/ (x1 x0)^2
= 2*(2*1.818 + 14.0909)/(10 0) 6*(130 30)/(10 0)^2
= -2.4545
d1 = 1/6 * (f"1(x1) - f"1(x0))/(x1 x0)
= 1/6 * (-2.4545 0)/(10 0)
= -0.0409
c1 = 1/2 * (x1*f"1(x0) x0*f"1(x1))/(x1 x0)
= 1/2 * (10*0 0*1.8181)/(10 0)
= 0
b1 = ((y1 y0) c1*(x21 x20) d1*( x31 x30))/(x1 x0)
= ((130 30) 0*(102 02) + 0.0409*(103 03))/(10 0)
= 14.09
a1 = y0 b1*x0 c1*x20 d1*x30
= 30
y1 = 30 + 14.09x - 0.0409x3 for 0 <= x <= 10
*/
/*
Second segment, i=2, for 10 <= x <= 30
f'2(x2) = 2/((x3 x2)/(y3 y2) + (x2 x1)/(y2 y1))
= 2/((50 30)/(150 150) + (30 10)/(150 130))
= 0
f'2(x1) = 2/((x2 x1)/(y2 y1) + (x1 x0)/(y1 y0))
= 1.8181
f"2(x1) = -2*(f'2(x2) + 2* f'2(x1))/(x2 x1) + 6*(y2 y1)/ (x2 x1)^2
= -2*(0 + 2*1.8181)/(30 10) + 6*(150 130)/(30 10)2
= -0.063636
f"2(x2) = 2*(2*f'2(x2) + f'2(x1))/(x2 x1) - 6*(y2 y1)/ (x2 x1)^2
= 2*(2*0 + 1.8181)/(30 10) 6*(150 130)/(30 10)^2
= -0.11818
d2 = 1/6 * (f"2(x2) - f"2(x1))/(x2 x1)
= 1/6 * (-0.11818 + 0.063636)/(30 10)
= -0.0004545
c2 = 1/2 * (x2*f"2(x1) x1*f"2(x2))/(x2 x1)
= 1/2 * (-30*0.063636 + 10*0.11818)/(30 10)
= -0.01818
b2 = ((y2 y1) c2*(x2^2 x1^2) d2*( x2^3 x1^3))/(x2 x1)
= ((150 130) + 0.01818*(302 102) + 0.0004545*(303 103))/(30 10)
= 2.31818
a2 = y1 b2*x1 c2*x1^2 d2*x1^3
= 130 2.31818*10 + 0.01818*102 + 0.0004545*103
= 109.09
y2 = 109.09 + 2.31818x - 0.01818x^2 - 0.0004545x^3 for 10 <= x <= 30
*/
int x;
long double a1, b1, c1, d1, a2, b2, c2, d2, fdx0, fddx0, fdx1, fdx2, fddx1, fddx2;
double x0 = data[0].x;
double y0 = data[0].y;
double x1 = data[1].x;
double y1 = data[1].y;
double x2 = data[2].x;
double y2 = data[2].y;
double x3 = data[3].x;
double y3 = data[3].y;
double dx0 = x1 - x0;
double dy0 = y1 - y0;
double dx1 = x2 - x1;
double dy1 = y2 - y1;
double dx2 = x3 - x2;
double dy2 = y3 - y2;
// First (leftmost) segment
fdx1 = 2.0 / ( dx1 / dy1 + dx0 / dy0 );
fdx0 = 3.0 / 2.0 * dy0 / dx0 - fdx1 / 2.0;
fddx0 = -2.0 * (fdx1 + 2.0 * fdx0) / dx0 + 6.0 * dy0 / (dx0*dx0);
fddx1 = 2.0 * (2.0 * fdx1 + fdx0) / dx0 - 6.0 * dy0 / (dx0*dx0);
d1 = 1.0 / 6.0 * (fddx1 - fddx0) / dx0;
c1 = 1.0 / 2.0 * (x1 * fddx0 - x0 * fddx1) / dx0;
b1 = (dy0 - c1 * (x1* x1 - x0*x0) - d1 * (x1*x1*x1 - x0*x0*x0)) / dx0;
a1 = y0 - b1*x0 - c1*x0*x0 - d1*x0*x0*x0;
// printf("dx0=%f, dy0=%f, dx1=%f, dy1=%f\n", dx0, dy0, dx1, dy1);
// printf("fdx0=%Lf, fdx1=%Lf, fddx0=%Lf, fddx1=%Lf\n", fdx0, fdx1, fddx0, fddx1);
// printf("a1=%Lf, b1=%Lf, c1=%Lf, d1=%Lf\n", a1, b1, c1, d1);
// values from worked example: deltas rather arbitrary, I'm afraid
CPPUNIT_ASSERT_DOUBLES_EQUAL(30.0, a1, 0);
CPPUNIT_ASSERT_DOUBLES_EQUAL(14.09, b1, 0.001);
CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, c1, 0);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.0409, d1, 0.00001);
for (x = 0; x <= 10; x++) {
double v = a1 + b1*x + c1*x*x + d1*x*x*x;
char msg[64];
snprintf(msg, 64, "interpolating %d: v=%f, x=%f...\n", x, v, g[x]);
CPPUNIT_ASSERT_DOUBLES_EQUAL(v, g[x], 0.000004);
}
// Second segment
fdx2 = 2.0 / ( dx2 / dy2 + dx1 / dy1 );
fddx1 = -2.0 * (fdx2 + 2.0 * fdx1) / dx1 + 6.0 * dy1 / (dx1*dx1);
fddx2 = 2.0 * (2.0 * fdx2 + fdx1) / dx1 - 6.0 * dy1 / (dx1*dx1);
d2 = 1.0 / 6.0 * (fddx2 - fddx1) / dx1;
c2 = 1.0 / 2.0 * (x2 * fddx1 - x1 * fddx2) / dx1;
b2 = (dy1 - c2 * (x2*x2 - x1*x1) - d2 * (x2*x2*x2 - x1*x1*x1)) / dx1;
a2 = y1 - b2*x1 - c2*x1*x1 - d2*x1*x1*x1;
// printf("dx0=%f, dy0=%f, dx1=%f, dy1=%f dx2=%f, dy2=%f\n", dx0, dy0, dx1, dy1, dx2, dy2);
// printf("fdx1=%Lf, fdx2=%Lf, fddx1=%Lf, fddx2=%Lf\n", fdx1, fdx2, fddx1, fddx2);
// printf("a2=%Lf, b2=%Lf, c2=%Lf, d2=%Lf\n", a2, b2, c2, d2);
// values from worked example: deltas rather arbitrary, I'm afraid
CPPUNIT_ASSERT_DOUBLES_EQUAL(109.09, a2, 0.001);
CPPUNIT_ASSERT_DOUBLES_EQUAL(2.31818, b2, 0.00001);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.01818, c2, 0.00001);
CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.0004545, d2, 0.0000001);
for (x = 10; x <= 30; x++) {
double v = a2 + b2*x + c2*x*x + d2*x*x*x;
char msg[64];
snprintf(msg, 64, "interpolating %d: v=%f, x=%f...\n", x, v, g[x]);
CPPUNIT_ASSERT_DOUBLES_EQUAL(v, g[x], 0.000008);
}
}