# Extend the cubic spline form 2D to 3D

I have here a cubic spline code written in 2d space can anyone show me in details how to extend this code to be 3D by adding z direction and draw the lines using gluCylinder instead of lines.

``````/*
*
* Algorithm derived from http://blog.ivank.net/interpolation-with-cubic-splines.html
*/

#include <Windows.h>
#include <iostream>
#include <vector>
#include <GL/glut.h>

#define WIDTH 400
#define HEIGHT 400
#define POINT_SIZE 10.0f

// Typedefs to make it easier to refer to our matrices
typedef std::vector<double> MatrixRow;
typedef std::vector<MatrixRow> Matrix;

// Struct for representing a point in 2D space
struct Point {
double x;
double y;
};

// Variables for dragging points arround
bool dragging = false;
Point lastPoint = Point{ -1, -1 };

// If the next click should delete a point instead of dragging or adding one
bool shouldDelete = false;

// The number of interpolation steps to do between points
int steps = 5;

// The list of spline points, sorted from lowest x value to highest
std::vector<Point> points;

{
INCREASE_STEPS,
DECREASE_STEPS,
QUIT
};

/**
* Helper function to a point while still maintaining the vector's sorting
* Also disallows two points having the same x coordinate
*/
void insertPoint(Point p)
{
// If the size is 0 just go ahead and insert
if (points.size() == 0) {
points.push_back(p);
}
else {
bool inserted = false;
for (auto it = points.begin(); it != points.end(); ++it) {
Point cur = *it;
if (cur.x == p.x) {
// We can't have overlapping x's for the calculation
return;
}
if (cur.x > p.x) {
points.insert(it, p);
inserted = true;
break;
}
}

if (!inserted) {
points.push_back(p);
}
}
}

/**************************************************************
* Spline functions
**************************************************************/

/**
* Solve the given Matrix and return the k values
*/
std::vector<double> solveMatrix(Matrix A)
{
std::vector<double> ks(points.size());

int m = A.size();

// Iterate through columns
for (int k = 0; k < m; k++) {
// Pivot for column
int i_max = 0;
int val_max = MININT;

for (int i = k; i < m; i++) {
if (A[i][k] > val_max) {
i_max = i;
val_max = A[i][k];
}
}

// Swap row k and the max row
MatrixRow temp = A[k];
A[k] = A[i_max];
A[i_max] = temp;

// Iterate on all the rows before the pivot
for (int i = k + 1; i < m; i++) {
for (int j = k + 1; j < m + 1; j++) {
// Make sure there's no division by 0
if (A[k][k] == 0) {
std::cerr << "Division by 0, can't calculate k values" << std::endl;
return ks;
}
A[i][j] = A[i][j] - A[k][j] * (A[i][k] / A[k][k]);
}

A[i][k] = 0;
}
}

// Iterate through columns
for (int i = m - 1; i >= 0; i--) {
// Make sure there's no division by 0
if (A[i][i] == 0) {
std::cerr << "Division by 0, can't calculate k values" << std::endl;
return ks;
}

double v = A[i][m] / A[i][i];
ks[i] = v;

// Iterate over rows
for (int j = i - 1; j >= 0; j--) {
A[j][m] -= A[j][i] * v;
A[j][i] = 0;
}
}

return ks;
}

