#pragma rtGlobals=3
#pragma version=4.03
#pragma IgorVersion=6
#pragma ModuleName=BaselineSpline

#include <Readback ModifyStr>

// written by tony withers, tony.withers@uwo.ca

// please tell me when you find bugs!

// How to use Baseline Spline Fit

// A 1d wave must be plotted in a graph window before selecting Spline 
// Baseline from the macros menu. The wave must be plotted on the default
// bottom and left axes. A control panel will be added to the left of the
// graph window. When 'edit mode' is selected you can drag the nodes. The
// spline should update as the nodes are repositioned. When in edit mode 
// you can add a new node at any position on the graph by positioning the
// mouse cursor and clicking while holding down the control key. Nodes 
// may be deleted by clicking them with the option key held down. 
// Clicking anywhere in the graph window with the shift key held down 
// will toggle between node editing mode and 'normal' mode where you can 
// interact with the graph as usual. execute 'DisplayHelpTopic 
// Interpolate2' for details of the methods of spline interpolation.

// Panel controls

// Data wave: select a trace for baseline correction. When switching 
// between traces the nodes will retain their X-positions but will be 
// placed on the selected data wave at those positions (unless they lie 
// outside of the range of the data wave).

// Akima: fit an Akima style spline through the nodes. Akima splines can 
// have sharper bends than cubic splines 

// Cubic: fit a cubic spline through the nodes.

// Ends: At the ends of the spline where the cubic function is not fully 
// constrained select a method to add a constraint. 1st deriv.: match 1st
// derivative; natural: match 2nd derivative. Affects only cubic spline 
// type.

// Linear: connect the nodes with linear segments.

// Smoothing: fit a smoothing spline, which does not have to pass exactly
// through the nodes. The degree of smoothing is adjustable. 

// Edit mode: toggle between node editing and normal modes

// Reset: set node positions by spacing a set number of nodes at regular 
// intervals along the section of the data wave displayed in the graph 
// window. No nodes will be set beyond the range of the graph axes.  

// Load: reload the node positions from a previous fit.

// Subtract: creates a copy of the baseline and a baseline-subtracted 
// output spectrum in the current data folder and appends them to the 
// plot; saves the node positions in a 2d wave in the current data 
// folder.


// Version history

// 4.03 11/11/17 Should be backward compatible for Igor 6.
// 4.02 Catch interpolate errors. 
// 4.01 Panel size decreased for compatibility with smaller or lower 
// resolution screens. Some bugs fixed.
// 4.00 Changed to panel interface. Major rewrite.
// 3.20 Removed option to prevent nodes from being inserted. Now the /NI 
// no insert flag is used by default, and we handle node insertion by 
// using a hook to collect ctrl-clicks.
// 3.10 Some changes to the global strings and variables in the package 
// folder. Tried and failed to add improved guesses for node positions.
// 3.01 11/4/16 Bug fix - waves with deltaX<0 produced nodes that were 
// reverse sorted in X, but Akima calculation cannot deal with reverse 
// sorted nodes. 
// 3.00 4/11/16 Added option for Akima style spline.
// 2.60 3/30/16 Added option to supress node addition. This is a toggle 
// that switches on the /NI option in GraphWaveEdit. Updated a lot of 
// other code. Changed some of the defaults (now saves nodes by default),
// so behaviour is not the same as previous version.
// 2.50 Log mode to save node positions to file. Can load node positions 
// from a file.
// 2.40 Added textbox reminder of how to toggle graph mode, renamed 
// procedures for consistency.
// 2.30 5/16/11 Bug fix - code was choking on non-existent wave in first 
// call to SBL_Clear(). Added panel to change fitting options.
// 2.10 Call SBL_Clear() from SBL_Subtract().
// 2.0 8/6/09 Incorporated suggestions from JJ Weimer that improve 
// functionality and clarity of code and compatibility with 6.10 ... and 
// then wrecked his nice programming with my clumsy code.
// 1.43b Allow refit function to switch interpolate flag when number of 
// nodes is small to allow linear fit.
// 1.42 6/8/08 Set y axis to manual range to avoid crazy rescaling if 
// spline shoots out of range. Show subtracted spectrum whilst fitting.
// 1.41 10/9/07 For XY data we now interpolate directly into XY baseline.
// Improved guesses for initial node positions for unevenly spaced XY 
// data.
// 1.40 10/9/07 Fixed bug that occurred with baseline subtraction of XY 
// data.
// 1.30 7/24/07 Handles XY data. Baseline is waveform, but data wave need
// not be (though it should be reasonably closely spaced).
// 1.20 7/3/07 Uses new graphwaveedit flag. Now requires version 6.0.1.

Menu "Analysis"
	Submenu "Packages"
		"Spline Baseline",/Q, SBL_init()
	end
end

Menu "Macros"	
	"Spline Baseline",/Q, SBL_init()
end

