How to draw a bezier curve through its control points

Hi there,

Really stupid misunderstanding on my part coming up, I’m sure, but I was wondering how to draw a curve through its control points. I’m really quite a newbie when it comes to beziers and NURBs etc., but my requirement is simple and I’m sure that it has been done many, many times - I just haven’t found a good enough resource for silly old me on the net.

Here’s my code using a Bezier curve - I’m actually wondering if NURBs will get me there. Perhaps someone can modify the code below to make the curve pass through the control points…

  
	glMap1f(
			GL_MAP1_VERTEX_3, 
			0.0, 1.0, 
			7, theNumControlPoints, 
			&(theControlPoints[0].mX)
			);
	glEnable(GL_MAP1_VERTEX_3);
	
	glBegin(GL_LINE_STRIP);
		GLfloat theStep = 0.0;
		for (unsigned i = 0; i <= 25; ++i, theStep += .04) 
			glEvalCoord1f(theStep);
	glEnd();

Thanks in advance.

Cheers,
-C

Hi Huntc,

A Bezier curve will only interpolate its first and last control points, but will be constrained to the convex hull of all control points. Likewise, a NURBS curve can be made to interpolate control points by replicating knot vector values for particular control point (if first and last are replicated and all weights = 1, it reduces to a Bezier curve).

You may want to look at Catmull-Rom splines, if you need to interpolate all points rather than fit an approximating curve through them. While not offered by OpenGL, Catmull-Rom splines are not at all difficult to implement. Note that this class of curve does not have the convex hull property, and certain classes of shapes, like circles, can’t be modelled with this curve (but can be with NURBS).

There are many other curve types, but I hope that this helps to get you started.

Cheers,

Bill

Thanks Bill - I’ll take a look - got any sample GL code for this? A quick search on the net hasn’t yielded a full GL example yet tho’ I will also look a bit harder when I get more of a chance.

Thanks again for your help.

im using Catmull-Rom splines precisly for this,
dont know where i got the info to implemtn them though, perhaps here
http://astronomy.swin.edu.au/%7Epbourke/index.html

Hi guys,

That’s a great website. I gave it a quick scan and saw a bunch of good info on curves.

I went ahead and threw together a quick sample demonstrating the basics of Catmull-Rom. The first function uniformly subdivides the control points into time intervals, while the second function simply evaluates the polynomial. The gist of it is to calculate the curve points only when the control points change, so as to avoid unnecessary calculations.

The code is C#, so I hope it is recognizable. I haven’t had much time to test this today, but I hope that it helps.

	public void Draw()
	{
	    // Precalculate curve points
	    Vector3[] control = new Vector3[My number of control points...];
            for (int i = 0; i < control.Length; i++)
                control[i] = My positional data...;
            
	    int numPoints = Something reasonable...;
            float dt = (float)(control.Length - 1) / (numPoints - 1);
            Vector3[] points = new Vector3[numPoints];
            for (int i = 0; i < points.Length; i++)
            {
                float t = i * dt;
                points[i] = CatmullRom(t, control);
            }

	    // Loop and draw ...
	    while(true)
		DrawPoints(points);
	}

	public override Vector3 CatmullRom(float t, Vector[] control)
        {    
            // Get the index of the control point corresponding
            // to the global time parameter, t, 0 <= t <= numcontrol.
            if (t <= 0)
                return control[0];
            if (t >= control.Length)
                return control[control.Length - 1];

            // Statr of local time interval.
            int u0 = (int)t;

            // Calc the local u relative to start of interval
            float u = t - u0;
            
            // Grab the 4 control points for the interval, 
            // replicating the endpoints as necessary.
            Vector3[] cp = new Vector[4];
            for (int i = 0; i < 4; i++)
            {
                int index = u0 - 1 + i;
                if (index < 0)
                    index = 0;
                else if (index >= control.Length)
                    index = control.Length - 1;
                cp[i] = control[index];
            }

            // Evaluate the polynomial at u.
            return 0.5F * ((2 * cp[1])
                   + u * ((-cp[0] + cp[2])
                   + u * ((2 * cp[0] - 5 * cp[1] + 4 * cp[2] - cp[3])
                   + u * (-cp[0] + 3 * cp[1] - 3 * cp[2] + cp[3]))));
        }

Cheers,

Bill

Thanks Bill and Zed. I’ll digest this!

Cheers,
-C

OK then, I’ve now got this working from a C++ perspective. Thanks so much for the C# code Bill. Most helpful.

In case others are interested, here’s my code. I’ve tried to make it as GL-like as I can so that it is easily substituted for Bezier curves.

Firstly to call it:

  
CatmullRomSpline theSpline(
                            3, 4,
                            &theCtrlpoints[0][0]
                            );
	
glBegin(GL_LINE_STRIP);
	GLfloat theVertexP[3];
	GLfloat theStep = 0.0;
	for (unsigned i = 0; i <= 25; ++i, theStep += .04) 
        {
	    theSpline.Eval(theStep, theVertexP);
	    glVertex3fv(theVertexP);
        }
glEnd();

Here’s my C++ header:

  
#include <boost/shared_array.hpp>

class CatmullRomSpline
{
  public:
    // ctor/dtor
    
