curve.f90:evaluate_curve_barycentric() broken when num_nodes exceeds 30

[saxenarohan97]

I want to use the module for fitting a curve through more than 30 points. However, I found that the interpolation "breaks" for >30 points. Minimal reproducible example:

def create_curve(npoints):

    x = sorted(np.random.randint(1000, size=npoints).astype('float'))
    y = sorted(np.random.randint(1000, size=npoints).astype('float'))

    nodes = np.asfortranarray([x, y])
    
    curve = bezier.Curve(nodes, degree=3)
    
    curve.plot(num_pts=100)

Note that the points I've passed to bezier.Curve are monotonically increasing, so the resulting approximation curve should also be monotonically increasing.

If I fit a curve through 30 points, then it works as expected:
>>> create_curve(30)

However, if I use >30 points:
>>> create_curve(31)

Then the curve loops in on itself - which it shouldn't - since all points are positive and monotonically increasing.

Any idea what might be the issue (or if there is a problem with my example code)?

[dhermes]

@saxenarohan97 Thanks for stopping by, it's great to hear from users! Sorry for my 12 day (:scream:) delay here (advent of code has been stealing all my free time).

However, I found that the interpolation "breaks" for >30 points

So, the Curve constructor does not do any interpolation. See in particular the note on the second argument degree:

degree (int) – The degree of the curve. This is assumed to correctly correspond to the number of nodes. Use from_nodes() if the degree has not yet been computed.

You should just use Curve.from_nodes(nodes) (the factory constructor). Using bezier.Curve(nodes, 3) is "undefined behavior" and I'm not actually sure what it's doing under the covers (will dig into that because I'm curious, but not because it's supported behavior).

Then the curve loops in on itself - which it shouldn't - since all points are positive and monotonically increasing.

I have observed this as well (see notebook), even in a case where I use degree elevation (which should only change the nodes, not the curve itself) and degree 29 (i.e. 30 nodes) plot

becomes

Thanks for surfacing this bug, hopefully something comes of it! (I have marked the issue as type: investigating for now.)

[dhermes]

So, one preliminary finding is that the fortran version of evaluate_multi seems to be what is causing the breakage. By using the pure Python (updated in the gist), plotting the same degree 30 (i.e. 31 nodes) curve we see

It does this by isolating the pure-Python components (from the Fortran speedup binary extension)

import bezier._curve_helpers


def _bezier_evaluate_multi_py(nodes, s_vals):
    one_less = 1.0 - s_vals
    return bezier._curve_helpers._evaluate_multi_barycentric(
        nodes, one_less, s_vals
    )

---

Another side note:

The reason this "worked" with a mismatch of degree and the number of nodes is that Curve.plot() uses Curve.evaluate_multi(), which passes nodes directly to a low-level helper without referencing Curve.degree.

[dhermes]

Regarding INTERPOLATION, you can do this with a little bit of math by using an "elevation matrix" starting from the identity matrix. For example, to produce the matrix in the Curve.reduce_() docs:

(1/3) * [3 1 0 0]
        [0 2 2 0]
        [0 0 1 3]

via

>>> import bezier
>>> import numpy as np
>>> degree3_identity = bezier.Curve.from_nodes(np.eye(3, order="F"))
>>> degree3_elevated = degree3_identity.elevate()
>>> degree3_elevated.nodes * 3                                                                                                                                                
array([[3., 1., 0., 0.],
       [0., 2., 2., 0.],
       [0., 0., 1., 3.]])

We can extend this to elevating multiple degrees, e.g. from degree 3 (4 nodes) to degree 30 (31 nodes):

>>> curve = bezier.Curve.from_nodes(np.eye(4, order="F"))
>>> for _ in range(31 - 4):
...     curve = curve.elevate()
... 
>>> curve.degree
30
>>> elevation_matrix = curve.nodes
>>> elevation_matrix.shape
(4, 31)

This "elevation matrix" E = E_{3,30} can be used to take nodes N_{3} of shape (d, 4) for a cubic and elevate them to new nodes N_{30} (defining the same curve) of shape (d, 31) via

N_{30} = N_{3} E

This is great if we know N_{3} but in your case you have N_{30} and want to compute N_{3}. In general, there is no unique solution to the system in that direction, but it's quite common to use the pseudo-inverse / least-squares solution:

N_{30} [ (E E^T)^{-1} E ]^T

which translates to

>>> interpolation_matrix = np.linalg.solve(
...     elevation_matrix.dot(elevation_matrix.T),
...     elevation_matrix,
... ).T

Doing this to our "broken" curve of degree 30 (31 nodes) from the gist produces a cubic that "looks right":

---

NOTE: When computing the elevation matrix, each intermediate stage can be computed exactly as rational numbers, so for some applications it may make sense to avoid roundoff as long as possible. This can be done with something like fractions.Fraction, something a bit more heavyweight like sympy.Matrix with sympy.Rational or by judicious tracking of a denominator along with integer matrices (inverse can be computed in a similar way by using denominators via Cramer's rule).

[saxenarohan97]

Hi @dhermes. Thanks a lot for your reply and detailed investigation of the issue. I too must apologize for my late response!

I suppose the main takeaway for users from your comment is that calling bezier.Curve directly is not recommended, rather bezier.Curve.from_nodes is the right method to be called.

However, if I replace bezier.Curve with bezier.Curve.from_nodes in my above code, it still gives me the same issue:

from bezier import Curve
import numpy as np

def create_curve(npoints):

    x = sorted(np.random.randint(1000, size=npoints).astype('float'))
    y = sorted(np.random.randint(1000, size=npoints).astype('float'))

    nodes = np.asfortranarray([x, y])
    curve = Curve.from_nodes(nodes)
    curve.plot(num_pts=100)

>>> create_curve(30)

>>> create_curve(31)

Am I still messing up something? :(

[dhermes]

Thanks for coming by @saxenarohan97.

Am I still messing up something? :(

I fixed the issue, but have not released the fix yet. You are doing nothing wrong, thanks for checking! I hope to make the 0.11.0 release today or over the weekend.

I suppose the main takeaway for users from your comment is that calling bezier.Curve directly is not recommended, rather bezier.Curve.from_nodes is the right method to be called.

Not necessarily, but rather that you should use the constructor with the correct degree. (But your accidental mis-usage of the constructor was great because it contributed the verify constructor argument feature added in #163.)

[dhermes]

@saxenarohan97 OK version 0.11.0 has been released