// initialize spline baseline package
function SBL_init()
	
	if (strlen(WinList("*",";","WIN:1"))==0)  // no graphs 
		doalert 0, "Spline baseline requires a trace plotted in a graph window."
		return 0
	endif
	
	// make sure the top graph is visible
	dowindow /F $WinName(0,1)
	
	// make sure graph has default axes
	getaxis /Q left
	if (v_flag==0)
		getaxis /Q bottom
	endif
	if (v_flag)
		doalert 0, "Trace must be plotted on bottom/left axes."
		return 0
	endif
	
	// clear any package detritus from the last used spline graph
	SBL_clear(1)
	
	// initialize (or reinitialize) the package data folder
	SBL_resetDFR()
	DFREF dfr = SBL_getDFREF()
	wave /SDFR=dfr W_spline_dependency
	
	SBL_makePanel()
		
	string s_graph=SBL_getGraph()
	controlinfo /W=$s_graph+"#SBL_panel" popTrace
	string s_trace=s_value
	
	variable success

	// place nodes within window
	success=SBL_resetNodes()
	if (!success)
		doalert 0, "Spline baseline requires a trace plotted in a graph window."
		SBL_clear(1)
		return 0
	endif
	
	// set the sub and base waves for this trace
	success=SBL_resetGraph(s_trace)
	if (!success)
		doalert 0, "Spline baseline requires a trace plotted in a graph window."
		SBL_Clear(1)
		return 0
	endif
	
	ModifyGraph zero(left)=1
	
	// set a dependency to trigger interpolation AND update graph when we adjust a node position
	W_spline_dependency:=SBL_update(root:Packages:SplineFit:W_nodesY)	
	// won't work if we change package folder
	
	GraphWaveEdit /W=$s_graph /NI/T=1 /M $"W_nodesY"
	setWindow $s_graph hook(SBL_hook)=SBL_graphHook
		
	return 1	
end

function SBL_graphHook(s)
	STRUCT WMWinHookStruct &s	
	
	if (s.eventcode!=5) // mouseup
		return 0
	endif
	
	variable v_key=GetKeyState(0) 
	if (v_key==0)
		return 0
	endif
	
	if(v_key&2^2) // shift click
		// toggle mode
		SBL_toggleEdit()
		return 0
	endif
	
	controlInfo /W=$s.winName+"#SBL_panel" check_edit
	if(v_value==0)
		return 0
	endif

	variable insert=0
	insert = (stringmatch(igorinfo(2),"Macintosh")) ? (v_key&2^4) : (v_key&2^0) // ctrl
	if (insert)
		DFREF dfr = SBL_getDFREF()
		wave /SDFR=dfr w_nodesY, w_nodesX
		variable v_X, v_Y
		v_X=AxisValFromPixel(s.winName, "bottom", s.mouseLoc.h)
		v_Y=AxisValFromPixel(s.winName, "left", s.mouseLoc.v)
		
		// deal with offset
		controlInfo /W=$s.winName+"#SBL_panel" popTrace
		string s_info=traceinfo(s.winName,s_value,0)
		v_Y-=GetNumFromModifyStr(s_info,"offset","{",1)
		
		w_nodesX[numpnts(w_nodesX)]={v_X}
		w_nodesY[numpnts(w_nodesY)]={v_Y}
		sort w_nodesX, w_nodesX, w_nodesY
		SBL_update({nan})
	endif

	return 0
end


function SBL_toggleEdit()
	
	string s_graph=SBL_getGraph()
	controlinfo /W=$s_graph+"#SBL_panel" check_edit
	if (v_value)
		GraphNormal /W=$s_graph
		ModifyGraph /W=$s_graph mode(W_nodesY)=2
		checkbox check_edit win=$s_graph+"#SBL_panel", value=0
	else
		// manually scale y axis
		getaxis /W=$s_graph /Q left
		setAxis /W=$s_graph left, V_min, V_max	
		
		ModifyGraph /W=$s_graph mode(W_nodesY)=3
		GraphWaveEdit /W=$s_graph /NI/T=1 /M $"W_nodesY"	
		checkbox check_edit win=$s_graph+"#SBL_panel", value=1	
	endif
end

function SBL_resetGraph(s_trace)
	string s_trace
	
	DFREF dfr = SBL_getDFREF()
	
	string s_graph=SBL_getGraph()
	wave /Z w_data=TraceNameToWaveRef(s_graph, s_trace)
	if (waveexists(w_data)==0)
		return 0
	endif
	wave /Z w_x=XWaveRefFromTrace(s_graph, s_trace)
	
	if (waveexists(w_x))
		// check for monotonic X wave
		make /free w_diff
		Differentiate w_x /D=w_diff
		Wavestats /Q/M=0 W_diff
		if (V_max>0 && V_min<0)
			doalert 0, "X wave is not monatonic! Quitting :("			
			SBL_clear(1)
			return 0
		endif	
	endif
	
//	string xAxisName= StringByKey("XAXIS", traceinfo(s_graph,s_trace,0))
//	string yAxisName= StringByKey("YAXIS", traceinfo(s_graph,s_trace,0))
		
	// manually scale y axis
	getaxis /Q left
	setAxis left, V_min, V_max
		
	duplicate /O w_data dfr:W_base, dfr:W_sub
	wave /SDFR=dfr w_nodesX, w_nodesY, w_spline_dependency, w_base, w_sub
	
	
	checkdisplayed /W=$s_graph W_spline_dependency
	if (v_flag==0)
		appendtograph /W=$s_graph W_spline_dependency
	endif
	
	// make sure base and sub waves are displayed properly
	removefromgraph /Z /W=$s_graph w_base, w_sub
	if (waveexists(w_x))	
		appendtograph /W=$s_graph W_base vs w_x
		appendtograph /W=$s_graph W_sub vs w_x
	else
		appendtograph /W=$s_graph W_base
		appendtograph /W=$s_graph W_sub
	endif

	ModifyGraph /W=$s_graph hideTrace(W_spline_dependency)=1
	
	// if data are offset, then apply offsets to baseline
	string s_info=traceinfo(s_graph,s_trace,0)
	variable trace_offset=GetNumFromModifyStr(s_info,"offset","{",1)
	
	make /n=3 /free w_c={0,10000,65535}
	SBL_chooseColour(w_c, s_graph, s_trace)
	
	ModifyGraph /W=$s_graph offset(W_base)={0,trace_offset}
	ModifyGraph /W=$s_graph rgb(W_base)=(w_c[0],w_c[1],w_c[2])
	modifygraph /W=$s_graph live(W_base)=1	
	ModifyGraph /Z /W=$s_graph offset(W_nodesY)={0,trace_offset}
	ModifyGraph /W=$s_graph rgb(W_sub)=(0,0,0)
	ReorderTraces /W=$s_graph $"w_nodesY",{$"w_base"}

	return 1
end

function SBL_resetNodes()
	
	string s_graph=SBL_getGraph()
	DFREF dfr = SBL_getDFREF()
	wave /SDFR=dfr W_nodesX, W_nodesY

	// define some local variables	
	variable i, p_low, p_high, delP, delX, numNodes=5
	string s_info
	variable trace_offset=0
	
	controlinfo /W=$s_graph+"#SBL_panel" popTrace
	string s_trace=s_value
	
	wave /Z w_data=TraceNameToWaveRef(s_graph, s_trace)
	wave /Z w_x=XWaveRefFromTrace(s_graph, s_trace)
	if (waveexists(w_data)==0)
		return 0
	endif	
	redimension/N=(numNodes) W_nodesX, W_nodesY
	
	// determine where to place nodes on the bottom axis
	
	// keep nodes within window 
	GetAxis/Q bottom
	make /n=2 /free xrange={v_min, v_max}
	sort xrange, xrange
	
	GetAxis/Q left
	make /n=2 /free yrange={v_min, v_max}
	sort yrange, yrange
	
	// keep nodes within extent of w_data	
	if (waveexists(w_x))
		xrange[0]=max(xrange[0],wavemin(w_x))
		xrange[1]=min(xrange[1],wavemax(w_x))
	else
		xrange[0]=max(xrange[0], min(pnt2x(w_data, 0), pnt2x(w_data, numpnts(w_data)-1)) )
		xrange[1]=min(xrange[1], max(pnt2x(w_data, 0), pnt2x(w_data, numpnts(w_data)-1)) )
	endif

	delX=(xrange[1]-xrange[0])/(numNodes-1)
	// initialize node positions		
	W_nodesX=xrange[0]+p*delX
	
	for (i=0;i<numNodes;i+=1)
		if (waveexists(w_x))
			// necessary for unevenly spaced data
			findlevel /Q/P w_x, W_nodesX[i]-delX/24
			if (v_flag)
				findlevel /Q/P w_x, W_nodesX[i] // this HAS to be within range
			endif
			p_low = V_LevelX
			
			findlevel /Q/P w_x, W_nodesX[i]+delX/24
			if (v_flag)
				findlevel /Q/P w_x, W_nodesX[i] // this HAS to be within range
			endif
			p_high = V_LevelX
			wavestats /Q/M=1/R=[p_low , p_high] /Z  w_data
		else
			wavestats /Q/M=1/R=(W_nodesX[i]-delX/24 , W_nodesX[i]+delX/24) /Z  w_data
		endif
		W_nodesY[i]=V_avg  // this handles a few NANs	
	endfor
	
	// akima calculation won't work if nodes are reverse sorted
	sort w_nodesX, w_nodesX, w_nodesY
	
	// if data are offset, then apply offsets to nodes
	s_info=traceinfo(s_graph,s_trace,0)
	trace_offset=GetNumFromModifyStr(s_info,"offset","{",1)
	
	// keep nodes within window
	W_nodesY=max(yrange[0]-trace_offset, w_nodesY)
	W_nodesY=min(yrange[1]-trace_offset, w_nodesY)
	
	checkdisplayed /W=$s_graph W_nodesY
	if (v_flag==0)
		appendtograph /W=$s_graph W_nodesY vs W_nodesX
	endif
	ModifyGraph /W=$s_graph mode(W_nodesY)=3,marker(W_nodesY)=11,rgb(W_nodesY)=(0,39168,0)
	ModifyGraph /W=$s_graph offset(W_nodesY)={0,trace_offset}, mode(W_nodesY)=3
		
	// manually scale y axis
	getaxis /W=$s_graph /Q left
	setAxis /W=$s_graph left, V_min, V_max
	
	GraphWaveEdit /W=$s_graph /NI/T=1 /M $"W_nodesY"
	Checkbox check_edit, win=$s_graph+"#SBL_panel", value=1
		
	return 1
end

// place nodes on w_data. doesn't plot waves.
// doesn't reset x positions unless they're out of range.
function SBL_setNodes()
	
	string s_graph=SBL_getGraph()
	controlInfo /W=$s_graph+"#SBL_panel" popTrace
	string s_trace=s_value	
	wave /Z w_data=TraceNameToWaveRef(s_graph, s_trace)
	wave /Z w_x=XWaveRefFromTrace(s_graph, s_trace)
	
	DFREF dfr = SBL_getDFREF()
	wave /SDFR=dfr W_nodesX, W_nodesY
	variable numNodes=numpnts(W_nodesY)
	variable i, v_pnt, v_x
	
	for (i=0;i<numNodes;i+=1)
		if (waveexists(w_x))
			// necessary for unevenly spaced data
			v_x=max(wavemin(w_x), W_nodesX[i])
			v_x=min(wavemax(w_x), v_x)
			
			findlevel /Q/P w_x, v_x
			if (V_flag)
				v_pnt=numpnts(w_x)/2 // shouldn't arrive here
			else
				v_pnt=V_LevelX
			endif
			W_nodesY[i]=w_data[v_pnt]
		else					
			// keep x values within range of w_data
			W_nodesX=max(min(pnt2x(w_data,0), pnt2x(w_data,numpnts(w_data)-1)),W_nodesX)
			W_nodesX=min(max(pnt2x(w_data,0), pnt2x(w_data,numpnts(w_data)-1)),W_nodesX)
			W_nodesY[i]=w_data(W_nodesX[i]) 
		endif
	endfor

#if IgorVersion() >= 7	
	// remove duplicate nodes, 
	// in case some were outside of the range of w_data
	duplicate /free w_nodesX w_duplicates
	findduplicates /SN=(nan) /SNDS=w_duplicates  w_nodesX
	i=0
	do
		if(numtype(w_duplicates[i])!=0)
			deletepoints i, 1, w_nodesX, w_nodesY, w_duplicates
		else
			i+=1
		endif
	while(i<numpnts(w_nodesX))
#else
// you're out of luck
#endif
	
	if(numpnts(W_nodesX)<2)
		SBL_resetNodes()
		doalert 0, "Nodes couldn't be placed on "+nameofwave(w_data)+" at the same x positions. Nodes have been reset."
	endif
	
end

function SBL_makePanel()
		
	DFREF dfr = SBL_getDFREF()
	
	// make sure the top graph is visible
	string s_graph=WinName(0,1)
	dowindow /F $s_graph
	
	// make sure there's space to the left of the graph for the panel
	getwindow /Z $s_graph wsize
	if (V_left<200)
#if IgorVersion() >= 7	
		movewindow /W=$s_graph 200, V_top, -1, -1
#else
		movewindow /W=$s_graph 200, V_top, 200+(V_right-V_left), v_bottom
#endif
	endif

	// make panel
	NewPanel /K=1/N=SBL_panel/W=(200,0,0,355)/HOST=$s_graph/EXT=1 as "Spline Controls"	
	ModifyPanel /W=$s_graph#SBL_panel, noEdit=1		
	
	string s_trace=stringfromlist(0, SBL_traces())
	wave /Z w=TraceNameToWaveRef(s_graph, s_trace)
	
	variable i=0, deltaY=25, font=12, groupw=180
	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=SBL_traces(), title="",pos={v_left,v_top+deltaY*i},size={130,20}
	PopupMenu popTrace, help={"select data wave" }, proc=SBL_popup, fSize=font
	i+=1.5	
	
	GroupBox group1 pos={v_left-20,v_top+deltaY*i},size={groupw,deltaY*5.3},title="Spline parameters"
	GroupBox group1 fSize=font
	i+=1	
	Checkbox check_akima, pos={v_left,v_top+deltaY*i}, fSize=font, value=0, proc=SBL_checkboxes, title="Akima"
	Checkbox check_akima, mode=1,help={"use Akima spline"}	
	i+=0.8
	Checkbox check_cubic, pos={v_left,v_top+deltaY*i}, fSize=font, value=1, proc=SBL_checkboxes, title="Cubic"
	Checkbox check_cubic, mode=1,help={"use cubic spline"}
	i+=0.8
	PopupMenu popE,pos={v_left+35,v_top+deltaY*i},size={46,21}, fSize=font, title="Ends:", proc=SBL_popup
	PopupMenu popE, value="1st deriv.;natural", mode=2
	i+=0.8
	Checkbox check_linear, pos={v_left,v_top+deltaY*i}, fSize=font, value=0, proc=SBL_checkboxes, title="Linear"
	Checkbox check_linear, mode=1,help={"lines between nodes"}
	i+=0.8
	Checkbox check_smoothing, pos={v_left,v_top+deltaY*i}, fSize=font, value=0, proc=SBL_checkboxes, title="Smoothing"
	Checkbox check_smoothing, mode=1,help={"use smoothing spline"}
	SetVariable setvarF, pos={v_left+100,v_top+deltaY*i},size={50,16},title="",bodyWidth=50
	SetVariable setvarF, value=_NUM:1, limits={0,inf,1}, fSize=font, proc=SBL_setvar
	i+=1.5
	
	GroupBox group2,pos={v_left-20,v_top+deltaY*i},size={groupw,deltaY*4.5},title="Node control"
	GroupBox group2,fSize=font
	i+=1	
	Checkbox check_edit, pos={v_left+30,v_top+deltaY*i}, fSize=font, value=1, proc=SBL_Checkboxes, title="Edit mode "
	Checkbox check_edit, mode=0,help={"Toggle to edit nodes"}, side=1
	i+=0.8
	SetDrawLayer ProgBack
	SetDrawEnv textyjust=2
	DrawText v_left-12,v_top+deltaY*i,"Shift-click to toggle mode, ctrl-\r/option-click to add/zap nodes"
	i+=1.5
	Button SBL_nodes,pos={v_left,v_top+deltaY*i},size={50,20},title="Reset", proc=SBL_Buttons
	Button SBL_nodes,help={"Distribute nodes over x-range of graph"}, fSize=font
	
	Button SBL_load,pos={v_left+80,v_top+deltaY*i},size={50,20},title="Load...", proc=SBL_Buttons
	Button SBL_load,help={"Load nodes used in a previous fit"} , fSize=font
	
	i+=1.5
	Button SBL_sub,pos={65,v_top+deltaY*i},size={70,20},title="Subtract", proc=SBL_Buttons
	Button SBL_sub,help={"Subtract baseline and save a copy"}, fSize=font
	
	return 1
