B-Spline Surface Implementation with OpenGL

In this implementation, instead of using any B-Spline built-in functions of OpenGL, I made my own functions calculate points on surface base on given knots and control points.

OK, let’s start!

Basic Information

First of all, I want to specify some notations in this implementation.

Here are two core expression of B-Spline curve.

and here is the linear combination of B-spline surface. They use the same expressions for the polynomial pieces ( N() function).

If you do not understand these expressions clearly yet, I recommend you read this tutorial below before going on. Knowledge in Unit 6 and 8 is enough for this implementation.

<CS3621 Introduction to Computing with Geometry Notes>

<You can download the complete runnable code and input file at the end of this page.>

Input Data (knots vector and control points)

Here is an example to help you have intuition view on input data. Since we are making a surface, we have U and V two directions.


6 6 //the number of control points in U and V directions.
2 2 //the degree of U and V direction
0 0 0 0.2 0.5 0.7 1 1 1 //knots in U direction
0 0 0 0.2 0.5 0.7 1 1 1 //knots in V direction
-30 -36 -30 //first control point’s 3D coordinates in U and V direction
-30 -26 -20 // secondcontrol point
....... //Other 36 control points in this example. control points’ coordinates in row-major order. The row is V direction, the column is U direction

From this example, we have these variables.


numContrU = 6, numContrV = 6 //numbers of control points

degreeU = 2, degreeV = 2 // degree
knotU[numContrU + degreeU + 1] = 0 0 0 0.2 0.5 0.7 1 1 1
knotV[numContrV + degreeV + 1] = 0 0 0 0.2 0.5 0.7 1 1 1


The introduction of the code for reading input data will be omitted here.

B-Spline Function (core code)

First, for the B-Spline basic functions(Cox-de Boor recursion formula)

I wrote a recursive function deBoor()


float deBoor(int i, int p, float u, float *knots){
float equationR, equationL;//right coefficient and left coefficient

if (p == 0){
if (u >= knots[i] && (u < knots[i + 1]))
return 1;
else return 0;
}
else{
if ((knots[i + p] - knots[i]) == 0) equationR = 0; //Prevent divide by zero
else equationR = ((u - knots[i]) / (knots[i + p] - knots[i]));
if ((knots[i + p + 1] - knots[i + 1]) == 0) equationL = 0; // Prevent divide by zero
else equationL = ((knots[i + p + 1] - u) / (knots[i + p + 1] - knots[i + 1]));
return equationR * deBoor(i, p - 1, u, knots) + equationL * deBoor(i + 1, p - 1, u, knots);
}
}


After calculating the deBoor(), it’s the time to get the points on the surface by Control Points. I wrote a function to do this job.


void bSpline(){

for (int u = 0; u < numSegment; u++) //for each point on the surface, 20*20 points for this surface
for (int v = 0; v < numSegment; v++){
for (int i = 0; i <= numContrU; i++)//for each control points
for (int j = 0; j <= numContrV; j++){
P[u][v] = P[u][v] + ctlPoints[i][j] * double(deBoor(i, degreeU, u*piece, knotU) * deBoor(j, degreeV, v*piece, knotV));
}
}
}


Normal Vectors

I wrote a function(vertexNormal) calculate the normal vector of each vertex(smooth shading) and face (flat shading) on the surface.


Vec3 faceNormal(Vec3 a, Vec3 b, Vec3 c){
Vec3 ab = a - b;
Vec3 ac = a - c;
Vec3 cross = ab.crossProduct(ac);
return cross.normalize();
}

void vertexNormal(){
for (int i = 0; i < numSegment-1; i++)
for (int j = 0; j < numSegment-1; j++){
facenormal[i][j] = faceNormal(P[i][j + 1], P[i + 1][j], P[i][j]);
}
//four corners
vertexnormal[0][0] = facenormal[0][0];
vertexnormal[numSegment][numSegment] = facenormal[numSegment - 1][numSegment - 1];
vertexnormal[0][numSegment] = facenormal[0][numSegment - 1];
vertexnormal[numSegment][0] = facenormal[numSegment - 1][0];
//four edge
for (int i = 1; i < numSegment; i++){
vertexnormal[i][0] = (facenormal[i - 1][0] + facenormal[i][0]) / 2;
vertexnormal[0][i] = (facenormal[0][i-1] + facenormal[0][i]) / 2;
vertexnormal[i][19] = (facenormal[i - 1][18] + facenormal[i][18]) / 2;
vertexnormal[19][i] = (facenormal[18][i-1] + facenormal[18][i]) / 2;
}
//others
for (int i = 1; i < numSegment-1; i++)
for (int j = 1; j < numSegment-1; j++){
vertexnormal[i][j] = (facenormal[i - 1][j - 1] + facenormal[i - 1][j] + facenormal[i][j - 1] + facenormal[i][j]) / 4;
}
}


Render the Sence

I used the GL_TRIANGLES to generate the surface in two models(smooth shading and flat shading).

If you use the GL_LIGHT_MODEL_TWO_SIDE model  for lighting, you have to be careful to arrange the order of three points of each triangle and their normal vectors.

Smooth Shading

Each loop generates two triangles, all these triangles will be connected together.