/**
* Calculates the k values for the points vector and returns them
*/
std::vector<double> getKs()
{
int n = points.size() - 1;
std::vector<double> ks;

// The n + 1 by n + 2 matrix that we will solve, represented as nested vectors
Matrix A(n + 1);
for (int i = 0; i < A.size(); i++) {
// Initialize the rows to be n + 2 long
A[i] = MatrixRow(n + 2);
}

// Iterate through the rows (skipping the first one)
for (int i = 1; i < n; i++) {
double prevXDiff = points[i].x - points[i - 1].x;
double nextXDiff = points[i + 1].x - points[i].x;
double prevYDiff = points[i].y - points[i - 1].y;
double nextYDiff = points[i + 1].y - points[i].y;

// If either diff is 0, we can't continue because division would fail
if (prevXDiff == 0 || nextXDiff == 0) {
std::cerr << "Two points are directly above each other at x=" << points[i].x
<< ", unable to calculate k values" << std::endl;
return ks; // Return exmpty vector
}

// Now calculate the columns
A[i][i - 1] = 1 / prevXDiff;
A[i][i] = 2 * (1 / prevXDiff + 1 / nextXDiff);
A[i][i + 1] = 1 / nextXDiff;
A[i][n + 1] = 3 * ((prevYDiff / (prevXDiff * prevXDiff)) + (nextYDiff / (nextXDiff * nextXDiff)));
}

double xdiff = points[1].x - points[0].x;
double ydiff = points[1].y - points[0].y;

// If xdiff is 0, we can't continue because division would fail
if (xdiff == 0) {
std::cerr << "Two points are directly above each other at x=" << points[0].x
<< ", unable to calculate k values" << std::endl;
return ks; // Return empty vector
}

// Calculate the first row
A[0][0] = 2 / xdiff;
A[0][1] = 1 / xdiff;
A[0][n + 1] = 3 * ydiff / (xdiff * xdiff);

xdiff = points[n].x - points[n - 1].x;
ydiff = points[n].y - points[n - 1].y;

// If xdiff is 0, we can't continue because division would fail
if (xdiff == 0) {
std::cerr << "Two points are directly above each other at x=" << points[0].x
<< ", unable to calculate k values";
return ks; // Return empty vector
}

// Calculate the last row
A[n][n - 1] = 1 / xdiff;
A[n][n] = 2 / xdiff;
A[n][n + 1] = 3 * ydiff / (xdiff * xdiff);

ks = solveMatrix(A);

return ks;
}

/**
* Given an x coordinate and a vector of k values, return
*  the y coordinate of the spline. The ks vector should be
*  the same length as the points vector
*/
double getSplineHeightAt(double x, std::vector<double> ks)
{
// If we don't have the right amount of K's we can't do anything
if (ks.size() != points.size()) {
std::cerr << "We don't have enough K's for the points we have, there might be two points directly above each other" << std::endl;
return -1;
}

// Find the index of the point we will use as one of the points in the spline
int index;
for (index = 1; index < points.size(); index++) {
if (points[index].x >= x) {
// We found the point we want
break;
}
}

// If we didn't get a valid index, give up
if (index == points.size()) {
return -1;
}

// If these are equal, the division after this would fail, so show an error and abort
if (points[index].x == points[index - 1].x) {
std::cerr << "Two points are directly above each other at x=" << points[index].x
<< ", unable to calculate spline for that point";
return -1;
}

// The difference between the current point's x and y coordinates and the previous one
double xdiff = points[index].x - points[index - 1].x;
double ydiff = points[index].y - points[index - 1].y;

// Perform the calculations needed to determine the y coordinate
double t = (x - points[index - 1].x) / xdiff;
double a = ks[index - 1] * xdiff - ydiff;
double b = -ks[index] * xdiff + ydiff;

// Interpolate to the point we want and then return it
double q = (1 - t) * points[index - 1].y + t * points[index].y + t * (1 - t) * (a * (1 - t) + b * t);
return q;
}

/**************************************************************
* Event handlers
**************************************************************/

void onDisplay(void)
{
// Clear the buffer
glClear(GL_COLOR_BUFFER_BIT);

glColor3f(1.0f, 1.0f, 1.0f);

glPushMatrix();

// Draw the spline (if we have at least 3 points)
if (points.size() >= 3) {
glBegin(GL_LINE_STRIP);
std::vector<double> ks = getKs();

// Loop over all the points (skipping the first
for (int i = 1; i < points.size(); ++i) {
Point cur = points[i];
Point prev = points[i - 1];

// Iterate over the X's between this point and the previous
double increment = (cur.x - prev.x) / steps;
for (double x = prev.x; x <= cur.x; x += increment) {
double y = getSplineHeightAt(x, ks);
if (y == -1) {
// There was an error
break;
}
glVertex2d(x, y);
}
}
glEnd();
}

// Draw the points
glPointSize(POINT_SIZE);
glBegin(GL_POINTS);
for (auto it = points.begin(); it != points.end(); ++it) {
Point point = *it;
glVertex2d(point.x, point.y);
}
glEnd();

glPopMatrix();

// Draw now
glutSwapBuffers();
}

