Algorithm to add Color in Bezier curves

Question

I'm playing with GD library for a while and more particuraly with Bezier curves atm.

I used some existant class which I modified a little (seriously eval()...). I found out it was a generic algorithm used in and convert for GD.

Now I want to take it to another level: I want some colors.

No problem for line color but with fill color it's harder.

My question is:

Is there any existant algorithm for that? I mean mathematical algorithm or any language doing it already so that I could transfer it to PHP + GD?

EDIT2 So, I tried @MizardX solution with a harder curve :

• 1st position : 50 - 50
• final position : 50 - 200
• 1st control point : 300 - 225
• 2nd control point : 300 - 25

Which should show this :

And gives this :

EDIT I already read about @MizardX solution. Using imagefilledpolygon to make it works.

But it doesn't work as expected. See the image below to see the problem. Top graph is what I expect (w/o the blackline for now, only the red part).

Coordinates used:

• first point is 100 - 100
• final point is 300 - 100
• first control point is 100 - 0
• final control point is 300 - 200

Bottom part is what I get with that kind of algorithm...

Answer 1

Convert the Bezier curve to a polyline/polygon, and fill that. If you evaluate the Bezier polynomial at close enough intervals (~1 pixel) it will be identical to an ideal Bezier curve.

I don't know how familiar you are with Bezier curves, but here is a crash course:

<?php

// Calculate the coordinate of the Bezier curve at $t = 0..1
function Bezier_eval($p1,$p2,$p3,$p4,$t) {
    // lines between successive pairs of points (degree 1)
    $q1  = array((1-$t) * $p1[0] + $t * $p2[0],(1-$t) * $p1[1] + $t * $p2[1]);
    $q2  = array((1-$t) * $p2[0] + $t * $p3[0],(1-$t) * $p2[1] + $t * $p3[1]);
    $q3  = array((1-$t) * $p3[0] + $t * $p4[0],(1-$t) * $p3[1] + $t * $p4[1]);
    // curves between successive pairs of lines. (degree 2)
    $r1  = array((1-$t) * $q1[0] + $t * $q2[0],(1-$t) * $q1[1] + $t * $q2[1]);
    $r2  = array((1-$t) * $q2[0] + $t * $q3[0],(1-$t) * $q2[1] + $t * $q3[1]);
    // final curve between the two 2-degree curves. (degree 3)
    return array((1-$t) * $r1[0] + $t * $r2[0],(1-$t) * $r1[1] + $t * $r2[1]);
}

// Calculate the squared distance between two points
function Point_distance2($p1,$p2) {
    $dx = $p2[0] - $p1[0];
    $dy = $p2[1] - $p1[1];
    return $dx * $dx + $dy * $dy;
}

// Convert the curve to a polyline
function Bezier_convert($p1,$p2,$p3,$p4,$tolerance) {
    $t1 = 0.0;
    $prev = $p1;
    $t2 = 0.1;
    $tol2 = $tolerance * $tolerance;
    $result []= $prev[0];
    $result []= $prev[1];
    while ($t1 < 1.0) {
        if ($t2 > 1.0) {
            $t2 = 1.0;
        }
        $next = Bezier_eval($p1,$p2,$p3,$p4,$t2);
        $dist = Point_distance2($prev,$next);
        while ($dist > $tol2) {
            // Halve the distance until small enough
            $t2 = $t1 + ($t2 - $t1) * 0.5;
            $next = Bezier_eval($p1,$p2,$p3,$p4,$t2);
            $dist = Point_distance2($prev,$next);
        }
        // the image*polygon functions expect a flattened array of coordiantes
        $result []= $next[0];
        $result []= $next[1];
        $t1 = $t2;
        $prev = $next;
        $t2 = $t1 + 0.1;
    }
    return $result;
}

// Draw a Bezier curve on an image
function Bezier_drawfilled($image,$p1,$p2,$p3,$p4,$color) {
    $polygon = Bezier_convert($p1,$p2,$p3,$p4,1.0);
    imagefilledpolygon($image,$polygon,count($polygon)/2,$color);
}

?>

Edit:

I forgot to test the routine. It is indeed as you said; It doesn't give a correct result. Now I have fixed two bugs:

1. I unintentionally re-used the variable names $p1 and $p2. I renamed them $prev and $next.
2. Wrong sign in the while-loop. Now it loops until the distance is small enough, instead of big enough.

Answer 2

I checked the algorithm for generating a Polygon ensuring a bounded distance between successive parameter-generated points, and seems to work well for all the curves I tested.

Code in Mathematica:

pts={{50,50},{300,225},{300,25},{50,200}};
f=BezierFunction[pts];
step=.1; (*initial step*)

While[ (*get the final step - Points no more than .01 appart*)
   Max[
     EuclideanDistance @@@ 
         Partition[Table[f[t],{t,0,1,step}],2,1]] > .01,
   step=step/2]

(*plot it*)
Graphics@Polygon@Table[f[t],{t,0,1,step}]  

.

.

The algorithm could be optimized (ie. generate less points) if you don't require the same parameter increment between points, meaning you can chose a parameter increment at each point that ensures a bounded distance to the next.

Random examples:

