subdivide a 2d cubic curve into specific number of points in uniform (by arc length) and non uniform (curve speed) using pure python and no external libraries

cubic-curve-subdivide.py

import logging
import sys
import time
"""
based on:
https://pomax.github.io/bezierinfo/#tracing
Tracing a curve at fixed distance intervals from "A Primer on Bézier Curves" by Pomax
"""
logger = logging.getLogger('my_timer')
handler = logging.StreamHandler(sys.stdout)
logger.addHandler(handler)

def timeit(method):
    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 cubic_curve(p0, p1, p2, p3, t):
    x = (1-t)**3 * p0[0] + 3*(1-t)**2 * t * p1[0] + 3*(1-t) * t**2 * p2[0] + t**3 * p3[0]
    y = (1-t)**3 * p0[1] + 3*(1-t)**2 * t * p1[1] + 3*(1-t) * t**2 * p2[1] + t**3 * p3[1]
    return (x, y)

def cubic_curve_derivative(p0, p1, p2, p3, t):
    x = 3*(1-t)**2 * (p1[0] - p0[0]) + 6*(1-t) * t * (p2[0] - p1[0]) + 3*t**2 * (p3[0] - p2[0])
    y = 3*(1-t)**2 * (p1[1] - p0[1]) + 6*(1-t) * t * (p2[1] - p1[1]) + 3*t**2 * (p3[1] - p2[1])
    return (x, y)

def magnitude_squared(vector):
    return vector[0]**2 + vector[1]**2

def arc_length_lookup_table(p0, p1, p2, p3, num_steps=20):
    # Define the function to integrate
    def integrand(t):
        speed_squared = magnitude_squared(cubic_curve_derivative(p0, p1, p2, p3, t))
        return speed_squared**0.5

    # Create the lookup table
    lookup_table = {0: 0}
    dt = 1 / num_steps
    length = 0
    for i in range(1, num_steps + 1):
        t = i / num_steps
        # Use the trapezoidal rule to approximate the integral
        length += 0.5 * (integrand(t - dt) + integrand(t)) * dt
        lookup_table[t] = length

    return lookup_table

def find_t_for_arc_length(s, lookup_table):
    t_values = list(lookup_table.keys())
    arc_lengths = list(lookup_table.values())
    idx = next((i for i, x in enumerate(arc_lengths) if x >= s), len(arc_lengths) - 1)
    t0 = t_values[idx - 1]
    t1 = t_values[idx]
    length0 = arc_lengths[idx - 1]
    length1 = arc_lengths[idx]
    return t0 + (s - length0) * (t1 - t0) / (length1 - length0)

@timeit
def subdivide_curve(num_segments, p0, p1, p2, p3, uniform=False):
    """
    Subdivides a cubic curve into a specified number of segments.

    This function takes in the control points of a cubic curve and divides it into
    a specified number of segments. The subdivision can be either uniform or based
    on the speed of the curve.

    Args:
        num_segments (int): The number of segments to divide the curve into.
        p0 (tuple): The first control point of the cubic curve.
        p1 (tuple): The second control point of the cubic curve.
        p2 (tuple): The third control point of the cubic curve.
        p3 (tuple): The fourth control point of the cubic curve.
        uniform (bool, optional): Whether to divide the curve uniformly. Defaults to False.

    Returns:
        list: A list of points representing the subdivided curve.
    """
    if uniform is False:
        t_values = [i / num_segments for i in range(num_segments + 1)]
        return [cubic_curve(p0, p1, p2, p3, t) for t in t_values]
    lookup_table = arc_length_lookup_table(p0, p1, p2, p3)
    total_length = list(lookup_table.values())[-1]
    target_arc_lengths = [(i / num_segments) * total_length for i in range(num_segments + 1)]
    interpolated_t_values = [find_t_for_arc_length(s, lookup_table) for s in target_arc_lengths]
    points = [cubic_curve(p0, p1, p2, p3, t) for t in interpolated_t_values]
    return points

def _test():
    pass
if __name__ == '__main__':
    size(1000, 2000)
    translate(200, 100)
    p0 = (0, 0)
    p1 = (660, 651)
    p2 = (-160, 670)
    p3 = (520, 80)
    num_segments = 30
    points_by_length = subdivide_curve(num_segments, p0, p1, p2, p3, False)
    fill(1, 0, 0)
    for p in points_by_length:
        oval(*p, 4, 4)
    fontSize(30)
    text('Curve Speed', (200, 0))

    fill(0)
    oval(*p0, 7, 7)
    oval(*p1, 7, 7)
    oval(*p2, 7, 7)
    oval(*p3, 7, 7)
    fill(None)
    stroke(.7)
    line(p0, p1)
    line(p2, p3)

    points_by_length = subdivide_curve(num_segments, p0, p1, p2, p3, True)
    translate(0, 1000)
    fill(0, 0, 1)    
    for p in points_by_length:
        oval(*p, 4, 4)

    fontSize(30)
    text('Uniform', (200, 0))

    fill(0)
    oval(*p0, 7, 7)
    oval(*p1, 7, 7)
    oval(*p2, 7, 7)
    oval(*p3, 7, 7)
    fill(None)
    stroke(.7)
    line(p0, p1)
    line(p2, p3)

Sample render: