libpappsomspp
Library for mass spectrometry
Loading...
Searching...
No Matches
pappso::LocalSignalToNoiseEstimator Class Reference

#include <localsignaltonoiseestimator.h>

Classes

struct  Parameters
struct  GaussianEstimateParams

Public Types

enum  IntensityThresholdCalculation { MANUAL = 0x0 , AUTOMAXBYSTDEV , AUTOMAXBYPERCENTILE }

Public Member Functions

 LocalSignalToNoiseEstimator (const Parameters &parameters)
virtual ~LocalSignalToNoiseEstimator ()
GaussianEstimateParams estimateGaussian (const TraceIterator &iter_first_data_point, const TraceIterator &iter_last_data_point) const
void initialize (const Trace &trace)
void computeSignaToNoiseRatio (const Trace &trace)
double getSignalToNoiseRatio (qsizetype index) const

Protected Attributes

Parameters m_parameters
double m_sparseWindowPercentage
double m_histogramOverflowPercentage
QList< double > m_signalToNoiseEstimates

Detailed Description

Definition at line 25 of file localsignaltonoiseestimator.h.

Member Enumeration Documentation

◆ IntensityThresholdCalculation

enum pappso::LocalSignalToNoiseEstimator::IntensityThresholdCalculation
Enumerator
MANUAL 
AUTOMAXBYSTDEV 
AUTOMAXBYPERCENTILE 

Definition at line 30 of file localsignaltonoiseestimator.h.

31 {
32 MANUAL = 0x0,
33 AUTOMAXBYSTDEV,
34 AUTOMAXBYPERCENTILE
35 };

Constructor & Destructor Documentation

◆ LocalSignalToNoiseEstimator()

pappso::LocalSignalToNoiseEstimator::LocalSignalToNoiseEstimator ( const Parameters & parameters)
explicit

Definition at line 21 of file localsignaltonoiseestimator.cpp.

22 : m_parameters(parameters)
23{
24}

References m_parameters.

◆ ~LocalSignalToNoiseEstimator()

pappso::LocalSignalToNoiseEstimator::~LocalSignalToNoiseEstimator ( )
virtual

Definition at line 26 of file localsignaltonoiseestimator.cpp.

27{
28}

Member Function Documentation

◆ computeSignaToNoiseRatio()

void pappso::LocalSignalToNoiseEstimator::computeSignaToNoiseRatio ( const Trace & trace)

Definition at line 73 of file localsignaltonoiseestimator.cpp.