/**
* Control dragging when the mouse moves with a button down
*/
void onActiveMouseMotion(int x, int y)
{
if (!dragging) {
return;
}

// Remove the previous point (if it exists)
if (!(lastPoint.x == -1 && lastPoint.y == -1)) {
for (auto it = points.begin(); it != points.end(); ++it) {
Point p = *it;
if (p.x == lastPoint.x && p.y == lastPoint.y) {
points.erase(it);
break;
}
}
}

// Make a new point for the current location and insert it
lastPoint = Point{ (double)x, (double)y };
insertPoint(lastPoint);

// Since we changed something we want to redraw
glutPostRedisplay();
}

/**
* If a mouse button is pressed, see if we need to drag or create a new point
*/
void onMouseButton(int button, int state, int x, int y)
{
// Do something with the mouse
std::cout << "Click: " << button << " " << state << " " << x << " " << y << std::endl;

if (button != GLUT_LEFT_BUTTON) {
// We only care about the left button
return;
}

if (state == 0) {
// See if they are trying to choose a point
for (auto it = points.begin(); it != points.end(); ++it) {
Point p = *it;

// See if the mouse is on the square for the point
if (x >= p.x - POINT_SIZE / 2 && x <= p.x + POINT_SIZE / 2 && y >= p.y - POINT_SIZE / 2 && y <= p.y + POINT_SIZE / 2) {
if (!shouldDelete) {
// We found a point to drag
dragging = true;
lastPoint = p;
}
else {
// We want to delete, so delete this point
points.erase(it);
}
break;
}
}
if (!dragging) {
lastPoint = Point{ -1, -1 };
}
}

if (state == 1) {
if (!dragging && !shouldDelete) {
// If we weren't dragging or deleting anything, insert a point on mouse up
insertPoint(Point{ (double)x,  (double)y });
}
else {
// Otherwise, say we're done dragging
dragging = false;
lastPoint = Point{ -1, -1 };
}
}

// Since we changed something we want to redraw
glutPostRedisplay();
}

void onKeyDown(unsigned char key, int x, int y)
{
// Do something with the keyboard
std::cout << "Key down: " << key << " " << x << " " << y << std::endl;
switch (key) {
case 'd':
case 'D':
shouldDelete = true;
break;
}
}

void onKeyUp(unsigned char key, int x, int y)
{
// Do something with the keyboard
std::cout << "Key up: " << key << " " << x << " " << y << std::endl;

switch (key) {
case 'q':
case 'Q':
exit(0);
break;
case 'd':
case 'D':
shouldDelete = false;
break;
case '+':
steps += 1;
if (steps > 15) steps = 15;
break;
case '-':
steps -= 1;
if (steps <= 0) steps = 1;
break;
}

// Since we changed something we want to redraw
glutPostRedisplay();
}

{
switch (option) {
case INCREASE_STEPS:
steps += 1;
if (steps > 15) steps = 15;
break;
case DECREASE_STEPS:
steps -= 1;
if (steps <= 0) steps = 1;
break;
case QUIT:
exit(0);
break;
}

// Since we changed something we want to redraw
glutPostRedisplay();
}

/**************************************************************
* Main/init functions
**************************************************************/

void initGL(void)
{
// Set clear color
glClearColor(0.0, 0.0, 0.0, 0.0);

// Initialize view
glMatrixMode(GL_PROJECTION);
gluOrtho2D(0, WIDTH, HEIGHT, 0);
}

/**
* Declare initial window size, position, and display mode
* Set up callback functions for glut events
*/
int main(int argc, char** argv)
{
// Init everything
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
glutInitWindowSize(WIDTH, HEIGHT);
glutCreateWindow("Cubic Spline");
initGL();

// Hook up callbacks
glutDisplayFunc(onDisplay);
glutMouseFunc(onMouseButton);
glutMotionFunc(onActiveMouseMotion);
glutKeyboardFunc(onKeyDown);
glutKeyboardUpFunc(onKeyUp);

// Go!
glutMainLoop();

return 0;
}
``````
###### Who is Participating?
I wear a lot of hats...

