How do I split a cubic bezier curve into equal segments?

Assume I have a cubic bezier curve (whose length I know). How do I split it into equal segments, let's say, three equal segments?

In other words, how do I calculate the t parameters that split it into three equal segments? Splitting at t = 0.33 and t = 0.66 will not produce equal segments. In fact, splitting at any t other than 0.5 on a symmetric curve will not produce equal segments.

Example of t = 0.33 and t = 0.66 on an actual curve:

I havn't any magic equation making the relationship between the length of a Bezier curve from X(0) and X(t).

The best I can propose is to compute the length by integration (by recursivity, dividing the Bezier curve in 2 until it is "tiny"), and afterward, compute the same integration stopping at the interesting distance (a fraction of the first computed distance).

In OCaml, this could be done easily. Note that bezier_stop takes 4 points - a,b,c,d - a current_length (the length of the Bezier segment before a) and the threshold. The idea, is if the desired point is in the Bezier, it will return Return t where t is what you want (relative to the current Bezier... the caller will have to scale this value). If the desired point is after the d point (your Bezier curve is too short), return Length l where l is the length of the current Bezier curve. (See Length and Result as tags which are added to the values).

You will find a Return 0.5 which will surprise you. At a very closed neighbour of the desired point (at a point we refused to dig further), we can't have a precise t. It is not important, since the imprecision is divided by 2 each time the recusive function returns. But an interpolation should be more efficient. ((threshold-current_len)/dist) typically). This should make the program usable with an higher epsilon.

Here, the program computes t for a length 2.0 and len-2.0. Since the proposed curve is symetrical, the sum must be 1.

Many addition of length have the same magnitude, this is a good way to improve the precision.

let epsilon = 0.001
let middle (x,y) (x',y') = ((x+.x')/.2.,(y+.y')/.2.)

let close (x,y) (x',y') = Float.abs(x-.x') +. Float.abs(y-.y') < epsilon

let square x = x *. x
let distance (x,y) (x',y') = Float.sqrt (square (x-.x')+.square (y-.y'))

let rec bezier_length a b c d =
  if close a d then
    distance a d
  else
    let e = middle a b in
    let f = middle b c in
    let g = middle c d in
    let h = middle e f in
    let i = middle f g in
    let j = middle h i in
    bezier_length a e h j +. bezier_length j i g d

type r = Length of float | Result of float

let rec bezier_stop a b c d  current_len threshold =
  if close a d then
    let dl = distance a d in
      if current_len +. dl < threshold then
        Length dl
      else
        Result 0.5
  else
    let e = middle a b in
    let f = middle b c in
    let g = middle c d in
    let h = middle e f in
    let i = middle f g in
    let j = middle h i in
    match bezier_stop a e h j current_len threshold with
    | Result r -> Result (r /. 2.)
    | Length l -> begin
      match bezier_stop j i g d (current_len +. l) threshold with
      | Result r' -> Result ((r' +. 1.)/. 2.)
      | Length l' -> Length (l +. l')
    end

let a = (0. , 0.)
let b = (0. , 1.)
let c = (10. , 1.)
let d = (10. , 0.)

let len = bezier_length a b c d

let () = 
  (match bezier_stop a b c d  0. len/.3 with 
   Result x -> print_float x
  | Length _ -> print_string "not found");
  print_string "\n";
  (match bezier_stop a b c d  0. (2.*.l/.3.) with 
   Result x -> print_float x
  | Length _ -> print_string "not found")

or in Python. The language makes it a little complex to return either a length, either a Returned t. (Here, with two classes).

import math

epsilon = 0.01
def middle(a,b):
    (x1,y1)=a
    (x2,y2)=b
    return ((x1+x2)/2, (y1+y2)/2)

def distance(a,b):
    (x1,y1)=a
    (x2,y2)=b
    return math.sqrt( (x1-x2)**2 + (y1-y2)**2)

def bezier_length(a, b, c, d):
    dist = distance(a,d)
    if dist < epsilon:
        return dist
    else:
        e = middle(a, b)
        f = middle(b, c)
        g = middle(c, d)
        h = middle(e, f)
        i = middle(f, g)
        j = middle(h, i)
        return bezier_length(a, e, h, j) + bezier_length(j, i, g, d)

class Length:
    def __init__(this,l):
        this.length = l

class Return:
    def __init__(this,t):
        this.t = t

def bezier_stop(a, b, c, d,  current_len, threshold):
    dist = distance(a,d)
    if dist < epsilon:
        if current_len + dist < threshold:
            return Length(dist)
        else:
            return Return((threshold-current_len)/dist)
    else:
        e = middle(a, b)
        f = middle(b, c)
        g = middle(c, d)
        h = middle(e, f)
        i = middle(f, g)
        j = middle(h, i)
        bs = bezier_stop(a, e, h, j, current_len, threshold)
        if isinstance(bs, Return):
            return Return(bs.t/2)
        else:
            l = bs.length
            bs2 = bezier_stop(j, i, g, d, current_len + l, threshold)
            if isinstance(bs2, Return):
                return Return(( bs2.t + 1) / 2.)
            else:
                return Length(l + bs2.length)

a = (0,0)
b = (0,1)
c = (10,1)
d = (10,0)

l = bezier_length(a, b, c, d)

bs1 = bezier_stop(a, b, c, d,  0., l/3)
print(bs1.t) 

bs2 = bezier_stop(a, b, c, d,  0., 2*l/3)
print(bs2.t)

You have the value of t in bs1.t and bs2.t. The sum is 1 with a very good precision. You can increase epsilon. The program will be faster be less precise.