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();

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.

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.

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]))));
}

#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));
}
}