libpappsomspp
Library for mass spectrometry
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 # pragma once
23 
24 //#include <boost/math/special_functions.hpp>
25 #include "QDebug"
26 #include "../resample/convert2dense.hpp"
27 #include "pwiz/utility/findmf/base/filter/filter.hpp"
28 #include "simplepicker.hpp"
29 #include "pwiz/utility/findmf/base/filter/gaussfilter.hpp"
30 #include "pwiz/utility/findmf/base/base/interpolate.hpp"
31 #include "pwiz/utility/findmf/base/resample/determinebinwidth.hpp"
32 #include "pwiz/utility/findmf/base/base/copyif.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  [[maybe_unused]] 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>( std::round( idx - integwith_ ) );
62  size_t end = static_cast<size_t>( std::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>( std::round( idx - integwith_ ) );
100  size_t end = static_cast<size_t>( std::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_;
171  ralab::base::resample::SamplingWith sw_;
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  {
189  c2d_.defBreak(massrange,ralab::base::resample::resolution2ppm(resolution));
191  ralab::base::filter::getGaussianFilterQuantile(filter_,width);
192  }
193 
194 
195  template<typename Tmass, typename Tintensity>
196  void operator()(Tmass begmz, Tmass endmz, Tintensity begint )
197  {
198  // get the first position of the peak begining (intensity > 0.1)
199  typename std::iterator_traits<Tintensity>::value_type minint = *std::upper_bound(begint,begint+std::distance(begmz,endmz),0.1);
200 
201  //determine sampling with
202  double a = sw_(begmz,endmz);
203 
204  //resmpale the spectrum
205  c2d_.am_ = a;
206  c2d_.convert2dense(begmz,endmz, begint, resampledintensity_);
207 
208  //smooth the resampled spectrum
209  ralab::base::filter::filter(resampledintensity_, filter_, smoothedintensity_, true);
210  //determine zero crossings
211  zerocross_.resize( smoothedintensity_.size()/2 );
212  size_t nrzerocross = simplepicker_( smoothedintensity_.begin( ), smoothedintensity_.end(), zerocross_.begin(), zerocross_.size());
213 
214  peakmass_.resize(nrzerocross);
215  //determine mass of zerocrossing
216  ralab::base::base::interpolate_linear( resampledmz_.begin(), resampledmz_.end(),
217  zerocross_.begin(), zerocross_.begin()+nrzerocross,
218  peakmass_.begin());
219 
220  //determine peak area
221  if(area_) {
222  peakarea_.resize(nrzerocross);
223  integrator_( zerocross_.begin(), zerocross_.begin() + nrzerocross,
224  smoothedintensity_.begin(),resampledintensity_.begin(), peakarea_.begin() );
225  } else {
226  //determine intensity
227  peakarea_.resize(nrzerocross);
228  ralab::base::base::interpolate_cubic( smoothedintensity_.begin(), smoothedintensity_.end(),
229  zerocross_.begin(), zerocross_.begin()+nrzerocross,
230  peakarea_.begin());
231  }
232 
233  TReal threshold = static_cast<TReal>(minint) * intensitythreshold_;
234 
235  if(maxnumbersofpeaks_ > 0) {
236  double threshmax = getNToppeaks();
237  if(threshmax > threshold)
238  threshold = threshmax;
239  }
240 
241 
242  if(threshold > 0.01) {
243  filter(threshold);
244  }
245  }
246 
247  /// get min instensity of peak to qualify for max-intensity;
248  TReal getNToppeaks() {
249  TReal intthres = 0.;
250  if(maxnumbersofpeaks_ < peakarea_.size())
251  {
252  std::vector<TReal> tmparea( peakarea_.begin(), peakarea_.end() );
253  std::nth_element(tmparea.begin(),tmparea.end() - maxnumbersofpeaks_, tmparea.end());
254  intthres = *(tmparea.end() - maxnumbersofpeaks_);
255  }
256  return intthres;
257  }
258 
259 
260  /// clean the masses using the threshold
261  void filter(TReal threshold) {
262  typename std::vector<TReal>::iterator a = ralab::base::utils::copy_if(peakarea_.begin(),peakarea_.end(),peakmass_.begin(),
263  peakmass_.begin(),boost::bind(std::greater<TReal>(),_1,threshold));
264  peakmass_.resize(std::distance(peakmass_.begin(),a));
265  typename std::vector<TReal>::iterator b = ralab::base::utils::copy_if(peakarea_.begin(),peakarea_.end(),
266  peakarea_.begin(),boost::bind(std::greater<TReal>(),_1,threshold));
267  peakarea_.resize(std::distance(peakarea_.begin(),b));
268  //int x = 1;
269  }
270 
271  const std::vector<TReal> & getPeakMass() {
272  return peakmass_;
273  }
274 
275  const std::vector<TReal> & getPeakArea() {
276  return peakarea_;
277  }
278 
279  const std::vector<TReal> & getResampledMZ() {
280  return resampledmz_;
281  }
282 
283  const std::vector<TReal> & getResampledIntensity() {
284  return resampledintensity_;
285  }
286 
287  const std::vector<TReal> & getSmoothedIntensity() {
288  return smoothedintensity_;
289  }
290 };
291 }//ms
292 }//base
293 }//ralab
ralab::base::ms::SimplePeakArea::integwith_
TReal integwith_
Definition: peakpickerqtof.hpp:45
ralab::base::ms::SimplePicker
computes first derivative of a sequence, looks for zero crossings
Definition: simplepicker.hpp:39
ralab::base::ms::PeakPicker::smoothwith_
TReal smoothwith_
Definition: peakpickerqtof.hpp:168
ralab::base::ms::LocalMinPeakArea::threshold_
TReal threshold_
Definition: peakpickerqtof.hpp:79
ralab::base::ms::PeakPicker::intensitythreshold_
TReal intensitythreshold_
Definition: peakpickerqtof.hpp:173
ralab::base::ms::PeakPicker::filter_
std::vector< TReal > filter_
Definition: peakpickerqtof.hpp:166
ralab::base::resample::Convert2Dense::defBreak
std::size_t defBreak(std::pair< double, double > &mzrange, double ppm)
computes split points of an map.
Definition: convert2dense.hpp:56
ralab::base::ms::PeakPicker::getResampledIntensity
const std::vector< TReal > & getResampledIntensity()
Definition: peakpickerqtof.hpp:283
ralab::base::ms::PeakPicker::PeakPicker
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)
Definition: peakpickerqtof.hpp:177
ralab::base::ms::PeakPicker::smoothedintensity_
std::vector< TReal > smoothedintensity_
Definition: peakpickerqtof.hpp:166
ralab::base::ms::PeakPicker::operator()
void operator()(Tmass begmz, Tmass endmz, Tintensity begint)
Definition: peakpickerqtof.hpp:196
ralab::base::ms::LocalMinPeakArea::integwith_
TReal integwith_
Definition: peakpickerqtof.hpp:78
ralab::base::ms::PeakPicker::resolution_
TReal resolution_
Definition: peakpickerqtof.hpp:163
ralab::base::ms::PeakPicker::maxnumbersofpeaks_
uint32_t maxnumbersofpeaks_
Definition: peakpickerqtof.hpp:175
ralab::base::ms::PeakPicker::getPeakArea
const std::vector< TReal > & getPeakArea()
Definition: peakpickerqtof.hpp:275
ralab::base::ms::LocalMinPeakArea::LocalMinPeakArea
LocalMinPeakArea(TReal integwith, TReal threshold=.1)
Definition: peakpickerqtof.hpp:81
ralab::base::resample::Convert2Dense::convert2dense
void convert2dense(Tmass beginMass, Tmass endMass, Tintens intens, Tout ass)
Converts a sparse spec to a dense spec.
Definition: convert2dense.hpp:69
ralab::base::ms::SimplePeakArea
Definition: peakpickerqtof.hpp:44
ralab::base::ms::PeakPicker::peakarea_
std::vector< TReal > peakarea_
Definition: peakpickerqtof.hpp:167
ralab::base::ms::PeakPicker::resampledmz_
std::vector< TReal > resampledmz_
Definition: peakpickerqtof.hpp:165
ralab::base::ms::PeakPicker::getSmoothedIntensity
const std::vector< TReal > & getSmoothedIntensity()
Definition: peakpickerqtof.hpp:287
ralab::base::ms::PeakPicker::c2d_
ralab::base::resample::Convert2Dense c2d_
Definition: peakpickerqtof.hpp:164
ralab::base::ms::LocalMinPeakArea::mextend
void mextend(TInt &start, TInt &end, TInt idx) const
exend peak to left and rigth
Definition: peakpickerqtof.hpp:119
ralab::base::ms::SimplePeakArea::SimplePeakArea
SimplePeakArea(TReal integwith)
Definition: peakpickerqtof.hpp:47
ralab::base::ms::PeakPicker::filter
void filter(TReal threshold)
clean the masses using the threshold
Definition: peakpickerqtof.hpp:261
ralab::base::ms::LocalMinPeakArea::operator()
void operator()(Tzerocross beginZ, Tzerocross endZ, Tintensity intensity, Tintensity resampled, Tout area) const
intagrates the peak intesnities
Definition: peakpickerqtof.hpp:89
ralab::base::ms::PeakPicker::value_type
TReal value_type
Definition: peakpickerqtof.hpp:160
ralab::base::ms::PeakPicker
Definition: peakpickerqtof.hpp:159
ralab::base::ms::PeakPicker::getPeakMass
const std::vector< TReal > & getPeakMass()
Definition: peakpickerqtof.hpp:271
ralab::base::ms::PeakPicker::getNToppeaks
TReal getNToppeaks()
get min instensity of peak to qualify for max-intensity;
Definition: peakpickerqtof.hpp:248
ralab::base::ms::PeakPicker::PeakIntegrator
TIntegrator< value_type > PeakIntegrator
Definition: peakpickerqtof.hpp:161
ralab::base::ms::PeakPicker::integrator_
PeakIntegrator integrator_
Definition: peakpickerqtof.hpp:172
ralab::base::ms::SimplePeakArea::operator()
void operator()(Tzerocross beginZ, Tzerocross endZ, [[maybe_unused]] Tintensity intensity, Tintensity resmpled, Tout area) const
intagrates the peak intesnities
Definition: peakpickerqtof.hpp:51
ralab::base::ms::LocalMinPeakArea
Definition: peakpickerqtof.hpp:76
ralab::base::ms::PeakPicker::getResampledMZ
const std::vector< TReal > & getResampledMZ()
Definition: peakpickerqtof.hpp:279
ralab::base::resample::Convert2Dense::getMids
void getMids(std::vector< double > &mids)
Definition: convert2dense.hpp:122
ralab::base::ms::PeakPicker::peakmass_
std::vector< TReal > peakmass_
Definition: peakpickerqtof.hpp:167
ralab::base::ms::LocalMinPeakArea::value_type
TReal value_type
Definition: peakpickerqtof.hpp:77
ralab::base::ms::PeakPicker::area_
bool area_
Definition: peakpickerqtof.hpp:174
simplepicker.hpp
ralab::base::resample::Convert2Dense
Definition: convert2dense.hpp:45
ralab::base::ms::PeakPicker::simplepicker_
ralab::base::ms::SimplePicker< TReal > simplepicker_
Definition: peakpickerqtof.hpp:170
ralab
Definition: peakpickerqtof.hpp:34
ralab::base::ms::PeakPicker::integrationWidth_
TReal integrationWidth_
Definition: peakpickerqtof.hpp:169
ralab::base::ms::PeakPicker::zerocross_
std::vector< TReal > zerocross_
Definition: peakpickerqtof.hpp:166
ralab::base::resample::Convert2Dense::am_
double am_
Definition: convert2dense.hpp:49
ralab::base::ms::PeakPicker::sw_
ralab::base::resample::SamplingWith sw_
Definition: peakpickerqtof.hpp:171
ralab::base::ms::PeakPicker::resampledintensity_
std::vector< TReal > resampledintensity_
Definition: peakpickerqtof.hpp:165