Skip to content

Instantly share code, notes, and snippets.

@ryanthejuggler
Created November 22, 2012 16:55
Show Gist options
  • Save ryanthejuggler/4132103 to your computer and use it in GitHub Desktop.
Save ryanthejuggler/4132103 to your computer and use it in GitHub Desktop.
C++ cubic spline interpolation
#include <std::vector>
#include <opencv2/core/core.hpp>
#include <iostream>
#include <stdio.h>
/*
This work by Ryan Muller released under the Creative Commons CC0 License
http://creativecommons.org/publicdomain/zero/1.0/
*/
/**
Spline interpolation of a parametric function.
INPUT: std::vector<double> x
A list of double values that represent sampled
points. The array index of each point will be
taken as the parameter t by which x will be
represented as a function.
OUTPUT: std::vector<cv::Vec4d> P
A list of cv::Vec4d representing polynomials. To
interpret segment [i]:
x(t) = P0*a + P1*b + P2*(a^3-a)/6 + P3*(b^3-b)/6
where a = t-i
b = i-t+1
*/
std::vector<cv::Vec4d> splinterp(std::vector<double> x){
std::vector<cv::Vec4d> out;
// spline size
int n=x.size();
// loop counter
int i;
// working variables
double p;
std::vector<double> u;
// second derivative
std::vector<double> z;
u.resize(n);
z.resize(n);
// set the second derivative to 0 at the ends
z[0] = u[0] = 0;
z[n-1] = 0;
// decomposition loop
for(i=1;i<n-1;i++){
p = 0.5*z[i-1] + 2.0;
z[i] = -0.5/p;
u[i] = x[i+1]+x[i-1]-2*x[i];
u[i] = (3*u[i]-0.5*u[i-1])/p;
}
// back-substitution loop
for(i=n-1;i>0;i--){
z[i] = z[i] * z[i+1] + u[i];
}
for(i=0;i<n-1;i++){
out.push_back(cv::Vec4d(
x[i+1],
x[i],
z[i+1],
z[i]
));
}
return out;
}
/**
A demo of Splinterp's capabilities.
Modify the function in the for loop to change
the sample points. When you run the program
it will generate a polynomial for each segment
and output it in text format. Copy and paste
the output into WxMaxima (andrejv.github.com/wxmaxima/)
to see a plot of your interpolated curves!
*/
void splinterpDemo(){
std::vector<double> samples;
std::vector<cv::Vec4d> spline;
cv::Vec4d p;
double x;
unsigned int i, imax=12;
for(i=0;i<imax;i++){
// modify this function
x = (cos(2*M_PI*i/(imax-1)));
samples.push_back(x);
}
spline = splinterp(samples);
for(i=0;i<spline.size();i++){
p = spline.at(i);
printf("x%d(t):= if t>= %d and t<=%d then %.8f * (t-%d) + %.8f * (%d-t) + "
"%.8f * ((t-%d)^3-(t-%d))/6 + %.8f * ((%d-t)^3-(%d-t))/6 $\n",
i,i,i+1,p[0],i,p[1],i+1,p[2],i,i,p[3],i+1,i+1);
}
printf("wxplot2d([");
for(i=0;i<spline.size();i++){
if(i>0) cout << ",";
printf("x%d(t)",i);
}
printf("],[t,0,%d],[gnuplot_preamble,\"set nokey;\"])$\n",spline.size());
}
@titsitits
Copy link

Hi,
I tried the algorithm, had the same bug, and arrived on your comment.
I'm not sure, but I think the second derivative at the last sample can just stay 0.
So no need for the line z[n-1] = u[n-1] in your proposition.
Besides I noticed that u[n-1] is never initialized, nor set in the first loop, so it would just give a completely random value.
Cheers

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment