OpenMS
ContinuousWaveletTransformNumIntegration.h
Go to the documentation of this file.
1 // Copyright (c) 2002-2023, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin
2 // SPDX-License-Identifier: BSD-3-Clause
3 //
4 // --------------------------------------------------------------------------
5 // $Maintainer: Timo Sachsenberg $
6 // $Authors: Eva Lange $
7 // --------------------------------------------------------------------------
8 //
9 
10 #pragma once
11 
12 #include <cmath>
13 
16 
17 
18 #ifdef DEBUG_PEAK_PICKING
19 #include <iostream>
20 #include <fstream>
21 #endif
22 
23 namespace OpenMS
24 {
32  {
33 public:
36 
44 
45 
49  {}
50 
53 
67  template <typename InputPeakIterator>
68  void transform(InputPeakIterator begin_input,
69  InputPeakIterator end_input,
70  float resolution)
71  {
72 #ifdef DEBUG_PEAK_PICKING
73  std::cout << "ContinuousWaveletTransformNumIntegration::transform: start " << begin_input->getMZ() << " until " << (end_input - 1)->getMZ() << std::endl;
74 #endif
75  if (fabs(resolution - 1) < 0.0001)
76  {
77  // resolution = 1 corresponds to the cwt at supporting points which have a distance corresponding to the minimal spacing in [begin_input,end_input)
78  SignedSize n = distance(begin_input, end_input);
79  signal_length_ = n;
80 
81  signal_.clear();
82  signal_.resize(n);
83 
84  // TODO avoid to compute the cwt for the zeros in signal
85 #ifdef DEBUG_PEAK_PICKING
86  std::cout << "---------START TRANSFORM---------- \n";
87 #endif
88  InputPeakIterator help = begin_input;
89  for (int i = 0; i < n; ++i)
90  {
91  signal_[i].setMZ(help->getMZ());
92  signal_[i].setIntensity((Peak1D::IntensityType)integrate_(help, begin_input, end_input));
93  ++help;
94  }
95 #ifdef DEBUG_PEAK_PICKING
96  std::cout << "---------END TRANSFORM----------" << std::endl;
97 #endif
98  // no zeropadding
99  begin_right_padding_ = n;
100  end_left_padding_ = -1;
101  }
102  else
103  {
104  SignedSize n = SignedSize(resolution * distance(begin_input, end_input));
105  double origin = begin_input->getMZ();
106  double spacing = ((end_input - 1)->getMZ() - origin) / (n - 1);
107 
108  std::vector<double> processed_input(n);
109  signal_.clear();
110  signal_.resize(n);
111 
112  InputPeakIterator it_help = begin_input;
113  processed_input[0] = it_help->getIntensity();
114 
115  double x;
116  for (SignedSize k = 1; k < n; ++k)
117  {
118  x = origin + k * spacing;
119  // go to the real data point next to x
120  while (((it_help + 1) < end_input) && ((it_help + 1)->getMZ() < x))
121  {
122  ++it_help;
123  }
124  processed_input[k] = getInterpolatedValue_(x, it_help);
125  }
126 
127  // TODO avoid to compute the cwt for the zeros in signal
128  for (Int i = 0; i < n; ++i)
129  {
130  signal_[i].setMZ(origin + i * spacing);
131  signal_[i].setIntensity((Peak1D::IntensityType)integrate_(processed_input, spacing, i));
132  }
133 
134  begin_right_padding_ = n;
135  end_left_padding_ = -1;
136  }
137  }
138 
149  void init(double scale, double spacing) override;
150 
151 protected:
152 
154  template <typename InputPeakIterator>
155  double integrate_(InputPeakIterator x, InputPeakIterator first, InputPeakIterator last)
156  {
157 #ifdef DEBUG_PEAK_PICKING
158  std::cout << "integrate_" << std::endl;
159 #endif
160 
161  double v = 0.;
162  double middle_spacing = wavelet_.size() * spacing_;
163 
164  double start_pos = ((x->getMZ() - middle_spacing) > first->getMZ()) ? (x->getMZ() - middle_spacing)
165  : first->getMZ();
166  double end_pos = ((x->getMZ() + middle_spacing) < (last - 1)->getMZ()) ? (x->getMZ() + middle_spacing)
167  : (last - 1)->getMZ();
168 
169  InputPeakIterator help = x;
170 
171 #ifdef DEBUG_PEAK_PICKING
172  std::cout << "integrate from middle to start_pos " << help->getMZ() << " until " << start_pos << std::endl;
173 #endif
174 
175  //integrate from middle to start_pos
176  while ((help != first) && ((help - 1)->getMZ() > start_pos))
177  {
178  // search for the corresponding data point of help in the wavelet (take the left most adjacent point)
179  double distance = fabs(x->getMZ() - help->getMZ());
180  Size index_w_r = (Size) Math::round(distance / spacing_);
181  if (index_w_r >= wavelet_.size())
182  {
183  index_w_r = wavelet_.size() - 1;
184  }
185  double wavelet_right = wavelet_[index_w_r];
186 
187 #ifdef DEBUG_PEAK_PICKING
188  std::cout << "distance x help " << distance << std::endl;
189  std::cout << "distance in wavelet_ " << index_w_r * spacing_ << std::endl;
190  std::cout << "wavelet_right " << wavelet_right << std::endl;
191 #endif
192 
193  // search for the corresponding datapoint for (help-1) in the wavelet (take the left most adjacent point)
194  distance = fabs(x->getMZ() - (help - 1)->getMZ());
195  Size index_w_l = (Size) Math::round(distance / spacing_);
196  if (index_w_l >= wavelet_.size())
197  {
198  index_w_l = wavelet_.size() - 1;
199  }
200  double wavelet_left = wavelet_[index_w_l];
201 
202  // start the interpolation for the true value in the wavelet
203 
204 #ifdef DEBUG_PEAK_PICKING
205  std::cout << " help-1 " << (help - 1)->getMZ() << " distance x, help-1" << distance << std::endl;
206  std::cout << "distance in wavelet_ " << index_w_l * spacing_ << std::endl;
207  std::cout << "wavelet_ at left " << wavelet_left << std::endl;
208 
209  std::cout << " intensity " << fabs((help - 1)->getMZ() - help->getMZ()) / 2. << " * " << (help - 1)->getIntensity() << " * " << wavelet_left << " + " << (help)->getIntensity() << "* " << wavelet_right
210  << std::endl;
211 #endif
212 
213  v += fabs((help - 1)->getMZ() - help->getMZ()) / 2. * ((help - 1)->getIntensity() * wavelet_left + help->getIntensity() * wavelet_right);
214  --help;
215  }
216 
217 
218  //integrate from middle to end_pos
219  help = x;
220 #ifdef DEBUG_PEAK_PICKING
221  std::cout << "integrate from middle to endpos " << (help)->getMZ() << " until " << end_pos << std::endl;
222 #endif
223  while ((help != (last - 1)) && ((help + 1)->getMZ() < end_pos))
224  {
225  // search for the corresponding datapoint for help in the wavelet (take the left most adjacent point)
226  double distance = fabs(x->getMZ() - help->getMZ());
227  Size index_w_l = (Size) Math::round(distance / spacing_);
228  if (index_w_l >= wavelet_.size())
229  {
230  index_w_l = wavelet_.size() - 1;
231  }
232  double wavelet_left = wavelet_[index_w_l];
233 
234 #ifdef DEBUG_PEAK_PICKING
235  std::cout << " help " << (help)->getMZ() << " distance x, help" << distance << std::endl;
236  std::cout << "distance in wavelet_ " << index_w_l * spacing_ << std::endl;
237  std::cout << "wavelet_ at left " << wavelet_left << std::endl;
238 #endif
239 
240  // search for the corresponding datapoint for (help+1) in the wavelet (take the left most adjacent point)
241  distance = fabs(x->getMZ() - (help + 1)->getMZ());
242  Size index_w_r = (Size) Math::round(distance / spacing_);
243  if (index_w_r >= wavelet_.size())
244  {
245  index_w_r = wavelet_.size() - 1;
246  }
247  double wavelet_right = wavelet_[index_w_r];
248 
249 #ifdef DEBUG_PEAK_PICKING
250  std::cout << " help+1 " << (help + 1)->getMZ() << " distance x, help+1" << distance << std::endl;
251  std::cout << "distance in wavelet_ " << index_w_r * spacing_ << std::endl;
252  std::cout << "wavelet_ at right " << wavelet_right << std::endl;
253 #endif
254 
255  v += fabs(help->getMZ() - (help + 1)->getMZ()) / 2. * (help->getIntensity() * wavelet_left + (help + 1)->getIntensity() * wavelet_right);
256  ++help;
257  }
258 
259 
260 #ifdef DEBUG_PEAK_PICKING
261  std::cout << "return" << (v / sqrt(scale_)) << std::endl;
262 #endif
263  return v / sqrt(scale_);
264  }
265 
267  double integrate_(const std::vector<double> & processed_input, double spacing_data, int index);
268 
270  inline double marr_(const double x) const
271  {
272  return (1 - x * x) * exp(-x * x / 2);
273  }
274 
275  };
276 } //namespace OpenMS
This class computes the continuous wavelet transformation using a marr wavelet.
Definition: ContinuousWaveletTransformNumIntegration.h:32
ContinuousWaveletTransform::PeakConstIterator PeakConstIterator
Profile data const iterator type.
Definition: ContinuousWaveletTransformNumIntegration.h:35
void transform(InputPeakIterator begin_input, InputPeakIterator end_input, float resolution)
Computes the wavelet transform of a given profile data interval [begin_input,end_input)
Definition: ContinuousWaveletTransformNumIntegration.h:68
double integrate_(InputPeakIterator x, InputPeakIterator first, InputPeakIterator last)
Computes the convolution of the wavelet and the profile data at position x with resolution = 1.
Definition: ContinuousWaveletTransformNumIntegration.h:155
ContinuousWaveletTransformNumIntegration()
Constructor.
Definition: ContinuousWaveletTransformNumIntegration.h:47
double integrate_(const std::vector< double > &processed_input, double spacing_data, int index)
Computes the convolution of the wavelet and the profile data at position x with resolution > 1.
void init(double scale, double spacing) override
Perform necessary preprocessing steps like tabulating the Wavelet.
double marr_(const double x) const
Computes the Marr wavelet at position x.
Definition: ContinuousWaveletTransformNumIntegration.h:270
~ContinuousWaveletTransformNumIntegration() override
Destructor.
Definition: ContinuousWaveletTransformNumIntegration.h:52
This class is the base class of the continuous wavelet transformation.
Definition: ContinuousWaveletTransform.h:21
double scale_
Spacing and scale of the wavelet and length of the signal.
Definition: ContinuousWaveletTransform.h:197
std::vector< double > wavelet_
The pre-tabulated wavelet used for the transform.
Definition: ContinuousWaveletTransform.h:194
SignedSize end_left_padding_
Definition: ContinuousWaveletTransform.h:204
double spacing_
Definition: ContinuousWaveletTransform.h:198
SignedSize begin_right_padding_
Definition: ContinuousWaveletTransform.h:205
std::vector< Peak1D >::const_iterator PeakConstIterator
Raw data const iterator type.
Definition: ContinuousWaveletTransform.h:24
SignedSize signal_length_
Definition: ContinuousWaveletTransform.h:199
std::vector< Peak1D > signal_
The transformed signal.
Definition: ContinuousWaveletTransform.h:191
float IntensityType
Intensity type.
Definition: Peak1D.h:36
int Int
Signed integer type.
Definition: Types.h:76
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:108
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:101
T round(T x)
Rounds the value.
Definition: MathFunctions.h:184
const double k
Definition: Constants.h:132
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:22