end

function SBL_setvar(s)
	STRUCT WMSetVariableAction &s
	
	if(s.eventCode==-1)
		return 0
	endif
	SBL_Update({nan})	
end

function SBL_popup(s)
	STRUCT WMPopupAction &s
	
	if(s.eventCode==-1)
		return 0
	endif
	
	if (stringmatch(s.ctrlName, "popTrace"))
		if ( SBL_ResetGraph(s.popStr) )
			SBL_setNodes()
			SBL_Update({nan})
		endif
		controlInfo /W=$s.win+"#SBL_panel" check_edit
		// not sure why this is needed, there's no reason why it should get out of sync
		if (v_value)
			string s_graph=parseFilepath(0, s.win, "#", 0, 0)
			GraphWaveEdit /W=$s_graph /NI/T=1 /M $"W_nodesY"
		endif
	else
		SBL_Update({nan})	
	endif	
end

function SBL_checkboxes(s)
	STRUCT WMCheckboxAction &s
	
	if(s.eventCode==-1)			
		SBL_clear(0)
		return 0
	endif

	strswitch(s.ctrlName)
		case "check_cubic":
			Checkbox check_akima, win=$s.win, value=0
			Checkbox check_smoothing, win=$s.win, value=0
			Checkbox check_linear, win=$s.win, value=0
			SBL_Update({nan})
			return 1	
		case "check_akima":
			Checkbox check_cubic, win=$s.win, value=0
			Checkbox check_smoothing, win=$s.win, value=0
			Checkbox check_linear, win=$s.win, value=0
			SBL_Update({nan})
			return 1
		case "check_smoothing":
			Checkbox check_cubic, win=$s.win, value=0
			Checkbox check_akima, win=$s.win, value=0
			Checkbox check_linear, win=$s.win, value=0
			SBL_Update({nan})
			return 1
		case "check_linear":
			Checkbox check_cubic, win=$s.win, value=0
			Checkbox check_akima, win=$s.win, value=0
			Checkbox check_smoothing, win=$s.win, value=0
			SBL_Update({nan})
			return 1
		case "check_edit":
			string s_graph=parseFilepath(0, s.win, "#", 0, 0)
			if (s.checked)
				// manually scale y axis
				getaxis /W=$s_graph /Q left
				setAxis /W=$s_graph left, V_min, V_max
				
				ModifyGraph /W=$s_graph mode(W_nodesY)=3
				GraphWaveEdit /W=$s_graph /NI/T=1 /M $"W_nodesY"
			else
				GraphNormal /W=$s_graph
				ModifyGraph /W=$s_graph mode(W_nodesY)=2		
			endif
	endswitch
	return 0
end

function SBL_buttons(s)
	STRUCT WMButtonAction &s
	
	if(s.eventCode!=2)			
		return 0
	endif
	
	if (stringmatch(s.ctrlname, "SBL_sub"))
		SBL_subtract()
	endif
	
	if (stringmatch(s.ctrlname, "SBL_nodes"))
		SBL_resetNodes()
	endif
	
	if (stringmatch(s.ctrlname, "SBL_load"))
		SBL_loadNodes()
	endif
	return 0
end

function /S SBL_traces()
	
	string s_graph=SBL_getGraph()
	string s_list=TraceNameList(s_graph,";",1+4)
	variable i=0
	do
		if(stringmatch(stringfromlist(i,s_list), "w_base"))
			s_list=RemoveListItem(i, s_list)
		elseif(stringmatch(stringfromlist(i,s_list), "*_sub"))
			s_list=RemoveListItem(i, s_list)	
		elseif(stringmatch(stringfromlist(i,s_list), "*_BL"))
			s_list=RemoveListItem(i, s_list)
		elseif(stringmatch(stringfromlist(i,s_list), "w_NodesY"))
			s_list=RemoveListItem(i, s_list)
		else
			i+=1
		endif
	while (i<itemsinlist(s_list))
	return s_list
end

// initialize (or reinitialize) the package data folder
function /DF SBL_resetDFR()
	
	NewDataFolder /O root:Packages
	NewDataFolder /O root:Packages:SplineFit	
	DFREF dfr = root:Packages:SplineFit	

	wave /Z W_nodesX=dfr:W_nodesX
	if (WaveExists(W_nodesX)==0)	
		make /O/n=5 dfr:W_nodesX, dfr:W_nodesY
	endif
	make /O/n=1  dfr:W_spline_dependency={nan}	
	return dfr
end

