diff --git a/Makefile b/Makefile
index 4a24c07a30161264ac5d90da4ec81964124e2b6f..40e5696794c08517424f92677213d4585cd3cf9b 100644
--- a/Makefile
+++ b/Makefile
@@ -7,7 +7,7 @@ PROJNAME = miscmaths
USRINCFLAGS = -I${INC_NEWMAT} -I${INC_PROB}
USRLDFLAGS = -L${LIB_NEWMAT} -L${LIB_PROB}
-OBJS = miscmaths.o optimise.o miscprob.o kernel.o histogram.o base2z.o t2z.o f2z.o volume.o volumeseries.o minimize.o
+OBJS = miscmaths.o optimise.o miscprob.o kernel.o histogram.o base2z.o t2z.o f2z.o volume.o volumeseries.o minimize.o cspline.o
LIBS = -lutils -lnewmat -lprob -lm
diff --git a/cspline.cc b/cspline.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f3a4e7efa4d9732cac7d10561076dfd5421814b0
--- /dev/null
+++ b/cspline.cc
@@ -0,0 +1,303 @@
+/* cspline
+
+ Cubic spline fitting and interpolation
+
+ Tim Behrens, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2000 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <unistd.h>
+
+#include "newmatap.h"
+#include "newmatio.h"
+#include "miscmaths/miscmaths.h"
+#include "cspline.h"
+
+#define WANT_STREAM
+#define WANT_MATH
+
+using namespace NEWMAT;
+using namespace std;
+
+///////////////////////////////////////////////////////
+
+
+
+namespace MISCMATHS{
+
+ // void Cspline::Cspline(){}
+
+ void Cspline::set(ColumnVector& pnodes,ColumnVector& pvals){
+ nodes=pnodes;vals=pvals;
+ fitted=false;
+ n=vals.Nrows();
+ }
+ void Cspline::set(ColumnVector& pnodes, Matrix& pcoefs){
+ nodes=pnodes;coefs=pcoefs;
+ fitted=false;
+ n=vals.Nrows();
+ }
+
+ void Cspline::diff(const ColumnVector& x, ColumnVector& dx ){
+ // dx should be of length length(x)-1
+ dx.ReSize(x.Nrows()-1);
+ for(int i=2;i<=x.Nrows();i++){
+ dx(i-1)=x(i)-x(i-1);
+ }
+ }
+
+
+ void Cspline::fit(){
+ if(vals.Nrows()<4){
+ cerr<<"Cspline::fit - You have less than 4 data pts for spline fitting."<<endl;
+ exit(-1);
+ }
+ if(nodes.Nrows()!=vals.Nrows()){
+ cerr<<"Nodes and VALS must be the same length in your spline"<<endl;
+ exit(-1);
+ }
+ int n=vals.Nrows();
+ ColumnVector s(n);
+ ColumnVector dx,dy,dydx(n-1);
+ diff(nodes,dx);
+ diff(vals,dy);
+
+ for(int i=1;i<=n-1;i++){
+ dydx(i)=dy(i)/dx(i);
+ }
+ ColumnVector b(n);
+ b=0;
+ for(int i=2;i<=b.Nrows()-1;i++){
+ b(i)=3*(dx(i)*dydx(i-1)+dx(i-1)*dydx(i));
+ }
+
+ float x31=nodes(3)-nodes(1),xn=nodes(n)-nodes(n-2);
+ b(1)=((dx(1)+2*x31)*dx(2)*dydx(1)+dx(1)*dx(1)*dydx(2))/x31;
+ b(n)=(dx(n-1)*dx(n-1)*dydx(n-2)+(2*xn+dx(n-1))*dx(n-2)*dydx(n-1))/xn;
+
+ Matrix tridiag(n,n);
+ tridiag=0;
+ ColumnVector y3(n);
+ for(int j=2;j<=n-1;j++){
+ tridiag(j,j-1)=dx(j);
+ tridiag(j,j)=2*(dx(j)+dx(j-1));
+ tridiag(j,j+1)=dx(j-1);
+ }
+ tridiag(1,1)=dx(2);tridiag(1,2)=x31;
+ tridiag(n,n-1)=xn;tridiag(n,n)=dx(n-2);
+ s=tridiag.i()*b;
+
+
+ ColumnVector d(n-1),c(n-1);
+ for(int j=1;j<n;j++){
+ d(j)=(s(j)+s(j+1)-2*dydx(j))/dx(j);
+ c(j)=(dydx(j)-s(j))/dx(j)-d(j);
+ }
+
+ coefs.ReSize(n-1,4);
+ for(int j=1;j<n;j++){
+ coefs(j,1)=vals(j);
+ coefs(j,2)=s(j);
+ coefs(j,3)=c(j);
+ coefs(j,4)=d(j)/dx(j);
+ }
+ cerr<<coefs<<endl;
+ }
+
+
+ float Cspline::interpolate(float xx) const{
+ // nodes must be monotonically increasing. I don't check this.
+ // On your head be it if you don't.
+ if(nodes.Nrows()!=vals.Nrows()){
+ cerr<<"Cspline:interpolate: Nodes and Vals should be the same length"<<endl;
+ exit(-1);
+ }
+
+ float ret;
+ if(!fitted){
+ cerr<<"Cspline::interpolate - Cspline has not been fitted"<<endl;
+ exit(-1);
+ }
+ else{
+
+ bool stop=false;
+ int ind=0;
+
+ if(xx<nodes(1)){
+ ind=1;
+ }
+ else if(xx>nodes(nodes.Nrows())){
+ ind=nodes.Nrows()-1;
+ }
+ else{
+ for(int i=1;i<nodes.Nrows();i++){
+ if(!stop){
+ if( (nodes(i)<=xx) && (nodes(i+1)>xx) ){
+ ind=i;
+ stop=true;
+ }
+ }
+ }
+ }
+
+ float a=coefs(ind,1);
+ float b=coefs(ind,2);
+ float c=coefs(ind,3);
+ float d=coefs(ind,4);
+ float t=xx-nodes(ind);
+ ret=a+b*t+c*t*t+d*t*t*t;
+
+ }
+ return ret;
+}
+
+ float Cspline::interpolate(float xx, int ind) const{
+ float ret;
+ if(!fitted){
+ cerr<<"Cspline::interpolate - Cspline has not been fitted"<<endl;
+ exit(-1);
+ }
+ else{
+ if(ind>n-1){
+ cerr<<"Cspline::interpolate - segment index is greater than number of segments - exiting"<<endl;
+ exit(-1);
+ }
+ else if(ind<1){
+ cerr<<"Cspline::interpolate - segment index is less than 1 - exiting"<<endl;
+ exit(-1);
+ }
+ float a=coefs(ind,1);
+ float b=coefs(ind,2);
+ float c=coefs(ind,3);
+ float d=coefs(ind,4);
+ float t=xx-nodes(ind);
+ ret=a+b*t+c*t*t+d*t*t*t;
+ }
+ return ret;
+ }
+
+
+ ColumnVector Cspline::interpolate(const ColumnVector& x) const{
+ // nodes must be monotonically increasing. I don't check this.
+ // On your head be it if you don't.
+
+ if(nodes.Nrows()!=vals.Nrows()){
+ cerr<<"Cspline::interpolate - Nodes and Vals should be the same length"<<endl;
+ exit(-1);
+ }
+
+ ColumnVector ret(x.Nrows());
+
+ if(!fitted){
+ cerr<<"Cspline::interpolate - Cspline has not been fitted"<<endl;
+ exit(-1);
+ }
+ else{
+
+ for(int xnum=1;xnum<=x.Nrows();xnum++){
+
+ float xx=x(xnum);
+
+ bool stop=false;
+ int ind=0;
+
+ if(xx<nodes(1)){
+ ind=1;
+ }
+ else if(xx>=nodes(nodes.Nrows())){
+ ind=nodes.Nrows()-1;
+ }
+ else{
+ for(int i=1;i<nodes.Nrows();i++){
+ if(!stop){
+ if( (nodes(i)<=xx) && (nodes(i+1)>xx) ){
+ ind=i;
+ stop=true;
+ }
+ }
+ }
+ }
+
+ float a=coefs(ind,1);
+ float b=coefs(ind,2);
+ float c=coefs(ind,3);
+ float d=coefs(ind,4);
+ float t=xx-nodes(ind);
+ ret(xnum)=a+b*t+c*t*t+d*t*t*t;
+ }
+ }
+ return ret;
+ }
+
+
+ ColumnVector Cspline::interpolate(const ColumnVector& x,const ColumnVector& indvec) const{
+ // nodes must be monotonically increasing. I don't check this.
+ // On your head be it if you don't.
+
+ if(nodes.Nrows()!=vals.Nrows()){
+ cerr<<"Cspline::interpolate - Nodes and Vals should be the same length"<<endl;
+ exit(-1);
+ }
+
+ ColumnVector ret(x.Nrows());
+
+ if(!fitted){
+ cerr<<"Cspline::interpolate - Cspline has not been fitted"<<endl;
+ exit(-1);
+ }
+ else{
+ for(int xnum=1;xnum<=x.Nrows();xnum++){
+
+ float xx=x(xnum);
+
+ int ind=int(indvec(xnum));
+
+ float a=coefs(ind,1);
+ float b=coefs(ind,2);
+ float c=coefs(ind,3);
+ float d=coefs(ind,4);
+ float t=xx-nodes(ind);
+ ret(xnum)=a+b*t+c*t*t+d*t*t*t;
+ }
+ }
+ return ret;
+ }
+
+
+
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/cspline.h b/cspline.h
new file mode 100644
index 0000000000000000000000000000000000000000..5c231f258feb94da4a26a58f95f1e9a3c782c94c
--- /dev/null
+++ b/cspline.h
@@ -0,0 +1,74 @@
+/* cspline
+
+ Cubic spline fitting and interpolation
+ Tim Behrens, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2000 University of Oxford */
+
+/* CCOPYRIGHT */
+#if !defined(__cspline_h)
+#define __cspline_h
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <unistd.h>
+
+#include "newmatap.h"
+#include "newmatio.h"
+#include "miscmaths/miscmaths.h"
+
+#define WANT_STREAM
+#define WANT_MATH
+
+using namespace NEWMAT;
+using namespace std;
+///////////////////////////////////////////////////////
+
+
+namespace MISCMATHS {
+ class Cspline{
+ public:
+ Cspline();
+ Cspline(ColumnVector& pnodes,ColumnVector& pvals):
+ nodes(pnodes),
+ vals(pvals),
+ n(nodes.Nrows())
+ {
+ fit();
+ fitted=true;
+ }
+
+ Cspline(ColumnVector& pnodes, Matrix& pcoefs) :
+ nodes(pnodes),
+ coefs(pcoefs),
+ n(nodes.Nrows())
+ { fitted=true;}
+
+ ~Cspline(){
+ fitted=false;
+ };
+
+ void set(ColumnVector& pnodes,ColumnVector& pvals);
+ void set(ColumnVector& pnodes, Matrix& pcoefs);
+
+ void fit();
+ float interpolate(float xx) const;
+ float interpolate(float xx,int ind) const;
+ ColumnVector interpolate(const ColumnVector& x) const;
+ ColumnVector interpolate(const ColumnVector& x, const ColumnVector& indvec) const;
+
+ protected:
+
+ bool fitted;
+ ColumnVector nodes;
+ ColumnVector vals;
+ Matrix coefs;
+ int n;
+ void diff(const ColumnVector& x, ColumnVector& dx );
+
+ };
+}
+
+
+#endif