PeakFinder

C++

Public Domain

Finds peaks and peak widths in sampled function data

Download (right click, save as, rename as appropriate)

Embed

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
#ifndef PEAK_FINDER_1D_H
#define PEAK_FINDER_1D_H

#define _USE_MATH_DEFINES
#include <functional>
#include <list>
#include <vector>
#include <cmath>
#include <float.h>

// Version 1.0 - 2009-05-08 - Ripped out of FunctionSampler. Reworked it
//                            to avoid copying into a vector. The code
//                            depends heavily on list iterators being able
//                            to go backwards.

// Uncomment the following line to print out (lots!) of debugging
// output that shows the subdivision process.
// #define DEBUG_PEAK_FINDER

#ifdef DEBUG_PEAK_FINDER
# include <iostream>
#endif

namespace PeakFinder1D{

template <class ResultType>
struct Peak{
	double x;
	ResultType value;
	
	// true if the peak is a maximum, false for a minimum
	bool is_local_max;
	
	// The full width at half maximum (fwhm) is defined as the x-extent away
	// from the peak until the function reaches the threshold value.
	// The threshold value for a local minimum is halfway from the peak's y
	// value to the maximum y value of any sample. For local maxima, it is
	// halfway from the peak's y value to the minimum function value.
	// The function does not have to be monotonically varying from the peak
	// to the threshold.
	// If the fwhm is not defined, this means the extent would reach beyond
	// the endpoints of the sample data.
	bool well_defined_fwhm;
	double fwhm;
	
	// The Q is defined only when the function value varies monotonically
	// from the peak to the fwhm threshold value on either side.
	// The Q, when it is defined, is fwhm/x.
	bool well_defined_Q;
	double Q;
	
