Home  · Classes  · Annotated Classes  · Modules  · Members  · Namespaces  · Related Pages
ContinuousWaveletTransformNumIntegration.h
Go to the documentation of this file.
1 // --------------------------------------------------------------------------
2 // OpenMS -- Open-Source Mass Spectrometry
3 // --------------------------------------------------------------------------
4 // Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
5 // ETH Zurich, and Freie Universitaet Berlin 2002-2017.
6 //
7 // This software is released under a three-clause BSD license:
8 // * Redistributions of source code must retain the above copyright
9 // notice, this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright
11 // notice, this list of conditions and the following disclaimer in the
12 // documentation and/or other materials provided with the distribution.
13 // * Neither the name of any author or any participating institution
14 // may be used to endorse or promote products derived from this software
15 // without specific prior written permission.
16 // For a full list of authors, refer to the file AUTHORS.
17 // --------------------------------------------------------------------------
18 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 // ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
22 // INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
23 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
24 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
25 // OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
26 // WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
27 // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
28 // ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 //
30 // --------------------------------------------------------------------------
31 // $Maintainer: Timo Sachsenberg $
32 // $Authors: Eva Lange $
33 // --------------------------------------------------------------------------
34 //
35 
36 #ifndef OPENMS_TRANSFORMATIONS_RAW2PEAK_CONTINUOUSWAVELETTRANSFORMNUMINTEGRATION_H
37 #define OPENMS_TRANSFORMATIONS_RAW2PEAK_CONTINUOUSWAVELETTRANSFORMNUMINTEGRATION_H
38 
39 #include <cmath>
40 
43 
44 
45 #ifdef DEBUG_PEAK_PICKING
46 #include <iostream>
47 #include <fstream>
48 #endif
49 
50 namespace OpenMS
51 {
59  {
60 public:
63 
71 
72 
76  {}
77 
80 
94  template <typename InputPeakIterator>
95  void transform(InputPeakIterator begin_input,
96  InputPeakIterator end_input,
97  float resolution)
98  {
99 #ifdef DEBUG_PEAK_PICKING
100  std::cout << "ContinuousWaveletTransformNumIntegration::transform: start " << begin_input->getMZ() << " until " << (end_input - 1)->getMZ() << std::endl;
101 #endif
102  if (fabs(resolution - 1) < 0.0001)
103  {
104  // resolution = 1 corresponds to the cwt at supporting points which have a distance corresponding to the minimal spacing in [begin_input,end_input)
105  SignedSize n = distance(begin_input, end_input);
106  signal_length_ = n;
107 
108  signal_.clear();
109  signal_.resize(n);
110 
111  // TODO avoid to compute the cwt for the zeros in signal
112 #ifdef DEBUG_PEAK_PICKING
113  std::cout << "---------START TRANSFORM---------- \n";
114 #endif
115  InputPeakIterator help = begin_input;
116  for (int i = 0; i < n; ++i)
117  {
118  signal_[i].setMZ(help->getMZ());
119  signal_[i].setIntensity((Peak1D::IntensityType)integrate_(help, begin_input, end_input));
120  ++help;
121  }
122 #ifdef DEBUG_PEAK_PICKING
123  std::cout << "---------END TRANSFORM----------" << std::endl;
124 #endif
125  // no zeropadding
126  begin_right_padding_ = n;
127  end_left_padding_ = -1;
128  }
129  else
130  {
131  SignedSize n = SignedSize(resolution * distance(begin_input, end_input));
132  double origin = begin_input->getMZ();
133  double spacing = ((end_input - 1)->getMZ() - origin) / (n - 1);
134 
135  std::vector<double> processed_input(n);
136  signal_.clear();
137  signal_.resize(n);
138 
139  InputPeakIterator it_help = begin_input;
140  processed_input[0] = it_help->getIntensity();
141 
142  double x;
143  for (SignedSize k = 1; k < n; ++k)
144  {
145  x = origin + k * spacing;
146  // go to the real data point next to x
147  while (((it_help + 1) < end_input) && ((it_help + 1)->getMZ() < x))
148  {
149  ++it_help;
150  }
151  processed_input[k] = getInterpolatedValue_(x, it_help);
152  }
153 
154  // TODO avoid to compute the cwt for the zeros in signal
155  for (Int i = 0; i < n; ++i)
156  {
157  signal_[i].setMZ(origin + i * spacing);
158  signal_[i].setIntensity((Peak1D::IntensityType)integrate_(processed_input, spacing, i));
159  }
160 
161  begin_right_padding_ = n;
162  end_left_padding_ = -1;
163  }
164  }
165 
176  virtual void init(double scale, double spacing);
177 
178 protected:
179 
181  template <typename InputPeakIterator>
182  double integrate_(InputPeakIterator x, InputPeakIterator first, InputPeakIterator last)
183  {
184 #ifdef DEBUG_PEAK_PICKING
185  std::cout << "integrate_" << std::endl;
186 #endif
187 
188  double v = 0.;
189  double middle_spacing = wavelet_.size() * spacing_;
190 
191  double start_pos = ((x->getMZ() - middle_spacing) > first->getMZ()) ? (x->getMZ() - middle_spacing)
192  : first->getMZ();
193  double end_pos = ((x->getMZ() + middle_spacing) < (last - 1)->getMZ()) ? (x->getMZ() + middle_spacing)
194  : (last - 1)->getMZ();
195 
196  InputPeakIterator help = x;
197 
198 #ifdef DEBUG_PEAK_PICKING
199  std::cout << "integrate from middle to start_pos " << help->getMZ() << " until " << start_pos << std::endl;
200 #endif
201 
202  //integrate from middle to start_pos
203  while ((help != first) && ((help - 1)->getMZ() > start_pos))
204  {
205  // search for the corresponding data point of help in the wavelet (take the left most adjacent point)
206  double distance = fabs(x->getMZ() - help->getMZ());
207  Size index_w_r = (Size) Math::round(distance / spacing_);
208  if (index_w_r >= wavelet_.size())
209  {
210  index_w_r = wavelet_.size() - 1;
211  }
212  double wavelet_right = wavelet_[index_w_r];
213 
214 #ifdef DEBUG_PEAK_PICKING
215  std::cout << "distance x help " << distance << std::endl;
216  std::cout << "distance in wavelet_ " << index_w_r * spacing_ << std::endl;
217  std::cout << "wavelet_right " << wavelet_right << std::endl;
218 #endif
219 
220  // search for the corresponding datapoint for (help-1) in the wavelet (take the left most adjacent point)
221  distance = fabs(x->getMZ() - (help - 1)->getMZ());
222  Size index_w_l = (Size) Math::round(distance / spacing_);
223  if (index_w_l >= wavelet_.size())
224  {
225  index_w_l = wavelet_.size() - 1;
226  }
227  double wavelet_left = wavelet_[index_w_l];
228 
229  // start the interpolation for the true value in the wavelet
230 
231 #ifdef DEBUG_PEAK_PICKING
232  std::cout << " help-1 " << (help - 1)->getMZ() << " distance x, help-1" << distance << std::endl;
233  std::cout << "distance in wavelet_ " << index_w_l * spacing_ << std::endl;
234  std::cout << "wavelet_ at left " << wavelet_left << std::endl;
235 
236  std::cout << " intensity " << fabs((help - 1)->getMZ() - help->getMZ()) / 2. << " * " << (help - 1)->getIntensity() << " * " << wavelet_left << " + " << (help)->getIntensity() << "* " << wavelet_right
237  << std::endl;
238 #endif
239 
240  v += fabs((help - 1)->getMZ() - help->getMZ()) / 2. * ((help - 1)->getIntensity() * wavelet_left + help->getIntensity() * wavelet_right);
241  --help;
242  }
243 
244 
245  //integrate from middle to end_pos
246  help = x;
247 #ifdef DEBUG_PEAK_PICKING
248  std::cout << "integrate from middle to endpos " << (help)->getMZ() << " until " << end_pos << std::endl;
249 #endif
250  while ((help != (last - 1)) && ((help + 1)->getMZ() < end_pos))
251  {
252  // search for the corresponding datapoint for help in the wavelet (take the left most adjacent point)
253  double distance = fabs(x->getMZ() - help->getMZ());
254  Size index_w_l = (Size) Math::round(distance / spacing_);
255  if (index_w_l >= wavelet_.size())
256  {
257  index_w_l = wavelet_.size() - 1;
258  }
259  double wavelet_left = wavelet_[index_w_l];
260 
261 #ifdef DEBUG_PEAK_PICKING
262  std::cout << " help " << (help)->getMZ() << " distance x, help" << distance << std::endl;
263  std::cout << "distance in wavelet_ " << index_w_l * spacing_ << std::endl;
264  std::cout << "wavelet_ at left " << wavelet_left << std::endl;
265 #endif
266 
267  // search for the corresponding datapoint for (help+1) in the wavelet (take the left most adjacent point)
268  distance = fabs(x->getMZ() - (help + 1)->getMZ());
269  Size index_w_r = (Size) Math::round(distance / spacing_);
270  if (index_w_r >= wavelet_.size())
271  {
272  index_w_r = wavelet_.size() - 1;
273  }
274  double wavelet_right = wavelet_[index_w_r];
275 
276 #ifdef DEBUG_PEAK_PICKING
277  std::cout << " help+1 " << (help + 1)->getMZ() << " distance x, help+1" << distance << std::endl;
278  std::cout << "distance in wavelet_ " << index_w_r * spacing_ << std::endl;
279  std::cout << "wavelet_ at right " << wavelet_right << std::endl;
280 #endif
281 
282  v += fabs(help->getMZ() - (help + 1)->getMZ()) / 2. * (help->getIntensity() * wavelet_left + (help + 1)->getIntensity() * wavelet_right);
283  ++help;
284  }
285 
286 
287 #ifdef DEBUG_PEAK_PICKING
288  std::cout << "return" << (v / sqrt(scale_)) << std::endl;
289 #endif
290  return v / sqrt(scale_);
291  }
292 
294  double integrate_(const std::vector<double> & processed_input, double spacing_data, int index);
295 
297  inline double marr_(const double x) const
298  {
299  return (1 - x * x) * exp(-x * x / 2);
300  }
301 
302  };
303 } //namespace OpenMS
304 #endif
const double k
double marr_(const double x) const
Computes the Marr wavelet at position x.
Definition: ContinuousWaveletTransformNumIntegration.h:297
ContinuousWaveletTransform::PeakConstIterator PeakConstIterator
Raw data const iterator type.
Definition: ContinuousWaveletTransformNumIntegration.h:62
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:135
double scale_
Spacing and scale of the wavelet and length of the signal.
Definition: ContinuousWaveletTransform.h:224
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:47
This class is the base class of the continuous wavelet transformation.
Definition: ContinuousWaveletTransform.h:47
This class computes the continuous wavelet transformation using a marr wavelet.
Definition: ContinuousWaveletTransformNumIntegration.h:57
T round(T x)
Rounds the value.
Definition: MathFunctions.h:138
double integrate_(InputPeakIterator x, InputPeakIterator first, InputPeakIterator last)
Computes the convolution of the wavelet and the raw data at position x with resolution = 1...
Definition: ContinuousWaveletTransformNumIntegration.h:182
void transform(InputPeakIterator begin_input, InputPeakIterator end_input, float resolution)
Computes the wavelet transform of a given raw data interval [begin_input,end_input) ...
Definition: ContinuousWaveletTransformNumIntegration.h:95
SignedSize begin_right_padding_
Definition: ContinuousWaveletTransform.h:232
ContinuousWaveletTransformNumIntegration()
Constructor.
Definition: ContinuousWaveletTransformNumIntegration.h:74
SignedSize end_left_padding_
Definition: ContinuousWaveletTransform.h:231
SignedSize signal_length_
Definition: ContinuousWaveletTransform.h:226
std::vector< Peak1D > signal_
The transformed signal.
Definition: ContinuousWaveletTransform.h:218
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:128
virtual ~ContinuousWaveletTransformNumIntegration()
Destructor.
Definition: ContinuousWaveletTransformNumIntegration.h:79
std::vector< double > wavelet_
The pre-tabulated wavelet used for the transform.
Definition: ContinuousWaveletTransform.h:221
std::vector< Peak1D >::const_iterator PeakConstIterator
Raw data const iterator type.
Definition: ContinuousWaveletTransform.h:51
double spacing_
Definition: ContinuousWaveletTransform.h:225
int Int
Signed integer type.
Definition: Types.h:103

OpenMS / TOPP release 2.3.0 Documentation generated on Tue Jan 9 2018 18:21:59 using doxygen 1.8.13