#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3
#pragma version=2.07
#pragma IgorVersion=6
#pragma moduleName=archull

#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

// This package adds a collection of algorithmically-defined baselines 
// for spectral data. It's named for the arc hull type baseline, in which
// a circular arc with adjustable depth is added to the spectrum (ie., 
// the spectrum is bent), 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. 

// How to use the Arc Hull package

// 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 from a popup menu which of the spectra in the 
// graph to use for the calculation. 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. 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. 

// Spline: baseline is a cubic spline using the convex hull as nodes. 
// When arc depth is non-zero the nodes are more like a subset of 
// alphashape vertices. Adjust Smooth and Depth parameters to produce an 
// optimal number of nodes.

// Iterative: the function selected in the Fit Func popup is fitted 
// iteratively. Any fit point with a negative residual (positive if the 
// negative peaks checkbox is selected) is replaced by the corresponding 
// point from the input data after each iteration, thus eliminating 
// points with a positive residual, so ultimately eliminating peaks from 
// the fit. The fit output is used as the input for the next iteration. 
// Lieber C.A. and Mahadevan-Jansen A. (2003) Applied Spectroscopy 57: 
// 1363-1367.

// 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, the smoothing is 
// applied only to a temporary copy that's 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. For convex baselines, use negative arc depth (you'll 
// have to type a negative value in the setvar).

// Negative Peaks checkbox: subtract the top part of the convex hull.
// Useful for dealing with downward facing peaks, or for when you have
// a reverse-scaled vertical axis.

// Use graph limits checkbox: calculate baseline within the horizontal 
// axis limits of the graph. Unless the entire wave can be seen in the 
// graph, the baseline-subtracted output wave will have a smaller range 
// than the data wave.

// 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 _BL and _sub suffixes) 
// using the current settings for arc depth, smoothing, base level and 
// checkbox values.

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

// For X-Y data, the X wave needs to be monatonic for Arc Hull calculation.

// Execute archull#demo() to demo.

// Version history
// 2.07 Bug fix: graph hook was updating fit needlessly.
// 2.06 Record fit coefficients for iterative fit in wave note
// 2.05 Made most functions static
// 2.04 Added headers for Project Updater.
// 2.03 11/11/17 Should be backward compatible for Igor 6.
// 2.02 decreased panel size for compatibility with smaller or lower 
// resolution screens. And a bug fix.
// 2.0 Introduces iterative fitting of CurveFit functions with 
// elimination of either positive or negative residual points
// Versions 1.3, 1.31, 1.32 Add option to create baseline within x-range 
// of graph and spline type baseline using convex hull as nodes
// Version 1.2, 14/4/17 Adds smoothing and base level options

// If "use graph limits" is selected, the baseline will update as the 
// plot area changes. You may want to disable this so that you can zoom 
// in or out and see the fit in a different context, then reenter a fit 
// parameter to update the fit. Set UPDATE_ON_ZOOM=0 to prevent this 
// behaviour:
static constant UPDATE_ON_ZOOM=1

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

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

static function CCH_init()
	
	if (CheckForUpdates())
		return 0 // this file has been updated; quit to allow user to reload
	endif
	
	if (strlen(WinList("*",";","WIN:1"))==0)  // no graphs 
		doalert 0, "Arc hull baseline requires a trace plotted in a graph window."
		return 0
	endif
	
	// 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:ArcHull
	
	// make control panel and append package waves to graph
	variable success	
	success=CCH_makePanel()	
	if (success)
		CCH_updatePlot()
	endif
end

