Split a cubic bezier curve at multiple points

I'm writing an algorithm that will split a cubic bezier curve into multiple curves (up to 4). I have the t values for each point I want to split at from the beginning. I also have an algorithm already for splitting the curve once:

SubdivPoints subdivideBezier(Vector2 p0, Vector2 p1, Vector2 p2, Vector2 p3, float t)
{
    Vector2 p11 = (p1 - p0) * t + p0;
    Vector2 p21 = (p2 - p1) * t + p1;
    Vector2 p31 = (p3 - p2) * t + p2;

    Vector2 p12 = (p21 - p11) * t + p11;
    Vector2 p22 = (p31 - p21) * t + p21;

    Vector2 p13 = (p22 - p12) * t + p12;

    return SubdivPoints(p11, p12, p22, p31, p13);
}

My question is, is there a simple way to expand this for splitting multiple times? I imagine after each split I want to recalculate the t values; one thing I'm wondering is if simple arithmetic will work here. E.g. Let's say I have t values 0.2 and 0.6. I split the curve at t = 0.2, giving me two curves. The second curve covers t values 0.2 < t < 1 from the original. If I recalculate the second t value of 0.6 by division: (0.6 - 0.2) / (1 - 0.2) = 0.5, then divide the second curve at t = 0.75, will that work? Or do I need to recalculate it some other way?

回答1:

very hard to say even if your current approach work as we do not see what is behind SubdivPoints. I can think of 2 approaches for this:

algebraic

if you look at your problem as a parameter t scaling of polynomial then it became clearer. For example you want the control points for part t=<0.25,0.5>. That means we need form new control points representing the same curve with parameter u=<0.0,1.0>. How to do it?
1. write BEZIER as polynomial
each axis can be done separately in the same manner s lets us focus on x-axis only. We have 4 BEZIER control points with x-coordinates (x0,x1,x2,x3). If we apply bernstein polynoms to form cubic polynomial we get:
```
x(t)=      (                           (    x0))
    +    t*(                  (3.0*x1)-(3.0*x0))
    +  t*t*(         (3.0*x2)-(6.0*x1)+(3.0*x0))
    +t*t*t*((    x3)-(3.0*x2)+(3.0*x1)-(    x0))
```
1. rescale parameter by substitution
Use linear interpolation for that so:
```
t0=0.25 -> u0=0.0
t1=0.50 -> u1=1.0
t=t0+(t1-t0)*(u-u0)/(u1-u0)
t=0.25+0.5*u
```
now rewrite the polynomial using u instead of t
```
x(t)=             (                           (    x0))
    +(0.25+u/2)  *(                  (3.0*x1)-(3.0*x0))
    +(0.25+u/2)^2*(         (3.0*x2)-(6.0*x1)+(3.0*x0))
    +(0.25+u/2)^3*((    x3)-(3.0*x2)+(3.0*x1)-(    x0))
```
1. change the polynomial form to match BEZIER equation again
Now you need to separate terms of polynomial for u^0,u^1,u^2,u^3 and reform each to match BEZIER style (from #1). From that extract new control points.

For example term u^3 is like this. The only possibility of getting u^3 is from
```
(1/4+u/2)^3= (1/8)*u^3 + (3/16)*u^2 + (3/32)*u + (1/64)
```
so separated u^3 term will be:
```
u*u*u*((    x3)-(3.0*x2)+(3.0*x1)-(    x0))/8
```
the others will be a bit more complicated as you need combine all lines together ... After simplifying you will need to separate the new coordinates. As you can see it is slight madness but there are algebraic solvers out there like Derive for Windows ...

Sorry I do not have time/mood/stomach to make full example of all the terms but you will see it will be quite a polynomial madness...
curve fitting

this is based on that you are finding the coordinates of your control points and checking how far it is from your desired curve. So test ""all possible" points and remember the closes match between original curve on target range and new one. That would be insolvable as there is infinite number of control points possible so we need to cut down those to manageable size by exploiting something. For example we now the control points will not be to far from the original control points so we can use that to limit area for each point. You can use approximation search for this or any other minimization technique.

you can also get much much better start point for this if you use interpolation cubic. See BEZIER vs Interpolation cubics. So:
1. compute 4 new interpolation cubic control points from your BEZIER curve
  
  so if we have the same ranges as before then
```
X0 = BEZIER(t0-(t1-t0))
X1 = BEZIER(t0)
X2 = BEZIER(t1)
X3 = BEZIER(t1+(t1-t0))
```
  these are interpolation cubic control points not BEZIER. they should be uniformly dispersed. X1,X2 are the curve endpoints and X0,X3 are before and after them to preserve local shape and linearity of parameter as much as possible
2. transfrom (X0,X1,X2,X3) back into BEZIER control points (x0,x1,x2,x3)
  
  You can use mine formula from link above:
```
// input: X0,Y0,..X3,Y3  ... interpolation control points
// output: x0,y0,..x3,y3 ... Bezier control points
    double x0,y0,x1,y1,x2,y2,x3,y3,m=1.0/9.0;
    x0=X1;             y0=Y1;
    x1=X1+((X1-X0)*m); y1=Y1+((Y1-Y0)*m);
    x2=X2+((X2-X3)*m); y2=Y2+((Y2-Y3)*m);
    x3=X2;             y3=Y2;
```
  As you can see each axis is computed the same way ...
3. apply approximation search for BEZIER control points
  
  the new (x0,x1,x2,x3) are not precise yet as by changing control points blindly we could miss some local extreme of curve distorting shape slightly. So we need to search for the real ones. Luckily the new control points (x0',x1',x2',x3') will be very close to these points so now we have to search each point only near its counter part with some reasonable radius around (like bounding box size/8 or whatever ... you will see the results and can tweak ...

[notes]

If you need the exact precision result then only the #1 approach is usable. Approach #2 will have always some error. If the shape does not need to be exact sometimes the interpolation cubic converted to BEZIER without the final fitting will suffice (if shape not too complex in/near cutted areas).

As I wrote before without knowing which principle is used in your SubdivPoints is impossible to answer what the result of repetitive use of it will be...

Also you did not specify the constrains of solution like accuracy of result, speed/runtime restrictions ... if the ranges are fixed (constant) or can vary during runtime (that will extremly complicate the #1 approach as you need to handle ranges as variable till the end)

回答2:

Since your subdivideBezier() function did follow the De Casteljau algorithm, I am assuming that it works to subdivide a cubic Bezier curve at parameter t. So, yes, to continue subdividing at a different parameter (say t2), all you need to do is to figure out which subdivided curve t2 falls on, compute the new t2* against that curve and subdivide. In your example where you want to subdivide at t=0.2 and 0.6, the new parameter for 0.6 should be (0.6-0.2)/(1-0.2) = 0.5. So, you can simply subdivide the 2nd curve at t=0.5.