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 "../resample/convert2dense.hpp"
26 #include "pwiz/utility/findmf/base/filter/filter.hpp"
27 #include "simplepicker.hpp"
28 #include "pwiz/utility/findmf/base/filter/gaussfilter.hpp"
29 #include "pwiz/utility/findmf/base/base/interpolate.hpp"
30 #include "pwiz/utility/findmf/base/resample/determinebinwidth.hpp"
31 #include "pwiz/utility/findmf/base/base/copyif.hpp"
32 
33 namespace ralab {
34 namespace base {
35 namespace ms {
36 
37 
38 /// resamples spectrum, apply smoothing,
39 /// determines zero crossings,
40 /// integrates peaks.
41 
42 template<typename TReal>
44  TReal integwith_;
45 
46  SimplePeakArea(TReal integwith):integwith_(integwith) {}
47 
48  /// intagrates the peak intesnities
49  template<typename Tzerocross, typename Tintensity, typename Tout>
50  void operator()( Tzerocross beginZ,
51  Tzerocross endZ,
52  [[maybe_unused]] Tintensity intensity,
53  Tintensity resmpled,
54  Tout area)const
55  {
56  typedef typename std::iterator_traits<Tout>::value_type AreaType;
57  for( ; beginZ != endZ ; ++beginZ, ++area )
58  {
59  size_t idx = static_cast<size_t>( *beginZ );
60  size_t start = static_cast<size_t>( std::round( idx - integwith_ ) );
61  size_t end = static_cast<size_t>( std::round( idx + integwith_ + 2.) );
62  AreaType aread = 0.;
63  for( ; start != end ; ++start )
64  {
65  aread += *(resmpled + start);
66  }
67  *area = aread;
68  }
69  }
70 };
71 
72 /// extends peak to the left and to the right to the next local minimum or a predefined threshol
73 /// or a maximum allowed extension.
74 template<typename TReal>
76  typedef TReal value_type;
77  TReal integwith_;
78  TReal threshold_;
79 
80  LocalMinPeakArea(TReal integwith,//!<maximal allowed peak width +- in pixel
81  TReal threshold = .1// minimum intensity
82  ):integwith_(integwith),threshold_(threshold) {}
83 
84 
85 
86  /// intagrates the peak intesnities
87  template< typename Tzerocross, typename Tintensity, typename Tout >
88  void operator()( Tzerocross beginZ,
89  Tzerocross endZ,
90  Tintensity intensity,
91  Tintensity resampled,
92  Tout area) const
93  {
94  typedef typename std::iterator_traits<Tout>::value_type AreaType;
95  for( ; beginZ != endZ ; ++beginZ, ++area )
96  {
97  size_t idx = static_cast<size_t>( *beginZ );
98  size_t start = static_cast<size_t>( std::round( idx - integwith_ ) );
99  size_t end = static_cast<size_t>( std::round( idx + integwith_ + 2) );
100 
101  Tintensity st = intensity + start;
102  Tintensity en = intensity + end;
103  Tintensity center = intensity + idx;
104  std::ptrdiff_t x1 = std::distance(st, center);
105  std::ptrdiff_t y1 = std::distance(center,en);
106  mextend(st, en, center);
107  std::ptrdiff_t x2 = std::distance(intensity,st);
108  std::ptrdiff_t y2 = std::distance(intensity,en);
109  std::ptrdiff_t pp = std::distance(st,en);
110  AreaType areav = std::accumulate(resampled+x2,resampled+y2,0.);
111  *area = areav;
112  }
113  }
114 
115 private:
116  ///exend peak to left and rigth
117  template<typename TInt >
118  void mextend( TInt &start, TInt &end, TInt idx) const
119  {
120  typedef typename std::iterator_traits<TInt>::value_type Intensitytype;
121  //
122  for(TInt intens = idx ; intens >= start; --intens) {
123  Intensitytype val1 = *intens;
124  Intensitytype val2 = *(intens-1);
125  if(val1 > threshold_) {
126  if(val1 < val2 ) {
127  start = intens;
128  break;
129  }
130  }
131  else {
132  start = intens;
133  break;
134  }
135  }
136 
137  for(TInt intens = idx ; intens <= end; ++intens) {
138  Intensitytype val1 = *intens;
139  Intensitytype val2 = *(intens+1);
140  if(val1 > threshold_) {
141  if(val1 < val2 ) {
142  end = intens;
143  break;
144  }
145  }
146  else {
147  end = intens;
148  break;
149  }
150  }
151  }
152 };
153 
154 /// resamples spectrum, apply smoothing,
155 /// determines zero crossings,
156 /// integrates peaks.
157 template<typename TReal, template <typename B> class TIntegrator >
158 struct PeakPicker {
159  typedef TReal value_type;
160  typedef TIntegrator<value_type> PeakIntegrator;
161 
162  TReal resolution_;
164  std::vector<TReal> resampledmz_, resampledintensity_; // keeps result of convert to dense
165  std::vector<TReal> filter_, zerocross_, smoothedintensity_; // working variables
166  std::vector<TReal> peakmass_, peakarea_; //results
167  TReal smoothwith_;
170  ralab::base::resample::SamplingWith sw_;
173  bool area_;
175 
176  PeakPicker(TReal resolution, //!< instrument resolution
177  std::pair<TReal, TReal> & massrange, //!< mass range of spectrum
178  TReal width = 2., //!< smooth width
179  TReal intwidth = 2., //!< integration width used for area compuation
180  TReal intensitythreshold = 10., // intensity threshold
181  bool area = true,//!< compute area or height? default - height.
182  uint32_t maxnumberofpeaks = 0, //!< maximum of peaks returned by picker
183  double c2d = 1e-5 //!< instrument resampling with small default dissables automatic determination
184  ): resolution_(resolution),c2d_( c2d ),smoothwith_(width),
186  intensitythreshold_(intensitythreshold),area_(area),maxnumbersofpeaks_(maxnumberofpeaks)
187  {
188  c2d_.defBreak(massrange,ralab::base::resample::resolution2ppm(resolution));
190  ralab::base::filter::getGaussianFilterQuantile(filter_,width);
191  }
192 
193 
194  template<typename Tmass, typename Tintensity>
195  void operator()(Tmass begmz, Tmass endmz, Tintensity begint )
196  {
197  typename std::iterator_traits<Tintensity>::value_type minint = *std::upper_bound(begint,begint+std::distance(begmz,endmz),0.1);
198 
199  //determine sampling with
200  double a = sw_(begmz,endmz);
201  //resmpale the spectrum
202  c2d_.am_ = a;
203  c2d_.convert2dense(begmz,endmz, begint, resampledintensity_);
204 
205  //smooth the resampled spectrum
206  ralab::base::filter::filter(resampledintensity_, filter_, smoothedintensity_, true);
207  //determine zero crossings
208  zerocross_.resize( smoothedintensity_.size()/2 );
209  size_t nrzerocross = simplepicker_( smoothedintensity_.begin( ), smoothedintensity_.end(), zerocross_.begin(), zerocross_.size());
210 
211  peakmass_.resize(nrzerocross);
212  //determine mass of zerocrossing
213  ralab::base::base::interpolate_linear( resampledmz_.begin(), resampledmz_.end(),
214  zerocross_.begin(), zerocross_.begin()+nrzerocross,
215  peakmass_.begin());
216 
217  //determine peak area
218  if(area_) {
219  peakarea_.resize(nrzerocross);
220  integrator_( zerocross_.begin(), zerocross_.begin() + nrzerocross,
221  smoothedintensity_.begin(),resampledintensity_.begin(), peakarea_.begin() );
222  } else {
223  //determine intensity
224  peakarea_.resize(nrzerocross);
225  ralab::base::base::interpolate_cubic( smoothedintensity_.begin(), smoothedintensity_.end(),
226  zerocross_.begin(), zerocross_.begin()+nrzerocross,
227  peakarea_.begin());
228  }
229 
230  TReal threshold = static_cast<TReal>(minint) * intensitythreshold_;
231 
232  if(maxnumbersofpeaks_ > 0) {
233  double threshmax = getNToppeaks();
234  if(threshmax > threshold)
235  threshold = threshmax;
236  }
237 
238 
239  if(threshold > 0.01) {
240  filter(threshold);
241  }
242  }
243 
244  /// get min instensity of peak to qualify for max-intensity;
245  TReal getNToppeaks() {
246  TReal intthres = 0.;
247  if(maxnumbersofpeaks_ < peakarea_.size())
248  {
249  std::vector<TReal> tmparea( peakarea_.begin(), peakarea_.end() );
250  std::nth_element(tmparea.begin(),tmparea.end() - maxnumbersofpeaks_, tmparea.end());
251  intthres = *(tmparea.end() - maxnumbersofpeaks_);
252  }
253  return intthres;
254  }
255 
256 
257  /// clean the masses using the threshold
258  void filter(TReal threshold) {
259  typename std::vector<TReal>::iterator a = ralab::base::utils::copy_if(peakarea_.begin(),peakarea_.end(),peakmass_.begin(),
260  peakmass_.begin(),boost::bind(std::greater<TReal>(),_1,threshold));
261  peakmass_.resize(std::distance(peakmass_.begin(),a));
262  typename std::vector<TReal>::iterator b = ralab::base::utils::copy_if(peakarea_.begin(),peakarea_.end(),
263  peakarea_.begin(),boost::bind(std::greater<TReal>(),_1,threshold));
264  peakarea_.resize(std::distance(peakarea_.begin(),b));
265  //int x = 1;
266  }
267 
268  const std::vector<TReal> & getPeakMass() {
269  return peakmass_;
270  }
271 
272  const std::vector<TReal> & getPeakArea() {
273  return peakarea_;
274  }
275 
276  const std::vector<TReal> & getResampledMZ() {
277  return resampledmz_;
278  }
279 
280  const std::vector<TReal> & getResampledIntensity() {
281  return resampledintensity_;
282  }
283 
284  const std::vector<TReal> & getSmoothedIntensity() {
285  return smoothedintensity_;
286  }
287 };
288 }//ms
289 }//base
290 }//ralab
ralab::base::ms::SimplePeakArea::integwith_
TReal integwith_
Definition: peakpickerqtof.hpp:44
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:167
ralab::base::ms::LocalMinPeakArea::threshold_
TReal threshold_
Definition: peakpickerqtof.hpp:78
ralab::base::ms::PeakPicker::intensitythreshold_
TReal intensitythreshold_
Definition: peakpickerqtof.hpp:172
ralab::base::ms::PeakPicker::filter_
std::vector< TReal > filter_
Definition: peakpickerqtof.hpp:165
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:280
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:176
ralab::base::ms::PeakPicker::smoothedintensity_
std::vector< TReal > smoothedintensity_
Definition: peakpickerqtof.hpp:165
ralab::base::ms::PeakPicker::operator()
void operator()(Tmass begmz, Tmass endmz, Tintensity begint)
Definition: peakpickerqtof.hpp:195
ralab::base::ms::LocalMinPeakArea::integwith_
TReal integwith_
Definition: peakpickerqtof.hpp:77
ralab::base::ms::PeakPicker::resolution_
TReal resolution_
Definition: peakpickerqtof.hpp:162
ralab::base::ms::PeakPicker::maxnumbersofpeaks_
uint32_t maxnumbersofpeaks_
Definition: peakpickerqtof.hpp:174
ralab::base::ms::PeakPicker::getPeakArea
const std::vector< TReal > & getPeakArea()
Definition: peakpickerqtof.hpp:272
ralab::base::ms::LocalMinPeakArea::LocalMinPeakArea
LocalMinPeakArea(TReal integwith, TReal threshold=.1)
Definition: peakpickerqtof.hpp:80
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:68
ralab::base::ms::SimplePeakArea
Definition: peakpickerqtof.hpp:43
ralab::base::ms::PeakPicker::peakarea_
std::vector< TReal > peakarea_
Definition: peakpickerqtof.hpp:166
ralab::base::ms::PeakPicker::resampledmz_
std::vector< TReal > resampledmz_
Definition: peakpickerqtof.hpp:164
ralab::base::ms::PeakPicker::getSmoothedIntensity
const std::vector< TReal > & getSmoothedIntensity()
Definition: peakpickerqtof.hpp:284
ralab::base::ms::PeakPicker::c2d_
ralab::base::resample::Convert2Dense c2d_
Definition: peakpickerqtof.hpp:163
ralab::base::ms::LocalMinPeakArea::mextend
void mextend(TInt &start, TInt &end, TInt idx) const
exend peak to left and rigth
Definition: peakpickerqtof.hpp:118
ralab::base::ms::SimplePeakArea::SimplePeakArea
SimplePeakArea(TReal integwith)
Definition: peakpickerqtof.hpp:46
ralab::base::ms::PeakPicker::filter
void filter(TReal threshold)
clean the masses using the threshold
Definition: peakpickerqtof.hpp:258
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:88
ralab::base::ms::PeakPicker::value_type
TReal value_type
Definition: peakpickerqtof.hpp:159
ralab::base::ms::PeakPicker
Definition: peakpickerqtof.hpp:158
ralab::base::ms::PeakPicker::getPeakMass
const std::vector< TReal > & getPeakMass()
Definition: peakpickerqtof.hpp:268
ralab::base::ms::PeakPicker::getNToppeaks
TReal getNToppeaks()
get min instensity of peak to qualify for max-intensity;
Definition: peakpickerqtof.hpp:245
ralab::base::ms::PeakPicker::PeakIntegrator
TIntegrator< value_type > PeakIntegrator
Definition: peakpickerqtof.hpp:160
ralab::base::ms::PeakPicker::integrator_
PeakIntegrator integrator_
Definition: peakpickerqtof.hpp:171
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:50
ralab::base::ms::LocalMinPeakArea
Definition: peakpickerqtof.hpp:75
ralab::base::ms::PeakPicker::getResampledMZ
const std::vector< TReal > & getResampledMZ()
Definition: peakpickerqtof.hpp:276
ralab::base::resample::Convert2Dense::getMids
void getMids(std::vector< double > &mids)
Definition: convert2dense.hpp:120
ralab::base::ms::PeakPicker::peakmass_
std::vector< TReal > peakmass_
Definition: peakpickerqtof.hpp:166
ralab::base::ms::LocalMinPeakArea::value_type
TReal value_type
Definition: peakpickerqtof.hpp:76
ralab::base::ms::PeakPicker::area_
bool area_
Definition: peakpickerqtof.hpp:173
simplepicker.hpp
ralab::base::resample::Convert2Dense
Definition: convert2dense.hpp:44
ralab::base::ms::PeakPicker::simplepicker_
ralab::base::ms::SimplePicker< TReal > simplepicker_
Definition: peakpickerqtof.hpp:169
ralab
Definition: peakpickerqtof.hpp:33
ralab::base::ms::PeakPicker::integrationWidth_
TReal integrationWidth_
Definition: peakpickerqtof.hpp:168
ralab::base::ms::PeakPicker::zerocross_
std::vector< TReal > zerocross_
Definition: peakpickerqtof.hpp:165
ralab::base::resample::Convert2Dense::am_
double am_
Definition: convert2dense.hpp:49
ralab::base::ms::PeakPicker::sw_
ralab::base::resample::SamplingWith sw_
Definition: peakpickerqtof.hpp:170
ralab::base::ms::PeakPicker::resampledintensity_
std::vector< TReal > resampledintensity_
Definition: peakpickerqtof.hpp:164