"The solutions and answers provided on Experts Exchange have been extremely helpful to me over the last few years. I wear a lot of hats - Developer, Database Administrator, Help Desk, etc., so I know a lot of things but not a lot about one thing. Experts Exchange gives me answers from people who do know a lot about one thing, in a easy to use platform." -Todd S.

Commented:
Why you need opengl to draw dots?
0
Author Commented:
not to draw dots, this draw cubic spline curve like in this page http://blog.ivank.net/interpolation-with-cubic-splines.html
0
Commented:
if the line is a 3d line and hasn't all points at a plane (cut) in the 3D space, the math is completely different and you can't use spline math which relies on 2D. you could see that by using a line that is on the surface of a 3D object, for example on a sphere. then all interpolation points are located at the spehre but any spline function between two base points has all points beside of begin and end within the the sphere what surely is not what you wanted. what you could do is to firstly find a 'manifold' function (mapping) which contains all interpolation points. you may think of a map made of a region of the hemisphere. any 3D point of the region was mapped homeomorphicly to a 2D point of the map. if you have found such a manifold you could map all 3D interpolation points to 2D, then calculate the spline, and remap spline points for drawing back to the 'manifold'. this requires high mathematics in complex analysis.

note, cubic splines are not 3D curves but polynoms of the form f(x) = a*x^3 + b*x^2 +c*x +d, where any point is a 2D point (x, y).

Sara
0
Author Commented:
but our professor said as long as you did the 2D part. it is easy to extend it to 3D!! how is going to be different math?
0
Commented:
can you exactly repeat what your professor said? or can you post  a full description of the task?

if you look at

http://www.uio.no/studier/emner/matnat/ifi/INF-MAT5340/v07/undervisningsmateriale/kap7.pdf

you easily can see that a function of type

z = f(x, y)

where the right side is a mixed 3-dim polynom like z = a*x^3*y^2 +b*x^2*y^2 + c*x^2 +d*y +e

can't be handled as simple as a "normal" spline function of the form a*x^3 + b*x^2 +c*x +d.

of course you could limit the spline functions to be a sum of an x^3 polynom and a y^3 polynom, which has 8 unknowns and where you easily can make partial derivations for both x and y and perhaps it is rather easy to make up an equations system for such kind of splines (though you need 8 equations for each spline).

but even if so, it is not easy to extend your current code to handle this, since the complexity was more than doubled.

Sara
0
Commented:
If you look into your code you draw DOTS, until they render a line.
That should be fine for 2048x2048 surface but now lets try 2048x2048x2048 space... and try to fill thet resonable... Now use opengl to turn our object... ouch.... we actually did not need dors, we needed lines (or triangles or curves) to start with.
0
Author Commented:

Extending 2D interpolation into 3D is actually very trivial. Your original interpolation was already done for 3D in fact, if you consider each control point <x,y> was just a very special case of <x,y,z> with z always == 0. The interpolation evaluator you used for 2D spline, inserted a number of new points between each pair of control points, by adding new interpolated values between Xs, and between Ys. We didn't add new Zs between each pair of Zs of control points, only because in our homework1, all Zs are set to be 0 (since it was 2D only) --- so that if we were doing interpolation then all added Zs would have been 0 anyways.

Now if you set non-zero Z values to the control points, and use exactly the same interpolation function to add Zs between the control points' Z values, you will be interpolating in 3D. In fact, cubic interpolation can be done in arbitrary dimensions --- you are using the same interpolation function to add more values between each pair of existing Xs, Ys, and Zs, ... etc.

and instead of drawing dots we will draw spheres he said ..
0
Commented:
ok. that makes things easier indeed.

what your professor said means that you simply use your current spline functionality twice, once to approximate

y = f(x)

and once for

z = g(z)

the only difference is that you use the z values instead of the y values for your second run (spline interpolation).

and instead of drawing dots we will draw spheres he said ..

not so quite clear what that means. can you post a sample?

a sphere is determined by center point and radius. so, given an x-y  diagram you could draw a circle with radius z instead of a dot for each 3d point, granted that all z values are positive.

if you have a x-y-z diagram (3D) the dots would be 3D. you could use spheres instead of dots but don't know what the radius is and what the advantage is over dots (which actually are spheres as well).

Sara
0

Experts Exchange Solution brought to you by