static function CCH_makePanel()

	// make sure the top graph is visible
	string graphStr=WinName(0,1)
	dowindow /F $graphStr
	
	// Create a data folder reference variable
	DFREF dfr = root:Packages:ArcHull
	
   // make sure there's space to the left of the graph for the panel
	getwindow /Z $graphStr wsize
	if (v_left<200)
		#if IgorVersion() >= 7	
			movewindow /W=$graphStr 200, v_top, -1, -1
		#else
			movewindow /W=$graphStr 200, v_top, 200+(v_right-v_left), v_bottom
		#endif
	endif
			
	// make panel
	NewPanel /K=1/N=CCH_panel/W=(200,0,0,355)/HOST=$graphStr/EXT=1 as "Arc Hull Controls"	
	ModifyPanel /W=$graphStr#CCH_panel, noEdit=1	
	
	string traceStr=stringfromlist(0, CCH_traces())
	wave /Z w=TraceNameToWaveRef(graphStr, traceStr)
	variable depth=(wavemax(w)-wavemin(w))*.15
	
	variable i=0, deltaY=25, groupw=180, font=12
	v_left=30
	v_top=5
	
	GroupBox group0,pos={v_left-20,v_top+deltaY*i},size={groupw,deltaY*2},title="Data wave"
	GroupBox group0,fSize=font
	i+=1	
	// store values internally in these controls
	PopupMenu popTrace, mode=1, Value=archull#CCH_traces(), title="",pos={v_left,v_top+deltaY*i},size={130,20}
	PopupMenu popTrace, help={"select data wave" }, proc=archull#CCH_popup, fSize=font
	i+=1.5	
	GroupBox group1 pos={v_left-20,v_top+deltaY*i},size={groupw,deltaY*9},title="Baseline parameters"
	GroupBox group1 fSize=font	
	i+=1	
	Checkbox check_CCH, pos={v_left,v_top+deltaY*i}, fSize=font, value=1, proc=archull#CCH_checkbox, title="Arc hull"
	Checkbox check_CCH, mode=1,help={"Arc hull type baseline"}
	i+=1
	Checkbox check_spline, pos={v_left,v_top+deltaY*i}, fSize=font, value=0, proc=archull#CCH_checkbox, title="Hull spline"
	Checkbox check_spline, mode=1,help={"cubic spline through nodes from convex hull"}
	i+=1
	Checkbox check_itFit, pos={v_left,v_top+deltaY*i}, fSize=font, value=0, proc=archull#CCH_checkbox, title="Iterative"
	Checkbox check_itFit, mode=1,help={"iterative fit with elimination"}
	i+=1
	SetVariable setvarDepth,pos={v_left+10,v_top+deltaY*i},size={130,16},title="Arc depth"
	SetVariable setvarDepth,limits={-inf,inf,abs(depth/10)},value=_NUM:depth, bodyWidth=70
	SetVariable setvarDepth,help={"depth of arc"}, fSize=font, proc=archull#CCH_SetVar	
	
	PopupMenu popFitFunc, mode=2, title="Fit Func",pos={v_left,v_top+deltaY*i},size={130,20}
	PopupMenu popFitFunc, help={"select fit function" }, proc=archull#CCH_popup, fSize=font, disable = 1
	PopupMenu popFitFunc, Value="line;poly 3;poly 4;gauss;lor;exp;dblexp;sin;hillequation;sigmoid;power;lognormal;"	
	i+=1
	SetVariable setvarSmooth,pos={v_left+10,v_top+deltaY*i},size={130,16},title="Smooth"
	SetVariable setvarSmooth,limits={0,32767,1},value=_NUM:0, bodyWidth=70
	SetVariable setvarSmooth,help={"binomial smoothing factor"}, fSize=font, proc=archull#CCH_SetVar
	i+=1
	SetVariable setvarOffset,pos={v_left+10,v_top+deltaY*i},size={130,16},title="Base level"
	SetVariable setvarOffset,limits={-inf,inf,1},value=_NUM:0, bodyWidth=70
	SetVariable setvarOffset,help={"Offset from zero"}, fSize=font, proc=archull#CCH_SetVar
	i+=1
	CheckBox check_neg, pos={v_left,v_top+deltaY*i}, fSize=font, value=0, proc=archull#CCH_checkbox, title="Negative peaks"
	CheckBox check_neg, help={"Fit baseline to top of convex hull"}
	i+=1
	CheckBox check_win, pos={v_left,v_top+deltaY*i}, fSize=font, value=0, proc=archull#CCH_checkbox, title="Use graph limits"
	CheckBox check_win, help={"Calculate baseline within limits of graph X-axis"}	
	i+=1.5
	Button CCH_sub,pos={65,v_top+deltaY*i},size={70,20},title="Subtract", proc=archull#CCH_buttons
	Button CCH_sub,help={"Subtract baseline and save a copy"}, fSize=font
	i+=1
	Button CCH_all,pos={65,v_top+deltaY*i},size={70,20},title="All in one", proc=archull#CCH_buttons
	Button CCH_all,help={"Subtract baseline from all traces"}, fSize=font
	
	if (strlen(traceStr)==0)
		return 0
	endif
	return 1	
end

// checks for change in range of x-axis 
static function CCH_graphHook(s)
	STRUCT WMWinHookStruct &s
		
	if ( (UPDATE_ON_ZOOM==0) || (s.eventCode!=8) )
		return 0
	endif
	
	string graphStr=CCH_getGraph()
	string panelStr=graphStr+"#CCH_panel"
	controlInfo /W=$panelStr check_win
	if (V_Value==0)
		return 0
	endif
	
	getaxis /Q bottom
	make /free/N=2 freeLimits={v_min, v_max}
	
	DFREF dfr = root:Packages:ArcHull
	wave /z w_limits=dfr:w_limits
	if (waveexists(w_limits)==0)
		duplicate /o freeLimits dfr:w_limits /Wave=w_limits
	endif
	
	w_limits-=freeLimits
	w_limits=abs(w_limits)
	if(sum(w_limits)>1e-12)
		CCH_doFit()
	endif
	w_limits=freeLimits
	
	return 0
end

// update baseline and ensure package waves are displayed
static function CCH_updatePlot()
	
	DFREF dfr = root:Packages:ArcHull
	
	string graphStr=CCH_getGraph()
	string panelStr=graphStr+"#CCH_panel"
		
	// *****
	// check traceinfo for axis flags, give up if not default
	// *****
		
	removefromgraph /W=$graphStr /Z Arc_Hull_Base, Arc_Hull_Sub, W_YHull
	
	controlInfo /W=$panelStr popTrace
	string traceStr=s_value
	
	// do fit to create baseline
	CCH_doFit()
	
	wave /SDFR=dfr w_base, w_sub	
	wave /Z w_dataX=XWaveRefFromTrace(graphStr, traceStr) 	
		
	if(waveexists(w_dataX))
		wave w_baseX=dfr:w_baseX
		appendtograph /W=$graphStr w_base/TN=Arc_Hull_Base vs w_baseX
		appendtograph /W=$graphStr w_sub/TN=Arc_Hull_Sub vs w_baseX
	else
		AppendToGraph /W=$graphStr w_base/TN=Arc_Hull_Base, w_sub/TN=Arc_Hull_Sub
	endif
	
	// figure out a contrasting colour for baseline and sub waves
	make /free w_RGB={0,10000,65535}
	CCH_chooseColour(w_RGB, graphStr, traceStr)	
	ModifyGraph /W=$graphStr rgb(Arc_Hull_Base)=(w_RGB[0],w_RGB[1],w_RGB[2])
	w_RGB={0,10000,30000}
	ModifyGraph /W=$graphStr rgb(Arc_Hull_Sub)=(w_RGB[0],w_RGB[1],w_RGB[2])
	
	
	// if data are offset in y, then apply offset to baseline
	string infoStr=traceinfo(graphStr,traceStr,0)
	variable trace_offset=GetNumFromModifyStr(infoStr,"offset","{",1)
	ModifyGraph offset(Arc_Hull_Base)={0,trace_offset}
	
	wave W_Xhull=W_Xhull, W_Yhull=W_Yhull
	controlInfo /W=$panelStr check_spline
	if(v_value)
		checkdisplayed /W=$graphStr W_Yhull
		if (v_flag==0)
			AppendToGraph /W=$graphStr W_Yhull vs W_XHull
			ModifyGraph /W=$graphStr mode(W_YHull)=3,marker(W_YHull)=42,rgb(W_YHull)=(0,65535,0)
		endif
		ModifyGraph offset(W_Yhull)={0,trace_offset}
	else
		removefromgraph /W=$graphStr /Z W_YHull
	endif
	setwindow $graphStr hook(ArcHullHook)=archull#CCH_graphHook
	return 1	
end

