Catmull Rom spline

Average rating
(0 votes)

// 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

Back to top