#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3
#pragma version=1.06
#pragma IgorVersion=7
#pragma moduleName=encircle

#include <Readback ModifyStr>

// --------------------- Project Updater header ----------------------
// If you're using Igor Pro 7 or later and have Project Updater 
// installed, this package can check periodically for new releases. 
// http: www.igorexchange.com/project/Updater
static constant kCHECK_FOR_UPDATES=1  // set to 1 to check for updates on init

static function CheckForUpdates()
	if(kCHECK_FOR_UPDATES)
		#if (exists("Updater#UpdateCheck")==6)
			return Updater#UpdateCheck("fileloc:"+FunctionPath(""))
		#endif
	endif
	return 0
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" to encicle all of the data in the wave.
// Right click within marquee and select a wavename to encircle a subset 
// of data within that wave. 
// 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
// It makes most sense to apply this to x-y data.

// Version history
// 1.06 allows selection of a subset of points using the marquee,
//      but still limited to one encircled group of points per wave
// 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

// Wishlist
// Allow user to select a region and encircle only the data that plot 
// within that region 
// could use GraphWaveDraw to select a group of points

Menu "TracePopup"
	"Encircle points", /Q, encircle#DrawBoundingSpline()
	"Remove encircled region", /Q, encircle#RemoveBoundingSpline()
end

Menu "GraphMarquee", dynamic
	submenu "Encircle points"
		encircle#TracesInMarquee(), /Q, encircle#DrawBoundingSpline()
	end
end

static function DrawBoundingSpline()
	
	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
	
	variable isMarquee=0
	
	// figure out graph and trace names
	GetLastUserMenuInfo
	if (cmpstr(s_value, "Encircle points"))
		isMarquee=1
		S_traceName=s_value				
	endif
	
	string s_info=traceinfo(S_graphName, S_traceName, 0)
	
	if(numtype(NumberByKey("AXISZ",s_info))==0)
		doalert 1, "you've selected a contour trace. Continue?"
		if (v_flag==2)
			return 0
		endif
	endif
	
	variable xOffset=GetNumFromModifyStr(s_info,"offset","{",0)
	variable yOffset=GetNumFromModifyStr(s_info,"offset","{",1)
	
	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)
	
	variable i
	if(isMarquee)
		make /free/n=0 xPoly, yPoly
		string StrXAxisName=stringbykey("XAXIS", TraceInfo("", S_traceName, 0 ))
		string StrYAxisName=stringbykey("YAXIS", TraceInfo("", S_traceName, 0 ))
		getmarquee /K $StrXAxisName, $StrYAxisName	
		xPoly={v_left,v_left,v_right,v_right}
		yPoly={v_top, v_bottom, v_bottom, v_top}
		xPoly-=xOffset
		yPoly-=yOffset
		findPointsInPoly xdata, ydata, xPoly, yPoly
		wave W_inPoly=W_inPoly
		make /free/n=0 x_subset, y_subset
		for(i=0;i<numpnts(ydata);i+=1)
			if(W_inPoly[i])
				x_subset[numpnts(x_subset)]={xdata[i]}
				y_subset[numpnts(y_subset)]={ydata[i]}
			endif
		endfor
		wave xdata=x_subset
		wave ydata=y_subset
	endif
	
	// 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+=xOffset
	W_YHull+=yOffset
	
	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 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=1.2 // expansion factor
	variable hscale=1+(v_expand-1)*(vpixels/hpixels)/(1+vpixels/hpixels)
	variable vscale=1+(v_expand-1)/(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

static 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.
static 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

static 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


// Traces plotted as x-y are returned when points lie within marquee.
// Waveform traces are returned when they pass through the marquee.
// For most uses the distinction is inconsequential.
// Designed for use in a dynamic menu - will run when menu is built and 
// after compilation 
static function /s TracesInMarquee()
	
	getmarquee /Z
	if(V_flag==0)
		return ""
	endif
	
	string inlist=TraceNameList(S_marqueeWin,";",5) // excludes hidden traces
	string outlist=""
	string s_trace, yAxisName, xAxisName, s_info
	variable xOffset, yOffset, xMult, yMult
	variable i, ntraces=itemsinlist(inlist)
	for(i=0;i<ntraces;i+=1)
		s_trace=stringfromlist(i,inlist)
		wave /Z w=TraceNameToWaveRef(S_marqueeWin, s_trace)
		s_info=traceinfo(S_marqueeWin,s_trace,0)
		xAxisName= StringByKey("XAXIS", s_info)
		yAxisName= StringByKey("YAXIS", s_info)
		
		getmarquee /Z $xAxisName, $yAxisName
		
		// deal with reversed axes by making v_top > v_bottom
		// and v_right > v_left
		variable swap=max(v_bottom, v_top)
		v_bottom=min(v_bottom, v_top)
		v_top=swap	
		swap=max(v_left, v_right)
		v_left=min(v_left, v_right)
		v_right=swap	
		
		// get trace offsets
		xOffset=GetNumFromModifyStr(s_info,"offset","{",0)
		yOffset=GetNumFromModifyStr(s_info,"offset","{",1)
		
		v_left-=xOffset
		v_right-=xOffset
		v_top-=yOffset
		v_bottom-=yOffset
		
		xMult=GetNumFromModifyStr(s_info,"muloffset","{",0)
		yMult=GetNumFromModifyStr(s_info,"muloffset","{",1)
		
		if(xMult)
			v_right/=xMult
			v_left/=xMult
		endif
		
		if(yMult)
			v_top/=yMult
			v_bottom/=yMult
		endif
		
		variable x_low, x_high, y_min, y_max 
			
		wave /z w_x=XWaveRefFromTrace(S_marqueeWin, s_trace)
		if (waveexists(w_x)) // deal with x-y data
			// be careful not to choke on category plot traces etc.
			if((wavetype(w_x)==0) || (numpnts(w_x)==0))
				continue
			endif
			make /free /n=4 xPoly, yPoly
			xPoly={v_left,v_left,v_right,v_right}
			yPoly={v_top,v_bottom,v_bottom,v_top}
			findPointsInPoly w_x, w, xPoly, yPoly
			
			wave W_inPoly
			if(wavemax(W_inPoly)==0)
				continue
			endif		
		else
			// check for waveform waves plotted left or right of marquee
			x_low=min(leftx(w), pnt2x(w,numpnts(w)-1))
			x_high=max(leftx(w), pnt2x(w,numpnts(w)-1))
			if (x_low>v_right || x_high<v_left)
				continue
			endif		
			// set y_max and y_min to max and min of points plotted  
			// within the x-axis range of marquee
			y_max=wavemax(w, v_left, v_right)
			y_min=wavemin(w, v_left, v_right)
			
			if (y_min>v_top || y_max<v_bottom)
				// wave is plotted above or below marquee
				continue
			endif
		endif

		outlist=addlistitem(stringfromlist(i,inlist), outlist)	
	endfor
	
	return outlist
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.
static 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