#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3		// Use modern global access method and strict wave access.
#pragma version=1.2

// Version 1.2, 14/4/17
// Adds smoothing and base level options

// Arc Hull is a baseline approximation for spectral data.
// A concave-upward arc is defined by a circle centered at the 
// midpoint of the spectrum. The depth of the arc is adjustable. 
// The arc is added to the spectrum, and the lower portion of a 
// convex hull is calculated for the resultant spectrum. 
// The baseline consists of the sum of the arc and the convex hull.

// Usage

// Arc Hull requires that a spectrum is plotted in the top graph.
// Selecting Macros - Arc Hull will add a panel on the left side 
// of the graph. You can choose which of the spectra in the graph
// to use for the calculation in a popup menu. A working copy of 
// the baseline and baseline-subracted spectrum are appended to the 
// plot; these are updated as you adjust the Arc Hull parameters.

// Adjust arc depth to find an optimal value for your data.
// The arc depth SetVar increments by 10% of its value.
// Setting arc depth to zero creates a convex hull (a.k.a. 'rubber band')
// type baseline. 

// Smooth: the baseline is calculated for a smoothed copy
// of the spectrum, then subtracted from the original. With the right 
// choice of smoothing factor, the baseline will pass through spectral 
// noise rather than underlying the noise. Note that the baseline-subtracted
// output is not smoothed with respect to the original data, it's only a 
// temporary copy of the spectrum that's smoothed and discarded after 
// calculating the baseline. 

// Base level: the output (baseline-subtracted) spectrum is offset from 
// zero by this amount. Useful for negative peaks that extend down from
// a non-zero value. Use negative arc depth for such spectra (you'll have
// to type a negative value in the setvar).

// Subtract: creates a copy of the baseline and a baseline-subtracted 
// output spectrum in the current data folder and appends them to the plot.

// All in one: creates a baseline-subtracted output spectrum for each 
// wave on the graph (other than waves with _base and _sub suffixes) using
// the current settings for arc depth, smoothing and base level.

// Closing the Arc Hull panel will remove temporary package-created
// waves from the graph and clean up after the package, as will 
// initializing on another graph window. 

Menu "Analysis"
	Submenu "Packages"
		"Arc Hull",/Q, CCH_init()
	end
end

Menu "Macros"	
	"Arc Hull",/Q, CCH_init()	
end

function CCH_clear()	
	
	GetWindow /Z CCH_panel  activeSW 
	string s_graph=parseFilepath(0, s_value, "#", 0, 0)
	
	KillWindow /Z CCH_panel
	if (strlen(winlist(s_graph,";","")))
		removefromgraph /W=$s_graph /Z Arc_Hull_Base, Arc_Hull_Sub
	endif
end

function CCH_init()
	
	if (strlen(WinList("*",";","WIN:1"))==0)  // no graphs 
		doalert 0, "Arc hull baseline requires a trace plotted in a graph window."
		return 0
	endif
	
	// make sure the top graph is visible
	string s_graph=WinName(0,1)
	dowindow /F $s_graph
	
	// clear any package detritus from the last used arc hull graph
	CCH_clear()
	
	// initialize (or reinitialize) the package data folder
	NewDataFolder /O root:Packages
	NewDataFolder /O root:Packages:CCHbaseline
	
	// Create a data folder reference variable
	DFREF dfr = root:Packages:CCHbaseline
	
	// make sure there's space to the left of the graph for the panel
	getwindow /Z $s_graph wsize
	if (V_left<200)
		movewindow /W=$s_graph 200, V_top, -1, -1
	endif
			
	// make panel
	NewPanel /K=1/N=CCH_panel/W=(200,0,0,0)/HOST=$s_graph/EXT=1 as "Arc Hull Controls"	
	ModifyPanel /W=$s_graph#CCH_panel, noEdit=1	
	
	string s_trace=stringfromlist(0, CCH_traces())
	wave w=TraceNameToWaveRef(s_graph, s_trace)
	variable depth=(wavemax(w)-wavemin(w))*.15
	
	// store values internally in these controls
	PopupMenu popTrace, mode=1, Value=CCH_traces(), title="",pos={30,20},size={130,20}
	PopupMenu popTrace, help={"select data wave" }, proc=CCH_popup
	SetVariable setvarDepth,pos={30,60},size={130,16},title="Arc depth"
	SetVariable setvarDepth,limits={-inf,inf,abs(depth/10)},value=_NUM:depth
	SetVariable setvarDepth,help={"depth of arc"}, fsize=14, proc=CCH_SetVar	
	SetVariable setvarSmooth,pos={30,90},size={130,16},title="Smooth"
	SetVariable setvarSmooth,limits={0,32767,1},value=_NUM:0
	SetVariable setvarSmooth,help={"binomial smoothing factor"}, fsize=14, proc=CCH_SetVar
	SetVariable setvarOffset,pos={30,120},size={130,16},title="Base level"
	SetVariable setvarOffset,limits={-inf,inf,1},value=_NUM:0
	SetVariable setvarOffset,help={"Offset from zero"}, fsize=14, proc=CCH_SetVar
	Button CCH_sub,pos={50,160},size={80,20},title="Subtract", proc=CCH_buttons
	Button CCH_sub,help={"subtract arc hull baseline"}
	Button CCH_all,pos={50,190},size={80,20},title="All in one", proc=CCH_buttons
	Button CCH_all,help={"subtract arc hull baseline from all traces"}	
	
	// set hook to intercept killing of panel 
	SetWindow $s_graph#CCH_panel, hook(CCH_Hook) = CCH_PanelHook
	
	if (strlen(s_trace)==0)
		return 0
	endif
	
	// do fit to create baseline
	CCH_doFit(s_trace)
	
	wave /SDFR=dfr w_base, w_sub	
	wave w_x=XWaveRefFromTrace(s_graph, s_trace) 
	removefromgraph /W=$s_graph /Z Arc_Hull_Base, Arc_Hull_Sub
	if(waveexists(w_x))
		appendtograph /W=$s_graph w_base/TN=Arc_Hull_Base vs w_x
		appendtograph /W=$s_graph w_sub/TN=Arc_Hull_Sub vs w_x
	else
		AppendToGraph /W=$s_graph w_base/TN=Arc_Hull_Base, w_sub/TN=Arc_Hull_Sub
	endif
	
end

Function CCH_SetVar(s) : SetVariableControl
	STRUCT WMSetVariableAction &s
	
	if(stringmatch(s.ctrlName,"setvarDepth"))	
		if (s.dval==0)
			SetVariable  setvarDepth win=$s.win,limits={-inf,inf,1}
		else
			// reset increment value to 10% of current value
			SetVariable  setvarDepth win=$s.win,limits={-inf,inf,abs(s.dval/10)}
		endif
	endif
	
	controlinfo /W= $s.win popTrace	
	string s_graph=parseFilepath(0, s.win, "#", 0, 0)
	CCH_doFit(s_value)
	
	return 0
End

Function CCH_popup(s)
	STRUCT WMPopupAction &s
	
	string s_graph=parseFilepath(0, s.win, "#", 0, 0)
	
	wave w=TraceNameToWaveRef(s_graph, s.popStr)
	if (waveexists(w)==0)
		return 0
	endif
	
	// set default depth each time a trace is selected
	variable depth=(wavemax(w)-wavemin(w))*0.15
	SetVariable setvarDepth,win=$s.win,limits={-inf,inf,abs(depth/10)},value=_NUM:depth
		
	CCH_doFit(s.popStr)		
	
	DFREF dfr = root:Packages:CCHbaseline	
	wave /SDFR=dfr w_base, w_sub
	
	removefromgraph /W=$s_graph /Z Arc_Hull_Base, Arc_Hull_Sub
	wave w_x=XWaveRefFromTrace(s_graph, s.popstr) 
	
	if(waveexists(w_x))
		appendtograph /W=$s_graph w_base/TN=Arc_Hull_Base vs w_x
		appendtograph /W=$s_graph w_sub/TN=Arc_Hull_Sub vs w_x
	else
		AppendToGraph /W=$s_graph w_base/TN=Arc_Hull_Base, w_sub/TN=Arc_Hull_Sub
	endif	
			
	return 0