static function CCH_checkbox(s)
	STRUCT WMCheckboxAction &s
	
	string graphStr=parseFilepath(0, s.win, "#", 0, 0)
	
	if(s.eventCode==-1)			
		removefromgraph /W=$graphStr /Z Arc_Hull_Base, Arc_Hull_Sub, W_YHull
		setwindow $graphStr hook(ArcHullHook)=$""
		return 0
	endif
	
	if(stringmatch(s.ctrlName, "check_CCH"))
		checkbox check_spline win=$s.win, value=0
		checkbox check_ItFit win=$s.win, value=0
		SetVariable setvarDepth win=$s.win,disable=0
		PopupMenu popFitFunc, win=$s.win,disable=1	
		
		controlinfo /W= $s.win setvarDepth
		if (v_value==0)
			// set default depth
			controlinfo /W= $s.win popTrace
			wave /Z w=TraceNameToWaveRef(graphStr, S_Value)
			if (waveexists(w)) // just in case
				variable depth=(wavemax(w)-wavemin(w))*0.15
				SetVariable setvarDepth,win=$s.win,limits={-inf,inf,abs(depth/10)},value=_NUM:depth
			endif
		endif
	endif
		
	if(stringmatch(s.ctrlName, "check_spline"))
		checkbox check_CCH win=$s.win, value=0
		checkbox check_ItFit win=$s.win, value=0
		PopupMenu popFitFunc, win=$s.win,disable=1	
		SetVariable  setvarSmooth win=$s.win, value=_NUM:1
		SetVariable setvarDepth win=$s.win,disable=0		
	endif
	
	if(stringmatch(s.ctrlName, "check_ItFit"))
		checkbox check_CCH win=$s.win, value=0
		checkbox check_spline win=$s.win, value=0
		PopupMenu popFitFunc, win=$s.win,disable=0	
		SetVariable setvarDepth win=$s.win, disable=1
		SetVariable setvarSmooth win=$s.win, value=_NUM:0//, disable=2
	endif
	
	CCH_updatePlot()
	
	return 0	
end

static function CCH_SetVar(s) : SetVariableControl
	STRUCT WMSetVariableAction &s
	
	if(s.eventCode==-1)
		return 0
	endif
	
	if(stringmatch(s.ctrlName,"setvarDepth"))	
		if (s.dval==0)
			
			string graphStr=parseFilepath(0, s.win, "#", 0, 0)
			controlinfo /W= $s.win popTrace
			wave /Z w=TraceNameToWaveRef(graphStr, s_value)
			if (waveexists(w)==0)
				return 0
			endif			
			SetVariable  setvarDepth win=$s.win,limits={-inf,inf,(wavemax(w)-wavemin(w))*.15}
		else
			// reset increment value to 10% of current value
			SetVariable  setvarDepth win=$s.win,limits={-inf,inf,abs(s.dval/10)}
		endif
	endif
	
	CCH_doFit()
	
	return 0
end

static function CCH_popup(s)
	STRUCT WMPopupAction &s
	
	if(s.eventCode==-1)
		return 0
	endif
		
	string graphStr=parseFilepath(0, s.win, "#", 0, 0)
	controlinfo /W= $s.win popTrace
	wave  /Z w=TraceNameToWaveRef(graphStr, s_value)
	if (waveexists(w)==0)
		return 0
	endif
	
	if(stringmatch(s.ctrlName, "popFitFunc"))
		CCH_doFit()
		return 1
	endif
		
	// set default depth each time a trace is selected
	controlinfo /W= $s.win check_CCH	
	if (V_Value)
		variable depth=(wavemax(w)-wavemin(w))*0.15
		SetVariable setvarDepth,win=$s.win,limits={-inf,inf,abs(depth/10)},value=_NUM:depth
	endif
	
	CCH_updatePlot() // does fit and ensures package waves are displayed
			
	return 0
end

static function CCH_buttons(s) : ButtonControl
	STRUCT WMButtonAction &s

	if(s.eventCode!=2)
		return 0
	endif
	
	string graphStr=parseFilepath(0, s.win, "#", 0, 0)	
	string traceStr
		
	strswitch(s.ctrlName)
		case "CCH_sub":	
			CCH_subtract(0)	
			break
		case "CCH_all":
			string msg="subtract baseline from all traces using current settings?\r"
			msg+="existing baselines will be overwritten!"
			doalert 1, msg
			if(v_flag==2)
				return 0
			endif
					
			string listStr=CCH_traces()
			PopupMenu popTrace, win=$s.win, Value=archull#CCH_traces()
			variable i
			for(i=0;i<itemsinlist(listStr); i+=1)				
				PopupMenu popTrace, win=$s.win, mode=i+1							
				CCH_doFit()
				CCH_subtract(1)
			endfor
			
			break
	endswitch

	return 0
end

