Dave Mitchell's DigitalMicrograph™ Scripting Website

// Parameterised (fast) cubic spline calculates the series of constatnts A, B, C and D

// in the equation y=A + Bx + cx^2 + dx^3. The array of constants are calculated in a single slow step by calculatesplineconstants().

// They are then passed to a second function (fastcubicspline()) which uses them in a loop for

// quick calculation a series of interpolated values.

// Derived from 'Numerical Analysis', 6th Edition, R. L. Burden and J. Douglas Faires, Brooks/Cole Publishing Co.

// version:20091007

// A series of x,y data points (with x increasing acroos the series) is passed into this function

// as an image. The function returns a series of four constants for each interval in the x data, for

// use by the fastcubicspline() function. This pair of functions is about 10x faster than the

// equivalent natural cubic spline function, which is a single function. This convenient to use but has to

// calculate all the constants at each step (slow).

// D. R. G. Mitchell, drmnospam@ansto.gov.au (remove the nospam)

// October 7th 2009, v1.1

// v1.1 Fixed a bug in the spine calculation function

// The array of data points (x0, y0; x1 y1; etc) is passed into the function as an image.

// The image contains the point values stored in pixels. The image is two pixels high

// and as many pixels wide as the data set. The format of the the data is:

// x0, x1, x2, x3, . . . . xn

// y0, y1, y2, y3, . . . . yn

// If the image is anything other than 2 pixels high, the function reports an error.

// If the x values are not sorted into increasing order - left to right, then the function reports an error.

// The function returns an image 5 pixels high and n pixels wide, where n is the number of data points.

// The format of the image is

// -, A1, A2, A3, . . . . An

// - ,B1, B2, B3, . . . .Bn

// -, C1, C2, C3, . . . . Cn

// -, D1, D2, D3, . . . . Dn

// -, X1, X2, X3, . . . . Xn

// The X values are sourced from the passed in array of data points.

// note the 0 positions (-) in the image are unused.

// The constants are used to calculate the interpolated y value as follows:

// y interpolated = A + Bxc + Cxc^2 + Dxc^3, where xc=xi-xvalue

// xi is the x value at the start of the interval in which the passed in x value (xvalue) lies.

// If the x data run 10, 25, 55, 72, 101 and the y value is required when x=65, then xc=55-65

// Function to calculate the cubic spline constants

image calculatesplineconstants(image dataset)