    CatmullRomSpline(unsigned long inStride, unsigned long inOrder, float* inPointsP) throw();
        // An array of a vertex expressed as x, y, z control points. The points array
        // will be copied by this object. Thus the caller is responsible
        // for any deallocation of the points array passed in.
        // Stride specifies the number of floats between the beginning
        // of  one  control point and the beginning of the next one in the
        // data structure referenced in inPointsP.
        // Order is the number of control points.
        
    // Accessors
    
    void Eval(float f, float* outVertexP) const throw();
        // For a given percentage, f (0.0 -> 1.0) return the vertex
  
  protected:
    // ctor/dtor

    CatmullRomSpline() throw();
  
  private:
    // Miscellaneous state
    
    struct Vertex
    {
        float x;
        float y;
        float z;
    };
    
    const unsigned long mOrder;
    boost::shared_array<Vertex> mPointsP;
};

  

…and finally here’s my implementation based on the work of Bill and others:

 
#include <limits>
#include <memory>

#include "CatmullRomSpline.h"

using namespace boost;
using namespace std;

CatmullRomSpline::CatmullRomSpline(unsigned long inStride, unsigned long inOrder, float* inPointsP) throw() :
    mOrder(inOrder)
{
    try
    {
        mPointsP = shared_array<Vertex>(new Vertex[mOrder]);
        
        // Copy in the points to our internal array
        
        Vertex* theNextPointP = mPointsP.get();
        const Vertex* theEndPointP = theNextPointP + mOrder;

        inStride -= 3;  // We'll be incrementing for 3 of the stride value
                        // so our stride can be decremented here

        while (theNextPointP < theEndPointP)
        {
            Vertex& thePoint = *theNextPointP++;

            thePoint.x = *inPointsP++;
            thePoint.y = *inPointsP++;
            thePoint.z = *inPointsP++;
            
            inPointsP += inStride;
        }
    }
    catch (const bad_alloc&)
    {}
}

void CatmullRomSpline::Eval(float f, float* outVertexP) const throw()
{
    // Determine which control point we're actually talking about
    
    const float t = f * mOrder;
    
    // Get the index of the mPointsP point corresponding
    // to the global time parameter
    
    if (t <= 0.0)           // Boundary condition
    {
        ::memmove(outVertexP, &(mPointsP[0].x), sizeof(Vertex));
    }
    else if (t >= mOrder)   // Boundary condition
    {
        ::memmove(outVertexP, &(mPointsP[mOrder - 1].x), sizeof(Vertex));
    }
    else                    // Common condition
    {
        // Start of local time interval.
        
        const unsigned u0 = unsigned(t);

        // Calc the local u relative to start of interval
        
        const float u = t - u0;

        // Grab the 4 mPointsP points for the interval 
        // replicating the endpoints as necessary
        
        const unsigned theIndexAtMinusOne = numeric_limits<unsigned>::max();
        
        Vertex* cpP[4];
        for (unsigned i = 0; i < 4; i++)
        {
            unsigned index = u0 - 1 + i;
            if (index == theIndexAtMinusOne) // Effectively == -1
                index = 0;
            else if (index >= mOrder)
                index = mOrder - 1;
            cpP[i] = mPointsP.get() + index;
        }

        const Vertex& cp0 = *cpP[0];
        const Vertex& cp1 = *cpP[1];
        const Vertex& cp2 = *cpP[2];
        const Vertex& cp3 = *cpP[3];

        // Evaluate the polynomial at u
        
        Vertex theVertex;
        
        theVertex.x = 
            0.5 * (
                (2 * cp1.x) +
                u * (
                    (-cp0.x + cp2.x) +
                    u * (
                        (2 * cp0.x - 5 * cp1.x + 4 * cp2.x - cp3.x) +
                        u * (-cp0.x + 3 * cp1.x - 3 * cp2.x + cp3.x)
                        )
                    )
                );
        
        theVertex.y = 
            0.5 * (
                (2 * cp1.y) +
                u * (
                    (-cp0.y + cp2.y) +
                    u * (
                        (2 * cp0.y - 5 * cp1.y + 4 * cp2.y - cp3.y) +
                        u * (-cp0.y + 3 * cp1.y - 3 * cp2.y + cp3.y)
                        )
                    )
                );
        
        theVertex.z = 
            0.5 * (
                (2 * cp1.z) +
                u * (
                    (-cp0.z + cp2.z) +
                    u * (
                        (2 * cp0.z - 5 * cp1.z + 4 * cp2.z - cp3.z) +
                        u * (-cp0.z + 3 * cp1.z - 3 * cp2.z + cp3.z)
                        )
                    )
                );

        ::memmove(outVertexP, &(theVertex.x), sizeof(Vertex));
    }
}
  

Enjoy.

Thanks again.

Cheers,
-C

One small improvement. In the case of:

catch (const bad_alloc&)
{}

resetting mOrder to zero (and removing its const declaration) is the correct thing to do.

One more small bug correction - the line in the implementation file that reads:

    const float t = f * mOrder;

should read:

    const float t = f * (mOrder - 1);

Overwise you’ll get the line drawing beyond what it should.

Cheers,
-C