Catmull Rom spline

// Create a spline that passes through a series of x-y points
// Catmull E. and Rom R. (1974) A class of local interpolating splines. 
// In Computer  Aided  Geometric Design, R.E. Barnhill and R.F. Reisenfeld, Eds. 
// Academic Press, New York, 1974, pp. 317–326.
function CatmullRomSpline(xdata, ydata [,segPoints, alpha])
    wave xdata, ydata   
    variable segPoints, alpha
    
    if( ParamIsDefault(segPoints) )
        segPoints=20
    endif
 
 
    if( ParamIsDefault(alpha) )
        alpha=0.5
    endif
    
    variable np=numpnts(ydata)  
 
    string newXname, newYname
    newXname=nameofwave(xdata)+"_CR"
    newYname=nameofwave(ydata)+"_CR"
    
    if(exists(newXname)||exists(newYname))
        doalert 1, "overwrite existing wave?"
        if (V_Flag==2)
            return 0
        endif
    endif
    
    make /o/n=0 $newXname /WAVE=CR_x, $newYname /WAVE=CR_y  
    make  /free /n=(segPoints) x_out, y_out
    make /free /n=4 x_in, y_in  
    
    variable i
    for(i=0;i<(np-1);i+=1)
        if (i==0) // first segment
            // reflect through end points
            x_in[0]=xdata[0]-(xdata[1]-xdata[0])
            x_in[1,3]=xdata[p-1]
            y_in[0]=ydata[0]-(ydata[1]-ydata[0])
            y_in[1,3]=ydata[p-1]
        elseif(i==(np-2)) // last segment
            x_in[0,2]=xdata[i-1+p]
            x_in[3]=xdata[np-1]+(xdata[np-1]-xdata[np-2])
            y_in[0,2]=ydata[i-1+p]
            y_in[3]=ydata[np-1]+(ydata[np-1]-ydata[np-2])
        else
            x_in=xdata[i-1+p]
            y_in=ydata[i-1+p]
        endif
        CatmullRomSegment(x_in, y_in, x_out, y_out, alpha)
        concatenate /NP {x_out}, CR_x
        concatenate /NP {y_out}, CR_y
    endfor
    // insert end point
    CR_x[numpnts(CR_x)]={xdata[np-1]}
    CR_y[numpnts(CR_y)]={ydata[np-1]}
end
 
function CatmullRomSegment(w_x, w_y, xout, yout, alpha)
    wave w_x, w_y // x and y values of four control points
    wave xout, yout // output waves
    variable alpha
    
    // alpha=0 for standard (uniform) Catmull-Rom spline
    // alpha=0.5 for centripetal Catmull-Rom spline
    // alpha=1 for chordal Catmull-Rom spline
    
    variable nPoints=numpnts(xout)
    variable t0, t1, t2, t3
 
    t0=0
    t1=sqrt( (w_x[1]-w_x[0])^2 + (w_y[1]-w_y[0])^2 )^alpha + t0 
    t2=sqrt( (w_x[2]-w_x[1])^2 + (w_y[2]-w_y[1])^2 )^alpha + t1
    t3=sqrt( (w_x[3]-w_x[2])^2 + (w_y[3]-w_y[2])^2 )^alpha + t2
    
    make /free/n=(nPoints) w_t,A1x,A1y,A2x,A2y,A3x,A3y,B1x,B1y,B2x,B2y
    
    w_t=t1+p/(nPoints)*(t2-t1)
    A1x=(t1-w_t)/(t1-t0)*w_x[0] + (w_t-t0)/(t1-t0)*w_x[1]
    A1y=(t1-w_t)/(t1-t0)*w_y[0] + (w_t-t0)/(t1-t0)*w_y[1]
    A2x=(t2-w_t)/(t2-t1)*w_x[1] + (w_t-t1)/(t2-t1)*w_x[2]
    A2y=(t2-w_t)/(t2-t1)*w_y[1] + (w_t-t1)/(t2-t1)*w_y[2]   
    A3x=(t3-w_t)/(t3-t2)*w_x[2] + (w_t-t2)/(t3-t2)*w_x[3]
    A3y=(t3-w_t)/(t3-t2)*w_y[2] + (w_t-t2)/(t3-t2)*w_y[3]   
    B1x=(t2-w_t)/(t2-t0)*A1x + (w_t-t0)/(t2-t0)*A2x
    B1y=(t2-w_t)/(t2-t0)*A1y + (w_t-t0)/(t2-t0)*A2y
    B2x=(t3-w_t)/(t3-t1)*A2x + (w_t-t1)/(t3-t1)*A3x
    B2y=(t3-w_t)/(t3-t1)*A2y + (w_t-t1)/(t3-t1)*A3y
    
    xout=(t2-w_t)/(t2-t1)*B1x + (w_t-t1)/(t2-t1)*B2x
    yout=(t2-w_t)/(t2-t1)*B1y + (w_t-t1)/(t2-t1)*B2y
end