#pragma rtGlobals=3		// Use modern global access method.
#pragma version=5.07

// Peak AutoFind
// Version 5.04, 12/14/2004, JP - Added AutomaticallyFindPeaks.
// Version 5.05, 10/31/2007, JW - Added right and left width estimates to the output wave. Code that doesn't know about them simply won't access those columns.
// Version 5.06, 9/16/2008, JW- Changed algorithm somewhat: 
// Original:
//		find minimum value in the smoothed 2nd derivative
//		require minimum to be negative
//		use FindLevel to find zero-crossing in the 2nd derivative in each direction
// New:
//		Use FindPeak to find local minimum
//		Use FindPeak in each direction from found minimum to find local maxima on either side
// The Original algorithm will often fail to find overlapped peaks because of the requirement that the minimum in the 2nd derivative be actually negative.
// Another tweak is that the New algorithm, when estimating peak height, uses the minimum level found on either side of the peak, rather than
// the level defined by the line connecting the minima. The sloping baseline also discriminates against overlapped peaks.

// You can choose the Original algorithm by creating a global variable in the root: data folder called V_PeakAutoFindClassic and setting its value to 1.

// Version 5.07, 8/12/2010, JW: Bug fix for rtGlobals=3; EstPeakNoiseAndSmfact could set imax to a number with fraction < .5, resulting in
// first loop running one too many times. This happened with a very smooth wave having 128 points.

Menu "Analysis"
	"Automatically Find Peaks", /Q, AutomaticallyFindPeaks()
EndMacro

Function AutomaticallyFindPeaks()

	String wname, xdata="_calculated_"
	Variable maxPeaks=100, minPeakPercent=5
	Prompt wname, "Peak Wave", popup, WaveList("*",";","DIMS:1,TEXT:0,CMPLX:0")+"_none_;"
	Prompt xdata, "X values", popup, "_calculated_;"+WaveList("*",";","DIMS:1,TEXT:0,CMPLX:0")
	Prompt maxPeaks, "Maximum Peaks"
	Prompt minPeakPercent, "Minimum Peak Amplitude (% max)"
	DoPrompt "Automatically Find Peaks", wname, xdata, maxPeaks, minPeakPercent
	if( V_Flag != 0 )
		return 0	// user cancelled
	endif
	
	WAVE/Z w=$wname
	WAVE/Z wx=$xdata
	Variable pBegin=0, pEnd= numpnts(w)-1
	Variable/C estimates= EstPeakNoiseAndSmfact(w,pBegin, pEnd)
	Variable noiselevel=real(estimates)
	Variable smoothingFactor=imag(estimates)
	
	Variable peaksFound= AutoFindPeaks(w,pBegin,pEnd,noiseLevel,smoothingFactor,maxPeaks)
	if( peaksFound > 0 )
		WAVE W_AutoPeakInfo
		// Remove too-small peaks
		peaksFound= TrimAmpAutoPeakInfo(W_AutoPeakInfo,minPeakPercent/100)
		if( peaksFound > 0 )
			// Make waves to display in a graph
			// The x values in W_AutoPeakInfo are still actually points, not X
			Make/O/N=(peaksFound) WA_PeakCentersY = w[W_AutoPeakInfo[p][0]]
			AdjustAutoPeakInfoForX(W_AutoPeakInfo,w,wx)
			Make/O/N=(peaksFound) WA_PeakCentersX = W_AutoPeakInfo[p][0]

			// Show W_AutoPeakInfo in a table, with dimension labels
			SetDimLabel 1, 0, center, W_AutoPeakInfo
			SetDimLabel 1, 1, width, W_AutoPeakInfo
			SetDimLabel 1, 2, height, W_AutoPeakInfo
			CheckDisplayed/A W_AutoPeakInfo
			if( V_Flag == 0 )
				Edit W_AutoPeakInfo.ld
			endif

			DoWindow ShowPeaks
			if( V_Flag == 0 )
				if( WaveExists(wx) )
					Display/N=ShowPeaks w vs wx				
				else
					Display/N=ShowPeaks w
				endif
				AppendToGraph/W=ShowPeaks WA_PeakCentersY vs WA_PeakCentersX
				ModifyGraph/W=ShowPeaks rgb(WA_PeakCentersY)=(0,0,65535)
				ModifyGraph/W=ShowPeaks mode(WA_PeakCentersY)=3
				ModifyGraph/W=ShowPeaks marker(WA_PeakCentersY)=19
			endif
		endif
	endif
	if( peaksFound < 1 )
		DoAlert 0, "No Peaks found!"
	endif
	return peaksFound
