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