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

#include <Readback ModifyStr>

// Written by tony withers, tony.withers@uwo.ca. I would be very happy to
// hear from you if you find this package useful, or if you have any
// suggestions for improvement. 

// 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 updates 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. Cubic and smoothing splines are 
// calculated by Igor's Interpolate2 operation; 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.

// PCHIP: fit a cubic Hermite spline through the nodes. This spline type 
// doesn't overshoot local maxima or minima in node sequence. First 
// derivative is continuous through the nodes, but second derivative may 
// be discontinuous.

// Cubic: fit a 'normal' cubic spline through the nodes. Second 
// derivative is continuous though 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.10 2/2/18 Introduces piecewise cubic Hermite spline
// 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_base, w_sub
	
// not sure why	the dependency wave was always plotted?
// 	wave /SDFR=dfr w_spline_dependency
//	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_RGB={0,10000,65535}
	SBL_chooseColour(w_RGB, s_graph, s_trace)
	
	ModifyGraph /W=$s_graph offset(W_base)={0,trace_offset}
	ModifyGraph /W=$s_graph rgb(W_base)=(w_RGB[0],w_RGB[1],w_RGB[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,380)/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*6.0},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_PCHIP, pos={v_left,v_top+deltaY*i}, fSize=font, value=0, proc=SBL_checkboxes, title="PCHIP"
	Checkbox check_PCHIP, mode=1,help={"use piecewise cubic Hermite 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
	
	string s_typeButtons="check_akima;check_PCHIP;check_cubic;check_linear;check_smoothing"
	variable i
	variable numTypes=itemsinlist(s_typeButtons)
	variable buttonNum=WhichListItem(s.ctrlName, s_typeButtons)

	if(buttonNum>=0)
		for (i=0;i<numTypes;i+=1)
			if(i!=buttonNum)
				Checkbox $StringFromList(i, s_typeButtons), win=$s.win, value=0
			endif
		endfor
		SBL_Update({nan})
		return 1
	endif
	
	// only remaining checkbox is 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

	return 1
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)
	string s_remove="w_base;w_NodesY;"
	s_remove+=listmatch(s_list, "*_sub")
	s_remove+=listmatch(s_list, "*_BL")
	s_list=RemoveFromList(s_remove, s_list, ";", 0)
	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
	// if type < 1, spline is calculated using a method other than Igor's interpolate2
	
	controlInfo /W=$s_graph+"#SBL_panel" check_akima
	if (v_value)
		type = 0
	endif
	
	controlInfo /W=$s_graph+"#SBL_panel" check_PCHIP
	if (v_value)
		type = -1
	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
		
	// check for small number of points
	if (type>0 && numpnts(W_nodesX)<4)
		type=1 // don't uncheck spline type, so will revert to desired type as nodes are added
	endif
		
	if (type>0) // types that use interpolate2
		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

	// do the spline calculation 
	if (type==-1) // PCHIP
		SBL_PCHIP(w_nodesX, w_nodesY, W_base, xwave=w_x)
	elseif (type==0) // Akima type spline
		SBL_Akima(W_nodesX, W_nodesY, w_base, xwave=w_x)	
	elseif(type > 0) // 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
		RemoveFromGraph /W=$s_graph/Z 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.
static 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 part based on code posted to IgorExchange by Michael Bongard
// edited 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
// after M. Bongard, 11/17/09