static function /s CCH_traces()
	
	string graphStr=CCH_getGraph()
	string listStr=TraceNameList(graphStr,";",1+4)
	string removeStr="W_YHull;"
	removeStr+=listmatch(listStr, "Arc_Hull_Base")
	removeStr+=listmatch(listStr, "*_sub")
	removeStr+=listmatch(listStr, "*_BL")
	listStr=RemoveFromList(removeStr, listStr, ";", 0)
	return listStr
end

static function CCH_doFit()
	
	string graphStr=CCH_getGraph()
	string panelStr=graphStr+"#CCH_panel"
	
	controlinfo  /W=$graphStr+"#CCH_panel" popTrace
	string traceStr=s_value

	
	DFREF dfr = root:Packages:ArcHull
	
	killwaves /Z dfr:w_baseX // won't die if it's used on plot....
	
	variable depth, varSm, varOffset, negPeaks, GraphLim
	controlinfo /W=$panelStr setvarDepth
	depth=v_value
	controlinfo /W=$panelStr setvarSmooth
	varSm=v_value
	controlinfo /W=$panelStr setvarOffset
	varOffset=v_value
	controlinfo /W=$panelStr check_neg
	negPeaks=v_value
	controlinfo /W=$panelStr check_win
	GraphLim=v_value
	
	wave /Z w=TraceNameToWaveRef(graphStr, traceStr)
	if (waveexists(w)==0 )
		return 0
	endif
	wave /Z w_x=XWaveRefFromTrace(graphStr, traceStr)
	
	make /free/n=2 w_pRange={0,numpnts(w)-1}, w_Xrange=nan
	variable p1=0, p2=numpnts(w)-1
		
	if(GraphLim)
		string xAxisName= StringByKey("XAXIS", traceinfo(graphStr,traceStr,0))
		getaxis /W=$graphStr /Q $xAxisName
		w_Xrange={v_min,V_max}
		sort w_Xrange, w_Xrange
		if(waveexists(w_x))
			if (w_x[1]<w_x[0])
				sort /R w_Xrange, w_Xrange
			endif
			findlevel /Q /P w_x, w_Xrange[0]
			if (v_flag==0)
				w_pRange[0]=V_LevelX
			endif
			findlevel /Q /P w_x, w_Xrange[1]
			if (v_flag==0)
				w_pRange[1]=V_LevelX
			endif		
		else
			w_pRange[0]=x2pnt(w,w_Xrange[0])
			w_pRange[1]=x2pnt(w,w_Xrange[1])
		endif		
	endif
	sort w_pRange,w_pRange
	p1=max(0, w_pRange[0])
	p2=min(p2, w_pRange[1])
	
	// create waves to display baseline and baseline subtraction
	duplicate /o /r=[p1,p2] w dfr:w_sub /WAVE=w_sub
	duplicate /o /r=[p1,p2] w dfr:w_base /WAVE=w_base
	if(waveexists(w_x))
		duplicate /o /r=[p1,p2] w_x dfr:w_baseX /WAVE=w_x
	endif
	
	note /k w_base
	string cmd=""
	
	if(numpnts(w)<3)
		w_sub=nan
		w_base=nan
		return 0
	endif
	
	controlinfo /W=$panelStr check_ItFit
	if(v_value)	
		duplicate /o w_sub dfr:w_ItFit /WAVE=w_ItFit
		if (varSm)
			smooth varSm, w_ItFit
		endif
		duplicate /free w_sub w_next
		variable i, conv=inf, success=0
		controlinfo /W=$panelStr popFitFunc
		string fName=s_value
		
		for (i=0;i<1e4;i+=1)		
			// fit w_ItFit with fName, output to w_base
			success=CCH_FitWrap(w_ItFit, w_x, w_base, fName)			
			if (success==0)
				break
			endif
			w_next =  (w_ItFit>w_base)  ? w_base : w_ItFit
			if (negPeaks)
				w_next =  (w_ItFit<w_base)  ? w_base : w_ItFit
			endif
			if(success && i<10) // make at least 10 iterations
				w_ItFit=w_next
				continue
			endif
			// convergence test - get a count of the number of points changed in this iteration
			w_ItFit=1/(w_next-w_ItFit)
			wavestats /Q/M=1 w_ItFit
			if (i>10 && v_npnts<=conv) // converged
				break
			else
				conv=v_npnts
				w_ItFit=w_next
			endif	
		endfor
		w_sub=w[p1+p]-w_base+varOffset	
		sprintf cmd, "Iterative Fit with %g iterations of %s" i+1, fName
		if (success)
			cmd+=", no fit errors"
		else
			cmd+=", fit error"
		endif
		if (varSm)
			cmd+=", smoothing="+num2str(varSm)
		endif
		cmd+=";"+note(w_base)
	else // do arc hull calculation			
		if(varSm>0)
			smooth varSm, w_sub
		endif
		// make a copy of the (possibly smoothed) data wave
		duplicate /free w_sub w_ref 
		// calculate arc hull based on (possibly smoothed) data
		ArcHull(w_sub, depth, w_x=w_x, top=negPeaks)
		w_base=w_ref-w_sub
		w_sub=w[p1+p]-w_base+varOffset
		sprintf cmd, "Arc hull baseline with arc depth = %g, smoothing = %g", depth, varSm
	endif
	
	controlinfo /W=$panelStr check_spline
	if (v_value)
		wave w_XHull=w_XHull, w_YHull=w_YHull
		
		if (depth!=0)
			w_YHull=getY(w_XHull, w_ref, w_x)
		endif
	
		if (numpnts(w_XHull)>3)
			if (waveexists(w_x))
				Interpolate2/T=2/E=2/Y=w_base/X=w_x/I=3 W_XHull,W_YHull
			else
				Interpolate2/T=2/E=2/Y=w_base/I=3 W_XHull,W_YHull
			endif
			w_sub=w[p1+p]-w_base+varOffset	
		endif
		sprintf cmd, "Cubic spline baseline with nodes from convex hull and smoothing = %g", varSm
	endif
	
	if(GraphLim)
		if(waveexists(w_x))
			v_min=w_x[0]
			v_max=w_x[numpnts(w_x)-1]
		endif
		sprintf cmd, "%s, xrange from %g to %g", cmd, v_min, v_max
	endif
	note /K w_base cmd
end

static function getY(v_x, w, w_x)
	variable v_x
	wave /Z w, w_x
	
	if(waveexists(w_x))
		findlevel /Q /P w_x, v_x
		if (v_flag==0)
			return w[V_LevelX]
		endif
	else
		return w(v_x)
	endif
	return nan
end

// subtracts arc hull baseline from w (w is overwritten)
function ArcHull(w, depth [, w_x, top])
	wave /Z w, w_x
	variable depth, top	
	
	variable useXwave=0
	if (paramisdefault(w_x)==0 && waveexists(w_x))
		useXwave=1
	endif	
	
	top = paramisdefault(top) ? 0 : top
	// top=1 for top hull
	
	// add concave function
	variable radius	
	if(useXwave)
		radius=(wavemax(w_x)-wavemin(w_x))/2
		w+=depth*(w_x[p]-w_x[0])^2/radius^2
	else
		variable x1 = leftx(w), x2=pnt2x(w,numpnts(w)-1)
		radius=abs(x2-x1)/2
		w+=depth*(x-x1)^2/radius^2
	endif
	
	duplicate /free w w_base
	
	if (useXwave==0)
		duplicate /free w w_x
		// redimension to make sure we don't get hull vertices 
		// outside the range of w_x owing to difference in precision
		// of output from ConvexHull
		redimension /D w_x
		w_x=x
	endif
	
	Convexhull  w_x,w
	wave w_XHull=w_XHull, w_YHull=w_YHull

	if (top==0)
		reverse /P w_XHull, w_YHull
	endif
	// negative depth will subtract top part of convex hull
	// 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

// a wrapper for built-in fit functions
static function CCH_FitWrap(w, w_x, w_out, fName)
	wave /Z w, w_x, w_out
	string fName
	
	variable V_FitError=0, V_fitOptions=4
	try
		strswitch (fName)
			case "line":
				CurveFit /Q/N line, w /D=w_out /X=w_x /NWOK; AbortOnRTE
				break
			case "poly 3":
				CurveFit /Q/N poly 3, w /D=w_out /X=w_x /NWOK; AbortOnRTE
				break
			case "poly 4":
				CurveFit /Q/N poly 4, w /D=w_out /X=w_x /NWOK; AbortOnRTE
				break
			case "gauss":
				CurveFit /Q/N gauss, w /D=w_out /X=w_x /NWOK; AbortOnRTE
				break
			case "lor":
				CurveFit /Q/N lor, w /D=w_out /X=w_x /NWOK; AbortOnRTE
				break
			case "exp":
				CurveFit /Q/N exp, w /D=w_out /X=w_x /NWOK; AbortOnRTE
				break
			case "dblexp":
				CurveFit /Q/N dblexp, w /D=w_out /X=w_x /NWOK; AbortOnRTE
				break
			case "sin":
				CurveFit /Q/N sin, w /D=w_out /X=w_x /NWOK; AbortOnRTE
				break
			case "hillequation":
				CurveFit /Q/N hillequation, w /D=w_out /X=w_x /NWOK; AbortOnRTE
				break
			case "sigmoid":
				CurveFit /Q/N sigmoid, w /D=w_out /X=w_x /NWOK; AbortOnRTE
				break
			case "power":
				CurveFit /Q/N power, w /D=w_out /X=w_x /NWOK; AbortOnRTE
				break
			case "lognormal":
				CurveFit /Q/N lognormal, w /D=w_out /X=w_x /NWOK; AbortOnRTE
				break				
		endswitch
	catch
		if (V_AbortCode == -4)
			Print "Error during curve fit:"
			Variable CFerror = GetRTError(1)	// 1 to clear the error
			Print GetErrMessage(CFerror)
		endif
		return 0
	endtry
	
	if (V_FitError==0) // record fit coefficients in baseline wavenote
		wave w_coef=w_coef
		string coefStr
		variable i
		sprintf coefstr, "w_coef={"
		for(i=0;i<numpnts(w_coef);i+=1)
			sprintf coefstr, "%s%g,", coefstr, w_coef[i]
		endfor
		coefstr=removeEnding(coefstr, ",")+"};"
		note /K w_out, coefstr
	endif
	
	return (V_FitError	==0)
end

static function CCH_subtract(overwrite)
	variable overwrite
	
	variable appendBL=0
	
	string graphStr=CCH_getGraph()
	string panelStr=graphStr+"#CCH_Panel"
	
	controlinfo  /W=$panelStr popTrace
	string traceStr=s_value
	
	wave /Z w=TraceNameToWaveRef(graphStr, traceStr)
	
	if (waveexists(w)==0)
		return 0
	endif
	
	wave /Z w_dataX=XWaveRefFromTrace(graphStr, traceStr)
	
	DFREF CCHdfr = root:Packages:ArcHull	
	wave /Z/SDFR=CCHdfr w_sub, w_base, w_baseX
	
	DFREF CurrentDfr = $getdatafolder(1)
		
	// save a copy of the baseline
	string strNewName=CleanupName(nameofwave(w)+"_BL",0)
	if (overwrite==0 && exists(strNewName))
		doalert 1, strNewName + " exists. Overwrite?"
		if(V_flag==2)
			return 0
		endif
	endif
	removefromgraph /W=$graphStr /Z $strNewName
	duplicate /o w_base CurrentDfr:$strNewName
	wave /Z newbase= $strNewName
	printf "Arc Hull created wave %s", strNewName
	
	// save a copy of baseline-subtracted wave
	strNewName=CleanupName(nameofwave(w)+"_sub",0)
	if (overwrite==0 && exists(strNewName))
		doalert 1, strNewName + " exists. Overwrite?"
		if(V_flag==2)
			return 0
		endif
	endif
	removefromgraph /W=$graphStr /Z $strNewName
	duplicate /o w_sub CurrentDfr:$strNewName
	wave subtracted=$strNewName
	
	note subtracted, note(w) // append note from data wave to output wave note
	note subtracted, note(w_base)
	printf ", %s", strNewName
	
	if (waveexists(w_dataX))				
		if (numpnts(w_baseX)<numpnts(w_dataX)) // using a subrange			
			strNewName = CleanupName(nameofwave(w_dataX),0)
			strNewName = UniqueName(strNewName, 1, 1)
			duplicate /o w_baseX CurrentDfr:$strNewName
			wave w_x = $strNewName
			printf ", %s", strNewName
			note w_x "X wave for " + nameofwave(subtracted)
		else
			wave w_x = w_dataX
		endif
	endif
	
	printf "\r"
	
	// add baseline-subtracted output wave to plot
	checkdisplayed /W=$graphStr subtracted
	if(v_flag==0)
		if(waveexists(w_x))
			appendtograph /W=$graphStr /C=(0,0,0) subtracted vs w_x
		else
			appendtograph /W=$graphStr /C=(0,0,0) subtracted
		endif
	endif
	
	// possibly plot baseline
	if(appendBL)
		checkdisplayed /W=$graphStr newbase
		if(v_flag==0)
			if(waveexists(w_x))
				appendtograph /W=$graphStr /C=(0,0,0) newbase vs w_x
			else
				appendtograph /W=$graphStr /C=(0,0,0) newbase
			endif
		endif
	endif	
