[Question] How can I find a cubic Bezier curve that passes through a specified point?

[miu200521358]

Thank you for publishing a very nice library!
I'm currently working on a problem with this library.

I have values measured over time.
Example)

ys = [0, 7.027885760013675e-07, 1.3352245924580508e-05, 7.185419521038572e-05, 0.000247954241316628, 0.0006740029696562511, 0.0015907714640571724, 0.0034294112039156, 0.006821667242492446, 0.013031496324572234, 0.024556098240622215, 0.04575246255862098, 0.08269242982725422, 0.1322613906866914, 0.17519720275554274, 0.20047545179963588]
xs = np.arange(0, len(ys))

fig, ax = plt.subplots(figsize=(7, 7))
ax.plot(xs, ys, 'bo-', ms=2)
plt.show()

I want to generate a cubic Bezier curve that passes through these values.
(Not the Catmull-Rom curve)

If you know how to solve this problem by using your library, please let me know.
Thank you.

[dhermes]

@miu200521358 I can likely help you but want some clarification. Your 16 points will overdetermine a cubic (4 control points). Do you mean a cubic that exactly fits or a cubic this is the "best fit" in some metric? Or did you just mean a Bézier curve (not not actually need to be cubic)?

[miu200521358]

What I want is an "exactly fits" cubic Bezier curve.
In the case of the example, I suspect that a Bezier curve like the one in the image below can be obtained.

I would like to know how to find p1 and p2 of this Bezier curve.
(P0 and p3 are known because they are the start and end points)
Thank you.

[dhermes]

So you have 16 data points, but only 4 free points for a cubic. One way to approximate the control points would be to use the 4x4 identity matrix as a symbolic representation of the 4 control points (1 per row).

For example, consider our "symbolic" control points p0, p1, p2, p3 and degree elevate to a curve with control points q0, q1, q2, q3, q4:

In [1]: import bezier

In [2]: import numpy as np

In [3]: representative = bezier.Curve.from_nodes(np.eye(4))

In [4]: representative.elevate().nodes
Out[4]: 
array([[1.  , 0.25, 0.  , 0.  , 0.  ],
       [0.  , 0.75, 0.5 , 0.  , 0.  ],
       [0.  , 0.  , 0.5 , 0.75, 0.  ],
       [0.  , 0.  , 0.  , 0.25, 1.  ]])

I.e. this says q0 = p0, q1 = 0.25 p0 + 0.75 p1, q2 = 0.5 p1 + 0.5 p2, etc.

In order to "project" our 16 points down to just 4, we can find a transform in the other direction by degree elevating from degree 3 to degree 15:

In [5]: elevated = representative

In [6]: for _ in range(4, 15 + 1):
   ...:     elevated = elevated.elevate()
   ...: 

In [7]: elevated.degree
Out[7]: 15

In [8]: transform = elevated.nodes.T

In [9]: transform.shape
Out[9]: (16, 4)

Using our 16 data points from above

In [10]: ys = np.array([0, 7.027885760013675e-07, 1.3352245924580508e-05, 7.185419521038572e-05, 0.000247954241316628, 0.0006740029696562511, 0.0015907714640571724, 0.00342941120391
    ...: 56, 0.006821667242492446, 0.013031496324572234, 0.024556098240622215, 0.04575246255862098, 0.08269242982725422, 0.1322613906866914, 0.17519720275554274, 0.20047545179963588
    ...: ])

In [11]: xs = np.arange(0, len(ys))

In [12]: nodes = np.vstack([xs, ys])

we can use a least squares projection to find some best fit

In [13]: reduced_t, residuals, rank, _ = np.linalg.lstsq(transform, nodes.T, rcond=None)

In [14]: residuals
Out[14]: array([5.47125882e-29, 8.07090066e-04])

In [15]: rank
Out[15]: 4

In [16]: reduced = reduced_t.T

This shows we have a very good fit for x and a pretty good fit for y. Now that we've got this curve, we've got our p1 and p2:

In [17]: reduced
Out[17]: 
array([[-9.02365312e-15,  5.00000000e+00,  1.00000000e+01,
         1.50000000e+01],
       [-2.26101770e-04,  8.68120557e-03, -5.17256724e-02,
         2.14974631e-01]])