{

// variables

number n, sizex, sizey, minx, maxx, yspline, i, prevval, thisx,m,j

// If there any problems with the passed in data image return a null image (1 x 1 pixels and set to 0)

// Use this for error checking

image nullimg=integerimage("",1,1,1,1)

// Check the passed in image dataset to ensure that the image 2 pixels high by n pixels wide.

getsize(dataset, sizex, sizey)

if(sizey!=2)

{

result("\nNatural Cubic Spline Error: Passed in dataset is the wrong size!")

return nullimg

}

// the number of data points

n=sizex-1

// Get the minimum and maximum values in the x data of the dataset (top row)

// note pixel position 0 IS used for the passed in data - only in the output constants, and x values

// (copied to the fastcubicspline() function) is pixel position 0 not used

minmax(dataset[0,0,1,sizex], minx, maxx)

// Sort through the data checking that x increases monotonically

prevval=0

for(i=0; i<sizex; i++)

{

thisx=getpixel(dataset, i,0)

if(i>0 && thisx<prevval)

{

result("\nConstrained Cubic Spline Error: Data for X must be in ascending order!")

return nullimg

}

prevval=thisx

}

// Arrays to store the data points

// note the data start at pixel position 1 - pixel position 0 is not used

image x=realimage("",4,sizex+1, 1)

image a=realimage("",4,sizex+1, 1)

image xa=realimage("",4,sizex+1, 1)

image h=realimage("",4,sizex+1, 1)

image xl=realimage("",4,sizex+1, 1)

image xu=realimage("",4,sizex+1, 1)

image xz=realimage("",4,sizex+1, 1)

image b=realimage("",4,sizex+1, 1)

image c=realimage("",4,sizex+1, 1)

image d=realimage("",4,sizex+1, 1)

image constantarray=realimage("",4,sizex+1,5) // stores all the constants and the x values

// set the x and y arrays to the respective values in the dataset passed in

// Note the passed in data has values running from pixel position 0 to xsize-1

// (remember an image 18 pixels wide has pixel positions 0 - 17. The x and y

// data are transferred to arrays x and a - both xisze+1 pixels wide. The data

// is shifted one pixel position to the left so that pixel position 0 is ignored.

// There are data in positions 1 to xisze+1.

x[0,1,1,sizex+1]=dataset[0,0,1,sizex]

a[0,1,1,sizex+1]=dataset[1,0,2,sizex]

m=n-1

// Step 1

for(i=0; i<m+1;i++)

{

h[0,i+1,1,i+2]=getpixel(x,i+2,0)-getpixel(x, i+1,0)

}

// Step 2

for(i=1; i<m+1;i++)

{

xa[0,i+1,1,i+2]=3*(getpixel(a,i+2,0)*getpixel(h,i,0)-getpixel(a,i+1,0)*(getpixel(x,i+2,0)-

getpixel(x,i,0))+getpixel(a,i,0)*getpixel(h,i+1,0))/(getpixel(h,i+1,0)*getpixel(h,i,0))

}

// Step 3

setpixel(xl,1,0,1)

setpixel(xu,1,0,0)

setpixel(xz,1,0,0)

// Step 4

for(i=1;i<m+1;i++)

{

xl[0,i+1,1,i+2]=2*(getpixel(x,i+2,0)-getpixel(x,i,0))-getpixel(h,i,0)*getpixel(xu,i,0)

xu[0,i+1,1,i+2]=getpixel(h,i+1,0)/getpixel(xl,i+1,0)

xz[0,i+1,1,i+2]=(getpixel(xa,i+1,0)-getpixel(h,i,0)*getpixel(xz,i,0))/getpixel(xl,i+1,0)

}

// Step 5

setpixel(xl,n+1,0,1)

setpixel(xz,n+1,0,0)

setpixel(c,n+1,0,getpixel(xz,n+1,0))

// step 6

for(i=0;i<m+1;i++)

{

j=m-i

c[0,j+1,1,j+2]=getpixel(xz,j+1,0)-getpixel(xu,j+1,0)*getpixel(c,j+2,0)

b[0,j+1,1,j+2]=(getpixel(a,j+2,0)-getpixel(a,j+1,0))/getpixel(h,j+1,0)-getpixel(h,j+1,0)*(getpixel(c,j+2,0)+2*getpixel(c,j+1,0))/3

d[0,j+1,1,j+2]=(getpixel(c,j+2,0)-getpixel(c,j+1,0))/(3*getpixel(h,j+1,0))

}

// Copy the a, b, c and d images to the array image and return it

constantarray[0,0,1,sizex+1]=a[0,0,1,sizex+1]

constantarray[1,0,2,sizex+1]=b[0,0,1,sizex+1]

constantarray[2,0,3,sizex+1]=c[0,0,1,sizex+1]

constantarray[3,0,4,sizex+1]=d[0,0,1,sizex+1]

constantarray[4,0,5,sizex+1]=x[0,0,1,sizex+1]

return constantarray

}

// Function to calculate the interpolated value of y from an array of cubic spline

// constants calculated from the function calculatesplineconstants(). See that function

// for more detail. The function takes the array of spline constants (an image), an x value

// from which to interpolate a yvalue, and a boolean (interpolate). If interpolate is set to 1 then

// interpolation outsid the range will be performed, otherwise values outside this range return a zero.

// the function returns the interpolated y value. As the constants are only calculated once in the

// other function, the computation in this function is very fast and suitable for use in a loop

// to calculate a range of y values.

number fastcubicspline(image constantarray, number xvalue, number extrapolate)

{

number minx, maxx, sizex, sizey, yspline, i, n

getsize(constantarray, sizex, sizey)

if(sizey!=5)

{

result("\nFast Cubic Spline Error: Passed in array of spline constants is the wrong size!")

return yspline

}

// Get the minimum and maximum values in the x data -bottom row of the array from position 1 to n

minmax(constantarray[4,1,5,sizex], minx, maxx) // ignore position 0

// Check that the passed in xvalue is within the range of the xvalues in the data set supplied

// If the extrapolate option is turned off (0) then the function returns a y value of zero

if(extrapolate==0) // do not extrapolate, any out of range values of x result a y value of zero

{

if(xvalue<minx || xvalue>maxx)

{

yspline=0

return yspline

}

}

// loop through the x data (row 5 - position 4 in the constantarray) to find which interval the passed in xvalue lies in

n=sizex-2 // note - 2 because pixel position 0 is unused.

for(i=1;i<n;i++)

{

if(xvalue<getpixel(constantarray,1,4)) break

if(xvalue>=getpixel(constantarray,i,4) && xvalue<getpixel(constantarray,i+1,4)) break

}

// Calculate the distance between the lower bound x data point for the interval and the passed in xvalue

number xcalc=xvalue-getpixel(constantarray,i,4)

// Compute the spline a=row 0, b=row 1, c=row 2 and d=row 3

// y=a +bxcalc+c*ccalc^2+dxdcalc^3

yspline=getpixel(constantarray,i,0)+getpixel(constantarray,i,1)*xcalc+getpixel(constantarray,i,2)*xcalc**2+getpixel(constantarray,i,3)*xcalc**3

return yspline

}

