Akima spline interpolation
andyfaff
Wed, 11/04/2015 - 03:48 pm
#pragma rtGlobals=3 // Use modern global access method and strict wave access.
// This code was translated into IGOR from the Python code available at
// http://www.lfd.uci.edu/~gohlke/code/akima.py.html
// Copyright (c) 2007-2015, Christoph Gohlke
// Copyright (c) 2007-2015, The Regents of the University of California
// Produced at the Laboratory for Fluorescence Dynamics
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
// * Neither the name of the copyright holders nor the names of any
// contributors may be used to endorse or promote products derived
// from this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//Akima's interpolation method uses a continuously differentiable sub-spline
//built from piecewise cubic polynomials. The resultant curve passes through
//the given data points and will appear smooth and natural.
//
//:Author:
// `Christoph Gohlke <http://www.lfd.uci.edu/~gohlke/>`_
//
//:Organization:
// Laboratory for Fluorescence Dynamics, University of California, Irvine
//
//References
//----------
//(1) A new method of interpolation and smooth curve fitting based
// on local procedures. Hiroshi Akima, J. ACM, October 1970, 17(4), 589-602.
Function akima(x, y, x_new, y_new)
// Return interpolated data using Akima's method.
//
// Parameters
// ----------
// x_new : wave
// 1D array of monotonically increasing real values.
// y : wave
// N-D array of real values. y's length along the interpolation
// axis must be equal to the length of x.
// x_new : wave
// New independent variables.
// y_new : wave
// Wave to receive results. The number of points must be the same as x_new
// make/d/n=2 y_new
// akima({0, 1, 2}, {0, 0, 1}, {0.5, 1.5}, y_new)
// print y_new
// y_new[0]= {-0.125,0.375}
Wave x, y, x_new, y_new
variable n, mm, mmm, mp, mpp, ii, wm
duplicate/free x_new, x_i
n = numpnts(x)
make/d/n=(numpnts(x) - 1)/free dx, dy, m, m1
dx = x[p + 1] - x[p]
dy = y[p + 1] - y[p]
m = dy / dx
mm = 2 * m[0] - m[1]
mmm = 2 * mm - m[0]
mp = 2 * m[n - 2] - m[n - 3]
mpp = 2 * mp - m[n - 2]
m1 = m
insertpoints 0, 2, m1
m1[0] = mmm
m1[1] = mm
insertpoints numpnts(m1), 2, m1
m1[numpnts(m1) - 2] = mp
m1[numpnts(m1) - 1] = mpp
make/d/free/n=(numpnts(m1) - 1) dm
dm = abs(m1[p + 1] - m1[p])
duplicate/r=[2, n + 1]/free dm, f1
duplicate/r=[0, n - 1]/free dm, f2
duplicate/free f1, f12
f12 = f1 + f2
duplicate/free f12, ids
ids = f12[p] > 1e-9 * wavemax(f12) ? p : -1
duplicate/free/r=[1, n] m1, b
b = m1[p + 1]
b = ids[p] >= 0 ? (f1[p] * m1[p + 1] + f2[p] * m1[p + 2]) / f12[p]: 0
duplicate/free m, c, d
c = (3.0 * m[p] - 2.0 * b[p] - b[p + 1]) / dx[p]
d = (b[p] + b[p + 1] - 2.0 * m) / dx ^ 2
duplicate/free x_i, bins
bins = binarysearch(x, x_i)
bins = min(bins[p], n - 2)
duplicate/free/r=[0, numpnts(x_i) - 1] bins, bb
duplicate/free x_i, wj
wj = x_i - x[bb]
duplicate/free bb, d_p, c_p, b_p
b_p = b[bb[p]]
c_p = c[bb[p]]
d_p = d[bb[p]]
y_new = ((wj * d_p + c_p) * wj + b_p) * wj + y[bb]
End
// This code was translated into IGOR from the Python code available at
// http://www.lfd.uci.edu/~gohlke/code/akima.py.html
// Copyright (c) 2007-2015, Christoph Gohlke
// Copyright (c) 2007-2015, The Regents of the University of California
// Produced at the Laboratory for Fluorescence Dynamics
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
// * Neither the name of the copyright holders nor the names of any
// contributors may be used to endorse or promote products derived
// from this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//Akima's interpolation method uses a continuously differentiable sub-spline
//built from piecewise cubic polynomials. The resultant curve passes through
//the given data points and will appear smooth and natural.
//
//:Author:
// `Christoph Gohlke <http://www.lfd.uci.edu/~gohlke/>`_
//
//:Organization:
// Laboratory for Fluorescence Dynamics, University of California, Irvine
//
//References
//----------
//(1) A new method of interpolation and smooth curve fitting based
// on local procedures. Hiroshi Akima, J. ACM, October 1970, 17(4), 589-602.
Function akima(x, y, x_new, y_new)
// Return interpolated data using Akima's method.
//
// Parameters
// ----------
// x_new : wave
// 1D array of monotonically increasing real values.
// y : wave
// N-D array of real values. y's length along the interpolation
// axis must be equal to the length of x.
// x_new : wave
// New independent variables.
// y_new : wave
// Wave to receive results. The number of points must be the same as x_new
// make/d/n=2 y_new
// akima({0, 1, 2}, {0, 0, 1}, {0.5, 1.5}, y_new)
// print y_new
// y_new[0]= {-0.125,0.375}
Wave x, y, x_new, y_new
variable n, mm, mmm, mp, mpp, ii, wm
duplicate/free x_new, x_i
n = numpnts(x)
make/d/n=(numpnts(x) - 1)/free dx, dy, m, m1
dx = x[p + 1] - x[p]
dy = y[p + 1] - y[p]
m = dy / dx
mm = 2 * m[0] - m[1]
mmm = 2 * mm - m[0]
mp = 2 * m[n - 2] - m[n - 3]
mpp = 2 * mp - m[n - 2]
m1 = m
insertpoints 0, 2, m1
m1[0] = mmm
m1[1] = mm
insertpoints numpnts(m1), 2, m1
m1[numpnts(m1) - 2] = mp
m1[numpnts(m1) - 1] = mpp
make/d/free/n=(numpnts(m1) - 1) dm
dm = abs(m1[p + 1] - m1[p])
duplicate/r=[2, n + 1]/free dm, f1
duplicate/r=[0, n - 1]/free dm, f2
duplicate/free f1, f12
f12 = f1 + f2
duplicate/free f12, ids
ids = f12[p] > 1e-9 * wavemax(f12) ? p : -1
duplicate/free/r=[1, n] m1, b
b = m1[p + 1]
b = ids[p] >= 0 ? (f1[p] * m1[p + 1] + f2[p] * m1[p + 2]) / f12[p]: 0
duplicate/free m, c, d
c = (3.0 * m[p] - 2.0 * b[p] - b[p + 1]) / dx[p]
d = (b[p] + b[p + 1] - 2.0 * m) / dx ^ 2
duplicate/free x_i, bins
bins = binarysearch(x, x_i)
bins = min(bins[p], n - 2)
duplicate/free/r=[0, numpnts(x_i) - 1] bins, bb
duplicate/free x_i, wj
wj = x_i - x[bb]
duplicate/free bb, d_p, c_p, b_p
b_p = b[bb[p]]
c_p = c[bb[p]]
d_p = d[bb[p]]
y_new = ((wj * d_p + c_p) * wj + b_p) * wj + y[bb]
End
Forum
Support
Gallery
Igor Pro 9
Learn More
Igor XOP Toolkit
Learn More
Igor NIDAQ Tools MX
Learn More
First, call calcIota to generate interpolation information; then you can interpolate using Akima's spline method with the akima() function.
// Akima.ipf: Routines to implement Akima-spline fitting, based on
// H. Akima, Journ. ACM, Vol 17, No 4, 1970 p 589-602
// M. Bongard, 11/17/09
#pragma IgorVersion=6.1 // Use of free waves
ThreadSafe Function calcIota(knotX, knotY[, dWave])
WAVE knotX // knot X locations
WAVE knotY // knot Y locations
WAVE dWave // Destination wave reference for iotas
if(!WaveExists(knotX) || !WaveExists(knotY))
Print "calcIota: ERROR -- requisite waves do not exist! Aborting..."
return -1
endif
if(numpnts(knotX) != numpnts(knotY))
Print "calcIota: ERROR -- knot waves must have same number of points! Aborting..."
return -1
endif
Variable numKnots = numpnts(knotX)
if(numKnots < 5)
Print "calcIota: ERROR -- Akima spline algorithm requires at least 5 knots. Aborting..."
return -1
endif
// Make intermediate ai, bi, mi arrays
Make/D/FREE/N=(numKnots + 4) kX, kY
Variable i, j
for(i = 2, j = 0; j < numKnots; i += 1, j += 1)
kX[i] = knotX[j]
kY[i] = knotY[j]
endfor
// Handle end-point extrapolation
Make/D/FREE/N=5 endX, endY
// RHS: end points are last three in dataset
Variable endStartPt
// RHS
endStartPt = numPnts(kX) - 5
endX = kX[p + endStartPt]
endY = kY[p + endStartPt]
calcEndPoints(endX, endY)
kX[numpnts(kX)-2] = endX[3]
kX[numpnts(kX)-1] = endX[4]
kY[numpnts(kX)-2] = endY[3]
kY[numpnts(kX)-1] = endY[4]
// LHS: end points are first three in dataset, but reversed in ordering
// (i.e. point 3 in Akima's notation == index 0)
endX = 0
endY = 0
for(i = 0, j = 2; i < 3; i += 1, j -= 1)
endX[j] = knotX[i]
endY[j] = knotY[i]
endfor
calcEndPoints(endX, endY)
kX[1] = endX[3]
kX[0] = endX[4]
kY[1] = endY[3]
kY[0] = endY[4]
// kX, kY are now properly populated, along with all necessary extrapolated endpoints
// computed as specified in Akima1970
Make/D/FREE/N=(numKnots + 4 - 1) mK
mK = (kY[p + 1] - kY[p]) / (kX[p + 1] - kX[p])
Make/O/N=(numKnots) knotIota
Variable denom, m1,m2,m3,m4
for(i = 2, j = 0; j < numKnots; i += 1, j += 1)
m1 = mK[i - 2]
m2 = mK[i - 1]
m3 = mK[i]
m4 = mK[i + 1]
denom = abs(m4 - m3) + abs(m2 - m1)
if(denom == 0)
knotIota[j] = 0.5 * (m2 + m3)
continue
endif
knotIota[j] = ( abs(m4 - m3)*m2 + abs(m2 - m1)*m3 ) / denom
endfor
if(!ParamIsDefault(dWave))
// Overwrite input destination wave with new iotas
Duplicate/O knotIota, dWave
Killwaves/Z knotIota
endif
End
// Given: 5-point knot wave knotX, knotY, with i=[0,2] representing the last three
// knot locations from data, compute end knots i=[3,4] appropriately.
ThreadSafe Function calcEndPoints(kX, kY)
WAVE kX, kY // knot X,Y coordinate locations, respectively
// Sanity checks
if( (numPnts(kX) != numpnts(kY)) || numpnts(kX) != 5)
Print "calcEndPoints: ERROR -- must have 5 points in knot wave! Aborting..."
return -1
endif
// First, compute X locations of knots, according to relations in eq. 8 of Akima1970:
kX[3] = kX[1] + kX[2] - kX[0]
kX[4] = 2*kX[2] - kX[0]
// Now all kX are known, so let's set up the line segment slope waves
// ai, bi, mi of eq. (12)-(14).
Make/N=4/FREE ai, bi, mi
ai = kX[p+1] - kX[p]
// ai is now determined completely
bi = kY[p + 1] - kY[p]
mi = bi/ai
// bi, mi determined on i=[0,1]
// Determine remainder of quantities by applying solutions of eq (9)
kY[3] = (2*mi[1] - mi[0])*(kX[3] - kX[2]) + kY[2]
mi[2] = (kY[3] - kY[2]) / (kX[3] - kX[2])
kY[4] = (2*mi[2] - mi[1])*(kX[4] - kX[3]) + kY[3]
mi[3] = (kY[4] - kY[3]) / (kX[4] - kX[3])
End
ThreadSafe Function akima(x, knotX, knotY, knotIota)
Variable x
WAVE knotX, knotY, knotIota
Variable i1, i2
// Find where x
Variable done = 0
i2 = -1
Variable numKnots = numpnts(knotX)
do
i2 += 1
if(knotX[i2] == x)
return knotY[i2]
endif
done = (knotX[i2] > x) ? 1 : 0
while(!done && (i2 < numKnots))
i1 = i2 - 1
Variable x1, x2, y1, y2, iota1, iota2
x1 = knotX[i1]
y1 = knotY[i1]
iota1 = knotIota[i1]
x2 = knotX[i2]
y2 = knotY[i2]
iota2 = knotIota[i2]
Variable p0, p1, p2, p3, tmp
p0 = y1
p1 = iota1
p2 = ( 3 * (y2 - y1)/(x2 - x1) - 2*iota1 - iota2 ) / (x2 - x1)
p3 = (iota1 + iota2 - 2 * (y2 - y1)/(x2 - x1)) / (x2 - x1)^2
tmp = x - x1
return p0 + p1 * tmp + p2 * tmp^2 + p3*tmp^3
End
November 4, 2015 at 07:11 pm - Permalink