End

Function/C EstPeakNoiseAndSmfact(w,pBegin,pEnd)
	Wave w
	Variable pBegin,pEnd
	
	if (abs(pBegin-pEnd) < 10)			// 10 is pretty arbitrary; this is intended to avoid trying to apply this test to unreasonably small waves. Even a 10-point wave is probably a mistake: a fit coefficient wave was selected by mistake or something.
		return cmplx(0,0)
	endif
		
	if( pBegin > pEnd )
		Variable tmp= pBegin
		pBegin= pEnd
		pEnd= tmp
	endif
	
	NewDataFolder/S/O afpTemp
	Duplicate/O/R=[pBegin,pEnd] w,wtmp
	Differentiate wtmp
	Duplicate/O wtmp, wtmpOrigDif

	Make/O/N=1000 hist
	Histogram/B=1 wtmp,hist
	Integrate hist

	FindLevel/Q hist,0.4*hist[999]
	Variable x0= V_LevelX
	FindLevel/Q hist,0.6*hist[999]
	Variable x1= V_LevelX
	Variable noiselevel= abs(2*(x1-x0)*deltax(w))
	Variable snr=(pnt2x(hist, 999)-pnt2x(hist, 0))/(x1-x0)
	
	if( pEnd>=numpnts(w) )
		pEnd= numpnts(w)-1
	endif
	if( pBegin<0 )
		pBegin= 0
	endif
	Variable maxSfact= (pEnd-pBegin+1)/20

	Variable nMaxSF=2* ceil(sqrt(maxSfact))
	Make/O/N=(nMaxSF) wsmdata=0,wsmFact= round((P/2)^2)


	Variable nLin=10
	Variable nSpaced=20
	//print maxSfact
	if( maxSfact<20 )
		nSpaced= 0
		nLin= maxSfact
	endif
	nMaxSF= nLin+nSpaced
	Variable a= (maxSfact-nLin)/nSpaced^2
	Make/O/N=(nMaxSF) wsmdata=0,wsmFact= p+1
	if( nSpaced>0 )
		wsmFact[nLin,*]= ceil(nLin+a*(p-nLin+1)^2)
	endif

	Variable i=1,imax= min(nMaxSF, numpnts(wsmFact))

	wsmdata[0]= snr
	do
		Duplicate/O wtmpOrigDif,wtmp
		Smooth/E=2/B=3 2*wsmFact[i]+1, wtmp
		Histogram/B=1 wtmp,hist
		Integrate hist
		FindLevel/Q hist,0.4*hist[999]
		x0= V_LevelX
		FindLevel/Q hist,0.6*hist[999]
		x1= V_LevelX
		snr= (pnt2x(hist, 999)-pnt2x(hist, 0))/(x1-x0)
		wsmdata[i]= snr
		i+=1
	while(i<imax)
//Duplicate/O wsmdata, root:snrData
//Duplicate/O wsmFact, root:smoothFactors
	WaveStats/Q/R=[2,] wsmdata
	Variable smFact= wsmFact[V_maxloc]

	// added heuristics
//	Variable smFactwpd=0
	do
		if( wsmdata[V_maxloc] < 100 )
			Variable didIt=0
			i=0
			do
				Variable findPeaksReturn = AutoFindPeaks(w,pBegin,pEnd,noiselevel*10,wsmFact[i],1)
//				Variable findPeaksReturnOriginal = AutoFindPeaksOriginal(w,pBegin,pEnd,noiselevel*10,wsmFact[i],1)
//				Variable findPeaksReturnNew = AutoFindPeaksNew(w,pBegin,pEnd,noiselevel*10,wsmFact[i],1)
//print "Smooth Factor:", wsmFact[i], "Original:", findPeaksReturnOriginal, "New:", findPeaksReturnNew
//				if( findPeaksReturnOriginal > 0 )
//				if( findPeaksReturnNew > 0 )
				if( findPeaksReturn > 0 )
					Wave wpd= W_AutoPeakInfo
