find inflection points on a cubic bezier curve using python

cubic-bezier-inflection-point.py

import logging
import sys
import time
logger = logging.getLogger('my_timer')
handler = logging.StreamHandler(sys.stdout)
logger.addHandler(handler)

def timeit(method):
    """
    A decorator that makes it possible to time functions.
    """
    def timed(*args, **kw):
        logger.setLevel(logging.DEBUG)
        ts = time.time()
        result = method(*args, **kw)
        te = time.time()
        if 'log_time' in kw:
            name = kw.get('log_name', method.__name__.upper())
            kw['log_time'][name] = int((te - ts) * 1000)
        else:
            logger.debug('%r  %2.2f ms' %(method.__name__, (te - ts) * 1000))
        logger.setLevel(logging.WARNING)
        return result
    return timed

def align_to_x(p0, p1, p2, p3):
    """
    moves p0 to (0, 0) and after rotate the points so y of p3 becomes 0, making the curve align to x axis.
    """
    # convert points to complex numbers
    p0 = complex(*p0)
    p1 = complex(*p1)
    p2 = complex(*p2)
    p3 = complex(*p3)

    # move p0 to (0, 0)
    p1 -= p0
    p2 -= p0
    p3 -= p0
    angle = p3.imag / p3.real

    # rotate using complex number multiplication
    rotation = complex(1, -angle) / (1 + angle**2)**0.5
    p1 *= rotation
    p2 *= rotation
    p3 *= rotation

    return ((0, 0), (p1.real, p1.imag), (p2.real, p2.imag), (p3.real, p3.imag))

@timeit
def calculate_inflection_points(p0, p1, p2, p3):
    # moving the curve to origin (0, 0) to remove p0 from equation
    _p0, _p1, _p2, _p3 = align_to_x(p0, p1, p2, p3)
    x2, y2 = _p1
    x3, y3 = _p2
    x4, _ = _p3

    # coefficients
    a = x3 * y2
    b = x4 * y2
    c = x2 * y3
    d = x4 * y3

    # coefficients of quadratic equation
    x = -3*a + 2*b + 3*c - d
    y = 3*a - b - 3*c
    z = c - a

    # the discriminator is negative or if x is zero
    if y**2 - 4*x*z < 0 or x == 0:
        return []

    # the roots of the quadratic equation
    t1 = (-y + (y**2 - 4*x*z)**0.5) / (2*x)
    t2 = (-y - (y**2 - 4*x*z)**0.5) / (2*x)

    # remove roots that are not in the Bézier interval [0,1]
    inflection_points = []
    for t in [t1, t2]:
        if 0 <= t <= 1:
            # Calculate the inflection point coordinates
            x_inflection = (1-t)**3 * p0[0] + 3*(1-t)**2*t * p1[0] + 3*(1-t)*t**2 * p2[0] + t**3 * p3[0]
            y_inflection = (1-t)**3 * p0[1] + 3*(1-t)**2*t * p1[1] + 3*(1-t)*t**2 * p2[1] + t**3 * p3[1]
            inflection_points.append((x_inflection, y_inflection))

    return inflection_points


size(500, 500)
fill(1)
rect(0, 0, 500, 500)
p0 = (60, 220)
p1 = (360, 450)
p2 = (80, 40)
p3 = (381, 170)

# p0, p1, p2, p3 = align_to_x(p0, p1, p2, p3)
# translate(200, 200)
stroke(0)
fill(None)
newPath()
moveTo(p0)
curveTo(p1, p2, p3)
drawPath()

ips = calculate_inflection_points(p0, p1, p2, p3)

if ips:
    for p in ips:
        oval(p[0]-5, p[1]-5, 10, 10)


# import os
# from datetime import datetime
# import subprocess

# fn = os.path.split(os.path.realpath(__file__))[1] + datetime.now().strftime("%Y-%m-%d-%H%M%S")
# saveImage(f"{fn}.png")
# subprocess.call(['open', "%s.png" %fn])

Sample render:

Based on Primer on Bezier Curves. This function still doesn't return inflection points where the curve creates a knot and self intersects.