We can turn this into a curve, and plot it, then compare that plot to the data points.

In [18]: import matplotlib.pyplot as plt

In [19]: import seaborn
 
In [20]: seaborn.set()

In [21]: curve = bezier.Curve.from_nodes(reduced)

In [21]: ax = curve.plot(256)

In [22]: _ = ax.plot(xs, ys, marker="o", linestyle="none")

As you can see, the fit is good, but far from perfect. This means other methods (this was essentially a least squares fit) would also fail to find an exact match.

[dhermes]

So I now realize my answer is "wrong" in that I treated the ys as control points vs. points on the curve. I'll re-do the answer correctly below (but am not yet done).

Backing up a bit, I'm assuming your points on the curve xs and ys are evenly spaced inputs, i.e. points B(s) on a curve for parameters s = 0, 1/n, ..., (n - 1)/n, 1. Here we can do a different kind of fit with our representative curve (that has quasi symbolic inputs). For example, if we set n = 4 so that there are 5 evenly spaced points s = 0, 0.25, 0.5, 0.75, 1.0 we see

In [1]: import bezier

In [2]: import numpy as np

In [3]: representative = bezier.Curve.from_nodes(np.eye(4))

In [4]: s_vals = np.linspace(0.0, 1.0, 5)

In [5]: s_vals
Out[5]: array([0.  , 0.25, 0.5 , 0.75, 1.  ])

In [6]: with np.printoptions(linewidth=1000):
   ...:     print(representative.evaluate_multi(s_vals))
   ...: 
[[1.       0.421875 0.125    0.015625 0.      ]
 [0.       0.421875 0.375    0.140625 0.      ]
 [0.       0.140625 0.375    0.421875 0.      ]
 [0.       0.015625 0.125    0.421875 1.      ]]

In other words, B(0) = p0, B(0.25) = 0.421875 p0 + 0.421875 p1 + 0.140625 p2 + 0.015625 p3, B(0.5) = 0.125 p0 + 0.375 p1 + 0.375 p2 + 0.125 p3 and so on. So if we had the five values B(0), B(0.25), ..., B(1) but did not have p0, p1, p2, p3 we could use this overdetermined 5 x 4 system to find a least squares solution.

This is exactly what we can do for our 16 points (n = 15 in the s = 0, 1/n, ... set up):

In [7]: s_vals = np.linspace(0.0, 1.0, 16)

In [8]: transform = representative.evaluate_multi(s_vals).T

In [9]: ys = np.array([0, 7.027885760013675e-07, 1.3352245924580508e-05, 7.185419521038572e-05, 0.000247954241316628, 0.0006740029696562511, 0.0015907714640571724, 0.0034294112039156, 0.006821667242492446, 0.0
   ...: 13031496324572234, 0.024556098240622215, 0.04575246255862098, 0.08269242982725422, 0.1322613906866914, 0.17519720275554274, 0.20047545179963588])

In [10]: xs = np.arange(0, len(ys))

In [11]: nodes = np.vstack([xs, ys])

In [12]: reduced_t, residuals, rank, _ = np.linalg.lstsq(transform, nodes.T, rcond=None)

In [13]: residuals
Out[13]: array([5.22252497e-29, 8.07090066e-04])

In [14]: rank
Out[14]: 4

In [15]: reduced = reduced_t.T

In [16]: with np.printoptions(linewidth=1000):
    ...:     print(reduced)
    ...: 
[[-1.76090753e-15  5.00000000e+00  1.00000000e+01  1.50000000e+01]
 [-2.26101770e-04  1.50843117e-02 -7.65425640e-02  2.14974631e-01]]

But as before the x residual is very tight and y residual is good but not great. Plotting confirms this:

In [17]: import matplotlib.pyplot as plt

In [18]: import seaborn

In [19]: seaborn.set()

In [20]: curve = bezier.Curve.from_nodes(reduced)

In [21]: ax = curve.plot(256)

In [22]: _ = ax.plot(xs, ys, marker="o", linestyle="none")

In [23]: plt.show()

[miu200521358]

Dear dhermes

Thank you for your consideration!
Certainly, it seems that an approximate curve can be obtained with this method.
I was very helpful.
I will continue to support your success.
Thank you very much!