25 #include <boost/math/special_functions.hpp> 43 template<
typename TReal>
50 template<
typename Tzerocross,
typename T
intensity,
typename Tout>
57 typedef typename std::iterator_traits<Tout>::value_type AreaType;
58 for( ; beginZ != endZ ; ++beginZ , ++area )
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.) );
64 for( ; start != end ; ++start )
66 aread += *(resmpled + start);
75 template<
typename TReal>
83 ):integwith_(integwith),threshold_(threshold){}
88 template<
typename Tzerocross,
typename T
intensity,
typename Tout >
95 typedef typename std::iterator_traits<Tout>::value_type AreaType;
96 for( ; beginZ != endZ ; ++beginZ , ++area )
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) );
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.);
118 template<
typename TInt >
119 void mextend( TInt &start, TInt &end, TInt idx)
const 121 typedef typename std::iterator_traits<TInt>::value_type Intensitytype;
123 for(TInt intens = idx ; intens >= start; --intens){
124 Intensitytype val1 = *intens;
125 Intensitytype val2 = *(intens-1);
126 if(val1 > threshold_){
138 for(TInt intens = idx ; intens <= end; ++intens){
139 Intensitytype val1 = *intens;
140 Intensitytype val2 = *(intens+1);
141 if(val1 > threshold_){
158 template<
typename TReal,
template <
typename B>
class TIntegrator >
178 std::pair<TReal, TReal> & massrange,
181 TReal intensitythreshold = 10.,
185 ): resolution_(resolution),c2d_( c2d ) ,smoothwith_(width),
186 integrationWidth_(intwidth),sw_(),integrator_(integrationWidth_),
187 intensitythreshold_(intensitythreshold),area_(area),maxnumbersofpeaks_(maxnumberofpeaks)
195 template<
typename Tmass,
typename T
intensity>
196 void operator()(Tmass begmz, Tmass endmz, Tintensity begint )
198 typename std::iterator_traits<Tintensity>::value_type minint = *std::upper_bound(begint,begint+std::distance(begmz,endmz),0.1);
201 double a = sw_(begmz,endmz);
204 c2d_.
convert2dense(begmz,endmz, begint, resampledintensity_);
209 zerocross_.resize( smoothedintensity_.size()/2 );
210 size_t nrzerocross = simplepicker_( smoothedintensity_.begin( ) , smoothedintensity_.end() , zerocross_.begin(), zerocross_.size());
212 peakmass_.resize(nrzerocross);
215 zerocross_.begin(), zerocross_.begin()+nrzerocross ,
220 peakarea_.resize(nrzerocross);
221 integrator_( zerocross_.begin(), zerocross_.begin() + nrzerocross ,
222 smoothedintensity_.begin(),resampledintensity_.begin(), peakarea_.begin() );
225 peakarea_.resize(nrzerocross);
227 zerocross_.begin(), zerocross_.begin()+nrzerocross ,
231 TReal threshold =
static_cast<TReal
>(minint) * intensitythreshold_;
233 if(maxnumbersofpeaks_ > 0){
234 double threshmax = getNToppeaks();
235 if(threshmax > threshold)
236 threshold = threshmax;
240 if(threshold > 0.01){
248 if(maxnumbersofpeaks_ < peakarea_.size())
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_);
261 peakmass_.begin(),boost::bind(std::greater<TReal>(),_1,threshold));
262 peakmass_.resize(std::distance(peakmass_.begin(),a));
264 peakarea_.begin(),boost::bind(std::greater<TReal>(),_1,threshold));
265 peakarea_.resize(std::distance(peakarea_.begin(),b));
282 return resampledintensity_;
286 return smoothedintensity_;
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)
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'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)
uint32_t maxnumbersofpeaks_
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) ...
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()
TReal intensitythreshold_
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.
SimplePeakArea(TReal integwith)
EQUISPACEINTERPOL Interpolation on a equidistantly spaced grid.
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)
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.
const std::vector< TReal > & getResampledMZ()
PeakIntegrator integrator_
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.
ralab::base::ms::SimplePicker< TReal > simplepicker_