Function /wave SBL_AkimaIota(knotX, knotY)
	wave knotX, knotY
 	
 	// removed some sanity checks, converted to wave function
	Variable numKnots = numpnts(knotX)
 
	if(numKnots < 5)
		Print "calcIota: ERROR -- Akima spline algorithm requires at least 5 knots. Aborting..."
		return $""
	endif
 	 
	Make/D/FREE/N=(numKnots)  kX, kY
	kX = knotX
	kY = knotY
 
	// Handle end-point extrapolation
	Make/D/FREE/N=5 endX, endY
	
	// RHS: end points are last three in dataset
	endX[0,2] = kX[p+numKnots-3]
	endY[0,2] = kY[p+numKnots-3]
 
	SBL_AkimaEndPoints(endX, endY)
 	
 	kX[numpnts(kX)]={endX[3],endX[4]}
 	kY[numpnts(kY)]={endY[3],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
	
	endX[0,2] = knotX[2-p]
	endY[0,2] = knotY[2-p]
 
	SBL_AkimaEndPoints(endX, endY)
	
	InsertPoints  0, 2, kX, kY
	kX[0,1] = endX[4-p]
	kY[0,1] = endY[4-p]
 
	// kX, kY are now properly populated, including extrapolated endpoints
 
	Make/D/FREE/N=(numKnots + 4 - 1)  mK
	mK = (kY[p+1]-kY)/(kX[p+1]-kX)	
 
	Make/free/N=(numKnots) w_Iota
	w_Iota= ( abs(mK[p+3]-mK[p+2])*mK[p+1] + abs(mK[p+1]-mK[p])*mK[p+2] ) / (abs(mK[p+3]-mK[p+2])+abs(mK[p+1]-mK))
	w_Iota= numtype(w_Iota)==0 ? w_Iota : 0.5*(mK[p+1]+mK[p+2])
	return w_Iota
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_AkimaEndPoints(kX, kY)
	WAVE kX, kY // node X,Y coordinate locations, respectively
 
	// compute X locations of knots, eq. 8 of Akima 1970:
	kX[3] = kX[1] + kX[2] - kX[0]
	kX[4] = 2*kX[2] - kX[0]
 
	// eq. (12)-(14), determine gradient on [0,1]
	Make/N=4/FREE  mi	
	mi = (kY[p+1]-kY)/(kX[p+1]-kX)
 
	// 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]
	// not sure why the next line was included
	//mi[3] = (kY[4] - kY[3]) / (kX[4] - kX[3])
End


function SBL_Akima(W_nodesX, W_nodesY, w_interp [, xwave])
	wave W_nodesX, W_nodesY, w_interp
	wave/Z xwave
	
	wave w_iota=SBL_AkimaIota(W_nodesX, W_nodesY)
	if (waveexists(w_iota))
		if(waveexists(xwave))
			MultiThread w_interp=SBL_AkimaThread(xwave[p], W_nodesX, W_nodesY, w_iota)
		else
			MultiThread w_interp=SBL_AkimaThread(x, W_nodesX, W_nodesY, w_iota)
		endif
		return 1
	endif	
	return 0
end


ThreadSafe Function SBL_AkimaThread(x, W_nodesX, W_nodesY, W_Iota)
	Variable x
	WAVE W_nodesX, W_nodesY, W_Iota
 	
	variable p0, p1
	if(x<wavemin(W_nodesX))
		p0=0
	elseif(x>=wavemax(W_nodesX))
		p0=numpnts(W_nodesX)-2
	else
		findlevel /P/EDGE=1/Q W_nodesX, x
		if(v_flag==1) // not found			
			return nan
		endif
		p0=floor(v_levelX)
	endif	
	p1=p0+1

	Variable x1, x2, y1, y2, iota1, iota2
	x1 = W_nodesX[p0]
	y1 = W_nodesY[p0]
	iota1 = W_Iota[p0]
 
	x2 = W_nodesX[p1]
	y2 = W_nodesY[p1]
	iota2 = W_Iota[p1]
 
	Variable coeff2, coeff3, delx
	coeff2 = ( 3 * (y2 - y1)/(x2 - x1) - 2*iota1 - iota2 ) / (x2 - x1)
	coeff3 = (iota1 + iota2 - 2 * (y2 - y1)/(x2 - x1)) / (x2 - x1)^2
	delx = x - x1
 
	return y1 + iota1*delx + coeff2*delx^2 + coeff3*delx^3
End

// ----------------------------- PCHIP ---------------------------

// populate w_interp with interpolated values from piecewise cubic 
// Hermite spline
function SBL_PCHIP(w_nodesX, w_nodesY, w_interp, [xwave])	
	wave w_nodesX, w_nodesY, w_interp
	wave /Z xwave
	
	if(numpnts(w_nodesX)==2) // linear interpolation
		if(waveexists(xwave))	
			w_interp=w_nodesY[0]+(xwave[p]-w_nodesX[0])/(w_nodesX[1]-w_nodesX[0])*(w_nodesY[1]-w_nodesY[0])
		else
			w_interp=w_nodesY[0]+(x-w_nodesX[0])/(w_nodesX[1]-w_nodesX[0])*(w_nodesY[1]-w_nodesY[0])
		endif
		return 1
	endif
	
	wave w_d=SBL_setDerivativesPCHIP(w_nodesX, w_nodesY)	
	if(waveexists(xwave))	
		w_interp=SBL_InterpolatePCHIP(w_nodesX, w_nodesY, w_d, xwave[p])
	else
		w_interp=SBL_InterpolatePCHIP(w_nodesX, w_nodesY, w_d, x)
	endif	
end

// The trick is to choose good values for the slope at each node. See:
// FN Fritsch and RE Carlson (1980) Monotone piecewise cubic 
// interpolation. SIAM J. Numer. Anal. 17: 238. DOI:10.1137/0717021 
// KW Brodlie (1980) A review of methods for curve and function drawing. 
// In Mathematical Methods in Computer Graphics and Design, KW Brodlie, 
// ed., Academic Press, London, pp. 1-37.
// FN Fritsch and J Butland (1984) A method for constructing local 
// monotone piecewise cubic interpolants, SIAM Journal on Scientific and 
// Statistical Computing 5: 300-304.
static function /WAVE SBL_setDerivativesPCHIP(w_nodesX, w_nodesY)
	wave w_nodesX, w_nodesY

	duplicate /free w_nodesX, w_d, w_a
	// w_d will be set to desired derivatives at node positions
	// w_a will be used as weighting for harmonic means
	w_d=0
	
	make /free/n=(numpnts(w_nodesX)-1), w_m, w_delx
	// w_m[i] is gradient between node i and i+1
	// w_delx[i] and w_delx[i-1] are x offsets of nodes i+1 and i-1	
	
	w_m = (w_nodesY[p+1]-w_nodesY)/(w_nodesX[p+1]-w_nodesX)
	w_delx = w_nodesX[p+1]-w_nodesX
	variable pmax=numpnts(w_d)-2  // pmax is numpnts(w_m)-1 = numpnts(w_d)-2
	
	w_a[1,pmax] = (1+(w_delx[p-1])/(w_delx[p-1]+w_delx))/3	
	
	// Brodlie modification of Butland formula (for same sign slopes)	
	w_d [1,pmax] = ( w_m!=0 && w_m[p-1]!=0 && (sign(w_m)==sign(w_m[p-1])) ) ? w_m[p-1]*w_m/(w_a*w_m[p-1]+(1-w_a)*w_m) : 0
		
	// deal with ends
	w_d[0]=SBL_EndDerivativePCHIP(w_delx[0], w_delx[1], w_m[0], w_m[1])
	w_d[pmax+1]=SBL_EndDerivativePCHIP(w_delx[pmax], w_delx[pmax-1], w_m[pmax], w_m[pmax-1])
	
	// clean up any division by zero errors
	w_d = numtype(w_d)==0 ? w_d : 0
			
	return w_d
end

// one-sided three-point estimate for the derivative adjusted to be shape
// preserving
static function SBL_EndDerivativePCHIP(delX0, delX1, m0, m1)
	variable delX0, delX1, m0, m1
	
	variable derivative
	
	derivative = ( m0*(2*delX0+delX1) - delX0*m1 ) / (delX0+delX1)
	if(sign(derivative)!=sign(m1))
		derivative=0
	elseif( (sign(m0)!=sign(m1)) && (abs(derivative)>3*abs(m0)) )
		derivative=3*m0
	endif
	return derivative
end

threadsafe function SBL_InterpolatePCHIP(w_nodesX, w_nodesY, w_nodesM, x)
	wave w_nodesX, w_nodesY, w_nodesM
	variable x
		
	variable p0, p1
	if(x<wavemin(w_nodesX))
		p0=0
	elseif(x>=wavemax(w_nodesX))
		p0=numpnts(w_nodesX)-2
	else
		findlevel /P/EDGE=1/Q w_nodesX, x
		if(v_flag==1) // not found			
			return nan
		endif
		p0=floor(v_levelX)
	endif
	
	p1=p0+1		
	return SBL_InterpolateHermite(x, w_nodesX[p0], w_nodesY[p0], w_nodesX[p1], w_nodesY[p1], w_nodesM[p0], w_nodesM[p1])
end

// https: en.wikipedia.org/wiki/Cubic_Hermite_spline
// calculate cubic function f(x) such that f(x1)=y1, f(x2)=y2, f'(x1)=m1 and f'(x2)=m2
threadsafe function SBL_InterpolateHermite(x, x0, y0, x1, y1, m0, m1)
	variable x, x0, y0, x1, y1, m0, m1
	
	variable t=(x-x0)/(x1-x0)
	m0*=(x1-x0)
	m1*=(x1-x0)
	return y0*(2*t^3-3*t^2+1) + m0*(t^3-2*t^2+t) + y1*(-2*t^3+3*t^2) + m1*(t^3-t^2)	
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