//					smFactwpd= floor(wpd[0][1]/3)
					smFact= round(wpd[0][1]/3)
					didIt= 1
//	print "TRIAL FIND",wpd[0][1],smFact,i,wsmFact[i],wsmdata[V_maxloc],wsmdata[i]
					break;
				endif
				i+=1
			while(i<imax)
			if( didIt )
				break
			endif
		endif

		// If really low snr and couldn't find a principal peak, force high smooth factors
		if( wsmdata[V_maxloc] < 20 )
			smFact= maxSfact
			break
		endif
		if( wsmdata[V_maxloc] < 30 )
			smFact= round(maxSfact/4)
			break
		endif
		if( wsmdata[V_maxloc] < 50 )
			smFact= round(maxSfact/6)
			break
		endif
	while(0)

	if( smFact < 2 )
		smFact= 2
	endif
	
//	smFact = max(smFact, smFactwpd)

	KillDataFolder :
//SetDataFolder ::
	return cmplx(noiselevel,smFact)
end

Function AutoFindPeaks(w,pBegin,pEnd,noiseEst,smFact,maxPeaks)
	Wave w
	Variable pBegin,pEnd
	Variable noiseEst,smFact
	Variable maxPeaks
	
	NVAR/Z doOriginal = root:V_PeakAutoFindClassic
	if ( (!NVAR_Exists(doOriginal)) || (doOriginal == 0) )
		return AutoFindPeaksNew(w,pBegin,pEnd,noiseEst,smFact,maxPeaks)
	else
		return AutoFindPeaksOriginal(w,pBegin,pEnd,noiseEst,smFact,maxPeaks)
	endif
end

static constant numPeakInfoColumns = 5		// JW 071031
Function AutoFindPeaksOriginal(w,pBegin,pEnd,noiseEst,smFact,maxPeaks)
	Wave w
	Variable pBegin,pEnd
	Variable noiseEst,smFact
	Variable maxPeaks
	
	if( pBegin > pEnd )
		Variable tmp= pBegin
		pBegin= pEnd
		pEnd= tmp
	endif
	
	Make/O/N=(0,numPeakInfoColumns) W_AutoPeakInfo
	
	NewDataFolder/S/O afpTemp2
	Duplicate/O/R=[pBegin,pEnd] w,wtmp1
	SetScale/P x,0,1,wtmp1					// we work in point numbers here

//	Duplicate/O wtmp1,wtmp2
	Smooth/B=3 smFact, wtmp1				// for peak amp determination
	
//Duplicate/O wtmp1, root:smooth1

	Duplicate/O wtmp1,wtmp2
//Duplicate/O wtmp1, root:debugDif

	Differentiate wtmp2
	Smooth/E=2/B=3 2*smFact, wtmp2
//Duplicate/O wtmp2, root:difsmooth2
	Differentiate wtmp2
	Smooth/E=2/B=3 2*smFact, wtmp2