End

Function CCH_buttons(s) : ButtonControl
	STRUCT WMButtonAction &s

	if(s.eventCode!=2)
		return 0
	endif
	
	string s_graph=parseFilepath(0, s.win, "#", 0, 0)	
	string s_trace
		
	strswitch(s.ctrlName)
		case "CCH_sub":	
			controlinfo /W=$s.win popTrace
			wave w=TraceNameToWaveRef(s_graph, s_value)
			if (waveexists(w)==0)
				return 0
			endif
			CCH_subtract(w)
			
			//CCH_clear()	
			break
		case "CCH_all":
			string msg="subtract arc hull baseline from all traces using current arc depth?\r"
			msg+="existing baselines will be overwritten!"
			doalert 1, msg
			if(v_flag==2)
				return 0
			endif
			string s_list=CCH_traces()
			variable i
			for(i=0;i<itemsinlist(s_list); i+=1)
				s_trace=stringfromlist(i, s_list)
				CCH_doFit(s_trace)
				wave w=TraceNameToWaveRef(s_graph, s_trace)
				CCH_subtract(w, overwrite=1)
			endfor
			
			break
	endswitch

	return 0
End

function /s CCH_traces()
	
	GetWindow /Z CCH_panel  activeSW 
	string s_graph=parseFilepath(0, s_value, "#", 0, 0)
	string s_list=TraceNameList(s_graph,";",1+4)
	variable i=0
	do
		if(stringmatch(stringfromlist(i,s_list), "*_base"))
			s_list=RemoveListItem(i, s_list)
		elseif(stringmatch(stringfromlist(i,s_list), "*_sub"))
			s_list=RemoveListItem(i, s_list)	
		else
			i+=1
		endif
	while (i<itemsinlist(s_list))
	return s_list
end

function CCH_doFit(s_trace)
	string s_trace
	
	GetWindow /Z CCH_panel activeSW 
	string s_graph=parseFilepath(0, s_value, "#", 0, 0)
	
	DFREF dfr = root:Packages:CCHbaseline
	
	variable depth, varSm, varOffset
	controlinfo /W=CCH_panel setvarDepth
	depth=v_value
	controlinfo /W=CCH_panel setvarSmooth
	varSm=v_value
	controlinfo /W=CCH_panel setvarOffset
	varOffset=v_value
	
	wave w=TraceNameToWaveRef(s_graph, s_trace)
	if (waveexists(w)==0)
		return 0
	endif
	wave w_x=XWaveRefFromTrace(s_graph, s_trace)
	
	// create waves to display baseline and baseline subtraction
	duplicate /o w dfr:w_sub /WAVE=w_sub
	duplicate /o w dfr:w_base /WAVE=w_base	
	
	if(varSm>0)
		smooth varSm, w_sub
		// make a copy of the smoothed data wave
		duplicate /free w_sub w_ref
		// calculate arc hull based on smoothed data
		ArcHull(w_sub, depth, w_x=w_x)
		w_base=w_ref-w_sub
		w_sub=w-w_base+varOffset
	else
		ArcHull(w_sub, depth, w_x=w_x)
		w_base=w-w_sub
		w_sub+=varOffset
	endif
	
	note /k w_base
	string cmd
	sprintf cmd, "Arc hull baseline with arc depth = %g, smoothing = %g", depth, varSm
	note w_base cmd
end

// subtracts arc hull baseline from w (w is overwritten)
function ArcHull(w, depth [, w_x])
	wave w, w_x
	variable depth

	variable useXwave=0
	if (paramisdefault(w_x)==0 && waveexists(w_x))
		useXwave=1
	endif
	
	// add concave function
	wavestats /Q w

	variable midX		
	if(useXwave)
		midX=(wavemax(w_x)-wavemin(w_x))/2
		w+=depth*(w_x[p]-midX)^2/(midX-w_x[0])^2
	else
		midX=leftx(w)+deltax(w)*(numpnts(w)/2-1)
		w+=depth*(x-midX)^2/(midX-leftx(w))^2
	endif
	
	duplicate /free w w_base
	
	if (useXwave==0)
		duplicate /free w w_x
		w_x=x
	endif
	
	Convexhull  w_x,w
	wave w_XHull=w_XHull, w_YHull=w_YHull
	if (depth>=0)
		reverse /P w_XHull, w_YHull
	endif
	// negative depth will subtract top part of convex hull
	// for negative peaks. An offset will likely need to be applied.
	
	wavestats /q w_XHull
	rotate -v_minloc, w_XHull, w_YHull
	setscale /p x, 0, 1, w_XHull, w_YHull
	wavestats /q w_XHull
	deletepoints v_maxloc+1, numpnts(w_XHull)-v_maxloc-1, w_XHull, w_YHull
	
	if(useXwave)
		Interpolate2/T=1/Y=w_base/X=w_x/I=3 W_XHull,W_YHull
	else
		Interpolate2 /T=1/Y=w_base/I=3 W_XHull,W_YHull
	endif

	w-=w_base
end

Function CCH_PanelHook(s)
	STRUCT WMWinHookStruct &s
	
	strswitch (s.eventName)
		case "kill":				
			string s_graph=parseFilepath(0, s.winName, "#", 0, 0)	
			removefromgraph /W=$s_graph /Z Arc_Hull_Base, Arc_Hull_Sub
			return 1
	endswitch
	return 0
end

function CCH_subtract(w [,overwrite])
	wave w
	variable overwrite
	
	if (paramisdefault(overwrite))
		overwrite=0
	endif
	
	GetWindow /Z CCH_panel  activeSW 
	string s_graph=parseFilepath(0, s_value, "#", 0, 0)
	
	
	DFREF dfr = root:Packages:CCHbaseline	
	wave /SDFR=dfr w_sub, w_base
	
	// save a copy of the baseline
	string strNewName=CleanupName(nameofwave(w)+"_base",0)
	if (overwrite==0 && exists(strNewName))
		doalert 1, strNewName+" exists. Overwrite?"
		if(V_flag==2)
			return 0
		endif
	endif
	duplicate /o w $strNewName
	wave newbase= $strNewName
	newbase=w_base
	
	// subtract baseline
	strNewName=CleanupName(nameofwave(w)+"_sub",0)
	if (overwrite==0 && exists(strNewName))
		doalert 1, strNewName+" exists. Overwrite?"
		if(V_flag==2)
			return 0
		endif
	endif
	duplicate /o w $strNewName
	wave subtracted= $strNewName
	subtracted=w_sub
	
	note subtracted, note(w_base) // append note from baseline wave to output wave note
	note newbase, note(w_base) 
	
	string s_list=WaveList("*", ";", "WIN:"+s_graph)
	wave w_x=XWaveRefFromTrace(s_graph, nameofwave(w))
	if(WhichListItem(nameofwave(newbase),s_list)==-1)	
		if(waveexists(w_x))
			appendtograph /W=$s_graph /C=(0,0,0) newbase vs w_x
		else
			appendtograph /W=$s_graph /C=(0,0,0) newbase
		endif
	endif
	if(WhichListItem(nameofwave(subtracted),s_list)==-1)
		if(waveexists(w_x))
			appendtograph /W=$s_graph /C=(0,0,0) subtracted vs w_x
		else
			appendtograph /W=$s_graph /C=(0,0,0) subtracted
		endif
	endif
end