ProteoWizard
peakpickerqtof.hpp
Go to the documentation of this file.
1 //
2 // $Id: peakpickerqtof.hpp 6337 2014-06-06 20:04:10Z witek96 $
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 
34 namespace 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>
45  TReal integwith_;
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;
78  TReal integwith_;
79  TReal threshold_;
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 >
159  struct PeakPicker{
160  typedef TReal value_type;
161  typedef TIntegrator<value_type> PeakIntegrator;
162 
163  TReal resolution_;
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
168  TReal smoothwith_;
172  PeakIntegrator integrator_;
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),
186  integrationWidth_(intwidth),sw_(),integrator_(integrationWidth_),
187  intensitythreshold_(intensitythreshold),area_(area),maxnumbersofpeaks_(maxnumberofpeaks)
188  {
189  c2d_.defBreak(massrange,ralab::base::resample::resolution2ppm(resolution));
190  c2d_.getMids(resampledmz_);
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
207  ralab::base::filter::filter(resampledintensity_ , filter_ , smoothedintensity_ , true);
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
214  ralab::base::base::interpolate_linear( resampledmz_.begin() , resampledmz_.end() ,
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);
226  ralab::base::base::interpolate_cubic( smoothedintensity_.begin() , smoothedintensity_.end() ,
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;
246  TReal getNToppeaks(){
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
TIntegrator< value_type > PeakIntegrator
resamples spectrum, apply smoothing, determines zero crossings, integrates peaks. ...
void operator()(Tmass begmz, Tmass endmz, Tintensity begint)
boost::uint32_t uint32_t
Definition: filter.hpp:47
const std::vector< TReal > & getPeakArea()
const std::vector< TReal > & getSmoothedIntensity()
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&#39;s.
TReal getNToppeaks()
get min instensity of peak to qualify for max-intensity;
const std::vector< TReal > & getPeakMass()
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)
ralab::base::resample::SamplingWith sw_
LocalMinPeakArea(TReal integwith, TReal threshold=.1)
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) ...
Definition: gaussfilter.hpp:59
void convert2dense(Tmass beginMass, Tmass endMass, Tintens intens, Tout ass)
Converts a sparse spec to a dense spec.
resamples spectrum, apply smoothing, determines zero crossings, integrates peaks. ...
const std::vector< TReal > & getResampledIntensity()
void operator()(Tzerocross beginZ, Tzerocross endZ, Tintensity intensity, Tintensity resampled, Tout area) const
intagrates the peak intesnities
ralab::base::resample::Convert2Dense c2d_
void filter(TReal threshold)
clean the masses using the threshold
std::vector< TReal > peakmass_
std::vector< TReal > zerocross_
std::vector< TReal > resampledmz_
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.
EQUISPACEINTERPOL Interpolation on a equidistantly spaced grid.
Definition: base.hpp:39
void getMids(std::vector< double > &mids)
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
double resolution2ppm(double resolution)
OutputIterator copy_if(InputIterator first, InputIterator last, InputIterator2 source, OutputIterator result, Predicate pred)
Definition: copyif.hpp:35
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.
Definition: interpolate.hpp:58
const std::vector< TReal > & getResampledMZ()
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
ralab::base::ms::SimplePicker< TReal > simplepicker_