function/DF SBL_getDFREF()
	DFREF dfr = root:Packages:SplineFit
	if (DataFolderRefStatus(dfr) != 1)	
		DFREF dfr = SBL_resetDFR()
	endif
	return dfr
end

// dependency function to define the spline curve
function SBL_update(w)
	wave w
	
	DFREF dfr = SBL_getDFREF()
	wave /Z/SDFR=dfr W_nodesX,W_nodesY,W_base, W_sub
	
	string s_graph=SBL_getGraph()
	if (strlen(s_graph)==0)
		return 0
	endif
		
	controlInfo /W=$s_graph+"#SBL_panel" popTrace
	string s_trace=s_value
	wave /Z w_data=TraceNameToWaveRef(s_graph, s_trace)
	wave /Z w_x=XWaveRefFromTrace(s_graph, s_trace)
	
	if (waveexists(w_data)==0)
		return nan
	endif
	
	if(numpnts(W_nodesX)<2)
		return nan
	endif
	
	// flags for interpolate2
	variable flagE, flagF, type=1 // set default to linear interpolation
	
	controlInfo /W=$s_graph+"#SBL_panel" check_akima
	if (v_value)
		type = 0
	endif
	
	if (type==0 && numpnts(W_nodesX)<5)
		doalert 0, "Akima spline requires at least 5 nodes\r - reverting to standard cubic spline"
		Checkbox check_akima, win=$s_graph+"#SBL_panel", value=0	
		Checkbox check_cubic, win=$s_graph+"#SBL_panel", value=1
		type=2
	endif
	
	if (type) // not an Akima spline
		do
			// collect info to set the flags for interpolate2
			controlInfo /W=$s_graph+"#SBL_panel" popE
			flagE=V_Value
			controlInfo /W=$s_graph+"#SBL_panel" check_cubic
			if (v_value)
				type=2
				break
			endif
			controlInfo /W=$s_graph+"#SBL_panel" check_smoothing
			if(v_value)
				type=3
				controlInfo /W=$s_graph+"#SBL_panel" setvarF
				flagF=V_value
				break
			endif
		while(0)
	endif
	
	// check for small number of points
	if (numpnts(W_nodesX)<4)
		type=1 // don't uncheck spline type, so will revert to desired type as nodes are added
	endif

	// do the spline calculation 
	if (type==0) // Akima type spline
		if (SBL_CalcIota(W_nodesX, W_nodesY))
			wave w_iota=dfr:w_iota
			if(waveexists(w_x))
				MultiThread W_base=SBL_Akima(w_x[p], W_nodesX, W_nodesY, w_iota)
			else
				MultiThread W_base=SBL_Akima(x, W_nodesX, W_nodesY, w_iota)
			endif
		endif
	else // cubic spline
		DebuggerOptions
		variable sav_debug=V_debugOnError 
		DebuggerOptions debugOnError=0 // switch this off in case interpolate fails
		try
			if(waveexists(w_x))
				Interpolate2 /E=(flagE) /F=(flagF) /T=(type) /I=3  /X=w_x /Y=W_base W_nodesX, W_nodesY; AbortOnRTE
			else
				Interpolate2 /E=(flagE) /F=(flagF) /T=(type) /I=3  /Y=W_base  W_nodesX, W_nodesY; AbortOnRTE
			endif
		catch
			Print "Error during interpolate:"
			Variable CFerror = GetRTError(1)
			Print GetErrMessage(CFerror)
		endtry
		DebuggerOptions debugOnError=sav_debug
	endif	
	
	w_sub=w_data-w_base
	
	return nan	
end

