#include <iostream>
#include <cmath>
using namespace std;
struct CubicPoly
{
float c0, c1, c2, c3;
float eval(float t)
{
float t2 = t*t;
float t3 = t2 * t;
return c0 + c1*t + c2*t2 + c3*t3;
}
};
/*
* Compute coefficients for a cubic polynomial
* p(s) = c0 + c1*s + c2*s^2 + c3*s^3
* such that
* p(0) = x0, p(1) = x1
* and
* p'(0) = t0, p'(1) = t1.
*/
void InitCubicPoly(float x0, float x1, float t0, float t1, CubicPoly &p)
{
p.c0 = x0;
p.c1 = t0;
p.c2 = -3*x0 + 3*x1 - 2*t0 - t1;
p.c3 = 2*x0 - 2*x1 + t0 + t1;
}
// standard Catmull-Rom spline: interpolate between x1 and x2 with previous/following points x1/x4
// (we don't need this here, but it's for illustration)
void InitCatmullRom(float x0, float x1, float x2, float x3, CubicPoly &p)
{
// Catmull-Rom with tension 0.5
InitCubicPoly(x1, x2, 0.5f*(x2-x0), 0.5f*(x3-x1), p);
}
// compute coefficients for a nonuniform Catmull-Rom spline
void InitNonuniformCatmullRom(float x0, float x1, float x2, float x3, float dt0, float dt1, float dt2, CubicPoly &p)
{
// compute tangents when parameterized in [t1,t2]
float t1 = (x1 - x0) / dt0 - (x2 - x0) / (dt0 + dt1) + (x2 - x1) / dt1;
float t2 = (x2 - x1) / dt1 - (x3 - x1) / (dt1 + dt2) + (x3 - x2) / dt2;
// rescale tangents for parametrization in [0,1]
t1 *= dt1;
t2 *= dt1;
InitCubicPoly(x1, x2, t1, t2, p);
}
struct Vec2D
{
Vec2D(float _x, float _y) : x(_x), y(_y) {}
float x, y;
};
float VecDistSquared(const Vec2D& p, const Vec2D& q)
{
float dx = q.x - p.x;
float dy = q.y - p.y;
return dx*dx + dy*dy;
}
void InitCentripetalCR(const Vec2D& p0, const Vec2D& p1, const Vec2D& p2, const Vec2D& p3,
CubicPoly &px, CubicPoly &py)
{
float dt0 = powf(VecDistSquared(p0, p1), 0.25f);
float dt1 = powf(VecDistSquared(p1, p2), 0.25f);
float dt2 = powf(VecDistSquared(p2, p3), 0.25f);
// safety check for repeated points
if (dt1 < 1e-4f) dt1 = 1.0f;
if (dt0 < 1e-4f) dt0 = dt1;
if (dt2 < 1e-4f) dt2 = dt1;
InitNonuniformCatmullRom(p0.x, p1.x, p2.x, p3.x, dt0, dt1, dt2, px);
InitNonuniformCatmullRom(p0.y, p1.y, p2.y, p3.y, dt0, dt1, dt2, py);
}
int main()
{
Vec2D p0(0,0), p1(1,1), p2(1.1,1), p3(2,0);
CubicPoly px, py;
InitCentripetalCR(p0, p1, p2, p3, px, py);
for (int i = 0; i <= 10; ++i)
cout << px.eval(0.1f*i) << " " << py.eval(0.1f*i) << endl;
}
I2luY2x1ZGUgPGlvc3RyZWFtPgojaW5jbHVkZSA8Y21hdGg+Cgp1c2luZyBuYW1lc3BhY2Ugc3RkOwoKc3RydWN0IEN1YmljUG9seQp7CglmbG9hdCBjMCwgYzEsIGMyLCBjMzsKCQoJZmxvYXQgZXZhbChmbG9hdCB0KQoJewoJCWZsb2F0IHQyID0gdCp0OwoJCWZsb2F0IHQzID0gdDIgKiB0OwoJCXJldHVybiBjMCArIGMxKnQgKyBjMip0MiArIGMzKnQzOwoJfQp9OwoKLyoKICogQ29tcHV0ZSBjb2VmZmljaWVudHMgZm9yIGEgY3ViaWMgcG9seW5vbWlhbAogKiAgIHAocykgPSBjMCArIGMxKnMgKyBjMipzXjIgKyBjMypzXjMKICogc3VjaCB0aGF0CiAqICAgcCgwKSA9IHgwLCBwKDEpID0geDEKICogIGFuZAogKiAgIHAnKDApID0gdDAsIHAnKDEpID0gdDEuCiAqLwp2b2lkIEluaXRDdWJpY1BvbHkoZmxvYXQgeDAsIGZsb2F0IHgxLCBmbG9hdCB0MCwgZmxvYXQgdDEsIEN1YmljUG9seSAmcCkKewogICAgcC5jMCA9IHgwOwogICAgcC5jMSA9IHQwOwogICAgcC5jMiA9IC0zKngwICsgMyp4MSAtIDIqdDAgLSB0MTsKICAgIHAuYzMgPSAyKngwIC0gMip4MSArIHQwICsgdDE7Cn0KCi8vIHN0YW5kYXJkIENhdG11bGwtUm9tIHNwbGluZTogaW50ZXJwb2xhdGUgYmV0d2VlbiB4MSBhbmQgeDIgd2l0aCBwcmV2aW91cy9mb2xsb3dpbmcgcG9pbnRzIHgxL3g0Ci8vICh3ZSBkb24ndCBuZWVkIHRoaXMgaGVyZSwgYnV0IGl0J3MgZm9yIGlsbHVzdHJhdGlvbikKdm9pZCBJbml0Q2F0bXVsbFJvbShmbG9hdCB4MCwgZmxvYXQgeDEsIGZsb2F0IHgyLCBmbG9hdCB4MywgQ3ViaWNQb2x5ICZwKQp7CgkvLyBDYXRtdWxsLVJvbSB3aXRoIHRlbnNpb24gMC41CiAgICBJbml0Q3ViaWNQb2x5KHgxLCB4MiwgMC41ZiooeDIteDApLCAwLjVmKih4My14MSksIHApOwp9CgovLyBjb21wdXRlIGNvZWZmaWNpZW50cyBmb3IgYSBub251bmlmb3JtIENhdG11bGwtUm9tIHNwbGluZQp2b2lkIEluaXROb251bmlmb3JtQ2F0bXVsbFJvbShmbG9hdCB4MCwgZmxvYXQgeDEsIGZsb2F0IHgyLCBmbG9hdCB4MywgZmxvYXQgZHQwLCBmbG9hdCBkdDEsIGZsb2F0IGR0MiwgQ3ViaWNQb2x5ICZwKQp7CiAgICAvLyBjb21wdXRlIHRhbmdlbnRzIHdoZW4gcGFyYW1ldGVyaXplZCBpbiBbdDEsdDJdCiAgICBmbG9hdCB0MSA9ICh4MSAtIHgwKSAvIGR0MCAtICh4MiAtIHgwKSAvIChkdDAgKyBkdDEpICsgKHgyIC0geDEpIC8gZHQxOwogICAgZmxvYXQgdDIgPSAoeDIgLSB4MSkgLyBkdDEgLSAoeDMgLSB4MSkgLyAoZHQxICsgZHQyKSArICh4MyAtIHgyKSAvIGR0MjsKCiAgICAvLyByZXNjYWxlIHRhbmdlbnRzIGZvciBwYXJhbWV0cml6YXRpb24gaW4gWzAsMV0KICAgIHQxICo9IGR0MTsKICAgIHQyICo9IGR0MTsKCiAgICBJbml0Q3ViaWNQb2x5KHgxLCB4MiwgdDEsIHQyLCBwKTsKfQoKc3RydWN0IFZlYzJECnsKCVZlYzJEKGZsb2F0IF94LCBmbG9hdCBfeSkgOiB4KF94KSwgeShfeSkge30KCWZsb2F0IHgsIHk7Cn07CgpmbG9hdCBWZWNEaXN0U3F1YXJlZChjb25zdCBWZWMyRCYgcCwgY29uc3QgVmVjMkQmIHEpCnsKCWZsb2F0IGR4ID0gcS54IC0gcC54OwoJZmxvYXQgZHkgPSBxLnkgLSBwLnk7CglyZXR1cm4gZHgqZHggKyBkeSpkeTsKfQoKdm9pZCBJbml0Q2VudHJpcGV0YWxDUihjb25zdCBWZWMyRCYgcDAsIGNvbnN0IFZlYzJEJiBwMSwgY29uc3QgVmVjMkQmIHAyLCBjb25zdCBWZWMyRCYgcDMsCglDdWJpY1BvbHkgJnB4LCBDdWJpY1BvbHkgJnB5KQp7CiAgICBmbG9hdCBkdDAgPSBwb3dmKFZlY0Rpc3RTcXVhcmVkKHAwLCBwMSksIDAuMjVmKTsKICAgIGZsb2F0IGR0MSA9IHBvd2YoVmVjRGlzdFNxdWFyZWQocDEsIHAyKSwgMC4yNWYpOwogICAgZmxvYXQgZHQyID0gcG93ZihWZWNEaXN0U3F1YXJlZChwMiwgcDMpLCAwLjI1Zik7CgoJLy8gc2FmZXR5IGNoZWNrIGZvciByZXBlYXRlZCBwb2ludHMKICAgIGlmIChkdDEgPCAxZS00ZikgICAgZHQxID0gMS4wZjsKICAgIGlmIChkdDAgPCAxZS00ZikgICAgZHQwID0gZHQxOwogICAgaWYgKGR0MiA8IDFlLTRmKSAgICBkdDIgPSBkdDE7CgoJSW5pdE5vbnVuaWZvcm1DYXRtdWxsUm9tKHAwLngsIHAxLngsIHAyLngsIHAzLngsIGR0MCwgZHQxLCBkdDIsIHB4KTsKCUluaXROb251bmlmb3JtQ2F0bXVsbFJvbShwMC55LCBwMS55LCBwMi55LCBwMy55LCBkdDAsIGR0MSwgZHQyLCBweSk7Cn0KCgppbnQgbWFpbigpCnsKCVZlYzJEIHAwKDAsMCksIHAxKDEsMSksIHAyKDEuMSwxKSwgcDMoMiwwKTsKCUN1YmljUG9seSBweCwgcHk7CglJbml0Q2VudHJpcGV0YWxDUihwMCwgcDEsIHAyLCBwMywgcHgsIHB5KTsKCWZvciAoaW50IGkgPSAwOyBpIDw9IDEwOyArK2kpCgkJY291dCA8PCBweC5ldmFsKDAuMWYqaSkgPDwgIiAiIDw8IHB5LmV2YWwoMC4xZippKSA8PCBlbmRsOwp9