Bézier curves, Loop and Blinn style

Question

A couple of days ago I started looking into efficiently drawing bezier curves, and I came across a this method developed by Charles Loop and Jim Blinn that seemed very interesting. How ever, after a lot of experimenting with their algorithm, I just can't seem to get it to be able to render the cubic curves. Quadratics are fine, no problem there.

The only resources I have found so far are as follows:

GPU Gems 3 Chapter 25

Curvy Blues

Resolution Independent Curve Rendering using Programmable Graphics Hardware

To get the testing up and running quickly, I'm doing this in XNA. Basically I'm passing texture coordinates with my vertices to the GPU, apply a perspective transform and use the formula mentioned in all the articles in a pixel shader to render the final result. How ever, the problem (I think) lays in how I calculate the texture coordinates. Check this code out:

public void Update()
{
    float a1 = Vector3.Dot(p1, Vector3.Cross(p4, p3));
    float a2 = Vector3.Dot(p2, Vector3.Cross(p1, p4));
    float a3 = Vector3.Dot(p3, Vector3.Cross(p2, p2));

    float d1 = a1 - 2 * a2 + 3 * a3;
    float d2 = -a2 + 3 * a3;
    float d3 = 3 * a3;

    float discr = d1 * d1 * (3 * d2 * d2 - 4 * d1 * d3);

    if (discr > 0)
    {
        Type = CurveTypes.Serpentine;

        float ls = 3 * d2 - (float)Math.Sqrt(9 * d2 * d2 - 12 * d1 * d3);
        float lt = 6 * d1;
        float ms = 3 * d2 + (float)Math.Sqrt(9 * d2 * d2 - 12 * d1 * d3);
        float mt = 6 * d1;

        TexCoord1 = new Vector3(ls * ms, (float)Math.Pow(ls, 3), (float)Math.Pow(ms, 3));
        TexCoord2 = new Vector3((3 * ls * ms - ls * mt - lt * ms) / 3, ls * ls * (ls - lt), ms * ms * (ms - mt));
        TexCoord3 = new Vector3((lt * (mt - 2 * ms) + ls * (3 * ms - 2 * mt)) / 3, (float)Math.Pow(lt - ls, 2) * ls, (float)Math.Pow(mt - ms, 2) * ms);
        TexCoord4 = new Vector3((lt - ls) * (mt - ms), -(float)Math.Pow(lt - ls, 3), -(float)Math.Pow(mt - ms, 3));
    }
    else if (discr == 0)
    {
        Type = CurveTypes.Cusp;
    }
    else if (discr < 0)
    {
        Type = CurveTypes.Loop;
    }
}

Excuse the mess, it's just some testing code. p1...p4 are the control points in world space, and TexCoord1...TexCoord4 are the corresponding texture coordinates. This is a replication of what is said in the GPU Gems article.

There are a few problems here, first when calculating a3, we use p2 for both parameters, which of course always results in a (0,0,0) vector, and taking the dot product of that and p3 will always give us 0. That't pretty useless, so why would they mention that in the article?

This will of course make discr incorrect, and we won't even be able to determine what type of curve it is.

After fiddling around with that code for a while, I decided to try to do it exactly the why they did in in the Loop and Blinn paper. From that I get something like this:

