bezier curve produces degenerated geometry

[Luispain]

Hi everyone,

I got through this unexpected error of the Geom.PyTree.bezier() function.

I tried to reproduce the bezier curve from this figure (from A. Dumont ISABE 2024):

for this, I used the following script:

import Converter.Internal as I
import Geom.PyTree as D

bezier_pts = [
    (-2.700357630392316,0.0,0.0),
    (-2.6976262207607693,0.5003924846883171,0.0),
    (-2.229549849879404,0.5809392845644833,0.0),
    (-1.6778309560438065,0.5808324846960423,0.0),
    (-1.4307118882078287,0.5816481614489555,0.0),
    (-1.502325057810188,0.551365667106162,0.0),
    (-1.2473574798293745,0.551365667106162,0.0),
    (-1.0967228817002745,0.551365667106162,0.0),
    (-0.8981130129034307,0.6206294118001975,0.0),
    (-0.6963471318324306,0.6206294118001975,0.0),
    (-0.3981910727598996,0.6206294118001975,0.0),
    (-0.451287357252268,0.6010038239297818,0.0),
    (-0.25278893984233664,0.6014126903437489,0.0),
    (-0.0568255764847736,0.6009770063552681,0.0),
    (0.09593117914870675,0.6517117720686091,0.0),
    (0.250115568012407,0.6511758414448765,0.0),
    (0.6516740694594296,0.6509959425706399,0.0),
    (0.9020416781260892,0.8222244732381354,0.0),
    (1.5001298098941005,0.8216327768603116,0.0),
    (2.504222073695696,0.822016667345834,0.0),
    (3.929648935342521,0.07179152820182805,0.0),
    (4.0,0.0,0.0)]

polyline = D.polyline(bezier_pts)
bezier = D.bezier(polyline,N=5000)

y = I.getNodeFromName(bezier,'CoordinateY')[1]

ymax = y.max()
if ymax < 0.1:
    raise ValueError(f"with this input, at least one point should be bigger than 0.1, but got ymax={ymax}")

but I got an all degenerated curve as a result (all y was approximately equal to 0).

I do not understand what I got wrong, since the polyline geometry seems OK.

I used dev version of Cassiopée at ONERA LD.

Best regards and thank you in advance

Luis

[Luispain]

Hi again,

Just in case it may find it useful for debugging, please beware of the expected result:

import Converter.Internal as I
import Geom.PyTree as D
import numpy as np

def bezier_curve_2D(points, N=1000):
    """
       https://stackoverflow.com/questions/12643079/b%C3%A9zier-curve-fitting-with-scipy

       Given a set of control points, return the
       bezier curve defined by the control points.

       points should be a list of lists, or list of tuples
       such as [ [1,1], 
                 [2,3], 
                 [4,5], ..[Xn, Yn] ]
        N is the number of evaluation points

        See http://processingjs.nihongoresources.com/bezierinfo/
    """
    from scipy.special import comb

    def bernstein_poly(i, n, t):
        """
        The Bernstein polynomial of n, i as a function of t
        """
        return comb(n, i) * ( t**(n-i) ) * (1 - t)**i

    nPoints = len(points)
    xPoints = np.array([p[0] for p in points])
    yPoints = np.array([p[1] for p in points])

    t = np.linspace(0.0, 1.0, N)

    polynomial_array = np.array([ bernstein_poly(i, nPoints-1, t) for i in range(0, nPoints)   ])

    xvals = np.dot(xPoints, polynomial_array)
    yvals = np.dot(yPoints, polynomial_array)

    return xvals, yvals

bezier_pts = [
    (-2.700357630392316,0.0,0.0),
    (-2.6976262207607693,0.5003924846883171,0.0),
    (-2.229549849879404,0.5809392845644833,0.0),
    (-1.6778309560438065,0.5808324846960423,0.0),
    (-1.4307118882078287,0.5816481614489555,0.0),
    (-1.502325057810188,0.551365667106162,0.0),
    (-1.2473574798293745,0.551365667106162,0.0),
    (-1.0967228817002745,0.551365667106162,0.0),
    (-0.8981130129034307,0.6206294118001975,0.0),
    (-0.6963471318324306,0.6206294118001975,0.0),
    (-0.3981910727598996,0.6206294118001975,0.0),
    (-0.451287357252268,0.6010038239297818,0.0),
    (-0.25278893984233664,0.6014126903437489,0.0),
    (-0.0568255764847736,0.6009770063552681,0.0),
    (0.09593117914870675,0.6517117720686091,0.0),
    (0.250115568012407,0.6511758414448765,0.0),
    (0.6516740694594296,0.6509959425706399,0.0),
    (0.9020416781260892,0.8222244732381354,0.0),
    (1.5001298098941005,0.8216327768603116,0.0),
    (2.504222073695696,0.822016667345834,0.0),
    (3.929648935342521,0.07179152820182805,0.0),
    (4.0,0.0,0.0)]

# # BUG
# polyline = D.polyline(bezier_pts)
# bezier = D.bezier(polyline,N=5000)
# y = I.getNodeFromName(bezier,'CoordinateY')[1]

# expected result:
xy_pts = [ [p[0], p[1]] for p in bezier_pts ]
x,y = bezier_curve_2D(xy_pts)

ymax = y.max()
if ymax < 0.1:
    raise ValueError(f"with this input, at least one point should be bigger than 0.1, but got ymax={ymax}")

Best regards

Luis

[vincentcasseau]

Hi Luis, It should be fixed now, see #212. Thanks

[Luispain]

thank you very much for the fast response and fix !! 👏🏻