// function to subtract the spline curve
function SBL_subtract()
	
	variable appendBL=0
	
	string s_graph=SBL_getGraph()
	controlInfo /W=$s_graph+"#SBL_panel" popTrace
	string s_trace=s_value
	wave /Z w_data=TraceNameToWaveRef(s_graph, s_trace)
	wave /Z w_x=XWaveRefFromTrace(s_graph, s_trace)
	
	if (waveexists(w_data)==0)
		return nan
	endif
	
	DFREF dfr = SBL_getDFREF()
	wave /SDFR=dfr W_nodesX,W_nodesY,W_base
	
	string strSubtracted, strBaseline, strNodes
	strSubtracted=CleanupName( nameofwave(w_data)+"_sub",0)
	strBaseline=CleanupName( nameofwave(w_data)+"_BL",0)
	strNodes=CleanupName( nameofwave(w_data)+"_SpNodes",0)
	
	DFREF currentDFR=GetDataFolderDFR()
	
	if (exists(strSubtracted))
		doalert 1, strSubtracted+" exists. Overwrite?"
		if(V_flag==2)
			return 0
		endif
	endif
	duplicate /o w_data currentDFR:$strSubtracted /WAVE=nob
	
	if (exists(strBaseline))
		doalert 1, strBaseline+" exists. Overwrite?"
		if(V_flag==2)
			return 0
		endif
	endif
	duplicate /o W_base currentDFR:$strBaseline /WAVE=bl
	
	if (exists(strNodes))
		doalert 1, strNodes+" exists. Overwrite?"
		if(V_flag==2)
			return 0
		endif
	endif
	duplicate /o W_nodesX currentDFR:$strNodes /WAVE=nd
	redimension /N=(-1, 2) nd
	nd[][0]=W_nodesX[p]
	nd[][1]=W_nodesY[p]
	
	nob=w_data-bl
	
	string s_note="Baseline Parameters\r"
	controlInfo /W=$s_graph+"#SBL_panel" check_akima
	if (v_value)
		s_note+="type=AkimaSpline;"
	endif
	controlInfo /W=$s_graph+"#SBL_panel" check_cubic
	if (v_value)
		s_note+="type=CubicSpline;"
		controlInfo /W=$s_graph+"#SBL_panel" popE
		s_note+="Ends="+s_value+";"
	endif
	controlInfo /W=$s_graph+"#SBL_panel" check_smoothing
	if(v_value)
		s_note+="type=SmoothingSpline;"
		controlInfo /W=$s_graph+"#SBL_panel" setvarF
		s_note+="SmoothingFactor="+num2str(v_value)+";"
	endif
	controlInfo /W=$s_graph+"#SBL_panel" check_linear
	if(v_value)
		s_note+="type=LinearSpline;"
	endif
	s_note+="numNodes="+num2str(numpnts(w_nodesX))+";"
	s_note+="data="+nameofwave(w_data)+";output="+nameofwave(nob)+","+nameofwave(bl)+","+nameofwave(nd)+";\r"
	note /k bl, s_note
	note nob, s_note // append note from baseline wave to output wave note
	printf s_note
		
	// add baseline-subtracted output wave to plot
	checkdisplayed /W=$s_graph nob
	if(v_flag==0)
		if(waveexists(w_x))
			appendtograph /W=$s_graph /C=(0,0,0) nob vs w_x
		else
			appendtograph /W=$s_graph /C=(0,0,0) nob
		endif
	endif
	
	// possibly plot baseline
	if(appendBL)
		checkdisplayed /W=$s_graph bl
		if(v_flag==0)
			if(waveexists(w_x))
				appendtograph /W=$s_graph /C=(0,0,0) bl vs w_x
			else
				appendtograph /W=$s_graph /C=(0,0,0) bl
			endif
		endif
	endif	
	
	return 1
end

// function to clear the spline curve
function SBL_clear(killPanel)
	variable killPanel
	
	string s_graph=SBL_getGraph()
	
	if (strlen(s_graph))
		GraphNormal /W=$s_graph
		RemoveFromGraph /W=$s_graph/Z W_spline_dependency,W_base,W_nodesY,w_sub
		setWindow $s_graph hook(SBL_hook)=$""
	endif
	
	DFREF dfr = SBL_getDFREF()
	wave /Z w=dfr:W_spline_dependency
	if (waveexists(w))
		w=nan // remove dependency
	endif
	
	if(killPanel && strlen(s_graph))
		killwindow $s_graph+"#SBL_panel"
	endif

	return 1
end

function SBL_loadNodes()
	DFREF dfr = SBL_getDFREF()
			
	string str_NodeWave
	variable LoadY
	
	prompt str_NodeWave, "Select saved nodes", popup, wavelist("*",";","DIMS:2") 
	prompt LoadY, " ", popup, "Load X positions only;Load X and Y positions;"
	DoPrompt "", str_NodeWave, LoadY
	if (V_Flag)
		return 0
	endif
	
	wave w_nodes=$str_NodeWave
	loadY-=1
	
	wave /SDFR=dfr W_nodesX, W_nodesY
	
	redimension /N=(dimsize(w_nodes, 0)), W_nodesX, W_nodesY
	W_nodesX=w_nodes[p][0]
	
	if (LoadY)
		W_nodesY=w_nodes[p][1]
		return 1
	endif

	SBL_setNodes()

	return 1
end

// set colour wave w_color for visibility against w_target.
function SBL_chooseColour(w_color, s_graph, s_trace)
	wave w_color
	string s_graph, s_trace
	
	// figure out a contrasting colour for baseline and sub waves
	variable c0,c1,c2
	string s_temp=traceinfo(s_graph, s_trace, 0)
	variable v_start=strsearch(s_temp, ";rgb(x)=", 0)
	sscanf s_temp[v_start,strlen(s_temp)-1], ";rgb(x)=(%d,%d,%d*", c0,c1,c2
	make /free w_cTrace={c0,c1,c2}, w_index={1,2,3}
	MakeIndex /R w_cTrace, w_index
	IndexSort  w_index, w_color 
	return 1	
end

// **********     Akima spline code     *************

// this procedure contains code posted to IgorExchange by Michael Bongard
// a couple of edits by tony.withers@uwo.ca to allow extrapolation beyond end nodes
 
// Akima.ipf: Routines to implement Akima-spline fitting, based on
// H. Akima, JACM, Vol 17, No 4, 1970 p 589-602
// M. Bongard, 11/17/09

// First, call calcIota to generate interpolation information; 
// then you can interpolate using Akima's spline method with the akima() function.
 