glShadeModel(GL_SMOOTH);
glBegin(GL_TRIANGLES);
for (int i = 0; i < numSegment-1; i++){
for (int j = 0; j < numSegment-1; j++){
glNormal3f(vertexnormal[i][j + 1].x, vertexnormal[i][j + 1].y, vertexnormal[i][j + 1].z);
glVertex3f(P[i][j + 1].x, P[i][j + 1].y, P[i][j + 1].z);
glNormal3f(vertexnormal[i + 1][j].x, vertexnormal[i + 1][j].y, vertexnormal[i + 1][j].z);
glVertex3f(P[i + 1][j].x, P[i + 1][j].y, P[i + 1][j].z);
glNormal3f(vertexnormal[i][j].x, vertexnormal[i][j].y, vertexnormal[i][j].z);
glVertex3f(P[i][j].x, P[i][j].y, P[i][j].z);

glNormal3f(vertexnormal[i + 1][j].x, vertexnormal[i + 1][j].y, vertexnormal[i + 1][j].z);
glVertex3f(P[i + 1][j].x, P[i + 1][j].y, P[i + 1][j].z);
glNormal3f(vertexnormal[i][j + 1].x, vertexnormal[i][j + 1].y, vertexnormal[i][j + 1].z);
glVertex3f(P[i][j + 1].x, P[i][j + 1].y, P[i][j + 1].z);
glNormal3f(vertexnormal[i + 1][j + 1].x, vertexnormal[i + 1][j + 1].y, vertexnormal[i + 1][j + 1].z);
glVertex3f(P[i + 1][j + 1].x, P[i + 1][j + 1].y, P[i + 1][j + 1].z);
}
}
glEnd();


Flat Shading

Different from the smooth shading, there is only one normal vector for each triangle.


Vec3 facenormal1, facenormal2;
glShadeModel(GL_FLAT);
glBegin(GL_TRIANGLES);
for (int i = 0; i < numSegment-1; i++){
for (int j = 0; j < numSegment-1; j++){
facenormal1 = faceNormal(P[i][j + 1], P[i + 1][j], P[i][j]);
facenormal2 = faceNormal(P[i + 1][j + 1], P[i + 1][j], P[i][j + 1]);

glNormal3f(facenormal1.x, facenormal1.y, facenormal1.z);
glVertex3f(P[i][j + 1].x, P[i][j + 1].y, P[i][j + 1].z);
glVertex3f(P[i + 1][j].x, P[i + 1][j].y, P[i + 1][j].z);
glVertex3f(P[i][j].x, P[i][j].y, P[i][j].z);

glNormal3f(facenormal2.x, facenormal2.y, facenormal2.z);
glVertex3f(P[i + 1][j].x, P[i + 1][j].y, P[i + 1][j].z);
glVertex3f(P[i][j + 1].x, P[i][j + 1].y, P[i][j + 1].z);
glVertex3f(P[i+1][j+1].x, P[i+1][j+1].y, P[i+1][j+1].z);
}
}
glEnd()


Now, you have a B-Spline surface!

I added some small features for this surface.

Other Features

Keyboard Contol


void processSpecialKeys(int key, int x, int y)
{
switch (key) {
case GLUT_KEY_F2: //Switch flat and smooth shading
flatSmooth = !flatSmooth;
break;
case GLUT_KEY_F3: //turn on/off the control polygon
showControlPolygon = !showControlPolygon;
break;
case GLUT_KEY_F4: //ture on/off the wire frame
showWireframe = !showWireframe;
glutPostRedisplay();
break;
//change the visualization view, there is some code works together in the renderScenc function.
case GLUT_KEY_LEFT:
yAngle += 5;
if (yAngle == 360) yAngle = 0;
break;
case GLUT_KEY_RIGHT:
yAngle -= 5;
if (yAngle == -360) yAngle = 0;
break;
case GLUT_KEY_UP:
xAngle += 5;
if (xAngle == 360) xAngle = 0;
break;
case GLUT_KEY_DOWN:
xAngle -= 5;
if (xAngle == -360) xAngle = 0;
break;
}
glutPostRedisplay();
}


Flat and Smooth Shading


if (flatSmooth){
/************ Smooth shading ****************/
//you can find this part of code above
}
else {
/**************Flat shading****************
//you can find this part of code above
}


Control Polygon


if (showControlPolygon) {
glDisable(GL_LIGHTING);
glColor3f(1.0, 1.0, 0.0);
for (int i = 0; i <= numContrU; i++) {
glBegin(GL_LINE_STRIP);
for (int j = 0; j <= numContrV; j++) {
glVertex3f(ctlPoints[i][j].x, ctlPoints[i][j].y, ctlPoints[i][j].z);
}
glEnd();
}
for (int i = 0; i <= numContrV; i++) {
glBegin(GL_LINE_STRIP);
for (int j = 0; j <= numContrU; j++) {
glVertex3f(ctlPoints[j][i].x, ctlPoints[j][i].y, ctlPoints[j][i].z);
}
glEnd();
}
glEnable(GL_LIGHTING);
}


Control Polygon


if (showWireframe) {
glDisable(GL_LIGHTING);
glColor3f(1.0, 1.0, 1.0);
for (int i = 0; i < numSegment; i++) {
glBegin(GL_LINE_STRIP);
for (int j = 0; j < numSegment; j++) {
glVertex3f(P[i][j].x, P[i][j].y, P[i][j].z);
}
glEnd();
glBegin(GL_LINE_STRIP);
for (int j = 0; j < numSegment; j++) {
glVertex3f(P[j][i].x, P[j][i].y, P[j][i].z);
}
glEnd();
}
glEnable(GL_LIGHTING);
}


Control Polygon


//Global declaration
float xAngle = 0, yAngle = 0;
//In the renderScence function
glRotatef(xAngle, 1.0f, 0.0f, 0.0f);
glRotatef(yAngle, 0.0f, 1.0f, 0.0f);
//Kyeboard contol
//In the kyeboard control function


Ok! That’s all!

Code and input file download

Thanks for reading. If you have any questions, please comment this articles.