#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3
#pragma version=1.05

#include <Readback ModifyStr>

static constant kCHECK_FOR_UPDATES=1  // set to 1 to check for updates on init
static constant kFREQUENCY=604800 // (seconds). updater won't check more frequently than this. weekly=604800.
static strconstant ksVERSION="1.05"
static strconstant ksPROJECT="Encircle"
static strconstant ksLOCATION="http://www.igorexchange.com/node/7895/release"
static strconstant ksSHORTNAME="Encircle"
static strconstant ksPREFSFILE="EncirclePrefs.bin"

static function CheckForUpdates()
#if (exists("Updater#CheckForUpdates")==6)
	return Updater#CheckForUpdates(kCHECK_FOR_UPDATES,kFREQUENCY,ksVERSION,ksPROJECT,ksLOCATION,ksSHORTNAME,ksPREFSFILE)
#endif
end

// by tony.withers@uwo.ca

// Description
// draws a polgon around a cloud of data points on a plot.

// How to use
// Right click on a data point to invoke the trace menu and select 
// "Encircle points".
// If you reshape the graph window you can repeat this procedure to 
// update the shape of the polygon.
// Selecting "Remove Encircled Region" deletes the shaded region.
// Polygons are drawn in the user back layer and can be edited with the 
// drawing tools.

// Details
// Polygon surrounds all of the data in the wave.
// It makes most sense to apply this to x-y data.

// Version history
// 1.03 if Package Updater is installed, will check IgorExchange for new releases 
// 1.02 deal with offsets
// 1.01 bug fix: remove gratuitous print statement

Menu "TracePopup"
	"Encircle points", /Q, DrawCatmullRomBoundingSpline()
	"Remove encircled region", /Q, RemoveBoundingSpline()
end

