2
0
Fork 0
OB-Xd/Modules/gin/utilities/spline.cpp
2020-05-21 09:14:57 +02:00

78 lines
2.3 KiB
C++
Executable file

/*==============================================================================
Copyright 2018 by Roland Rabien, Devin Lane
For more information visit www.rabiensoftware.com
==============================================================================*/
/* "THE BEER-WARE LICENSE" (Revision 42): Devin Lane wrote this file. As long as you retain
* this notice you can do whatever you want with this stuff. If we meet some day, and you
* think this stuff is worth it, you can buy me a beer in return. */
Spline::Spline (const Array<juce::Point<double>>& points)
{
jassert (points.size() >= 3); // "Must have at least three points for interpolation"
points.size ();
int n = points.size() - 1;
Array<double> b, d, a, c, l, u, z, h;
a.resize (n);
b.resize (n);
c.resize (n + 1);
d.resize (n);
h.resize (n + 1);
l.resize (n + 1);
u.resize (n + 1);
z.resize (n + 1);
l.set (0, 1.0);
u.set (0, 0.0);
z.set (0, 0.0);
h.set (0, points[1].getX() - points[0].getX());
for (int i = 1; i < n; i++)
{
h.set (i, points[i+1].getX() - points[i].getX());
l.set (i, (2 * (points[i+1].getX() - points[i-1].getX())) - (h[i-1]) * u[i-1]);
u.set (i, l[i] == 0 ? 0 :h[i] / l[i]);
a.set (i, h[i] == 0 || h[i-1] == 0 ? 0 : (3.0 / (h[i])) * (points[i+1].getY() - points[i].getY()) - (3.0 / (h[i-1])) * (points[i].getY() - points[i-1].getY()));
z.set (i, l[i] == 0 ? 0 : (a[i] - (h[i-1]) * z[i-1]) / l[i]);
}
l.set (n, 1.0);
z.set (n, 0.0);
c.set (n, 0.0);
for (int j = n - 1; j >= 0; j--)
{
c.set (j, z[j] - u[j] * c[j+1]);
b.set (j, h[j] == 0 ? 0 : (points[j+1].getY() - points[j].getY()) / (h[j]) - ((h[j]) * (c[j+1] + 2.0 * c[j])) / 3.0);
d.set (j, h[j] == 0 ? 0 : (c[j+1] - c[j]) / (3.0 * h[j]));
}
for (int i = 0; i < n; i++)
elements.add (Element (points[i].getX(), points[i].getY(), b[i], c[i], d[i]));
}
double Spline::operator[] (double x) const
{
return interpolate (x);
}
double Spline::interpolate (double x) const
{
if (elements.size() == 0)
return 0.0;
int i;
for (i = 0; i < elements.size(); i++)
{
if (! (elements[i] < x))
break;
}
if (i != 0) i--;
return elements[i].eval (x);
}