@micycle Thanks for sharing! I am terribly sorry for such a long silence. This task fell completely out of my focus for a long time for a serious reason. But I am getting back to it.
I’ve tried your algorithms. They work, it’s amazing! But how can I adjust the accuracy level? The result I get is quite rough. What I am trying to get must be much more precise.
I somehow ported the algorithm by Philip J. Schneider from JS and Python to Processing. It works but not perfectly. Sometimes it makes too big inaccuracies and next to them — unnecessary short straight segments. I guess it is my fault, but I can not find the cause 
I would be happy if someone can help me to find mistakes in the code. I want to publish it as a library in the end.
Here is the ported code. It uses PVector_D — it is a clone of PVector class but with doubles instead of floats. With floats it was glitching sometimes.
ArrayList <PVector_D[]> fitCurve(ArrayList <PVector_D> points, double maxError) {
//Remove duplicate points
for (int i = points.size()-1; i>0; i--) {
PVector_D p = points.get(i);
PVector_D prevP = points.get(i-1);
if (p.x == prevP.x && p.y == prevP.y) points.remove(i);
}
if (points.size() < 2) {
return null;
}
int len = points.size();
PVector_D leftTangent = createTangent(points.get(1), points.get(0));
PVector_D rightTangent = createTangent(points.get(len - 2), points.get(len - 1));
return fitCubic(points, leftTangent, rightTangent, maxError);
}
/**
* Fit a Bezier curve to a (sub)set of digitized points.
* Your code should not call this function directly. Use {@link fitCurve} instead.
*
* @param {Array<Array<Number>>} points - Array of digitized points, e.g. [[5,5],[5,50],[110,140],[210,160],[320,110]]
* @param {Array<Number>} leftTangent - Unit tangent vector at start point
* @param {Array<Number>} rightTangent - Unit tangent vector at end point
* @param {Number} error - Tolerance, squared error between points and fitted curve
* @returns {Array<Array<Array<Number>>>} Array of Bezier curves, where each element is [first-point, control-point-1, control-point-2, second-point] and points are [x, y]
*/
ArrayList <PVector_D[]> fitCubic(ArrayList <PVector_D> points, PVector_D leftTangent, PVector_D rightTangent, double error) {
int MaxIterations = 100; //Max times to try iterating (to find an acceptable curve)
PVector_D centerVector, toCenterTangent, fromCenterTangent; //Unit tangent vector(s) at splitPoint
PVector_D [] bezCurve = new PVector_D [4]; //Control points of fitted Bezier curve
double [] u = new double [points.size()]; //Parameter values for point
double [] uPrime = new double [points.size()]; //Improved parameter values
double maxError ; //Maximum fitting error
double prevErr;
int splitPointIndex; //Point to split point set at if we need more than one curve
int prevSplitIndex;
ArrayList <PVector_D[]> beziers = new ArrayList <PVector_D []> (); //Array of fitted Bezier curves if we need more than one curve
double dist;
//Use heuristic if region only has two points in it
if (points.size() == 2) {
PVector_D tempVector = PVector_D.sub(points.get(0), points.get(1));
dist = tempVector.mag() / 3.0;
bezCurve [0] = points.get(0);
bezCurve [1] = PVector_D.add(points.get(0), PVector_D.mult(leftTangent, dist));
bezCurve [2] = PVector_D.add(points.get(1), PVector_D.mult(rightTangent, dist));
bezCurve [3] = points.get(1);
beziers.add(bezCurve);
return beziers;
}
//Parameterize points, and attempt to fit curve
u = chordLengthParameterize(points);
bezCurve = generateBezier(points, u, leftTangent, rightTangent);
//println("first attempt:");
//println(bezCurve[0], bezCurve[1], bezCurve[2], bezCurve[3]);
ArrayList temp = computeMaxError(points, bezCurve, u);
maxError = (double)temp.get(0);
splitPointIndex = (int)temp.get(1);
if (maxError == 0 || maxError < error) {
beziers.add(bezCurve);
return beziers;
}
//If error not too large, try some reparameterization and iteration
if (maxError < error * error) {
uPrime = u;
prevErr = maxError;
prevSplitIndex = splitPointIndex;
for (int i = 0; i < MaxIterations; i++) {
uPrime = reparameterize(bezCurve, points, uPrime);
//println(str(uPrime));
bezCurve = generateBezier(points, uPrime, leftTangent, rightTangent);
//println("second attempt:");
//println(bezCurve[0], bezCurve[1], bezCurve[2], bezCurve[3]);
ArrayList temp2 = computeMaxError(points, bezCurve, u);
maxError = (double)temp2.get(0);
splitPointIndex = (int)temp2.get(1);
if (maxError < error) {
beziers.add(bezCurve);
return beziers;
} else if (splitPointIndex == prevSplitIndex) {
double errChange = maxError / prevErr;
if (errChange > 0.9999 && errChange < 1.0001) {
break;
}
}
prevErr = maxError;
prevSplitIndex = splitPointIndex;
}
}
//Fitting failed -- split at max error point and fit recursively
beziers = new ArrayList <PVector_D []> ();
/*
//To create a smooth transition from one curve segment to the next, we
//calculate the line between the points directly before and after the
//center, and use that as the tangent both to and from the center point.
//println(splitPointIndex);
centerVector = PVector_D.sub(points.get(splitPointIndex - 1), points.get(splitPointIndex + 1));
//However, this won't work if they're the same point, because the line we
//want to use as a tangent would be 0. Instead, we calculate the line from
//that "double-point" to the center point, and use its tangent.
if (centerVector.x == 0 && centerVector.y == 0) {
//[x,y] -> [-y,x]: http://stackoverflow.com/a/4780141/1869660
centerVector = PVector_D.sub(points.get(splitPointIndex - 1), points.get(splitPointIndex));
PVector_D temp3 = new PVector_D (-centerVector.y, centerVector.x);
centerVector.x = temp3.x;
centerVector.y = temp3.y;
}
toCenterTangent = centerVector.get();
toCenterTangent.normalize();
//To and from need to point in opposite directions:
fromCenterTangent = PVector_D.mult(toCenterTangent, -1);
*/
toCenterTangent = createTangent(points.get(splitPointIndex - 1), points.get(splitPointIndex));
fromCenterTangent = createTangent(points.get(splitPointIndex + 1), points.get(splitPointIndex));
/*
Note: An alternative to this "divide and conquer" recursion could be to always
let new curve segments start by trying to go all the way to the end,
instead of only to the end of the current subdivided polyline.
That might let many segments fit a few points more, reducing the number of total segments.
However, a few tests have shown that the segment reduction is insignificant
(240 pts, 100 err: 25 curves vs 27 curves. 140 pts, 100 err: 17 curves on both),
and the results take twice as many steps and milliseconds to finish,
without looking any better than what we already have.
*/
ArrayList <PVector_D> before = new ArrayList <PVector_D> ();
for (int i = 0; i < splitPointIndex + 1; i++) {
PVector_D p = points.get(i);
before.add(p);
}
ArrayList <PVector_D> after = new ArrayList <PVector_D> ();
for (int i = splitPointIndex; i < points.size(); i++) {
PVector_D p = points.get(i);
after.add(p);
}
beziers.addAll(fitCubic(before, leftTangent, toCenterTangent, error));
beziers.addAll(fitCubic(after, fromCenterTangent, rightTangent, error));
return beziers;
}
/**
* Use least-squares method to find Bezier control points for region.
*
* @param {Array<Array<Number>>} points - Array of digitized points
* @param {Array<Number>} parameters - Parameter values for region
* @param {Array<Number>} leftTangent - Unit tangent vector at start point
* @param {Array<Number>} rightTangent - Unit tangent vector at end point
* @returns {Array<Array<Number>>} Approximated Bezier curve: [first-point, control-point-1, control-point-2, second-point] where points are [x, y]
*/
PVector_D [] generateBezier(ArrayList <PVector_D> points, double [] parameters, PVector_D leftTangent, PVector_D rightTangent) {
PVector_D [] bezCurve = new PVector_D [4];
PVector_D firstPoint = points.get(0);
PVector_D lastPoint = points.get(points.size()-1);
bezCurve [0] = firstPoint;
bezCurve [3] = lastPoint;
//Precomputed rhs for eqn
PVector_D [][] A;
//Compute the A's
A = makeZeroVectors(parameters.length);
for (int i = 0; i < parameters.length; i++) {
double u = parameters[i];
double ux = 1 - u;
A[i][0] = PVector_D.mult(leftTangent, (3 * u * (ux * ux)));
A[i][1] = PVector_D.mult(rightTangent, (3 * ux * (u * u)));
}
//Matrices C & X
double [][] C = {{0, 0}, {0, 0}};
double [] X = {0, 0};
for (int i = 0; i < points.size(); i++) {
double u = parameters[i];
C[0][0] += PVector_D.dot(A[i][0], A[i][0]);
C[0][1] += PVector_D.dot(A[i][0], A[i][1]);
C[1][0] += PVector_D.dot(A[i][0], A[i][1]);
C[1][1] += PVector_D.dot(A[i][1], A[i][1]);
PVector_D [] tmp1 = {firstPoint, firstPoint, lastPoint, lastPoint};
//println(parameters[i]);
PVector_D tmp2 = PVector_D.sub(points.get(i), bezierQ(tmp1, u));
//println(tmp2);
X[0] += PVector_D.dot(A[i][0], tmp2);
X[1] += PVector_D.dot(A[i][1], tmp2);
}
//Determinants of matrices C and X
double det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
double det_C0_X = C[0][0] * X[1] - C[1][0] * X[0];
double det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1];
//Alpha values, left and right
double alpha_l, alpha_r;
// Result = test ? expression1 : expression2
alpha_l = det_C0_C1 == 0 ? 0 : det_X_C1 / det_C0_C1;
alpha_r = det_C0_C1 == 0 ? 0 : det_C0_X / det_C0_C1;
//If alpha negative, use the Wu/Barsky heuristic (see text).
//If alpha is 0, you get coincident control points that lead to
//divide by zero in any subsequent NewtonRaphsonRootFind() call.
PVector_D tempVector = PVector_D.sub(firstPoint, lastPoint);
double segLength = tempVector.mag();
double epsilon = 1.0e-6 * segLength;
if (alpha_l < epsilon || alpha_r < epsilon) {
//Fall back on standard (probably inaccurate) formula, and subdivide further if needed.
bezCurve[1] = PVector_D.add(firstPoint, PVector_D.mult(leftTangent, segLength / 3.0));
bezCurve[2] = PVector_D.add(lastPoint, PVector_D.mult(rightTangent, segLength / 3.0));
//println("the first pair of points: ");
//println(bezCurve[1], bezCurve[2]);
} else {
//First and last control points of the Bezier curve are
//positioned exactly at the first and last data points
//Control points 1 and 2 are positioned an alpha distance out
//on the tangent vectors, left and right, respectively
//println("parameters are: alfa_l = "+alpha_l+", alfa_l = "+alpha_r+", epsilon = "+epsilon);
//println(det_C0_X);
bezCurve[1] = PVector_D.add(firstPoint, PVector_D.mult(leftTangent, alpha_l));
bezCurve[2] = PVector_D.add(lastPoint, PVector_D.mult(rightTangent, alpha_r));
//println("the second pair of points: ");
//println(bezCurve[1], bezCurve[2]);
}
return bezCurve;
}
/**
* Given set of points and their parameterization, try to find a better parameterization.
*
* @param {Array<Array<Number>>} bezier - Current fitted curve
* @param {Array<Array<Number>>} points - Array of digitized points
* @param {Array<Number>} parameters - Current parameter values
* @returns {Array<Number>} New parameter values
*/
double [] reparameterize(PVector_D [] bezier, ArrayList <PVector_D> points, double [] parameters) {
/*
var j, len, point, results, u;
results = [];
for (j = 0, len = points.length; j < len; j++) {
point = points[j], u = parameters[j];
results.push(newtonRaphsonRootFind(bezier, point, u));
}
return results;
//*/
//println("reparametrized");
for (int i = 0; i < parameters.length; i++) {
parameters[i] = newtonRaphsonRootFind(bezier, points.get(i), parameters[i]);
}
return parameters;
}
/**
* Use Newton-Raphson iteration to find better root.
*
* @param {Array<Array<Number>>} bez - Current fitted curve
* @param {Array<Number>} point - Digitized point
* @param {Number} u - Parameter value for "P"
* @returns {Number} New u
*/
double newtonRaphsonRootFind(PVector_D [] bez, PVector_D point, double u) {
/*
Newton's root finding algorithm calculates f(x)=0 by reiterating
x_n+1 = x_n - f(x_n)/f'(x_n)
We are trying to find curve parameter u for some point p that minimizes
the distance from that point to the curve. Distance point to curve is d=q(u)-p.
At minimum distance the point is perpendicular to the curve.
We are solving
f = q(u)-p * q'(u) = 0
with
f' = q'(u) * q'(u) + q(u)-p * q''(u)
gives
u_n+1 = u_n - |q(u_n)-p * q'(u_n)| / |q'(u_n)**2 + q(u_n)-p * q''(u_n)|
*/
//println("here");
PVector_D d = PVector_D.sub(bezierQ(bez, u), point);
PVector_D qprime = bezierQPrime(bez, u);
double numerator = PVector_D.dot(d, qprime); // not sure
double denominator = summaXY(squareVector(qprime)) + 2*PVector_D.dot(d, bezierQPrimePrime(bez, u));
if (denominator == 0) {
return u;
} else {
return u - numerator / denominator;
}
}
/**
* Assign parameter values to digitized points using relative distances between points.
*
* @param {Array<Array<Number>>} points - Array of digitized points
* @returns {Array<Number>} Parameter values
*/
double [] chordLengthParameterize(ArrayList <PVector_D> points) {
double [] u = new double [points.size()];
u[0] = 1;
for (int i = 1; i < points.size(); i ++) {
PVector_D p = points.get(i);
PVector_D prevP = points.get(i-1);
PVector_D tempVector = PVector_D.sub(p, prevP);
u[i] = u[i-1] + tempVector.mag();
if (u[i] == 0) println("ALARM 1");
}
for (int i = 1; i < u.length; i ++) {
u[i] = u[i] / u[i-1];
//println(u[i-1]);
}
return u;
}
/**
* Find the maximum squared distance of digitized points to fitted curve.
*
* @param {Array<Array<Number>>} points - Array of digitized points
* @param {Array<Array<Number>>} bez - Fitted curve
* @param {Array<Number>} parameters - Parameterization of points
* @returns {Array<Number>} Maximum error (squared) and point of max error
*/
ArrayList computeMaxError(ArrayList <PVector_D> points, PVector_D [] bez, double [] parameters) {
//println(bez);
double maxDist = 0; //Maximum error
int splitPointIndex = floor(points.size() / 2); //Point of maximum error. Let it be in the middle for the beginning
double [] t_distMap = mapTtoRelativeDistances(bez, 10);
for (int i = 1; i < points.size()-1; i++) {
PVector_D point = points.get(i);
//Find 't' for a point on the bez curve that's as close to 'point' as possible:
double t = find_t(bez, parameters[i], t_distMap, 10);
//double t = parameters[i];
//println("the parameters are:");
//println("bez: " + bez[0] + bez[1] + bez[2] + bez[3] + " & t = " + t);
PVector_D somePoint = bezierQ(bez, t);
//println("the points are: " + somePoint + " & " + point);
PVector_D v = PVector_D.sub(somePoint, point); //Vector from point to curve
double dist = v.magSq(); //squared dist
//println("The actual error is: " + dist);
if (dist > maxDist) {
maxDist = dist;
//println("boom");
splitPointIndex = i;
}
}
ArrayList outputs = new ArrayList();
outputs.add(maxDist);
outputs.add(splitPointIndex);
return outputs;
}
//Sample 't's and map them to relative distances along the curve:
double [] mapTtoRelativeDistances(PVector_D [] bez, int B_parts) {
PVector_D B_t_curr;
double [] B_t_dist = new double [B_parts];
PVector_D B_t_prev = bez[0];
double sumLen = 0;
for (int i = 0; i < B_parts; i++) {
B_t_curr = bezierQ(bez, i / B_parts);
PVector_D tempVector = PVector_D.sub(B_t_curr, B_t_prev);
sumLen += tempVector.mag();
B_t_dist[i] = sumLen;
B_t_prev = B_t_curr;
}
//Normalize B_length to the same interval as the parameter distances; 0 to 1:
for (int i = 0; i < B_t_dist.length; i ++) {
double f = B_t_dist[i];
f /= sumLen;
B_t_dist[i] = f;
}
return B_t_dist;
}
double find_t(PVector_D [] bez, double param, double [] t_distMap, int B_parts) {
if (param < 0) return 0;
if (param > 1) return 1;
/*
'param' is a value between 0 and 1 telling us the relative position
of a point on the source polyline (linearly from the start (0) to the end (1)).
To see if a given curve - 'bez' - is a close approximation of the polyline,
we compare such a poly-point to the point on the curve that's the same
relative distance along the curve's length.
But finding that curve-point takes a little work:
There is a function "B(t)" to find points along a curve from the parametric parameter 't'
(also relative from 0 to 1: http://stackoverflow.com/a/32841764/1869660
http://pomax.github.io/bezierinfo/#explanation),
but 't' isn't linear by length (http://gamedev.stackexchange.com/questions/105230).
So, we sample some points along the curve using a handful of values for 't'.
Then, we calculate the length between those samples via plain euclidean distance;
B(t) concentrates the points around sharp turns, so this should give us a good-enough outline of the curve.
Thus, for a given relative distance ('param'), we can now find an upper and lower value
for the corresponding 't' by searching through those sampled distances.
Finally, we just use linear interpolation to find a better value for the exact 't'.
More info:
http://gamedev.stackexchange.com/questions/105230/points-evenly-spaced-along-a-bezier-curve
http://stackoverflow.com/questions/29438398/cheap-way-of-calculating-cubic-bezier-length
http://steve.hollasch.net/cgindex/curves/cbezarclen.html
https://github.com/retuxx/tinyspline
*/
double lenMax, lenMin, tMax, tMin;
double t = 0;
//Find the two t-s that the current param distance lies between,
//and then interpolate a somewhat accurate value for the exact t:
for (int i = 0; i < B_parts; i++) {
if (param <= t_distMap[i]) {
tMin = (i - 1) / B_parts;
tMax = i / B_parts;
lenMin = t_distMap[i - 1];
lenMax = t_distMap[i];
t = (param - lenMin) / (lenMax - lenMin) * (tMax - tMin) + tMin;
break;
}
}
return t;
}
//Creates a vector of length 1 which shows the direction from B to A
PVector_D createTangent(PVector_D pointA, PVector_D pointB) {
PVector_D tangent = PVector_D.sub(pointA, pointB);
tangent.normalize();
return tangent;
}
PVector_D [][] makeZeroVectors(int x) {
PVector_D [][] zs = new PVector_D[x][2];
for (int i = 0; i < x; i++) {
zs[i][0] = new PVector_D(0, 0);
zs[i][1] = new PVector_D(0, 0);
}
return zs;
}
PVector_D squareVector(PVector_D p) {
return new PVector_D (p.x*p.x, p.y*p.y);
}
double summaXY(PVector_D p) {
return p.x+p.y;
}
PVector_D bezierQ(PVector_D [] ctrlPoly, double t) {
double tx = 1.0 - t;
PVector_D pA = PVector_D.mult(ctrlPoly[0], tx * tx * tx);
PVector_D pB = PVector_D.mult(ctrlPoly[1], 3 * tx * tx * t);
PVector_D pC = PVector_D.mult(ctrlPoly[2], 3 * tx * t * t);
PVector_D pD = PVector_D.mult(ctrlPoly[3], t * t * t);
return PVector_D.add(PVector_D.add(pA, pB), PVector_D.add(pC, pD));
}
PVector_D bezierQPrime(PVector_D [] ctrlPoly, double t) {
double tx = 1.0 - t;
PVector_D pA = PVector_D.mult(PVector_D.sub(ctrlPoly[1], ctrlPoly[0]), 3 * tx * tx);
PVector_D pB = PVector_D.mult(PVector_D.sub(ctrlPoly[2], ctrlPoly[1]), 6 * tx * t);
PVector_D pC = PVector_D.mult(PVector_D.sub(ctrlPoly[3], ctrlPoly[2]), 3 * t * t);
return PVector_D.add(PVector_D.add(pA, pB), pC);
}
PVector_D bezierQPrimePrime(PVector_D [] ctrlPoly, double t) {
return PVector_D.add(PVector_D.mult(PVector_D.add(PVector_D.sub(ctrlPoly[2], PVector_D.mult(ctrlPoly[1], 2)), ctrlPoly[0]), 6 * (1.0 - t)), PVector_D.mult(PVector_D.add(PVector_D.sub(ctrlPoly[3], PVector_D.mult(ctrlPoly[2], 2)), ctrlPoly[1]), 6 * t));
}