// Main program starts here

// Most of this code is associated with creating some example data and displaying it

// and drawing the spline through it. The actual spline part is a small loop at the end.

number i, dataxval, datayval

// Create an image to store the data in two rows of x and y x=5,12,52 etc.

image temp=integerimage("",2,1,18,2)

// These are the 18 x values stored in the top row of the image

setpixel(temp, 0,0,5)

setpixel(temp, 1,0,12)

setpixel(temp, 2,0,52)

setpixel(temp, 3,0,70)

setpixel(temp, 4,0,73)

setpixel(temp, 5,0,76)

setpixel(temp, 6,0,98)

setpixel(temp, 7,0,99)

setpixel(temp, 8,0,100)

setpixel(temp, 9,0,105)

setpixel(temp, 10,0,127)

setpixel(temp, 11,0,146)

setpixel(temp, 12,0,168)

setpixel(temp, 13,0,180)

setpixel(temp, 14,0,188)

setpixel(temp, 15,0,200)

setpixel(temp, 16,0,220)

setpixel(temp, 17,0,240)

// These are the 18 y values stored in the second row of the image

setpixel(temp, 0,1,8)

setpixel(temp, 1,1,3)

setpixel(temp, 2,1,5)

setpixel(temp, 3,1,6)

setpixel(temp, 4,1,5)

setpixel(temp, 5,1,6)

setpixel(temp, 6,1,12)

setpixel(temp, 7,1,8)

setpixel(temp, 8,1,7)

setpixel(temp, 9,1,4)

setpixel(temp, 10,1,-2)

setpixel(temp, 11,1,8)

setpixel(temp, 12,1,3)

setpixel(temp, 13,1,5)

setpixel(temp, 14,1,6)

setpixel(temp, 15,1,5 )

setpixel(temp, 16,1,6)

setpixel(temp, 17,1,5)

// Create a real image to display the interpolated values in (splineimg)

// and also a copy of the image (tempimg) to store the data values. This

// is added to the plot as a separate layer, so that both the data and

// the spline plot can be displayed simultaneously.

image splineimg:=realimage("",4,256,1)

setname(splineimg,"Fast Cubic Spline Interpolation")

image tempimg=imageclone(splineimg)

// Add the 18 x/y value pairs specified above to the tempimg image

for(i=0; i<18; i++)

{

dataxval=getpixel(temp, i, 0)

datayval=getpixel(temp, i,1)

setpixel(tempimg,dataxval, 0,datayval)

}

// Calculate the array of spline constants for the data stored in temp

image arrayimg=calculatesplineconstants(temp)

// Display the splineimg and get its lineplotimagedisplay()

showimage(splineimg)

lineplotimagedisplay lindisp=splineimg.imagegetimagedisplay(0)

// Add the image containing all the data points to the splineimg - with the name "Data"

// which appears in the legend and set the drawing style of slice 1 - the data - to 2 - filled.

// and set its colour to blue (0,0,255)

lindisp.imagedisplayaddimage(tempimg, "Data")

lineplotimagedisplaysetslicedrawingstyle(lindisp,1,2)

lineplotimagedisplaysetslicecomponentcolor(lindisp, 1, 1,0,0,255)

updateimage(tempimg)

// Get the zeroth slice which will hold the spline plot - set its name to Fast Spline

// and turn the plot legend on

object sliceid=lindisp.ImageDisplayGetSliceIDByIndex(0)

lindisp.ImageDisplaySetSliceLabelByID(sliceid, "Fast Spline")

lindisp.LinePlotImageDisplaySetLegendShown(1)

// Set the spline plot so that it appears as a red line

lineplotimagedisplaysetslicedrawingstyle(lindisp,0,1)

lineplotimagedisplaysetslicecomponentcolor(lindisp, 0, 0,255,0,0)

// This is the important part of the code which calls the spline function

// passing in to it the array of constants (and the x data) as well as an x value (i) and a boolean - extrapolate.

for(i=0; i<256; i=i+1)

{

number yspline=fastcubicspline(arrayimg, i,1)

setpixel(splineimg,i,0,yspline)

}