74{
75 // First element in the Trace
76 TraceIterator iter_first_data_point = trace.begin();
77 // Last element in the scan
78 TraceIterator iter_last_data_point = trace.end();
79
80 // Reset counter for sparse windows
81 m_sparseWindowPercentage = 0;
82 // Reset counter for histogram overflow
83 m_histogramOverflowPercentage = 0;
84
85 // Reset the results
86 m_signalToNoiseEstimates.clear();
87 m_signalToNoiseEstimates.resize(trace.size());
88
89 // Maximal range of histogram needs to be calculated first
90 if(m_parameters.maxIntMode == AUTOMAXBYSTDEV)
91 {
92 // Use MEAN+auto_m_parameters.maxIntensity*STDEV as threshold
93 GaussianEstimateParams gaussian_estimate_global =
94 estimateGaussian(iter_first_data_point, iter_last_data_point);
95 m_parameters.maxIntensity =
96 gaussian_estimate_global.mean +
97 std::sqrt(gaussian_estimate_global.variance) * m_parameters.maxIntStdDevFactor;
98
99 // qDebug() << "Computed max intensity:" << m_parameters.maxIntensity;
100 }
101 else if(m_parameters.maxIntMode == AUTOMAXBYPERCENTILE)
102 {
103 // get value at "m_parameters.maxIntPercentile"th percentile
104 // we use a histogram approach here as well.
105 if((m_parameters.maxIntPercentile < 0) || (m_parameters.maxIntPercentile > 100))
106 {
107 qFatal() << "The m_parameters.maxIntPercentile value is not in the [0, 100] range:"
108 << m_parameters.maxIntPercentile;
109 }
110
111 // Histogram with 100 percent values all initialized to 0.
112 QList<int> histogram_of_intensity_distribution(100, 0);
113
114 // Find maximum of current scan
115 auto max_intensity_data_point_iterator = std::max_element(
116 trace.begin(), trace.end(), [](const DataPoint &a, const DataPoint &b) {
117 return a.y > b.y;
118 });
119
120 double max_intensity = max_intensity_data_point_iterator->y;
121
122 // qDebug() << "Max intensity value in the trace:" << max_intensity;
123
124 double histogram_bin_size = max_intensity / 100;
125
126 // Fill histogram
127 for(const DataPoint &data_point : trace)
128 {
129 // Increment count at histogram bin that is the bin corresponding
130 // to the iterated data point intensity in the percentage scale above.
131 ++histogram_of_intensity_distribution[(int)((data_point.y - 1) / histogram_bin_size)];
132 }
133
134 // Add up element counts in histogram until ?th percentile is reached
135 int elements_below_max_int_percentile =
136 (int)(m_parameters.maxIntPercentile * trace.size() / 100);
137
138 int elements_seen = 0;
139 int i = -1;
140
141 TraceIterator the_iterator = iter_first_data_point;
142
143 while(the_iterator != iter_last_data_point &&
144 elements_seen < elements_below_max_int_percentile)
145 {
146 ++i;
147 elements_seen += histogram_of_intensity_distribution[i];
148 ++the_iterator;
149 }
150
151 // The max intensity is the median of the histogram.
152 m_parameters.maxIntensity = (((double)i) + 0.5) * histogram_bin_size;
153
154 // qDebug() << "Computed max intensity:" << m_parameters.maxIntensity;
155 }
156 else // if (m_parameters.maxIntMode == MANUAL)
157 {
158 if(m_parameters.maxIntensity <= 0)
159 {
160 qFatal() << "The method to determine the maximal intensity is set to MANUAL"
161 << "However, the m_parameters.maxIntensity value is <= 0.";
162 }
163
164 // qDebug() << "The max intensity was set manually:" << m_parameters.maxIntensity;
165 }
166
167 if(m_parameters.maxIntensity < 0)
168 {
169 qFatal() << "The m_parameters.maxIntensity "
170 "value should be positive! ";
171 }
172
173 TraceIterator window_pos_center = iter_first_data_point;
174 TraceIterator window_pos_border_left = iter_first_data_point;
175 TraceIterator window_pos_border_right = iter_first_data_point;
176
177 double window_half_size = m_parameters.windowSize / 2;
178 // At least size of 1 for intensity bins
179 double bin_size = std::max(1.0, m_parameters.maxIntensity / m_parameters.binCount);
180 int bin_count_minus_one = m_parameters.binCount - 1;
181
182 std::vector<int> histogram(m_parameters.binCount, 0);
183 std::vector<double> bin_value(m_parameters.binCount, 0);
184 // calculate average intensity that is represented by a bin
185 for(int bin = 0; bin < m_parameters.binCount; bin++)
186 {
187 histogram[bin] = 0;
188 bin_value[bin] = (bin + 0.5) * bin_size;
189 }
190 // bin in which a datapoint would fall
191 int to_bin = 0;
192
193 // index of bin where the median is located
194 int median_bin = 0;
195 // additive number of elements from left to x in histogram
196 int element_inc_count = 0;
197
198 // tracks elements in current window, which may vary because of unevenly
199 // spaced data
200 int elements_in_window = 0;
201 // number of windows
202 int window_count = 0;
203
204 // number of elements where we find the median
205 int element_in_window_half = 0;
206
207 double noise; // noise value of a datapoint
208
209 while(window_pos_center != iter_last_data_point)
210 {
211 // Erase all elements from histogram that will leave the window on the
212 // LEFT side
213 while((*window_pos_border_left).x <
214 (*window_pos_center).x - window_half_size)
215 {
216 to_bin = std::max(
217 std::min<int>((int)((*window_pos_border_left).y / bin_size),
218 bin_count_minus_one),
219 0);
220 --histogram[to_bin];
221 --elements_in_window;
222 ++window_pos_border_left;
223 }
224
225 // Add all elements to histogram that will enter the window on the RIGHT
226 // side
227 while((window_pos_border_right != iter_last_data_point) &&
228 ((*window_pos_border_right).x <=
229 (*window_pos_center).x + window_half_size))
230 {
231 // std::cerr << (*window_pos_border_right).y << " " <<
232 // bin_size << " " << m_parameters.binCountminus_1 << std::endl;
233 to_bin = std::max(
234 std::min<int>((int)((*window_pos_border_right).y / bin_size),
235 bin_count_minus_one),
236 0);
237 ++histogram[to_bin];
238 ++elements_in_window;
239 ++window_pos_border_right;
240 }
241
242 if(elements_in_window < m_parameters.minRequiredElements)
243 {
244 noise = m_parameters.noiseValueForEmptyWindow;
245 ++m_sparseWindowPercentage;
246 }
247 else
248 {
249 // Find bin i where ceil[elements_in_window/2] <= sum_c(0..i){
250 // histogram[c] }
251 median_bin = -1;
252 element_inc_count = 0;
253 element_in_window_half = (elements_in_window + 1) / 2;
254 while(median_bin < bin_count_minus_one &&
255 element_inc_count < element_in_window_half)
256 {
257 ++median_bin;
258 element_inc_count += histogram[median_bin];
259 }
260
261 // increase the error count
262 if(median_bin == bin_count_minus_one)
263 {
264 ++m_histogramOverflowPercentage;
265 }
266
267 // just avoid division by 0
268 noise = std::max(1.0, bin_value[median_bin]);
269 }
270
271 // Store result
272 m_signalToNoiseEstimates[window_count] = (*window_pos_center).y / noise;
273 // qDebug() << "Added new signal to noise estimate:"
274 // << m_signalToNoiseEstimates.at(window_count);
275
276 // advance the window center by one datapoint
277 ++window_pos_center;
278 ++window_count;
279
280 } // end while
281
282 // qDebug() << "Done filling in the signal to noise estimates, the are"
283 // << m_signalToNoiseEstimates.size() << "estimates.";
284
285 m_sparseWindowPercentage = m_sparseWindowPercentage * 100 / window_count;
286 m_histogramOverflowPercentage =
287 m_histogramOverflowPercentage * 100 / window_count;
288
289 // warn if percentage of sparse windows is above 20%
290 if(m_sparseWindowPercentage > 20)
291 {
292 qWarning() << m_sparseWindowPercentage
293 << "% of all windows were sparse. You should consider increasing"
294 << "'m_parameters.windowSize' or decrease 'm_parameters.minRequiredElements'";
295 }
296
297 // warn if percentage of possibly wrong median estimates is above 1%
298 if(m_histogramOverflowPercentage > 1)
299 {
300 qWarning() << m_histogramOverflowPercentage
301 << "% of all Signal-to-Noise estimates are too high, because the "
302 "median was found in the rightmost histogram-bin. "
303 << "You should consider increasing 'm_parameters.maxIntensityy' (and maybe "
304 "'m_parameters.binCount' with it, to keep bin width reasonable)";
305 }
306}
pappso::LocalSignalToNoiseEstimator::estimateGaussian
GaussianEstimateParams estimateGaussian(const TraceIterator &iter_first_data_point, const TraceIterator &iter_last_data_point) const
Definition localsignaltonoiseestimator.cpp:32
pappso::TraceIterator
std::vector< DataPoint >::const_iterator TraceIterator
Definition localsignaltonoiseestimator.h:23

