calculating cubic bezier curve offset using normal vector doesnt produce correct offset

Question

I'm trying to offset a cubic bezier curve, by doing the following

find the x and y derivative to get the tangent vector (dx, dy)
rotate the vector 90deg by (-dy, dx) to get the normal
making it a unit vector by dividing by the magnitude
multiplying by the desired offset, and adding it the the point on the curve

def cubic_offset(P0, P1, P2, P3, t, dist): 
    initx, inity = cubic(P0, P1, P2, P3, t) 
    dx, dy = cubic_dt(P0, P1, P2, P3, t) # this is the tangent vector
    
    normx, normy = -dy, dx
    mag = sqrt(normx**2 + normy**2)
    normx, normy = normx/mag, normy/mag
    
    return initx + dist*normx, inity + dist*normy

however, this isnt working, esspecialy near extrema

I also saw this blog post https://observablehq.com/@s-silva/bezier-curve-offsets that confirms this process actually works, and provides a preview that doesnt produce errors.

Our code matches almost perfectly, however his works and mine doesnt.

my derivative is also correct as Im using it to get the extrema and the bounding box which is working perfectly fine.

"Our code matches almost perfectly" except it's also _wildly_ different because their code uses the Chebyshev approximation of the curve, while your code doesn't. Which means your derivatives are going to potentially wildly differ too (as Bezier derivatives don't give you the geometric derivative, they give you the _time_ derivative) — Mike 'Pomax' Kamermans, Jul 01 '23 at 17:52
see [express "t" variable from Cubic Bezier Curve equation](https://stackoverflow.com/questions/60073813/is-it-possible-to-express-t-variable-from-cubic-bezier-curve-equation/60113617?r=SearchResults&s=1%7C48.2554#60113617) and also [draw outline for some connected lines](https://stackoverflow.com/a/22068534/2521214) — Spektre, Jul 02 '23 at 10:10
Mike 'Poxmax' Kamermans, I figured that an approximation is still the same thing, sorry. I don't know what you mean by bezier derivatives, but treating the curve as a parametric equation (f(t), g(t)) the vector (f'(t), g'(t)) is the vector tangent to the curve at t and it is in that regard a "Geometric derivative" as well as a time derivative. A "normal" derivative dy/dx would be calculated by g'(t)/f'(t) which could be what you meant. Anyway, I figured the error and posted an Answer. — 12ksins, Jul 02 '23 at 19:47

12ksins · Answer 1 · 2023-07-02T20:15:57.123

Solved:

Because the offset was only wrong around extrema, where the derivative is close to 0, I figured it could have something to do with floating-point arethmatic.

I copied the derivative thats on Wikipedia

3(1-t)^2(P1 - P0) + 6t(1-t)(P2-P1) + 3t^2(P3-P2)

and it worked perfectly fine, even though it is mathematically equivalent to my derivative

3t^2(P3 - 3P2 + 3P1 - P0) + 6t(P2 - 2P1 + P0) + 3(P1 - P0)

Apparently, when calculated by Python, the average difference between them for t between 0 and 1 and P's between 0 and 700 is 0.049 (Calculated 2,000,000 times) and for t between 0 and 100 the average difference is a staggering 43.2.

Any corrections to how this error is produced would be most appreciated, as I'm not that knowledgeable in such things.

EDIT:

Turns out it's not floating-points to blame, but tuples! or maybe something else I dont know, any help would be appreciated.

from random import randint, uniform

def dt1(P0, P1, P2, P3, t):
    return (3*(1-t)**2*(P1[0] - P0[0]) + 6*(1-t)*t*(P2[0] - P1[0]) + 3*t**2*(P3[0] - P2[0]),
            3*(1-t)**2*(P1[1] - P0[1]) + 6*(1-t)*t*(P2[1] - P1[1]) + 3*t**2*(P3[1] - P2[1]))

def dt2(P0, P1, P2, P3, t):
    return (t**2*3*(P3[0]-3*P2[0]+3*P1[0]-P0[0])+t*6*(P2[0]-2*P1[0]+P0[0])+3*(P1[0]- P0[0]),
            t**2*3*(P3[1]-3*P2[1]+3*P1[1]-P0[1])+t*6*(P2[1]-2*P1[1]+P0[0])+3*(P1[1]- P0[0]))
            

def dt1_(P0, P1, P2, P3, t):
    return 3*(1-t)**2*(P1 - P0) + 6*(1-t)*t*(P2 - P1) + 3*t**2*(P3 - P2)

def dt2_(P0, P1, P2, P3, t):
    return t**2*3*(P3-3*P2+3*P1-P0)+t*6*(P2-2*P1+P0)+3*(P1- P0)

reps = 1_000_000

# with tuples t between 0 and 1
sum = 0
for i in range(reps):
    Points = [(randint(0, 700), randint(0, 700)) for _ in range(4)]
    t = uniform(0, 1)
    dt_1 = dt1(*Points, t)
    dt_2 = dt2(*Points, t)
    sum += dt_1[0] - dt_2[0]
    sum += dt_1[1] - dt_2[1]
    
print(f"for t between 0 and 1, with tuples, the average diffrence is {sum/(reps*2)}")

# with tuples, t between 0 and 100
sum = 0
for i in range(reps):
    Points = [(randint(0, 700), randint(0, 700)) for _ in range(4)]
    t = uniform(0, 100)
    dt_1 = dt1(*Points, t)
    dt_2 = dt2(*Points, t)
    sum += dt_1[0] - dt_2[0]
    sum += dt_1[1] - dt_2[1]
    
print(f"for t between 0 and 100, the average diffrence is {sum/(reps*2)}")

# without tuples, t between 0 and 1
sum = 0
for i in range(reps*2):
    Points = [randint(0, 700) for _ in range(4)]
    t = uniform(0, 1)
    dt_1 = dt1_(*Points, t)
    dt_2 = dt2_(*Points, t)
    sum += dt_1 - dt_2

print(f"for t between 0 and 1, without tuples, the average diffrence is {sum/(reps*2)}")

# without tuples, t between 0 and 1
sum = 0
for i in range(reps*2):
    Points = [randint(0, 700) for _ in range(4)]
    t = uniform(0, 100)
    dt_1 = dt1_(*Points, t)
    dt_2 = dt2_(*Points, t)
    sum += dt_1 - dt_2

print(f"for t between 0 and 100, without tuples, the average diffrence is {sum/(reps*2)}")

Results:

for t between 0 and 1, with tuples, the average diffrence is 0.10319855785147072
for t between 0 and 100, with tuples, the average diffrence is -21.841299912204903
for t between 0 and 1, without tuples, the average diffrence is -3.370475946951057e-17
for t between 0 and 100, without tuples, the average diffrence is -1.1170571903237891e-12

`"Turns out it's not floating-points to blame"`, no, it very much is. Remove the tuples, i.e. `def dt1(P0, P1, P2, P3, t): return 3*(1-t)**2*(P1-P0)+...` so that you're only evaluating one dimension, then do the same for dt2, and compare them: we can trivially see these yield completely different results for the same inputs. E.g. `dt1(0,1,1,0, 0.2) ->1.8000000000000003` but `dt2(0,1,1,0, 0.2) -> 1.7999999999999998`. This is a *classic* example of floating point rounding. For numerically unstable math like this, use [mpmath](https://mpmath.org) instead. — Mike 'Pomax' Kamermans, Jul 19 '23 at 17:55

calculating cubic bezier curve offset using normal vector doesnt produce correct offset

1 Answers1