	Peak(ResultType P):value(P){}
};

template <class ResultType>
void FindPeaks(
	const std::list<std::pair<double, ResultType> > &data,
	std::list<Peak<ResultType> > &peaks
){
	typedef typename std::list<std::pair<double, ResultType> > input_list_t;
	typedef typename input_list_t::const_iterator list_iter_t;
	
	peaks.erase(peaks.begin(), peaks.end());
	if(data.size() < 3){ return; }
	
	// First pass, find maximum and minimum values
	double MinY = DBL_MAX, MaxY = -DBL_MAX;
	for(list_iter_t i = data.begin(); i != data.end(); ++i){
		double yi = static_cast<double>(i->second);
		if(yi < MinY){ MinY = yi; }
		if(yi > MaxY){ MaxY = yi; }
	}
#ifdef DEBUG_PEAK_FINDER
	std::cerr << "* MinY = " << MinY
	          << ", MaxY = " << MaxY << std::endl;
#endif
	
	// Find x-coords of peaks
	// Scan and record last point that was less than current value
	// We maintain pointers to three points in the data such that
	// the values at these points are not equal
	list_iter_t ii = data.begin();
	list_iter_t ji = ii; ++ji;
	list_iter_t ki = ji; ++ki;

	double iival = static_cast<double>(ii->second);
	double jival = static_cast<double>(ji->second);
	// Move ji forward until value at ji != value at ii
	while(iival == jival && ki != data.end()){
		++ji;
		jival = static_cast<double>(ji->second);
		++ki;
	}
	if(ki == data.end()){ return; }
#ifdef DEBUG_PEAK_FINDER
	std::cerr << "* Starting xi = " << ii->first
	                   << ", xj = " << ji->first
	                   << ", xk = " << ki->first << std::endl;
#endif
	do{
		// Move ki forward until value at ki != value at ji
		double kival = static_cast<double>(ki->second);
		while(jival == kival && ki != data.end()){
			++ki;
			kival = static_cast<double>(ki->second);
		}
		if(ki == data.end()){ return; }
		// At this point, val(ii) != val(i+1) and val(ji) != val(ki)
		// ji >= ii+1, val(ii+1) = ... = val(ji), ki == ji+1
		
#ifdef DEBUG_PEAK_FINDER
		std::cerr << "* Current xi = " << ii->first
						  << ", xj = " << ji->first
						  << ", xk = " << ki->first << std::endl;
		std::cerr << "*         yi = " << iival
						  << ", yj = " << jival
						  << ", yk = " << kival << std::endl;
#endif

		// A local extremum is when the middle point (at ji) is less than
		// or greater than both the other two values
		bool is_max = ((iival < jival) && (jival > kival));
		bool is_min = ((iival > jival) && (jival < kival));
		if(is_max || is_min){
			Peak<ResultType> lpeak(ji->second);
			lpeak.is_local_max = is_max;

#ifdef DEBUG_PEAK_FINDER
			std::cerr << "*  Found local ";
			if(is_max){
				std::cerr << "max";
			}else{
				std::cerr << "min";
			}
			std::cerr << std::endl;
#endif

			list_iter_t ip1i = ii; ++ip1i;
			if(ip1i == ji){
				lpeak.x = ji->first;
			}else{
				lpeak.x = 0.5*(ip1i->first + ji->first);
			}
			
			// Crawl away from the extremum on either side
			// until we reach the threshold.
			// The threshold is defined to be halfway from the extremum
			// value to the global min or max value (depending on if this
			// peak is a local max or min).
			double threshold;
			if(is_max){
				threshold = MinY + 0.5*(jival - MinY);
			}else{
				threshold = MaxY - 0.5*(MaxY - jival);
			}
			
			// Maintain pointers iibound, and kibound which are used
			// to crawl outward (iibound crawls left, kibound goes right)
			list_iter_t iibound = ii, kibound = ki;
			// Also maintain points ip1ibound, km1ibound:
			list_iter_t ip1ibound; // should always be iibound+1
			list_iter_t km1ibound; // should always be kibound-1
			
			// We want to see if we reach the threshold. These are flags
			// indicating if we reached it while the function was still
			// monotone.
			bool ibound_found = false, kbound_found = false;
			bool monotonic_until_ibound = true;
			bool monotonic_until_kbound = true;
			// Keep the value at iibound and ip1ibound around for use
			double iiboundval = iival;
			ip1ibound = iibound; ++ip1ibound;
			double ip1iboundval = static_cast<double>(ip1ibound->second);
			while(iibound != data.begin()){
				if( (is_max && (iiboundval < threshold)) ||
				    (is_min && (iiboundval > threshold)) ){
					ibound_found = true;
					break;
				}
				ip1ibound = iibound;
				ip1iboundval = iiboundval;
				--iibound;
				if(data.begin() == iibound){ break; }
				iiboundval = static_cast<double>(iibound->second);
				
				if( (is_max && (iiboundval > ip1iboundval)) ||
				    (is_min && (iiboundval < ip1iboundval)) ){
					monotonic_until_ibound = false;
				}
			}
			// Do the same thing for kibound, etc.
			double kiboundval = kival;
			km1ibound = ki; --km1ibound;
			double km1iboundval = static_cast<double>(km1ibound->second);
			while(kibound != data.end()){
				if( (is_max && (kiboundval < threshold)) ||
				    (is_min && (kiboundval > threshold)) ){
					kbound_found = true;
					break;
				}
				km1ibound = kibound;
				km1iboundval = kiboundval;
				++kibound;
				if(data.end() == kibound){ break; }
				kiboundval = static_cast<double>(kibound->second);
				
				if( (is_max && (kiboundval > km1iboundval)) ||
				    (is_min && (kiboundval < km1iboundval)) ){
					monotonic_until_kbound = false;
				}
			}
			// Now determine if the peak has computable properties
			lpeak.well_defined_fwhm = false;
			lpeak.well_defined_Q = false;
			if(ibound_found && kbound_found){
				lpeak.well_defined_fwhm = true;
				
				// Find the treshold x/y values using interpolation
				double xi, xk, ti, tk;
				ti = (threshold - iiboundval) / (ip1iboundval-iiboundval);
				xi = iibound->first + ti*(ip1ibound->first-iibound->first);
				tk = (threshold - km1iboundval) / (kiboundval-km1iboundval);
				xk = km1ibound->first + tk*(ip1ibound->first-iibound->first);
				lpeak.fwhm = xk - xi;
				
				if(monotonic_until_ibound && monotonic_until_kbound){
					lpeak.well_defined_Q = true;
					lpeak.Q = lpeak.x / (xk - xi);
				}
			}
			peaks.push_back(lpeak);
		}
		// Move all the pointers forwards
		ii = ji; iival = jival;
		ji = ki; jival = kival;
		++ki;
	}while(ki != data.end());
}

}; // end namespace PeakFinder1D

#endif // PEAK_FINDER_1D_H