end

static function CCH_clear()	
	
	string graphStr=CCH_getGraph()
	
	if (strlen(graphStr))
		removefromgraph /W=$graphStr /Z Arc_Hull_Base, Arc_Hull_Sub, W_YHull
		setwindow $graphStr hook(ArcHullHook)=$""
		KillWindow $graphStr+"#CCH_panel"
	endif
end

// set colour wave w_color for visibility against trace.
static function CCH_chooseColour(w_RGB, graphStr, traceStr)
	wave w_RGB
	string graphStr, traceStr
	
	// figure out a contrasting colour for baseline and sub waves
	variable c0,c1,c2
	string infoStr=traceinfo(graphStr, traceStr, 0)
	variable startIndex=strsearch(infoStr, ";rgb(x)=", 0)
	sscanf infoStr[startIndex,strlen(infoStr)-1], ";rgb(x)=(%d,%d,%d*", c0,c1,c2
	make /free w_TraceRGB={c0,c1,c2}, w_index={1,2,3}
	MakeIndex /R w_TraceRGB, w_index
	IndexSort  w_index, w_RGB 
	return 1	
end

// for Igor 6 compatibility
static function /T CCH_getGraph()

	#if IgorVersion() >= 7	
		GetWindow /Z CCH_panel  activeSW
		return parseFilepath(0, s_value, "#", 0, 0)
	#endif	
	
	string listStr=WinList("*", ";","WIN:1")
	string panelStr, graphStr
	variable i
	do
		graphStr=stringfromlist(i, listStr)
		panelStr=graphStr+"#CCH_panel"
		if (wintype(panelStr))
			return graphStr
		endif
		i+=1
	while (i<itemsinlist(listStr))
	return ""
end


// Shen et al (2018) Automatic baseline correction method for the 
// open-path Fourier transform infrared spectra by using simple iterative
// averaging. Optics Express 26, A609-A614. 
// https:doi.org/10.1364/OE.26.00A609
// This is slow, so I haven't incorporated it into the archull panel.
static function CCH_ItAv(w)
	wave w
	duplicate/o w w_bl
	duplicate /free w, w_previous, w_diff
	variable Pmax=numpnts(w)-1, i=0
	variable Pstart, Pend, Sabs, previousSabs, deltaSabs
	
	variable threshold=0.005, maxits=1e3
	variable timeout=datetime+3
	do
		previousSabs=Sabs
		w_previous=w_bl
		Pstart=1
		Pend=Pmax-1	
		do
			w_bl[Pstart,Pend]=min(w_bl, (w_bl[p-1]+w_bl[p+1])/2)
			Pstart+=1
			Pend-=1
		while(Pstart<Pmax/2)
		
		w_diff=abs(w_bl-w_previous)
		Sabs=sum(w_diff)
		deltaSabs=abs(Sabs-previousSabs)		
		i+=1
		if (datetime>timeout)
			print "baseline calculation timed out"
			break
		endif
	while((deltaSabs/Sabs)>threshold && i<maxits)
	print i, deltaSabs/Sabs
end

static function makeSpectrum()
	
	string strName=uniquename("foo", 1, 0)
	make /n=1000 $(uniquename("foo", 1, 0)) /WAVE=foo
	variable a=enoise(5), b=enoise(1e-3), c=enoise(1e-5), d=enoise(400)
	foo=gnoise(0.1)+a+b*x+c*(x-d)^2
	variable i, height, position, width
	for(i=0;i<5;i+=1)
		height=150+enoise(50)
		position=500+enoise(400)
		width=20+enoise(10)
		foo+=height*gauss(x, position, width/sqrt(2))
	endfor
	display /K=1 foo
end
	
// execute archull#demo() to create a fake spectrum for demo
static function demo()
	makespectrum()
	CCH_init()
end