ProteoWizard
peakpickerqtof.hpp
Go to the documentation of this file.
1//
2// $Id$
3//
4//
5// Original author: Witold Wolski <wewolski@gmail.com>
6//
7// Copyright : ETH Zurich
8//
9// Licensed under the Apache License, Version 2.0 (the "License");
10// you may not use this file except in compliance with the License.
11// You may obtain a copy of the License at
12//
13// http://www.apache.org/licenses/LICENSE-2.0
14//
15// Unless required by applicable law or agreed to in writing, software
16// distributed under the License is distributed on an "AS IS" BASIS,
17// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18// See the License for the specific language governing permissions and
19// limitations under the License.
20//
21
22#ifndef PEAKPICKER_H
23#define PEAKPICKER_H
24
25#include <boost/math/special_functions.hpp>
33
34namespace ralab{
35 namespace base{
36 namespace ms{
37
38
39 /// resamples spectrum, apply smoothing,
40 /// determines zero crossings,
41 /// integrates peaks.
42
43 template<typename TReal>
46
47 SimplePeakArea(TReal integwith):integwith_(integwith){}
48
49 /// intagrates the peak intesnities
50 template<typename Tzerocross, typename Tintensity, typename Tout>
51 void operator()( Tzerocross beginZ,
52 Tzerocross endZ,
53 Tintensity intensity,
54 Tintensity resmpled,
55 Tout area)const
56 {
57 typedef typename std::iterator_traits<Tout>::value_type AreaType;
58 for( ; beginZ != endZ ; ++beginZ , ++area )
59 {
60 size_t idx = static_cast<size_t>( *beginZ );
61 size_t start = static_cast<size_t>( boost::math::round( idx - integwith_ ) );
62 size_t end = static_cast<size_t>( boost::math::round( idx + integwith_ + 2.) );
63 AreaType aread = 0.;
64 for( ; start != end ; ++start )
65 {
66 aread += *(resmpled + start);
67 }
68 *area = aread;
69 }
70 }
71 };
72
73 /// extends peak to the left and to the right to the next local minimum or a predefined threshol
74 /// or a maximum allowed extension.
75 template<typename TReal>
77 typedef TReal value_type;
80
81 LocalMinPeakArea(TReal integwith,//!<maximal allowed peak width +- in pixel
82 TReal threshold = .1// minimum intensity
83 ):integwith_(integwith),threshold_(threshold){}
84
85
86
87 /// intagrates the peak intesnities
88 template< typename Tzerocross, typename Tintensity, typename Tout >
89 void operator()( Tzerocross beginZ,
90 Tzerocross endZ,
91 Tintensity intensity,
92 Tintensity resampled,
93 Tout area) const
94 {
95 typedef typename std::iterator_traits<Tout>::value_type AreaType;
96 for( ; beginZ != endZ ; ++beginZ , ++area )
97 {
98 size_t idx = static_cast<size_t>( *beginZ );
99 size_t start = static_cast<size_t>( boost::math::round( idx - integwith_ ) );
100 size_t end = static_cast<size_t>( boost::math::round( idx + integwith_ + 2) );
101
102 Tintensity st = intensity + start;
103 Tintensity en = intensity + end;
104 Tintensity center = intensity + idx;
105 std::ptrdiff_t x1 = std::distance(st, center);
106 std::ptrdiff_t y1 = std::distance(center,en);
107 mextend(st , en , center);
108 std::ptrdiff_t x2 = std::distance(intensity,st);
109 std::ptrdiff_t y2 = std::distance(intensity,en);
110 std::ptrdiff_t pp = std::distance(st,en);
111 AreaType areav = std::accumulate(resampled+x2,resampled+y2,0.);
112 *area = areav;
113 }
114 }
115
116 private:
117 ///exend peak to left and rigth
118 template<typename TInt >
119 void mextend( TInt &start, TInt &end, TInt idx) const
120 {
121 typedef typename std::iterator_traits<TInt>::value_type Intensitytype;
122 //
123 for(TInt intens = idx ; intens >= start; --intens){
124 Intensitytype val1 = *intens;
125 Intensitytype val2 = *(intens-1);
126 if(val1 > threshold_){
127 if(val1 < val2 ){
128 start = intens;
129 break;
130 }
131 }
132 else{
133 start = intens;
134 break;
135 }
136 }
137
138 for(TInt intens = idx ; intens <= end; ++intens){
139 Intensitytype val1 = *intens;
140 Intensitytype val2 = *(intens+1);
141 if(val1 > threshold_){
142 if(val1 < val2 ){
143 end = intens;
144 break;
145 }
146 }
147 else{
148 end = intens;
149 break;
150 }
151 }
152 }
153 };
154
155 /// resamples spectrum, apply smoothing,
156 /// determines zero crossings,
157 /// integrates peaks.
158 template<typename TReal, template <typename B> class TIntegrator >
160 typedef TReal value_type;
161 typedef TIntegrator<value_type> PeakIntegrator;
162
165 std::vector<TReal> resampledmz_, resampledintensity_; // keeps result of convert to dense
166 std::vector<TReal> filter_, zerocross_, smoothedintensity_; // working variables
167 std::vector<TReal> peakmass_, peakarea_; //results
174 bool area_;
176
177 PeakPicker(TReal resolution, //!< instrument resolution
178 std::pair<TReal, TReal> & massrange, //!< mass range of spectrum
179 TReal width = 2., //!< smooth width
180 TReal intwidth = 2., //!< integration width used for area compuation
181 TReal intensitythreshold = 10., // intensity threshold
182 bool area = true,//!< compute area or height? default - height.
183 uint32_t maxnumberofpeaks = 0, //!< maximum of peaks returned by picker
184 double c2d = 1e-5 //!< instrument resampling with small default dissables automatic determination
185 ): resolution_(resolution),c2d_( c2d ) ,smoothwith_(width),
187 intensitythreshold_(intensitythreshold),area_(area),maxnumbersofpeaks_(maxnumberofpeaks)
188 {
192 }
193
194
195 template<typename Tmass, typename Tintensity>
196 void operator()(Tmass begmz, Tmass endmz, Tintensity begint )
197 {
198 typename std::iterator_traits<Tintensity>::value_type minint = *std::upper_bound(begint,begint+std::distance(begmz,endmz),0.1);
199
200 //determine sampling with
201 double a = sw_(begmz,endmz);
202 //resmpale the spectrum
203 c2d_.am_ = a;
204 c2d_.convert2dense(begmz,endmz, begint, resampledintensity_);
205
206 //smooth the resampled spectrum
208 //determine zero crossings
209 zerocross_.resize( smoothedintensity_.size()/2 );
210 size_t nrzerocross = simplepicker_( smoothedintensity_.begin( ) , smoothedintensity_.end() , zerocross_.begin(), zerocross_.size());
211
212 peakmass_.resize(nrzerocross);
213 //determine mass of zerocrossing
215 zerocross_.begin(), zerocross_.begin()+nrzerocross ,
216 peakmass_.begin());
217
218 //determine peak area
219 if(area_){
220 peakarea_.resize(nrzerocross);
221 integrator_( zerocross_.begin(), zerocross_.begin() + nrzerocross ,
222 smoothedintensity_.begin(),resampledintensity_.begin(), peakarea_.begin() );
223 }else{
224 //determine intensity
225 peakarea_.resize(nrzerocross);
227 zerocross_.begin(), zerocross_.begin()+nrzerocross ,
228 peakarea_.begin());
229 }
230
231 TReal threshold = static_cast<TReal>(minint) * intensitythreshold_;
232
233 if(maxnumbersofpeaks_ > 0){
234 double threshmax = getNToppeaks();
235 if(threshmax > threshold)
236 threshold = threshmax;
237 }
238
239
240 if(threshold > 0.01){
241 filter(threshold);
242 }
243 }
244
245 /// get min instensity of peak to qualify for max-intensity;
247 TReal intthres = 0.;
248 if(maxnumbersofpeaks_ < peakarea_.size())
249 {
250 std::vector<TReal> tmparea( peakarea_.begin() , peakarea_.end() );
251 std::nth_element(tmparea.begin(),tmparea.end() - maxnumbersofpeaks_ , tmparea.end());
252 intthres = *(tmparea.end() - maxnumbersofpeaks_);
253 }
254 return intthres;
255 }
256
257
258 /// clean the masses using the threshold
259 void filter(TReal threshold){
260 typename std::vector<TReal>::iterator a = ralab::base::utils::copy_if(peakarea_.begin(),peakarea_.end(),peakmass_.begin(),
261 peakmass_.begin(),boost::bind(std::greater<TReal>(),_1,threshold));
262 peakmass_.resize(std::distance(peakmass_.begin(),a));
263 typename std::vector<TReal>::iterator b = ralab::base::utils::copy_if(peakarea_.begin(),peakarea_.end(),
264 peakarea_.begin(),boost::bind(std::greater<TReal>(),_1,threshold));
265 peakarea_.resize(std::distance(peakarea_.begin(),b));
266 //int x = 1;
267 }
268
269 const std::vector<TReal> & getPeakMass(){
270 return peakmass_;
271 }
272
273 const std::vector<TReal> & getPeakArea(){
274 return peakarea_;
275 }
276
277 const std::vector<TReal> & getResampledMZ(){
278 return resampledmz_;
279 }
280
281 const std::vector<TReal> & getResampledIntensity(){
282 return resampledintensity_;
283 }
284
285 const std::vector<TReal> & getSmoothedIntensity(){
286 return smoothedintensity_;
287 }
288 };
289 }//ms
290 }//base
291}//ralab
292
293
294
295#endif // PEAKPICKER_H
void interpolate_cubic(YInputIterator begY, YInputIterator endY, XInputIterator begX, XInputIterator endX, OutputIterator out, int start_index=0, typename std::iterator_traits< OutputIterator >::value_type epsilon=std::numeric_limits< typename std::iterator_traits< OutputIterator >::value_type >::epsilon())
cubic interpolation on equidistantly spaced y's.
void interpolate_linear(YInputIterator begY, YInputIterator endY, XInputIterator begX, XInputIterator endX, OutputIterator out, int start_index=0, typename std::iterator_traits< OutputIterator >::value_type epsilon=std::numeric_limits< typename std::iterator_traits< OutputIterator >::value_type >::epsilon())
affine interpolation on equidistantly spaced y.
void filter(const TContainer &data, const TContainer &filter, TContainer &result, bool circular=false, uint32_t sides=2)
Applies linear convolution (filtering) to a univariate time series.
Definition filter.hpp:112
TReal getGaussianFilterQuantile(std::vector< TReal > &gauss, TReal fwhm=20, TReal quantile=0.01)
generate the gauss filter function for filtering of peaks with fwhm (full width at half max)
double resolution2ppm(double resolution)
OutputIterator copy_if(InputIterator first, InputIterator last, InputIterator2 source, OutputIterator result, Predicate pred)
Definition copyif.hpp:35
EQUISPACEINTERPOL Interpolation on a equidistantly spaced grid.
Definition base.hpp:40
extends peak to the left and to the right to the next local minimum or a predefined threshol or a max...
void mextend(TInt &start, TInt &end, TInt idx) const
exend peak to left and rigth
void operator()(Tzerocross beginZ, Tzerocross endZ, Tintensity intensity, Tintensity resampled, Tout area) const
intagrates the peak intesnities
LocalMinPeakArea(TReal integwith, TReal threshold=.1)
resamples spectrum, apply smoothing, determines zero crossings, integrates peaks.
std::vector< TReal > resampledmz_
std::vector< TReal > smoothedintensity_
std::vector< TReal > filter_
TIntegrator< value_type > PeakIntegrator
TReal getNToppeaks()
get min instensity of peak to qualify for max-intensity;
ralab::base::resample::Convert2Dense c2d_
const std::vector< TReal > & getSmoothedIntensity()
std::vector< TReal > peakarea_
const std::vector< TReal > & getResampledMZ()
ralab::base::ms::SimplePicker< TReal > simplepicker_
void filter(TReal threshold)
clean the masses using the threshold
PeakPicker(TReal resolution, std::pair< TReal, TReal > &massrange, TReal width=2., TReal intwidth=2., TReal intensitythreshold=10., bool area=true, uint32_t maxnumberofpeaks=0, double c2d=1e-5)
std::vector< TReal > resampledintensity_
ralab::base::resample::SamplingWith sw_
void operator()(Tmass begmz, Tmass endmz, Tintensity begint)
const std::vector< TReal > & getPeakMass()
const std::vector< TReal > & getPeakArea()
std::vector< TReal > peakmass_
const std::vector< TReal > & getResampledIntensity()
std::vector< TReal > zerocross_
resamples spectrum, apply smoothing, determines zero crossings, integrates peaks.
void operator()(Tzerocross beginZ, Tzerocross endZ, Tintensity intensity, Tintensity resmpled, Tout area) const
intagrates the peak intesnities
std::size_t defBreak(std::pair< double, double > &mzrange, double ppm)
computes split points of an map.
void getMids(std::vector< double > &mids)
void convert2dense(Tmass beginMass, Tmass endMass, Tintens intens, Tout ass)
Converts a sparse spec to a dense spec.