//Duplicate/O wtmp2, root:dif2smooth3
	Duplicate/O wtmp2,wtmp3				// we mung one copy and need an unmunged version also
	
	Variable avgWidth=0					// for width not too far from average width criteria
	
	Variable i=0,peakNum=0,numBadPeaks=0
	do
		WaveStats/Q wtmp2
		Variable x0= V_minloc		// really, need to determine if + or - peaks and use right one

		if( V_min>=0 )
			break
		endif

		FindLevel/Q/R=[x0,] wtmp3,0
		Variable xr= V_LevelX
		if( V_Flag!=0 )
			xr= numpnts(wtmp3)-1
		endif
		FindLevel/Q/R=[x0,0] wtmp3,0		// note search is from right to left
		Variable xl= V_LevelX
		if( V_Flag!=0 )
			xl= 0
		endif
		wtmp2[xl-1,xr+1]=NaN						// don't find this peak again

		Variable widthEst
		Variable rightWidthEst, leftWidthEst, leftWidthFraction	// JW 071031
		do
			if( (x0-xl) < 1 )
				widthEst= xr-x0				// if up against the left edge, use right width
				break
			endif
			if( (xr-x0) < 1 )
				widthEst= x0-xl				// similar for right edge
				break
			endif
			Variable ratio= (xr-x0)/ (x0-xl)	// right width/left width
			if( (ratio < 0.5) || (ratio>2) )
				widthEst=2* min(xr-x0, x0-xl)	// take smaller of widths if one is much larger
				break
			endif
			widthEst= xr-xl
		while(0)

		rightWidthEst = xr-x0						// JW 071031
		leftWidthEst = x0-xl						// JW 071031
		leftWidthFraction = leftWidthEst/widthEst	// JW 071031
		
		if( !(widthEst>3) )						// this probably will neverhappen but if it did, we are probably out of real peaks
			break
		endif

		Variable impulseWidth= 2*(2*smFact+1)
		if( widthEst > 1.3*impulseWidth )
			widthEst= sqrt(widthEst^2 - impulseWidth^2)
		else
			widthEst= widthEst/2
		endif
		leftWidthEst = widthEst*leftWidthFraction
		rightWidthEst = widthEst - leftWidthEst
		
		Variable yl= wtmp1[xl], y0= wtmp1[x0], yr= wtmp1[xr]
		Variable bl0= ((yr-yl)/(xr-xl)*(x0-xl)+yl		// y at x0 for line between left and right inflection points
		
		Variable heightEst= 2*(y0-bl0)
		
		Variable avgNoiseEst= noiseEst/(1.35*sqrt(2*smFact+1))
		Variable minH= avgNoiseEst*8
		
		// throw in an additional penalty if width is far away from the average
		if( avgWidth>0 )
			minH *= sqrt( (widthEst/avgWidth)^2 + (avgWidth/widthEst)^2 )
		endif
		if( heightEst >  minH )
			Redimension/N=(peakNum+1,numPeakInfoColumns) W_AutoPeakInfo
			avgWidth= (avgWidth*peakNum+widthEst)/(peakNum+1)
			W_AutoPeakInfo[peakNum]={{x0+pBegin},{widthEst},{heightEst},{leftWidthEst},{rightWidthEst}}
			peakNum+=1
		else
			if( peakNum == 0 )
				break					// if very first peak is bad, then give up
			endif
			numBadPeaks += 1
			if( numBadPeaks > 3 )
				break
			endif
		endif
	while(peakNum<maxPeaks)
	KillDataFolder :

	return 	peakNum
end

// The AutoFindPeaks functions make a temporary data folder and then kill it. This function depends on that to get rid of some temporary waves.
Function findPeaksIn2ndDeriv(win)
	Wave win
	
	Make/N=(100,7)/O results
	Variable startP = 0
	Variable numPeak=0
	Variable lowBad = 0
	Variable highBad = 0

	findpeak/B=1/N/P/R=[startP,]/Q win
	results[numPeak][0]=V_PeakLoc
	results[numPeak][1]=V_PeakVal
	if (V_flag)
		// couldn't find even one low point???
		Redimension/N=0 results
		return -1
	endif
	findpeak/B=1/P/R=[V_PeakLoc,0]/Q win
	results[numPeak][2]=V_PeakLoc
	results[numPeak][3]=V_PeakVal
	if (V_flag)
		lowBad = 1
	endif
	findpeak/B=1/P/R=[results[numPeak][0],]/Q win
	results[numPeak][4]=V_PeakLoc
	results[numPeak][5]=V_PeakVal
	if (V_flag)
		// Couldn't find the other side of the first peak
		Redimension/N=0 results
		return -1
	endif
	if (lowBad)
		startP = V_PeakLoc
	endif
	
	Variable leftPeakLoc, leftPeakVal
	Variable midPeakLoc, midPeakVal
	Variable rightPeakLoc, rightPeakVal
	do
		findpeak/B=1/N/P/R=[startP,]/Q win
		midPeakLoc=V_PeakLoc
		midPeakVal=V_PeakVal
		if (V_flag)
			break;
		endif
		findpeak/B=1/P/R=[V_PeakLoc,0]/Q win
		leftPeakLoc=V_PeakLoc
		leftPeakVal=V_PeakVal
		if (V_flag)
			break;
		endif
		findpeak/B=1/P/R=[midPeakLoc,]/Q win
		rightPeakLoc=V_PeakLoc
		rightPeakVal=V_PeakVal
		if (V_flag)
			break;
		endif
		
		results[numPeak][0]=midPeakLoc
		results[numPeak][1]=midPeakVal
		results[numPeak][2]=leftPeakLoc
		results[numPeak][3]=leftPeakVal
		results[numPeak][4]=rightPeakLoc
		results[numPeak][5]=rightPeakVal
		results[numPeak][6]=min(abs(midPeakVal - leftPeakVal), abs(midPeakVal - rightPeakVal))
	
		startP=V_PeakLoc
		numPeak += 1
		if (numPeak >= DimSize(results, 0))
			Redimension/N=(numPeak+100, -1) results
		endif
	while(1)
	
	Redimension/N=(numPeak, -1) results
	Make/O/N=(numPeak) sortwave
	sortwave = results[p][6]
	MakeIndex/R sortwave, sortwave
	Duplicate/O results, resultscopy
//Duplicate/O results, root:resultscopy
//Wave resultscopy = root:resultscopy
	results = resultscopy[sortwave[p]][q]
//Duplicate/O results, root:resultsSorted
end

// for debugging
//Function PutLinesOnGraph(win, graph, numpeaks)
//	Wave win
//	String Graph
//	Variable numpeaks
//	
//	Variable numRows = DimSize(win, 0)
//	Variable i
//	
//	SetDrawLayer/W=$Graph/K progBack
//	for (i = 0; i < numpeaks; i += 1)
//		SetDrawEnv/W=$Graph fillfgc=(50000,50000,50000),linethick=1,ycoord=prel,xcoord=bottom
//		DrawRect/W=$Graph win[i][2], 1-i*.01, win[i][4], 0+(numpeaks-i)*.01
//		SetDrawEnv/W=$Graph linefgc=(30000,30000,30000),ycoord=prel,xcoord=bottom
//		DrawLine/W=$Graph win[i][0], 0, win[i][0], 1
//		SetDrawEnv/W=$Graph textxjust=1,ycoord=prel,xcoord=bottom
//		DrawText/W=$Graph win[i][0], 1-i*.01, num2str(i)
//	endfor
//end

Function AutoFindPeaksNew(w,pBegin,pEnd,noiseEst,smFact,maxPeaks)
	Wave w
	Variable pBegin,pEnd
	Variable noiseEst,smFact
	Variable maxPeaks
	
	if( pBegin > pEnd )
		Variable tmp= pBegin
		pBegin= pEnd
		pEnd= tmp
	endif
	
	Make/O/N=(0,numPeakInfoColumns) W_AutoPeakInfo
	
	NewDataFolder/S/O afpTemp2
	Duplicate/O/R=[pBegin,pEnd] w,wtmp1
	SetScale/P x,0,1,wtmp1					// we work in point numbers here

	Smooth/B=3 smFact, wtmp1				// for peak amp determination
	
//Duplicate/O wtmp1, root:smooth1

	Duplicate/O wtmp1,wtmp2
//Duplicate/O wtmp1, root:debugDif

	Differentiate wtmp2
	Smooth/E=2/B=3 2*smFact, wtmp2
//Duplicate/O wtmp2, root:difsmooth2
	Differentiate wtmp2
	Smooth/E=2/B=3 2*smFact, wtmp2
//Duplicate/O wtmp2, root:dif2smooth3
	Duplicate/O wtmp2,wtmp3				// we mung one copy and need an unmunged version also
	
	findPeaksIn2ndDeriv(wtmp2)
//PutLinesOnGraph(results, "Graph7", DimSize(results, 0))
	Wave Results
//print GetWavesDataFolder(results, 2)
	maxPeaks = min(maxPeaks, DimSize(results, 0))
	Variable nRows = DimSize(results, 0)
	
	Variable avgWidth=0					// for width not too far from average width criteria
	
	Variable i=0,peakNum=0,numBadPeaks=0
	for (i = 0; i < nRows; i += 1)
		Variable x0= Results[i][0]
		Variable xr= Results[i][4]
		Variable xl= Results[i][2]

		Variable widthEst
		Variable rightWidthEst, leftWidthEst, leftWidthFraction	// JW 071031
		do
			if( (x0-xl) < 1 )
				widthEst= xr-x0				// if up against the left edge, use right width
				break
			endif
			if( (xr-x0) < 1 )
				widthEst= x0-xl				// similar for right edge
				break
			endif
			Variable ratio= (xr-x0)/ (x0-xl)	// right width/left width
			if( (ratio < 0.5) || (ratio>2) )
				widthEst=2* min(xr-x0, x0-xl)	// take smaller of widths if one is much larger
				break
			endif
			widthEst= xr-xl
		while(0)

		rightWidthEst = xr-x0						// JW 071031
		leftWidthEst = x0-xl						// JW 071031
		leftWidthFraction = leftWidthEst/widthEst	// JW 071031
		
		if( !(widthEst>3) )						// this probably will neverhappen but if it did, we are probably out of real peaks
			break
		endif

		Variable impulseWidth= 2*(2*smFact+1)
		if( widthEst > 1.3*impulseWidth )
			widthEst= sqrt(widthEst^2 - impulseWidth^2)
		else
			widthEst= widthEst/2
		endif
		widthEst /= sqrt(6)
		leftWidthEst = widthEst*leftWidthFraction
		rightWidthEst = widthEst - leftWidthEst
		
		Variable yl= wtmp1[xl], y0= wtmp1[x0], yr= wtmp1[xr]
		Variable bl0 = min(yr, yl)
		
		Variable heightEst= 1.3*(y0-bl0)
		if (heightEst < 0)
			continue
		endif
		
		Variable avgNoiseEst= noiseEst/(1.35*sqrt(2*smFact+1))
		Variable minH= avgNoiseEst*8
		Variable saveMinH = minH
		
		// throw in an additional penalty if width is far away from the average
		if( avgWidth>0 )
			minH *= sqrt( (widthEst/avgWidth)^2 + (avgWidth/widthEst)^2 )
		endif
//print "i=",i,"; point=", x0, "; original minH=",saveMinH,"; minH=",minH,"; heightEst=",heightEst, "; widthEst=", widthEst
		if( heightEst >  minH )
			Redimension/N=(peakNum+1,numPeakInfoColumns) W_AutoPeakInfo
			avgWidth= (avgWidth*peakNum+widthEst)/(peakNum+1)
			W_AutoPeakInfo[peakNum]={{x0+pBegin},{widthEst},{heightEst},{leftWidthEst},{rightWidthEst}}
			peakNum+=1
		else
			if( peakNum == 0 )
				break					// if very first peak is bad, then give up
			endif
			numBadPeaks += 1
			if( numBadPeaks > 3 )
				break
			endif
		endif
		if(peakNum>=maxPeaks)
			break;
		endif
	endfor
	KillDataFolder :

	return 	peakNum
end


Function AdjustAutoPeakInfoForX(wpi,yData,xData)
	Wave wpi,yData
	WAVE/Z xData
	
	Variable imax= DimSize(wpi,0),i=0
	do
		if( WaveExists(xData) )
			Variable p0= wpi[i][0]
			Variable pw= wpi[i][1]/2
			wpi[i][0]=xData[p0]
			wpi[i][1]= abs(xData[p0+pw] - xData[p0-pw])
			Variable pLw = wpi[i][3]
			wpi[i][3]= abs(xData[p0] - xData[p0-pLw])
			Variable pRw = wpi[i][4]
			wpi[i][4]= abs(xData[p0-pRw] - xData[p0])
		else
			wpi[i][0]=pnt2x(yData,wpi[i][0])
			wpi[i][1]=  abs(wpi[i][1]*deltax(yData))
			wpi[i][3]=  abs(wpi[i][3]*deltax(yData))
			wpi[i][4]=  abs(wpi[i][4]*deltax(yData))
		endif
		i+=1
	while(i<imax)
end

Function TrimAmpAutoPeakInfo(wpi,gMinPeakFraction)
	Wave wpi
	Variable gMinPeakFraction
	
	Variable i= DimSize(wpi,0)-1					// index of last row
	Variable ymin= wpi[0][2]*gMinPeakFraction	// user want peaks to be bigger than this
	
	do
		if( wpi[i][2] < ymin )
			DeletePoints i,i,wpi
		endif
		i -= 1
	while(i>0)

	return DimSize(wpi,0)
end