34 #ifndef OPENMS_TRANSFORMATIONS_RAW2PEAK_PEAKPICKERHIRES_H
35 #define OPENMS_TRANSFORMATIONS_RAW2PEAK_PEAKPICKERHIRES_H
44 #include <gsl/gsl_spline.h>
45 #include <gsl/gsl_interp.h>
50 #define DEBUG_PEAK_PICKING
51 #undef DEBUG_PEAK_PICKING
82 template <
typename PeakType>
87 output.SpectrumSettings::operator=(input);
88 output.MetaInfoInterface::operator=(input);
95 if (input.size() < 5)
return;
100 if (signal_to_noise_ > 0.0)
106 for (
Size i = 2; i < input.size() - 2; ++i)
108 double central_peak_mz = input[i].getMZ(), central_peak_int = input[i].getIntensity();
109 double left_neighbor_mz = input[i - 1].getMZ(), left_neighbor_int = input[i - 1].getIntensity();
110 double right_neighbor_mz = input[i + 1].getMZ(), right_neighbor_int = input[i + 1].getIntensity();
113 double left_to_central = std::fabs(central_peak_mz - left_neighbor_mz);
114 double central_to_right = std::fabs(right_neighbor_mz - central_peak_mz);
115 double min_spacing = (left_to_central < central_to_right) ? left_to_central : central_to_right;
117 double act_snt = 0.0, act_snt_l1 = 0.0, act_snt_r1 = 0.0;
119 if (signal_to_noise_ > 0.0)
127 if (act_snt >= signal_to_noise_
128 && left_to_central < 1.5 * min_spacing
129 && central_peak_int > left_neighbor_int
130 && act_snt_l1 >= signal_to_noise_
131 && central_to_right < 1.5 * min_spacing
132 && central_peak_int > right_neighbor_int
133 && act_snt_r1 >= signal_to_noise_)
139 double act_snt_l2 = 0.0, act_snt_r2 = 0.0;
141 if (signal_to_noise_ > 0.0)
148 && std::fabs(left_neighbor_mz - input[i - 2].getMZ()) < 1.5 * min_spacing
149 && left_neighbor_int < input[i - 2].getIntensity()
150 && act_snt_l2 >= signal_to_noise_)
152 ((i + 2) < input.size()
153 && std::fabs(input[i + 2].getMZ() - right_neighbor_mz) < 1.5 * min_spacing
154 && right_neighbor_int < input[i + 2].getIntensity()
155 && act_snt_r2 >= signal_to_noise_)
163 std::map<double, double> peak_raw_data;
165 peak_raw_data[central_peak_mz] = central_peak_int;
166 peak_raw_data[left_neighbor_mz] = left_neighbor_int;
167 peak_raw_data[right_neighbor_mz] = right_neighbor_int;
174 Size missing_left(0);
175 Size missing_right(0);
177 while ((i - k + 1) > 0
178 && (missing_left < 2)
179 && input[i - k].getIntensity() <= peak_raw_data.begin()->second)
182 double act_snt_lk = 0.0;
184 if (signal_to_noise_ > 0.0)
190 if (act_snt_lk >= signal_to_noise_ && std::fabs(input[i - k].getMZ() - peak_raw_data.begin()->first) < 1.5 * min_spacing)
192 peak_raw_data[input[i -
k].getMZ()] = input[i -
k].getIntensity();
196 peak_raw_data[input[i -
k].getMZ()] = input[i -
k].getIntensity();
206 while ((i + k) < input.size()
207 && (missing_right < 2)
208 && input[i + k].getIntensity() <= peak_raw_data.rbegin()->second)
211 double act_snt_rk = 0.0;
213 if (signal_to_noise_ > 0.0)
218 if (act_snt_rk >= signal_to_noise_ && std::fabs(input[i + k].getMZ() - peak_raw_data.rbegin()->first) < 1.5 * min_spacing)
220 peak_raw_data[input[i +
k].getMZ()] = input[i +
k].getIntensity();
224 peak_raw_data[input[i +
k].getMZ()] = input[i +
k].getIntensity();
243 const Size num_raw_points = peak_raw_data.size();
245 std::vector<double> raw_mz_values;
246 std::vector<double> raw_int_values;
248 for (std::map<double, double>::const_iterator map_it = peak_raw_data.begin(); map_it != peak_raw_data.end(); ++map_it)
250 raw_mz_values.push_back(map_it->first);
251 raw_int_values.push_back(map_it->second);
255 gsl_interp_accel * spline_acc = gsl_interp_accel_alloc();
256 gsl_interp_accel * first_deriv_acc = gsl_interp_accel_alloc();
257 gsl_spline * peak_spline = gsl_spline_alloc(gsl_interp_cspline, num_raw_points);
258 gsl_spline_init(peak_spline, &(*raw_mz_values.begin()), &(*raw_int_values.begin()), num_raw_points);
263 double max_peak_mz = central_peak_mz, max_peak_int = central_peak_int;
264 double threshold = 0.000001;
265 double lefthand = left_neighbor_mz;
266 double righthand = right_neighbor_mz;
268 bool lefthand_sign = 1;
269 double eps = std::numeric_limits<double>::epsilon();
275 double mid = (lefthand + righthand) / 2;
277 double midpoint_deriv_val = gsl_spline_eval_deriv(peak_spline, mid, first_deriv_acc);
280 if (!(std::fabs(midpoint_deriv_val) > eps))
285 bool midpoint_sign = (midpoint_deriv_val < 0.0) ? 0 : 1;
287 if (lefthand_sign ^ midpoint_sign)
303 while (std::fabs(lefthand - righthand) > threshold);
306 max_peak_mz = (lefthand + righthand) / 2;
307 max_peak_int = gsl_spline_eval(peak_spline, max_peak_mz, spline_acc);
311 peak.
setMZ(max_peak_mz);
313 output.push_back(peak);
316 gsl_spline_free(peak_spline);
317 gsl_interp_accel_free(spline_acc);
318 gsl_interp_accel_free(first_deriv_acc);
328 template <
typename PeakType>
333 output.ChromatogramSettings::operator=(input);
334 output.MetaInfoInterface::operator=(input);
341 input_spectrum.push_back(*it);
343 pick(input_spectrum, output_spectrum);
346 output.push_back(*it);
354 template <
typename PeakType,
typename ChromatogramPeakT>
366 bool ms1_only = param_.getValue(
"ms1_only").toBool();
370 for (
Size scan_idx = 0; scan_idx != input.
size(); ++scan_idx)
372 if (ms1_only && (input[scan_idx].getMSLevel() != 1))
374 output[scan_idx] = input[scan_idx];
378 pick(input[scan_idx], output[scan_idx]);
380 setProgress(++progress);
387 setProgress(++progress);
400 void updateMembers_();
virtual double getSignalToNoise(const PeakIterator &data_point)
Definition: SignalToNoiseEstimator.h:130
const String & getName() const
Definition: MSChromatogram.h:194
A 2-dimensional raw data point or peak.
Definition: Peak2D.h:55
UInt getMSLevel() const
Returns the MS level.
Definition: MSSpectrum.h:231
Size size() const
Definition: MSExperiment.h:117
void setMZ(CoordinateType coordinate)
Mutable access to the m/z coordinate (index 1)
Definition: Peak2D.h:197
Peak data (also called centroided data or stick data)
Definition: SpectrumSettings.h:74
The representation of a chromatogram.
Definition: MSChromatogram.h:53
const String & getName() const
Returns the name.
Definition: MSSpectrum.h:243
void setName(const String &name)
Sets the name.
Definition: MSSpectrum.h:249
void resize(Size s)
Definition: MSExperiment.h:122
void setName(const String &name)
Sets the name.
Definition: MSChromatogram.h:200
void addChromatogram(const MSChromatogram< ChromatogramPeakType > &chromatogram)
adds a chromatogram to the list
Definition: MSExperiment.h:762
void setIntensity(IntensityType intensity)
Non-mutable access to the data point intensity (height)
Definition: Peak2D.h:167
virtual void init(const PeakIterator &it_begin, const PeakIterator &it_end)
Set the start and endpoint of the raw data interval, for which signal to noise ratios will be estimat...
Definition: SignalToNoiseEstimator.h:111
void pick(const MSSpectrum< PeakType > &input, MSSpectrum< PeakType > &output) const
Applies the peak-picking algorithm to a single spectrum (MSSpectrum). The resulting picked peaks are ...
Definition: PeakPickerHiRes.h:83
void setMSLevel(UInt ms_level)
Sets the MS level.
Definition: MSSpectrum.h:237
void clear(bool clear_meta_data)
Clears all data and meta data.
Definition: MSSpectrum.h:605
Representation of a mass spectrometry experiment.
Definition: MSExperiment.h:68
double signal_to_noise_
Definition: PeakPickerHiRes.h:397
void setType(SpectrumType type)
sets the spectrum type
void setRT(DoubleReal rt)
Sets the absolute retention time (is seconds)
Definition: MSSpectrum.h:221
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
void clear(bool clear_meta_data)
Clears all data and meta data.
Definition: MSExperiment.h:809
Base class for all classes that want to report their progess.
Definition: ProgressLogger.h:56
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:90
DoubleReal getRT() const
Definition: MSSpectrum.h:215
void pickExperiment(const MSExperiment< PeakType, ChromatogramPeakT > &input, MSExperiment< PeakType, ChromatogramPeakT > &output) const
Applies the peak-picking algorithm to a map (MSExperiment). This method picks peaks for each scan in ...
Definition: PeakPickerHiRes.h:355
void clear(bool clear_meta_data)
Clears all data and meta data.
Definition: MSChromatogram.h:564
This class implements a fast peak-picking algorithm best suited for high resolution MS data (FT-ICR-M...
Definition: PeakPickerHiRes.h:68
Description of the experimental settings.
Definition: ExperimentalSettings.h:59
const std::vector< MSChromatogram< ChromatogramPeakType > > & getChromatograms() const
returns the chromatogram list
Definition: MSExperiment.h:768
void pick(const MSChromatogram< PeakType > &input, MSChromatogram< PeakType > &output) const
Definition: PeakPickerHiRes.h:329