function DrawCatmullRomBoundingSpline()
	
	if (CheckForUpdates())
		return 0 // this file has been updated; quit to allow user to reload
	endif
	
	// initialize (or reinitialize) the package data folder
	NewDataFolder /O root:Packages
	NewDataFolder /O root:Packages:Encircle
	
	// Create a data folder reference variable
	DFREF dfr = root:Packages:Encircle
	
	// figure out graph and trace names
	GetLastUserMenuInfo
	string s_info=traceinfo(S_graphName, S_traceName, 0)
	//print s_info
	if(numtype(NumberByKey("AXISZ",s_info))==0)
		doalert 1, "you've selected a contour trace. Continue?"
		if (v_flag==2)
			return 0
		endif
	endif	
	
	if(strlen(stringbykey("XWAVE",s_info))==0)
		duplicate /free tracenametowaveref(S_graphName, S_traceName), xdata
		xdata=x
	else
		wave xdata=XWaveRefFromTrace(S_graphName, S_traceName)
	endif	
	wave ydata=tracenametowaveref(S_graphName, S_traceName)
	
	// number of points to calculate between nodes
	variable nPoints=20	
	variable alpha=0.5
	// alpha=0 for standard (uniform) Catmull-Rom spline
	// alpha=0.5 for centripetal Catmull-Rom spline
	// alpha=1 for chordal Catmull-Rom spline
	
	if (numpnts(ydata)<3)
		doalert 0, "must have more than two data points"
		return 0
	endif
	
	Convexhull xdata,ydata
	wave W_XHull=W_XHull, W_YHull=W_YHull
	
	// deal with offsets
	W_XHull+=GetNumFromModifyStr(s_info,"offset","{",0)
	W_YHull+=GetNumFromModifyStr(s_info,"offset","{",1)
	
	variable vertices=numpnts(W_XHull)	
	
	make /free/n=(vertices+3) x_hull, y_hull
	x_hull=W_XHull[p-(vertices)*(p>(vertices-1))]
	y_hull=W_YHull[p-(vertices)*(p>(vertices-1))]	

	// find the centroid of the hull
	variable x0=0, y0=0, varA=0
	variable i, factor
	for(i=0;i<vertices; i+=1)
		factor=x_hull[i]*y_hull[i+1] - x_hull[i+1]*y_hull[i]
		x0+=(x_hull[i]+x_hull[i+1])*factor
		y0+=(y_hull[i]+y_hull[i+1])*factor
		varA+=factor
	endfor
	varA*=0.5
	x0/=6*varA
	y0/=6*varA
	
	string newXname, newYname
	newXname=cleanupname(S_traceName, 0)+"_xCR"
	newYname=cleanupname(S_traceName, 0)+"_yCR"

	// could end up overwriting a different trace's shaded region
	// if you use long or liberal wave names	
	make /o/n=0 dfr:$newXname /WAVE=CR_x, dfr:$newYname /WAVE=CR_y	 
	make /free /n=(nPoints) x_out, y_out
	make /free /n=4 x_in, y_in	
	
	for(i=1;i<=vertices;i+=1)
		x_in=x_hull[i-1+p]
		y_in=y_hull[i-1+p]
		CatmullRomSegment(x_in, y_in, x_out, y_out, alpha)
		concatenate /NP {x_out}, CR_x
		concatenate /NP {y_out}, CR_y
	endfor
	// wrap around all the way back to first point
	CR_x[numpnts(CR_x)]={CR_x[0]}
	CR_y[numpnts(CR_y)]={CR_y[0]}
	
	// add a center point at the start of polygon
	// so that we can scale outwards
	insertpoints 0, 2, CR_x, CR_y
	CR_x[0]={x0,nan}
	CR_y[0]={y0,nan}	
		
	// need to know the aspect ratio of the window...
	getwindow /Z $S_graphName psize
	variable hpixels=V_right-V_left
	variable vpixels=V_bottom-V_top
	
	string s_Xax=stringbykey("XAXIS",s_info)
	string s_Yax=stringbykey("YAXIS",s_info)
	getaxis /Q $s_Xax
	variable xaxrange=abs(V_max-V_min)
	getaxis /Q $s_Yax
	variable yaxrange=abs(V_max-V_min)
	variable xrange=(wavemax(x_hull)-wavemin(x_hull))
	variable yrange=(wavemax(y_hull)-wavemin(y_hull))
	
	// ... to figure out the aspect ratio of the bounding shape 
	hpixels=xrange*hpixels/xaxrange // horizontal spread of data in pixels
	vpixels=yrange*vpixels/yaxrange // vertical spread of data in pixels
	
	// expand a bit to inset symbols
	variable v_expand=0.2 // expansion factor
	variable hscale=1+v_expand*(vpixels/hpixels)/(1+vpixels/hpixels)
	variable vscale=1+v_expand/(1+vpixels/hpixels)
		
	// prepare to draw the polygon
	string objName=cleanupname(S_traceName, 0)+"loop"	
	DrawAction /W=$S_graphName/L=UserBack getgroup=$objName, delete
	
	SetDrawLayer /W=$S_graphName UserBack
	SetDrawEnv /W=$S_graphName gstart,gname=$objName	
	SetDrawEnv /W=$S_graphName xcoord=$s_Xax,ycoord=$s_Yax
	
	variable v_opacity=4000
	variable v_thickness=2
	
	string s_color=""
	Variable v1,v2
	v1= strsearch(s_info, "rgb(x)=", 0)
	if( v1 >= 0 )
		v1+=7 // characters in key
		v2=strsearch(s_info, ")", v1)	
		s_color=s_info[v1,v2]
		variable red, green, blue
		sscanf s_color, "(%d,%d,%d)", red, green, blue
		SetDrawEnv /W=$S_graphName linefgc= (red, green, blue),fillfgc= (red, green, blue, v_opacity)
	endif
	SetDrawEnv /W=$S_graphName fillpat= 1,linethick= v_thickness
	drawpoly /W=$S_graphName x0, y0, hscale, vscale, CR_x, CR_y
	SetDrawEnv /W=$S_graphName gstop,gname=$objName
end

function RemoveBoundingSpline()
	GetLastUserMenuInfo
	string objName=cleanupname(S_traceName, 0)+"loop"	
	
	DrawAction /L=UserBack /W=$S_graphName getgroup=$objName, delete
	DFREF dfr = root:Packages:Encircle
	string basename=cleanupname(S_traceName, 0)
	killwaves /Z dfr:$basename+"_xCR", dfr:$basename+"_yCR"	
end

// Create a spline that passes through a series of x-y points,
// looping back to the start point.
// Can be used with ConvexHull to make a loop around a cloud of points.
function CatmullRomSplineLoop(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=(np+2) x_temp, y_temp
	x_temp=xdata[p-(np-1)*(p>(np-1))]
	y_temp=ydata[p-(np-1)*(p>(np-1))]	
	make  /free /n=(segPoints) x_out, y_out
	make /free /n=4 x_in, y_in	
	
	variable i
	for(i=1;i<np;i+=1)
		x_in=x_temp[i-1+p]
		y_in=y_temp[i-1+p]
		CatmullRomSegment(x_in, y_in, x_out, y_out, alpha)
		concatenate /NP {x_out}, CR_x
		concatenate /NP {y_out}, CR_y
	endfor
	// wrap around all the way back to first point
	CR_x[numpnts(CR_x)]={CR_x[0]}
	CR_y[numpnts(CR_y)]={CR_y[0]}
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

// this function isn't needed for drawing shaded regions.

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