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 matricestypedef std::vector<double> MatrixRow;typedef std::vector<MatrixRow> Matrix;// Struct for representing a point in 2D spacestruct Point { double x; double y;};// Variables for dragging points arroundbool dragging = false;Point lastPoint = Point{ -1, -1 };// If the next click should delete a point instead of dragging or adding onebool shouldDelete = false;// The number of interpolation steps to do between pointsint steps = 5;// The list of spline points, sorted from lowest x value to higheststd::vector<Point> points;// Options for the menuenum MENU_OPTIONS{ 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) { // We only care about this if we're 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();}void onMenu(int option){ 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); glLoadIdentity(); 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); // Set up the menu glutCreateMenu(onMenu); glutAddMenuEntry("Increase # of interpolation steps", INCREASE_STEPS); glutAddMenuEntry("Decrease # of interpolation steps", DECREASE_STEPS); glutAddMenuEntry("Quit", QUIT); glutAttachMenu(GLUT_RIGHT_BUTTON); // Go! glutMainLoop(); return 0;}

"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.

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

The most successful MSPs rely on metrics – known as key performance indicators (KPIs) – for making informed decisions that help their businesses thrive, rather than just survive. This eBook provides an overview of the most important KPIs used by top MSPs.

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.

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.

AaeshahAuthor Commented:

/this the professor reply:

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 ..

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

Experts Exchange Solution brought to you by

Your issues matter to us.

Facing a tech roadblock? Get the help and guidance you need from experienced professionals who care. Ask your question anytime, anywhere, with no hassle.

It's more than this solution.Get answers and train to solve all your tech problems - anytime, anywhere.Try it for freeEdge Out The Competitionfor your dream job with proven skills and certifications.Get started todayStand Outas the employee with proven skills.Start learning today for freeMove Your Career Forwardwith certification training in the latest technologies.Start your trial today