Cubic Bezier Curve Flattening using Forward Differencing

/*
The MIT License (MIT)

Copyright (c) 2014-2015 Russel Lindsay

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
THE SOFTWARE.
*/

// flatten (approximate) a cubic bezier curve using a forward differencing algorithm. 
public void FlattenForwardDifferencing(
    Point p1,   // start point
    Point p2,   // control 1
    Point p3,   // control 2
    Point p4,   // end point
    int steps,  // how many segments to flatten to
    out Point[] points)
{
    points = new Point[steps + 1];

    /*
        * coefficients for control points
        * 
        * | 1  0  0  0|   |A|   |       A        |
        * |-3  3  0  0| * |B| = |   -3A + 3B     |
        * | 3 -6  3  0|   |C|   |  3A - 6B + 3C  |
        * |-1  3 -3  1|   |D|   |-A + 3B - 3C + D|
        *                     
        * http://www.wolframalpha.com/input/?i=%7B%7B1%2C+0%2C+0%2C+0%7D%2C%7B-3%2C+3%2C+0%2C+0%7D%2C+%7B3%2C+-6%2C+3%2C+0%7D%2C+%7B-1%2C+3%2C+-3%2C+1%7D%7D*%7B%7Bb_0%7D%2C+%7Bb_1%7D%2C+%7Bb_2%7D%2C+%7Bb_3%7D%7D
        * 
        */
    // same as notes above, but factored for speed
    float 
        cx = 3 * (p2.X - p1.X),
        cy = 3 * (p2.Y - p1.Y),

        bx = 3 * (p3.X - p2.X) - cx,
        by = 3 * (p3.Y - p2.Y) - cy,

        ax = p4.X - p1.X - cx - bx,
        ay = p4.Y - p1.Y - cy - by;

    /**
        * Formula for difference between t and t+h for cubic, quadratic, and linear polynomials
        * 
        * --------------------
        * Cubic polynomial (4 points, a, b, c, d)
        * C(t) = a*t^3 + b*t^2 + c*t + d
        * 
        * Cubic polynomial difference
        * C(t + h) - C(t)                                  http://www.wolframalpha.com/input/?i=%28a*%28t+%2B+h%29%5E3+%2B+b*%28t+%2B+h%29%5E2+%2B+c*%28t+%2B+h%29+%2B+d%29+-+%28a*t%5E3+%2B+b*t%5E2+%2B+c*t+%2B+d%29+
        * = ah^3 + 3ah^2 t + 3aht^2 + bh^2 + 2bht + ch     
     i)   = (3ah)t^2 + (3ah^2 + 2bh)t + ah^3 + bh^2 + ch     
        * 
        * --------------------
        * quadratic polynomial
        * Q(t) = at^2 + bt + c
        * 
        * quadratic polynomial difference
        * Q(t + h) - Q(t)                                  http://www.wolframalpha.com/input/?i=%28a%28t%2Bh%29%5E2+%2B+b%28t%2Bh%29+%2B+c%29+-+%28at%5E2+%2B+bt+%2B+c%29
    ii)   = (2ah)t + ah^2 + bh
        * 
        * ------------------
        * linear polynomial
        * L(t) = at + b
        * 
        * linear polynomial difference
        * L(t + h) - L(t)                                  http://www.wolframalpha.com/input/?i=%28a%28t+%2B+h%29+%2B+b%29+-+%28at+%2B+b%29
   iii)   = ah
        *
        * ----------------------------------------------------------------------------------------------------------------
        * 
        * Notice that the difference between t and t+h for a cubic polynomial is quadratic ((3ak)t^2 + (3ak^2 + 2bk)t + ak^3 + bk^2 + ck)
        * Likewise the difference in the quadratic case is linear, and in the linear case is a constant
        * 
        * The technique of forward differencing allows us to compute the polynomial at t=0,
        * and then incrementally calculate points by adding (and updating) the differences.
        * 
        * First Difference:
        * differance between t and t+h. Since the curve is cubic, difference is
        * i) (3ah)t^2 + (3ah^2 + 2bh)t + ah^3 + bh^2 + ch
        * Since we are calculating from start of curve, t = 0.
        * Substituting 
        *      t = 0 to calculate difference at start of curve,
    D1   = ah^3 + bh^2 + ch                                    http://www.wolframalpha.com/input/?i=%28a*%280+%2B+h%29%5E3+%2B+b*%280+%2B+h%29%5E2+%2B+c*%280+%2B+h%29+%2B+d%29+-+%28a*0%5E3+%2B+b*0%5E2+%2B+c*0+%2B+d%29
        * 
        * Second Difference:   
        * difference between two successive first differences. Since the form of a first difference is quadratic, the difference
        * between two of them is
        * ii) (2ah)t + ah^2 + bh
        * Substituting 
        *      t = 0 to calculate difference at start of curve,
        *      a = 3ah, 
        *      b = 3ah^2 + 2bh    (see D1 calculation)
        * we get
     D2   = (6ah^2)t + 6ah^3 + 2bh^2                            http://www.wolframalpha.com/input/?i=%282*3*a*h*h%29t+%2B+3*a*h*h%5E2+%2B+%283*a*h%5E2+%2B+2*b*h%29h
        * = 6ah^3 + 2bh^2
        * 
        * Third difference:
        * difference between two successive second differences. since the form is linear
        * iii) ah
        * substituting
        *      a = 6ah^2  (see D2 calculation)
     D3   = 6ah^3                                               http://www.wolframalpha.com/input/?i=6*a*h%5E2+*+h
        * 
        */

    // step size
    float 
        h = 1.0f / steps,
        hh = h * h,
        hhh = hh * h;

    // Note these multiplications can/should be optimised. 
    // They are left in this more readable form for testing purposes
    // First Difference, ah^3 + bh^2 + ch 
    float 
        d1x = ax * hhh + bx * hh + cx * h,
        d1y = ay * hhh + by * hh + cy * h,

        //Second Difference, 6ah^3 + 2bh^2
        d2x = 6 * ax * hhh + 2 * bx * hh,
        d2y = 6 * ay * hhh + 2 * by * hh,

        // Third Difference, 6ah^3
        d3x = 6 * ax * hhh,
        d3y = 6 * ay * hhh;

    // start from the first point (t=0)
    // and calculate each successive point, saving it to currentX, currentY
    float 
        currentX = p1.X,
        currentY = p1.Y;
    points[0].X = currentX;
    points[0].Y = currentY;

    // Forward Differencing
    // float errors will accumulate over time. Storing differences as double can help, 
    // but the minor loss in precision is probably not worth the performance loss 
    // due to the float casts that would be needed when assigning to result
    for (int i = 1; i < steps; i++) // our loop skips the first and last points. Since we already know them we dont need to calculate.
    {
        // calculate & save next point
        currentX += d1x;
        currentY += d1y;

        // increment differences
        d1x += d2x;
        d1y += d2y;

        d2x += d3x;
        d2y += d3y;

        points[i].X = currentX;
        points[i].Y = currentY;

        /*
            * Notes for anyone who would like to try and make an adaptive 
            * forward differencing algorithm.
            * You will need to change the step size according to a flattness criteria (that you implement yourself).
            * 
            * We can double or half the step size as follows:
            * 
            * Halving (given in the text)
            * ---------------------------
            * D3' = D3/8
            * D2' = D2/4 - D3'
            * D1' = (D1-D2') / 2
            * 
            * Doubling (inverse of above)
            * ---------------------------
            * D1' = 2*D1 + D2          http://www.wolframalpha.com/input/?i=inverse%28a+%3D+%28A+-+b%29+%2F+2%29
            * D2' = 4(D2 + D3)         http://www.wolframalpha.com/input/?i=inverse%28b+%3D+B%2F4+-+c%29
            * D3' = 8D3                http://www.wolframalpha.com/input/?i=inverse%28c+%3D+C+%2F+8%29
            * 
            */

    }

    // no need to calcualte the last point, we already know it. This is also more accurate since we avoid accumulation errors.
    points[steps].X = p4.X;
    points[steps].Y = p4.Y;
}