Answer 3

Generate a list of successive points which lie along the curve (p_list)).

You create a line between the two end points of the curve (l1).

Then you are going to find the normal of the line (n1). Using this normal find the distance between the two furthest points (p_max1, and p_max2) along this normal (d1). Divide this distance into n discrete units (delta).

Now shift l1 along n1 by delta, and solve for the points of intersection (start with brute force and check for a solution between all the line segments in p_list). You should be able to get two points of intersection for each shift of l1, excepting boundaries and self intersection where you may have only have a single point. Hopefully the quad routine can have two points of the quad be at the same location (a triangle) and fill without complaint otherwise you'll need triangles in this case.

Sorry I didn't provide pseudo code but the idea is pretty simple. It's just like taking the two end points and joining them with a ruler and then keeping that ruler parallel to the original line start at one end and with successive very close pencil marks fill in the whole figure. You'll see that when you create your little pencil mark (a fine rectangle) that the rectangle it highly unlikely to use the points on the curve. Even if you force it to use a point on one side of the curve it would be quite the coincidence for it to exactly match a point on the other side, for this reason it is better to just calculate new points. At the time of calculating new points it would probably be a good idea to regenerate the curves p_list in terms of these points so you can fill it more quickly (if the curve is to stay static of course otherwise it wouldn't make any sense).

Answer 4

This answer is very similar to @MizardX's, but uses a different method to find suitable points along the Bezier for a polygonal approximation.

function split_cubic($p, $t)
{
    $a_x = $p[0] + ($t * ($p[2] - $p[0]));
    $a_y = $p[1] + ($t * ($p[3] - $p[1]));
    $b_x = $p[2] + ($t * ($p[4] - $p[2]));
    $b_y = $p[3] + ($t * ($p[5] - $p[3]));
    $c_x = $p[4] + ($t * ($p[6] - $p[4]));
    $c_y = $p[5] + ($t * ($p[7] - $p[5]));
    $d_x = $a_x + ($t * ($b_x - $a_x));
    $d_y = $a_y + ($t * ($b_y - $a_y));
    $e_x = $b_x + ($t * ($c_x - $b_x));
    $e_y = $b_y + ($t * ($c_y - $b_y));
    $f_x = $d_x + ($t * ($e_x - $d_x));
    $f_y = $d_y + ($t * ($e_y - $d_y));

    return array(
        array($p[0], $p[1], $a_x, $a_y, $d_x, $d_y, $f_x, $f_y),
        array($f_x, $f_y, $e_x, $e_y, $c_x, $c_y, $p[6], $p[7]));
}

$flatness_sq = 0.25; /* flatness = 0.5 */
function cubic_ok($p)
{
    global $flatness_sq;

    /* test is essentially:
     * perpendicular distance of control points from line < flatness */

    $a_x = $p[6] - $p[0];  $a_y = $p[7] - $p[1];
    $b_x = $p[2] - $p[0];  $b_y = $p[3] - $p[1];
    $c_x = $p[4] - $p[6];  $c_y = $p[5] - $p[7];
    $a_cross_b = ($a_x * $b_y) - ($a_y * $b_x);
    $a_cross_c = ($a_x * $c_y) - ($a_y * $c_x);
    $d_sq = ($a_x * $a_x) + ($a_y * $a_y);
    return max($a_cross_b * $a_cross_b, $a_cross_c * $a_cross_c) < ($flatness_sq * $d_sq);
}

$max_level = 8;
function subdivide_cubic($p, $level)
{
    global $max_level;

    if (($level == $max_level) || cubic_ok($p)) {
        return array();
    }

    list($q, $r) = split_cubic($p, 0.5);
    $v = subdivide_cubic($q, $level + 1);
    $v[] = $r[0]; /* add a point where we split the cubic */
    $v[] = $r[1];
    $v = array_merge($v, subdivide_cubic($r, $level + 1));
    return $v;
}

function get_cubic_points($p)
{
    $v[] = $p[0];
    $v[] = $p[1];
    $v = array_merge($v, subdivide_cubic($p, 0));
    $v[] = $p[6];
    $v[] = $p[7];
    return $v;
}

function imagefilledcubic($img, $p, $color)
{
    $v = get_cubic_points($p);
    imagefilledpolygon($img, $v, count($v) / 2, $color);
}

The basic idea is to recursively split the cubic in half until the bits we're left with are almost flat. Everywhere we split the cubic, we stick a polygon point.

split_cubic splits the cubic in two at parameter $t. cubic_ok is the "are we flat enough?" test. subdivide_cubic is the recursive function. Note that we stick a limit on the recursion depth to avoid nasty cases really screwing us up.

Your self-intersecting test case:

$img = imagecreatetruecolor(256, 256);

imagefilledcubic($img, array(
    50.0, 50.0, /* first point */
    300.0, 225.0, /* first control point */
    300.0, 25.0, /* second control point */
    50.0, 200.0), /* last point */
    imagecolorallocate($img, 255, 255, 255));

imagepng($img, 'out.png');

imagedestroy($img);

Gives this output:

I can't figure out how to make PHP nicely anti-alias this; imageantialias($img, TRUE); didn't seem to work.