Function SBL_CalcIota(knotX, knotY)
	WAVE knotX // knot X locations
	WAVE knotY // knot Y locations
 	
 	// removed some sanity checks
 	DFREF dfr = SBL_getDFREF()
 	
	Variable numKnots = numpnts(knotX)
 
	if(numKnots < 5)
		Print "calcIota: ERROR -- Akima spline algorithm requires at least 5 knots. Aborting..."
		return -1
	endif
 
	// Make intermediate ai, bi, mi arrays
	Make/D/FREE/N=(numKnots + 4)  kX, kY
 
	Variable i, j
	for(i = 2, j = 0; j < numKnots; i += 1, j += 1)
		kX[i] = knotX[j]
		kY[i] = knotY[j]
	endfor
 
	// Handle end-point extrapolation
	Make/D/FREE/N=5 endX, endY
	// RHS: end points are last three in dataset
	Variable endStartPt
	// RHS
	endStartPt = numPnts(kX) - 5
	endX = kX[p + endStartPt]
	endY = kY[p + endStartPt]
 
	SBL_CalcEndPoints(endX, endY)
 
	kX[numpnts(kX)-2] = endX[3]
	kX[numpnts(kX)-1] = endX[4]
	kY[numpnts(kX)-2] = endY[3]
	kY[numpnts(kX)-1] = endY[4]
 
	// LHS: end points are first three in dataset, but reversed in ordering
	// (i.e. point 3 in Akima's notation == index 0)
	endX = 0
	endY = 0
 
	for(i = 0, j = 2; i < 3; i += 1, j -= 1)
		endX[j] = knotX[i]
		endY[j] = knotY[i]
	endfor
 
	SBL_CalcEndPoints(endX, endY)
 
	kX[1] = endX[3]
	kX[0] = endX[4]
	kY[1] = endY[3]
	kY[0] = endY[4]
 
	// kX, kY are now properly populated, along with all necessary extrapolated endpoints
	// computed as specified in Akima1970
 
	Make/D/FREE/N=(numKnots + 4 - 1)  mK
	mK = (kY[p + 1] - kY[p]) / (kX[p + 1] - kX[p])	
 
	Make/O/N=(numKnots) dfr:w_Iota /wave=w_iota
 
	Variable denom, m1,m2,m3,m4
	for(i = 2, j = 0; j < numKnots; i += 1, j += 1)
		m1 = mK[i - 2]
		m2 = mK[i - 1]
		m3 = mK[i]
		m4 = mK[i + 1]
 
		denom = abs(m4 - m3) + abs(m2 - m1)
		if(denom == 0)
			w_Iota[j] = 0.5 * (m2 + m3)
			continue
		endif
 
		w_Iota[j] = ( abs(m4 - m3)*m2 + abs(m2 - m1)*m3 ) / denom
	endfor
	return 1
End
 
// Given: 5-point knot wave knotX, knotY, with i=[0,2] representing the last three
// knot locations from data, compute end knots i=[3,4] appropriately.
ThreadSafe Function SBL_CalcEndPoints(kX, kY)
	WAVE kX, kY // knot X,Y coordinate locations, respectively
 
	// Sanity checks
	if( (numPnts(kX) != numpnts(kY)) || numpnts(kX) != 5)
		Print "calcEndPoints: ERROR -- must have 5 points in knot wave! Aborting..."
		return -1
	endif
 
	// First, compute X locations of knots, according to relations in eq. 8 of Akima1970:
	kX[3] = kX[1] + kX[2] - kX[0]
	kX[4] = 2*kX[2] - kX[0]
 
	// Now all kX are known, so let's set up the line segment slope waves
	// ai, bi, mi of eq. (12)-(14).
	Make/N=4/FREE ai, bi, mi
 
	ai = kX[p+1] - kX[p]
	// ai is now determined completely
 
	bi = kY[p + 1] - kY[p]
	mi = bi/ai
	// bi, mi determined on i=[0,1]	
 
	// Determine remainder of quantities by applying solutions of eq (9)
	kY[3] = (2*mi[1] - mi[0])*(kX[3] - kX[2]) + kY[2]
	mi[2] = (kY[3] - kY[2]) / (kX[3] - kX[2])
 
	kY[4] = (2*mi[2] - mi[1])*(kX[4] - kX[3]) + kY[3]
	mi[3] = (kY[4] - kY[3]) / (kX[4] - kX[3])
End
 
ThreadSafe Function SBL_akima(x, knotX, knotY, knotIota)
	Variable x
	WAVE knotX, knotY, knotIota
 
	Variable i1, i2
 
	// Find where x
	Variable done = 0
	i2 = -1
	Variable numKnots = numpnts(knotX)
	do
		i2 += 1
 
		if(knotX[i2] == x)
			return knotY[i2]
		endif
 
		done = (knotX[i2] > x) ? 1 : 0
 
	while(!done && (i2 < numKnots-1)) // edited this to allow extrapolation beyond end knot
	i1 = i2 - 1
	
	// this edit to allow extrapolation beyond start knot
	if (i1<0)
		i1+=1
		i2+=1
	endif
 
	Variable x1, x2, y1, y2, iota1, iota2
	x1 = knotX[i1]
	y1 = knotY[i1]
	iota1 = knotIota[i1]
 
	x2 = knotX[i2]
	y2 = knotY[i2]
	iota2 = knotIota[i2]
 
	Variable p0, p1, p2, p3, tmp
	p0 = y1
	p1 = iota1
	p2 = ( 3 * (y2 - y1)/(x2 - x1) - 2*iota1 - iota2 ) / (x2 - x1)
	p3 = (iota1 + iota2 - 2 * (y2 - y1)/(x2 - x1)) / (x2 - x1)^2
 
	tmp = x - x1
 
	return p0 + p1 * tmp + p2 * tmp^2 + p3*tmp^3
End


// for Igor 6 compatibility
function /T SBL_getGraph()

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