# 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))
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

Forum

Support

Gallery