I'm developing a script to connect point with curved lines.
The points are in ascendent order in x-axis like this:
I'm studying Bezier Curves, but I don't think it's the best solution (see https://www.geogebra.org/m/qcnExXbn):
But, instead of straight lines, I just would like to connect these points in a smooth way.
Could anyone help me with a formula or algorithm?
Cubic Bezier spline is perfectly suitable to smoothly connect the points.
Given an ordered sequence of $n$ points $p_0,\dots p_{n-1}$, we need to define $n$ cubic Bezier segments. The points on the $k$-th segment are defined parametrically as
\begin{align} s_k(t)=&a_k(1-t)^3+3b_k(1-t)^2t+3c_k(1-t)t^2+d_k t^3 ,\quad t\in[0,1] \tag{1}\label{1} . \end{align}
The first and second derivatives of \eqref{1} are given by \begin{align} s_k'(t)=&3((b_k-a_k)(1-t)^2+2(c_k-b_k)(1-t)t+(d_k-c_k)t^2) \tag{2}\label{2} ,\\ s_k''(t)=&6((a_k-2b_k+c_k)(1-t)+(b_k-2c_k+d_k)t) \tag{3}\label{3} . \end{align}
So, we need to define $a_k,b_k,c_k,d_k$ that satisfy (all indices here are taken $\mod\,n$): \begin{align} s_{k}'(0)=s_{k-1}'(1)&\Rightarrow &\quad b_k-a_k =&d_{k-1}-c_{k-1} \\ s_{k+1}'(0)=s_{k}'(1)&\Rightarrow &\quad b_{k+1}-a_{k+1}=& d_k-c_k \\ s_{k}''(0)=s_{k-1}''(1)&\Rightarrow &\quad c_k-2b_k+a_k=& d_{k-1}-2c_{k-1}+b_{k-1} \end{align} Excluding $c_{k-1}$ and using $d_k=a_{k+1}$, we arrive at $n\times n$ linear system for $b_k$: \begin{align} b_{k-1}+4b_k+b_{k+1}=4a_k+2a_{k+1} \tag{4}\label{4} \\ \text{for }k=0,\dots,n-1 . \end{align}
Then $c_k$ can be easily found:
\begin{align} c_k&=2a_{k+1}-b_{k+1} . \end{align}
Example:
In case if the curve is not closed, you can try this answer.
I've actually done a lot of curve drawing "by hand" - all the curves in the figures in Complex Made Simple were drawn using PostScript's curveto function, which is Bezier curves. If you try polynomial interpolation I predict you won't like the results - the curve passing though $p_1,\dots p_n$ has funny wiggles near $p_1$ due to the exact location of the other points.
I don't know what you're talking about when you say Bezier curves don't touch the points. A bezier curve is specified by four points; it passes through two of those points exactly, and the tangent vectors at those two endpoints are determined by the two "control points". The interface in terms of "control points" makes no sense to me - I wrote code to convert endpoints and tangent vectors at endpoints to endpoints and control points.
You don't say in what sense you want to "draw" this curve. If you want to write PostScript, here's what to do.
Say you want a curve $c:[0,1]\to\Bbb R^2$ with $c(0)=(x,y)$, $c(1)=(xx,yy)$, $c'(0)=(dx,dy)$ and $c'(1)=(dxx,dyy)$. The following two lines of PostScript give you exactly that:
x y moveto
x+dx/3 y+dy/3 xx-dxx/3 y-dyy/3 xx yy curveto
Not that text literally - you want strings giving the numeric value of expressions. For exammple if $(x,y)=(0,0)$, $(xx,yy)=(1,1)$, $(dx,dy)=(1,1)$ and $(dxx,dyy)=(2,3)$ you'd say
0 0 moveto
0.3333 0.3333 0.3333 0 1 1 curveto
I actually did that to get a curve passing through those points. Then I found this thing doesn't recognize eps as an image format. Here's a screen shot. (This was tedious enough - your points $A,B,\dots$ are just marked with X's.):
A path built from a series of cubic Bézier curves is fine for this application, you just need to add control points between your main points.
A single cubic Bézier curve is defined by its 2 endpoints and 2 control points. The control points lie on the tangents which pass through the endpoints; the distance between an endpoint and its associated control point affects the curvature.
To make a path passing through multiple points, we use a series of Bézier curves, connected at the endpoints, with control points chosen so that the tangents match where the curves connect (otherwise, we get cusps at the connections).
For the example points in the question, we need 5 curves, one from A to B, one from B to C, etc. Since we just have a set of points, with no tangent (or curvature) data, we construct a "default" Bézier path, using each main point's immediate neighbours to determine the tangents. Like this:
The tangent at B is parallel to AC, the tangent at C is parallel to BD, etc. The point A only has 1 neighbour, so we use AB for its tangent. Similarly, the tangent at F is EF.
The distance from a main point to its associated control point is one third of the (straight line) distance between that curve's endpoints. Eg, the distance from B to its left control point is AB/3, and the distance from B to its right control point is BC/3.
See How to find a bezier curve between two points and two tangents without anchor points for an explanation of the "1/3 rule".
Here's the Sage / Python script I used to create that diagram. SageMath is a vast mathematics system built on top of Python. It has nice plotting capabilities (courtesy of matplotlib), and its built-in vector type is handy for this kind of task. If you can read plain Python, (hopefully) you shouldn't have problems reading my script.
from itertools import chain
def bezier_fd(pos):
"Bezier path plot, using finite differences"
# Make the normalized tangent vectors
tgt = [v - u for u, v in zip(pos, pos[2:])]
if pos[0] == pos[-1]:
end = [pos[1] - pos[-2]]
tgt = end + tgt + end
else:
tgt = [pos[1] - pos[0]] + tgt + [pos[-1] - pos[-2]]
tgt = [u.normalized() for u in tgt]
# Build Bezier curves
pos_tgt = zip(pos, tgt)
p0, t0 = next(pos_tgt)
row = [p0]
curves = []
for p, t in pos_tgt:
s = abs(p - p0) / 3
row.extend([p0 + t0*s, p - t*s, p])
curves.append(row)
row = []
p0, t0 = p, t
return curves
pos = [
(0.5, 0.5),
(1, -0.5),
(1.5, 1),
(2.25, 1.1),
(2.6, -0.5),
(3, 0.5),
]
# Plot the main points
pos = [vector(p) for p in pos]
P = points(pos, color="blue")
# Label the points
ne = vector((0.1, 0.1))
nw = vector((-0.1, 0.1))
offsets = (nw, nw, nw, ne, ne, ne)
P += sum(text(c, p + delta)
for c, p, delta in zip("ABCDEF", pos, offsets))
# Compute & plot the bezier curves
curves = bezier_fd(pos)
P += bezier_path(curves, color="red")
# The main points & the control points in a flat list
pts = list(chain.from_iterable(curves))
# The tangents
P += sum(line(t, color="#333", linestyle="--")
for t in chain([pts[:2]], zip(pts[2::3], pts[4::3]), [pts[-2:]]))
# The rest of the hull
P += sum(line(t, color="#333", linestyle=":")
for t in zip(pts[1::3], pts[2::3]))
# The control points
P += points([p for i, p in enumerate(pts) if i%3], color="magenta")
P.show(xmin=0, aspect_ratio=1, ticks=(0.5, 0.5), gridlines="minor")
Here's a live version of the script, running on the SageMathCell server.
And here's an SVG version of the basic path, without visible control points, etc.
Just for fun, here's a single interactive Bézier curve. You can drag points around and watch what happens. (It should work with a mouse or touchscreen). You can run it on the server, or download it to run it locally.
Given an arbitrary set of $n$ points, it's possible to find the equation of a unique polynomial of degree at most $n-1$ that passes through all $n$ points.
Let's consider a very simple example; say we we wanted to find the parabola that passes through the points $(-4,2)$, $(-2,-1)$, and $(1,1)$.
The general equation of a parabola is $$y=ax^2+bx+c$$ so we can substitute the $x$- and $y$-coordinates of the three points we have to create a system of three equations with three unknowns, like so: $$\begin{align} 2&=16a-4b+c \\ -1&=4a-2b+c \\ 1&=a+b+c \end{align}$$ Solving this system, we find that $a=\frac{13}{30}$, $b=\frac{11}{10}$, and $c=\frac{-8}{15}$. The parabola $y=\frac{13}{30}x^2+\frac{11}{10}x-\frac{8}{15}$ does indeed pass cleanly through all three points, so we're done.
And of course, you could do this with a system of six equations with six unknowns as well, but that would just take a little longer than this example.
Edit:
Mathematica reveals that the interpolating polynomial for those six points is approximately $$y=1.51771x^5-10.5198x^4+23.6478x^3-17.2688x^2-0.696339x+2.81946$$ The graph of that function looks like this:
I created my own code to generate what you're looking for using the math in this video: https://www.youtube.com/watch?v=WCGKqJrf4N4
I required this code for a flutter project using Dart.
class LinearInterpolationPainter extends CustomPainter {
final List<double> data;
final double spacing;
LinearInterpolationPainter({
this.data = const [],
this.spacing,
});
@override
void paint(Canvas canvas, Size size) {
//Used to determine the pixels between each data point
double _spacing = spacing ?? size.width / (data.length - 1);
var paint = Paint();
// TODO: Set properties to paint
paint.color = Colors.green.shade800;
paint.style = PaintingStyle.stroke;
paint.strokeWidth = 2.0;
var path = Path();
/// Set the path to initial data point
path.moveTo(0, size.height - data[0]);
// Draw path
for (int x = 0; x < data.length - 1; x++) {
for (int space = 0; space < _spacing; space++) {
path.lineTo(
(x * _spacing) + space,
size.height - interpolate(x + (space / _spacing)),
);
}
}
/// Draw path to final data point
path.lineTo(
(data.length - 1) * _spacing,
size.height - data.last,
);
canvas.drawPath(path, paint);
}
@override
bool shouldRepaint(CustomPainter oldDelegate) {
return true;
}
/// Calculates the y-value for an x-value given a set of points
double interpolate(double a) {
double result = 0;
for (int j = 0; j < data.length; j++) {
double product = data[j];
for (int i = 0; i < data.length; i++) {
if (i != j) product *= (a - i) / (j - i);
}
result += product;
}
return result;
}
}