References pappso::a, AUTOMAXBYPERCENTILE, AUTOMAXBYSTDEV, pappso::b, estimateGaussian(), m_histogramOverflowPercentage, m_parameters, m_signalToNoiseEstimates, m_sparseWindowPercentage, pappso::LocalSignalToNoiseEstimator::GaussianEstimateParams::mean, and pappso::LocalSignalToNoiseEstimator::GaussianEstimateParams::variance.

Referenced by initialize().

◆ estimateGaussian()

LocalSignalToNoiseEstimator::GaussianEstimateParams pappso::LocalSignalToNoiseEstimator::estimateGaussian ( const TraceIterator & iter_first_data_point,
const TraceIterator & iter_last_data_point ) const

Definition at line 32 of file localsignaltonoiseestimator.cpp.

35{
36 int size = 0;
37 // add up
38 double v = 0;
39 double m = 0;
40 TraceIterator run = iter_first_data_point;
41 while(run != iter_last_data_point)
42 {
43 m += (*run).y;
44 ++size;
45 ++run;
46 }
47 // average
48 m = m / size;
49
50 // determine variance
51 run = iter_first_data_point;
52 while(run != iter_last_data_point)
53 {
54 double tmp(m - (*run).y);
55 v += tmp * tmp;
56 ++run;
57 }
58 v = v / ((double)size); // divide by n
59
60 GaussianEstimateParams estimate_params = {m, v};
61 return estimate_params;
62}

Referenced by computeSignaToNoiseRatio().

◆ getSignalToNoiseRatio()

double pappso::LocalSignalToNoiseEstimator::getSignalToNoiseRatio ( qsizetype index) const

Definition at line 312 of file localsignaltonoiseestimator.cpp.

313{
314 Q_ASSERT(index < m_signalToNoiseEstimates.size());
315 return m_signalToNoiseEstimates[index];
316}

References m_signalToNoiseEstimates.

Referenced by pappso::HighResPeakPicker::pick().

◆ initialize()

void pappso::LocalSignalToNoiseEstimator::initialize ( const Trace & trace)

Definition at line 65 of file localsignaltonoiseestimator.cpp.

66{
67 computeSignaToNoiseRatio(trace);
68}

References computeSignaToNoiseRatio().

Referenced by pappso::HighResPeakPicker::pick().

Member Data Documentation

◆ m_histogramOverflowPercentage

double pappso::LocalSignalToNoiseEstimator::m_histogramOverflowPercentage
protected

Definition at line 128 of file localsignaltonoiseestimator.h.

Referenced by computeSignaToNoiseRatio().

◆ m_parameters

Parameters pappso::LocalSignalToNoiseEstimator::m_parameters
protected

Definition at line 123 of file localsignaltonoiseestimator.h.

Referenced by LocalSignalToNoiseEstimator(), and computeSignaToNoiseRatio().

◆ m_signalToNoiseEstimates

QList<double> pappso::LocalSignalToNoiseEstimator::m_signalToNoiseEstimates
protected

Definition at line 130 of file localsignaltonoiseestimator.h.

Referenced by computeSignaToNoiseRatio(), and getSignalToNoiseRatio().

◆ m_sparseWindowPercentage

double pappso::LocalSignalToNoiseEstimator::m_sparseWindowPercentage
protected

Definition at line 126 of file localsignaltonoiseestimator.h.

Referenced by computeSignaToNoiseRatio().

The documentation for this class was generated from the following files: