Akima spline interpolation

 
#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
Here is an alternate implementation based on the same reference. It has been tested against the published test cases for the algorithm.

First, call calcIota to generate interpolation information; then you can interpolate using Akima's spline method with the akima() function.

#pragma rtGlobals=3 
 
// 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

Forum

Support

Gallery

Igor Pro 9

Learn More

Igor XOP Toolkit

Learn More

Igor NIDAQ Tools MX

Learn More