I need a simple C++ code snippet to do B-spline interpolation

[JoshuaMostafa]

Hello,

I am working on a project, and part of it needs smoothing a set of points. Say i have a set of points in 3d space and want to generate the best smooth line that represents them.

For now, I am just connecting them together, but I really need to do B-Spline. The points are saved in a text file. Something like this:

-129.87 9.87 50.71
-19.27 1.36 21.55
1.44 2.23 45.63
2.85 5.78 0.79

First row represents the first point, second row represents the second point, ..etc.

I am looking for the simplest code that can do b-spline for these control points. The number of control points is not fixed. Different set have different number of points, but that should not be a problem.

Any help is appreciated.

Thank you,

[48T7M4Gy]

Google c++ b-spine is a possibility. maybe boost has something and probably one of the graphics languages does too. if you intend programming it yourself, and from memory the maths is not all that hard, start with Wikipedia or math lab.

[JoshuaMostafa]

I have been doing some search, but unfortunately I have not found anything that I am able to implement easily or even understand. I think I should go and brush up on B-spline as a concept first.

One of the professors gave me this code to help out, but it is written in C. I don't know much of what's going on there. I think this is for 2d and I need to expand it to 3d and convert it to C++.

#include <stdio.h>
#include <math.h>

#define PI 3.14159265

void spline(float *, float *, int, float, float, float *);
void splint(float *, float *, float *, int, float, float *);

void main(void)
{
  int i, nfunc, np;
  float f, x, y, yp1, ypn, *xa, *ya, *y2;
  
  np = 10;       /* number of points */
  
  xa = (float *) malloc((np+1)*sizeof(float));
  ya = (float *) malloc((np+1)*sizeof(float));
  y2 = (float *) malloc((np+1)*sizeof(float));
  
  printf("\nsine function from 0 to pi\n");
  for (i=1; i<=np; i++) {
    xa[i] = i*PI/np;
    ya[i] = sin(xa[i]);
  }
  yp1 = cos(xa[1]);
  ypn = cos(xa[np]);
  
  /* call spline to get second derivatives */
  spline(xa, ya, np, yp1, ypn, y2);

  /* call splint for interpolations */
  printf("\n%9s %13s %17s\n", "x", "f(x)", "interpolation");
  for (i=1; i<=10; i++) {
    x = (-0.05+i/10.0)*PI;
    f = sin(x);
    splint(xa, ya, y2, np, x, &y);
    printf("%12.6f %12.6f %12.6f\n", x, f, y);
  }
}


void spline(float *x, float *y, int n, float yp1, float ypn, float *y2)
{
  int i, k;
  float p, qn, sig, un, *u;
  
  u = (float *) malloc(n*sizeof(float));
  
  if (yp1 > 0.99e30)
    y2[1] = u[1] = 0.0;
  else {
    y2[1] = -0.5;
    u[1] = (3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
  }
  
  for (i=2; i<=n-1; i++) {
    sig = (x[i]-x[i-1])/(x[i+1]-x[i-1]);
    p = sig*y2[i-1]+2.0;
    y2[i] = (sig-1.0)/p;
    u[i] = (y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
    u[i] = (6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
  }
  
  if (ypn > 0.99e30)
    qn = un = 0.0;
  else {
    qn = 0.5;
    un = (3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
  }
  
  y2[n] = (un-qn*u[n-1])/(qn*y2[n-1]+1.0);
  
  for (k=n-1; k>=1; k--)
    y2[k] = y2[k]*y2[k+1]+u[k];
}


void splint(float *xa, float *ya, float *y2a, int n, float x, float *y)
{
  int klo, khi, k;
  float h, b, a;
  
  klo=1; khi=n;
  while (khi-klo > 1) {
    k=(khi+klo) >> 1;
    if (xa[k] > x) khi=k;
    else klo=k;
  }
  
  h = xa[khi]-xa[klo];
  if (h == 0.0)
    printf("Error: Bad xa input to routine splint.\n");
  
  a = (xa[khi]-x)/h;
  b = (x-xa[klo])/h;
  *y = a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
}

[48T7M4Gy]

Looks very 'interesting'.

This also might be useful. Beyond that I can't really say any more. I roughly know that splines are about but I have never used or studied the technology other than in passing
http://eigen.tuxfamily.org/dox/unsupported/classEigen_1_1Spline.html