public void Update()
{
    Matrix m1 = new Matrix(
        p4.X, p4.Y, 1, 0,
        p3.X, p3.Y, 1, 0,
        p2.X, p2.Y, 1, 0,
        0, 0, 0, 1);
    Matrix m2 = new Matrix(
        p4.X, p4.Y, 1, 0,
        p3.X, p3.Y, 1, 0,
        p1.X, p1.Y, 1, 0,
        0, 0, 0, 1);
    Matrix m3 = new Matrix(
        p4.X, p4.Y, 1, 0,
        p2.X, p2.Y, 1, 0,
        p1.X, p1.Y, 1, 0,
        0, 0, 0, 1);
    Matrix m4 = new Matrix(
        p3.X, p3.Y, 1, 0,
        p2.X, p2.Y, 1, 0,
        p1.X, p1.Y, 1, 0,
        0, 0, 0, 1);

    float det1 = m1.Determinant();
    float det2 = -m2.Determinant();
    float det3 = m3.Determinant();
    float det4 = -m4.Determinant();

    float tet1 = det1 * det3 - det2 * det2;
    float tet2 = det2 * det3 - det1 * det4;
    float tet3 = det2 * det4 - det3 * det3;

    float discr = 4 * tet1 * tet3 - tet2 * tet2;

    if (discr > 0)
    {
        Type = CurveTypes.Serpentine;

        float ls = 2 * det2;
        float lt = det3 + (float)((1 / Math.Sqrt(3)) * Math.Sqrt(3 * det3 * det3 - 4 * det2 * det4));
        float ms = 2 * det2;
        float mt = det3 - (float)((1 / Math.Sqrt(3)) * Math.Sqrt(3 * det3 * det3 - 4 * det2 * det4));

        TexCoord1 = new Vector3(lt * mt, (float)Math.Pow(lt, 3), (float)Math.Pow(mt, 3));
        TexCoord2 = new Vector3(-ms * lt - ls * mt, -3 * ls * lt * lt, -3 * ms * mt * mt);
        TexCoord3 = new Vector3(ls * ms, 3 * ls * ls * lt, 3 * ms * ms * mt);
        TexCoord4 = new Vector3(0, -ls * ls * ls, -ms * ms * ms);
    }
    else if (discr == 0)
    {
        Type = CurveTypes.Cusp;
    }
    else if (discr < 0)
    {
        Type = CurveTypes.Loop;
    }
}

Guess what, that didn't work either. How ever, discr seem to be at least a little more correct now. At least it has the right sign, and it is zero when the control points are arranged to form a cusp. I still get the same visual result though, except the curve disappears randomly for a while (the pixel shader formula is always greater than zero) and returns after I move the control point back to more of a square shape. Here is the pixel shader code by the way:

 PixelToFrame PixelShader(VertexToPixel PSIn)
 {
     PixelToFrame Output = (PixelToFrame)0;

     if(pow(PSIn.TexCoords.x, 3) - PSIn.TexCoords.y * PSIn.TexCoords.z > 0)
     {
     Output.Color = float4(0,0,0,0.1);
     }
     else
     {
     Output.Color = float4(0,1,0,1);
     }

    return Output;
}

That is about all useful information I can think of right now. Does anyone have any idea what's going on? Because I'm running out of them.

Answer 1

I was looking at the paper and your code, and it seams you're missing the multiplication to the M3 matrix.

Your p1, p2, p3 and p4 coordinates should be placed in a matrix and multiplied by the M3 matrix, before using it to calculate the determinants. eg.

Matrix M3 = Matrix(
    1, 0, 0, 0,
    -3, 3, 0, 0,
    3, -6, 3, 0,
    -1, 3, -3, 1);
Matrix B = Matrix(
    p1.X, p1.Y, 0, 1,
    p2.X, p2.Y, 0, 1,
    p3.X, p3.Y, 0, 1,
    p4.X, p4.Y, 0, 1);
Matrix C = M3*B;

Then you use each row of the C matrix as the coordinates for the m1 to m4 matrices in your code. Where the first and second values of the row are the x,y coordinates and the last is the w coordinate.

Finally the matrix of texture coordinates needs to be mutiplied by the inverse of M3 eg.

Matrix invM3 = Matrix(
    1, 0, 0, 0,
    1, 0.3333333, 0, 0,
    1, 0.6666667, 0.333333, 0,
    1, 1, 1, 1);
Matrix F = Matrix(
    TexCoord1,
    TexCoord2,
    TexCoord3,
    TexCoord4);
Matrix result = invM3*F;

Each row of the resulting matrix corresponds to the texture coordinates needed for the shader.

I haven't implemented it myself yet, so cannot guarantee that that it will solve your problem. It simply is what I noticed to be missing from your implementation after reading the paper.

I hope this helps, if I'm wrong please tell me cause I will be trying this out soon.