#pragma TextEncoding = "UTF-8" #pragma rtGlobals=3 #pragma version=2.02 #pragma IgorVersion=7 #pragma moduleName=archull #include // 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. // Usage // Arc Hull requires that a spectrum is plotted in the top graph. // Selecting Macros - Arc Hull will add a panel on the left side of the // graph. You can choose 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. // 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 _base 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. // Version history // Version 2.02 decreased panel size for compatibility with smaller or // lower resolution screens. And a bug fix. // Version 2.0 // Introduces iterative fitting of CurveFit functions with elimination // of either positive or negative residual points // Version 1.3, 1.31, 1.32 // Adds 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. Set UPDATE_ON_ZOOM=0 to prevent // this behaviour: constant UPDATE_ON_ZOOM=1 Menu "Analysis" Submenu "Packages" "Arc Hull",/Q, CCH_init() end end Menu "Macros" "Arc Hull",/Q, CCH_init() end function CCH_init() if (strlen(WinList("*",";","WIN:1"))==0) // no graphs doalert 0, "Arc hull baseline requires a trace plotted in a graph window." return 0 endif // 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 function CCH_makePanel() // make sure the top graph is visible string s_graph=WinName(0,1) dowindow /F $s_graph // 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 $s_graph wsize if (V_left<200) movewindow /W=$s_graph 200, V_top, -1, -1 endif // make panel NewPanel /K=1/N=CCH_panel/W=(200,0,0,355)/HOST=$s_graph/EXT=1 as "Arc Hull Controls" ModifyPanel /W=$s_graph#CCH_panel, noEdit=1 string s_trace=stringfromlist(0, CCH_traces()) wave /Z w=TraceNameToWaveRef(s_graph, s_trace) 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++ // store values internally in these controls PopupMenu popTrace, mode=1, Value=CCH_traces(), title="",pos={v_left,v_top+deltaY*i},size={130,20} PopupMenu popTrace, help={"select data wave" }, proc=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++ Checkbox check_CCH, pos={v_left,v_top+deltaY*i}, fSize=font, value=1, proc=CCH_checkbox, title="Arc hull" Checkbox check_CCH, mode=1,help={"Arc hull type baseline"} i++ Checkbox check_spline, pos={v_left,v_top+deltaY*i}, fSize=font, value=0, proc=CCH_checkbox, title="Hull spline" Checkbox check_spline, mode=1,help={"cubic spline through nodes from convex hull"} i++ Checkbox check_itFit, pos={v_left,v_top+deltaY*i}, fSize=font, value=0, proc=CCH_checkbox, title="Iterative" Checkbox check_itFit, mode=1,help={"iterative fit with elimination"} i++ 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=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=CCH_popup, fSize=font, disable = 1 PopupMenu popFitFunc, Value="line;poly 3;poly 4;gauss;lor;exp;dblexp;sin;hillequation;sigmoid;power;lognormal" i++ 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=CCH_SetVar i++ 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=CCH_SetVar i++ CheckBox check_neg, pos={v_left,v_top+deltaY*i}, fSize=font, value=0, proc=CCH_checkbox, title="Negative peaks" CheckBox check_neg, help={"Fit baseline to top of convex hull"} i++ CheckBox check_win, pos={v_left,v_top+deltaY*i}, fSize=font, value=0, proc=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=CCH_buttons Button CCH_sub,help={"Subtract baseline and save a copy"}, fSize=font i++ Button CCH_all,pos={65,v_top+deltaY*i},size={70,20},title="All in one", proc=CCH_buttons Button CCH_all,help={"Subtract baseline from all traces"}, fSize=font if (strlen(s_trace)==0) return 0 endif return 1 // CCH_updatePlot() // setwindow $s_graph hook(ArcHullHook)=CCH_graphHook end Function CCH_graphHook(s) STRUCT WMWinHookStruct &s Variable hookResult = 0 switch(s.eventCode) case 8: if (UPDATE_ON_ZOOM) CCH_doFit() endif break endswitch return hookResult // 0 if nothing done, else 1 End // update baseline and ensure package waves are displayed function CCH_updatePlot() DFREF dfr = root:Packages:ArcHull GetWindow /Z CCH_panel activeSW string s_graph=parseFilepath(0, s_value, "#", 0, 0) // ***** // check traceinfo for axis flags, give up if not default // ***** removefromgraph /W=$s_graph /Z Arc_Hull_Base, Arc_Hull_Sub, W_YHull controlInfo /W=CCH_panel popTrace string s_trace=s_value // do fit to create baseline CCH_doFit() wave /SDFR=dfr w_base, w_sub wave /Z w_dataX=XWaveRefFromTrace(s_graph, s_trace) if(waveexists(w_dataX)) wave w_baseX=dfr:w_baseX appendtograph /W=$s_graph w_base/TN=Arc_Hull_Base vs w_baseX appendtograph /W=$s_graph w_sub/TN=Arc_Hull_Sub vs w_baseX else AppendToGraph /W=$s_graph w_base/TN=Arc_Hull_Base, w_sub/TN=Arc_Hull_Sub endif // if data are offset in y, then apply offset to baseline string s_info=traceinfo(s_graph,s_trace,0) variable trace_offset=GetNumFromModifyStr(s_info,"offset","{",1) ModifyGraph offset(Arc_Hull_Base)={0,trace_offset} wave W_Xhull=W_Xhull, W_Yhull=W_Yhull controlInfo /W=CCH_panel check_spline if(v_value) checkdisplayed /W=$s_graph W_Yhull if (v_flag==0) AppendToGraph /W=$s_graph W_Yhull vs W_XHull ModifyGraph /W=$s_graph 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=$s_graph /Z W_YHull endif setwindow $s_graph hook(ArcHullHook)=CCH_graphHook return 1 end function CCH_checkbox(s) STRUCT WMCheckboxAction &s string s_graph=parseFilepath(0, s.win, "#", 0, 0) if(s.eventCode==-1) removefromgraph /W=$s_graph /Z Arc_Hull_Base, Arc_Hull_Sub, W_YHull setwindow $s_graph 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(s_graph, 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 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 s_graph=parseFilepath(0, s.win, "#", 0, 0) controlinfo /W= $s.win popTrace wave /Z w=TraceNameToWaveRef(s_graph, 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 Function CCH_popup(s) STRUCT WMPopupAction &s if(s.eventCode==-1) return 0 endif string s_graph=parseFilepath(0, s.win, "#", 0, 0) controlinfo /W= $s.win popTrace wave /Z w=TraceNameToWaveRef(s_graph, 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 Function CCH_buttons(s) : ButtonControl STRUCT WMButtonAction &s if(s.eventCode!=2) return 0 endif string s_graph=parseFilepath(0, s.win, "#", 0, 0) string s_trace strswitch(s.ctrlName) case "CCH_sub": 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 s_list=CCH_traces() PopupMenu popTrace, win=$s.win, Value=CCH_traces() variable i for(i=0;iw_base) ? w_base : w_ItFit if (negPeaks) w_next = (w_ItFit10 && v_npnts<=v_conv) // converged break else v_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 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, v_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=CCH_panel 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 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, v_top]) wave /Z w, w_x variable depth, v_top variable useXwave=0 if (paramisdefault(w_x)==0 && waveexists(w_x)) useXwave=1 endif v_top = paramisdefault(v_top) ? 0 : v_top // v_top=1 for top hull // add concave function wavestats /Q w 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 (v_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 function CCH_FitWrap(w, w_x, w_out, fName) wave /Z w, w_x, w_out string fName variable V_FitError=0 try strswitch (fName) case "line": CurveFit /Q line, w /D=w_out /X=w_x /NWOK; AbortOnRTE break case "poly 3": CurveFit /Q poly 3, w /D=w_out /X=w_x /NWOK; AbortOnRTE break case "poly 4": CurveFit /Q poly 4, w /D=w_out /X=w_x /NWOK; AbortOnRTE break case "gauss": CurveFit /Q gauss, w /D=w_out /X=w_x /NWOK; AbortOnRTE break case "lor": CurveFit /Q lor, w /D=w_out /X=w_x /NWOK; AbortOnRTE break case "exp": CurveFit /Q exp, w /D=w_out /X=w_x /NWOK; AbortOnRTE break case "dblexp": CurveFit /Q dblexp, w /D=w_out /X=w_x /NWOK; AbortOnRTE break case "sin": CurveFit /Q sin, w /D=w_out /X=w_x /NWOK; AbortOnRTE break case "hillequation": CurveFit /Q hillequation, w /D=w_out /X=w_x /NWOK; AbortOnRTE break case "sigmoid": CurveFit /Q sigmoid, w /D=w_out /X=w_x /NWOK; AbortOnRTE break case "power": CurveFit /Q power, w /D=w_out /X=w_x /NWOK; AbortOnRTE break case "lognormal": CurveFit /Q 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 return (V_FitError ==0) end function CCH_subtract(overwrite) variable overwrite variable appendBL=0 controlinfo /W=CCH_panel popTrace string s_trace=s_value GetWindow /Z CCH_panel activeSW string s_graph=parseFilepath(0, s_value, "#", 0, 0) wave /Z w=TraceNameToWaveRef(s_graph, s_trace) if (waveexists(w)==0) return 0 endif wave /Z w_dataX=XWaveRefFromTrace(s_graph, s_trace) 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=$s_graph /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=$s_graph /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)