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-2015.
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: Alexandra Zerck $
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  unsigned int zeros = 0)
99  {
100 
101 #ifdef DEBUG_PEAK_PICKING
102  std::cout << "ContinuousWaveletTransformNumIntegration::transform: start " << begin_input->getMZ() << " until " << (end_input - 1)->getMZ() << std::endl;
103 #endif
104  if (fabs(resolution - 1) < 0.0001)
105  {
106  // resolution = 1 corresponds to the cwt at supporting points which have a distance corresponding to the minimal spacing in [begin_input,end_input)
107  SignedSize n = distance(begin_input, end_input);
108  signal_length_ = n;
109 
110  signal_.clear();
111  signal_.resize(n);
112 
113  // TODO avoid to compute the cwt for the zeros in signal
114 #ifdef DEBUG_PEAK_PICKING
115  std::cout << "---------START TRANSFORM---------- \n";
116 #endif
117  InputPeakIterator help = begin_input;
118  for (int i = 0; i < n; ++i)
119  {
120  signal_[i].setMZ(help->getMZ());
121  signal_[i].setIntensity((Peak1D::IntensityType)integrate_(help, begin_input, end_input));
122  ++help;
123  }
124 #ifdef DEBUG_PEAK_PICKING
125  std::cout << "---------END TRANSFORM----------" << std::endl;
126 #endif
127  // no zeropadding
128  begin_right_padding_ = n;
129  end_left_padding_ = -1;
130  }
131  else
132  {
133  SignedSize n = SignedSize(resolution * distance(begin_input, end_input));
134  double origin = begin_input->getMZ();
135  double spacing = ((end_input - 1)->getMZ() - origin) / (n - 1);
136 
137  // zero-padding at the ends?
138  if (zeros > 0)
139  {
140  n += (2 * zeros);
141  }
142 
143 
144  std::vector<double> processed_input(n);
145  signal_.clear();
146  signal_.resize(n);
147 
148  InputPeakIterator it_help = begin_input;
149  if (zeros > 0)
150  {
151  processed_input[0] = it_help->getMZ() - zeros * spacing;
152  for (unsigned int i = 0; i < zeros; ++i) processed_input[i] = 0;
153  }
154  else processed_input[0] = it_help->getIntensity();
155 
156  double x;
157  for (SignedSize k = 1; k < n - (int)zeros; ++k)
158  {
159  x = origin + k * spacing;
160  // go to the real data point next to x
161  while (((it_help + 1) < end_input) && ((it_help + 1)->getMZ() < x))
162  {
163  ++it_help;
164  }
165  processed_input[k] = getInterpolatedValue_(x, it_help);
166  }
167  if (zeros > 0)
168  {
169  for (unsigned int i = 0; i < zeros; ++i) processed_input[n - zeros + i] = 0;
170  }
171 
172  // TODO avoid to compute the cwt for the zeros in signal
173  for (Int i = 0; i < n; ++i)
174  {
175  signal_[i].setMZ(origin + i * spacing);
176  signal_[i].setIntensity((Peak1D::IntensityType)integrate_(processed_input, spacing, i));
177  }
178 
179  if (zeros == 0)
180  {
181  begin_right_padding_ = n;
182  end_left_padding_ = -1;
183  }
184  else
185  {
186  begin_right_padding_ = n - zeros;
187  end_left_padding_ = zeros - 1;
188  }
189  }
190  }
191 
202  virtual void init(double scale, double spacing);
203 
204 protected:
205 
207  template <typename InputPeakIterator>
208  double integrate_(InputPeakIterator x, InputPeakIterator first, InputPeakIterator last)
209  {
210 #ifdef DEBUG_PEAK_PICKING
211  std::cout << "integrate_" << std::endl;
212 #endif
213 
214  double v = 0.;
215  Size middle = wavelet_.size();
216 
217  double start_pos = ((x->getMZ() - (middle * spacing_)) > first->getMZ()) ? (x->getMZ() - (middle * spacing_))
218  : first->getMZ();
219  double end_pos = ((x->getMZ() + (middle * spacing_)) < (last - 1)->getMZ()) ? (x->getMZ() + (middle * spacing_))
220  : (last - 1)->getMZ();
221 
222  InputPeakIterator help = x;
223 
224 #ifdef DEBUG_PEAK_PICKING
225  std::cout << "integrate from middle to start_pos " << help->getMZ() << " until " << start_pos << std::endl;
226 #endif
227 
228  //integrate from middle to start_pos
229  while ((help != first) && ((help - 1)->getMZ() > start_pos))
230  {
231  // search for the corresponding datapoint of help in the wavelet (take the left most adjacent point)
232  double distance = fabs(x->getMZ() - help->getMZ());
233  Size index_w_r = (Size) Math::round(distance / spacing_);
234  if (index_w_r >= wavelet_.size())
235  {
236  index_w_r = wavelet_.size() - 1;
237  }
238  double wavelet_right = wavelet_[index_w_r];
239 
240 #ifdef DEBUG_PEAK_PICKING
241  std::cout << "distance x help " << distance << std::endl;
242  std::cout << "distance in wavelet_ " << index_w_r * spacing_ << std::endl;
243  std::cout << "wavelet_right " << wavelet_right << std::endl;
244 #endif
245 
246  // search for the corresponding datapoint for (help-1) in the wavelet (take the left most adjacent point)
247  distance = fabs(x->getMZ() - (help - 1)->getMZ());
248  Size index_w_l = (Size) Math::round(distance / spacing_);
249  if (index_w_l >= wavelet_.size())
250  {
251  index_w_l = wavelet_.size() - 1;
252  }
253  double wavelet_left = wavelet_[index_w_l];
254 
255  // start the interpolation for the true value in the wavelet
256 
257 #ifdef DEBUG_PEAK_PICKING
258  std::cout << " help-1 " << (help - 1)->getMZ() << " distance x, help-1" << distance << std::endl;
259  std::cout << "distance in wavelet_ " << index_w_l * spacing_ << std::endl;
260  std::cout << "wavelet_ at left " << wavelet_left << std::endl;
261 
262  std::cout << " intensity " << fabs((help - 1)->getMZ() - help->getMZ()) / 2. << " * " << (help - 1)->getIntensity() << " * " << wavelet_left << " + " << (help)->getIntensity() << "* " << wavelet_right
263  << std::endl;
264 #endif
265 
266  v += fabs((help - 1)->getMZ() - help->getMZ()) / 2. * ((help - 1)->getIntensity() * wavelet_left + help->getIntensity() * wavelet_right);
267  --help;
268  }
269 
270 
271  //integrate from middle to end_pos
272  help = x;
273 #ifdef DEBUG_PEAK_PICKING
274  std::cout << "integrate from middle to endpos " << (help)->getMZ() << " until " << end_pos << std::endl;
275 #endif
276  while ((help != (last - 1)) && ((help + 1)->getMZ() < end_pos))
277  {
278  // search for the corresponding datapoint for help in the wavelet (take the left most adjacent point)
279  double distance = fabs(x->getMZ() - help->getMZ());
280  Size index_w_l = (Size) Math::round(distance / spacing_);
281  if (index_w_l >= wavelet_.size())
282  {
283  index_w_l = wavelet_.size() - 1;
284  }
285  double wavelet_left = wavelet_[index_w_l];
286 
287 #ifdef DEBUG_PEAK_PICKING
288  std::cout << " help " << (help)->getMZ() << " distance x, help" << distance << std::endl;
289  std::cout << "distance in wavelet_ " << index_w_l * spacing_ << std::endl;
290  std::cout << "wavelet_ at left " << wavelet_left << std::endl;
291 #endif
292 
293  // search for the corresponding datapoint for (help+1) in the wavelet (take the left most adjacent point)
294  distance = fabs(x->getMZ() - (help + 1)->getMZ());
295  Size index_w_r = (Size) Math::round(distance / spacing_);
296  if (index_w_r >= wavelet_.size())
297  {
298  index_w_r = wavelet_.size() - 1;
299  }
300  double wavelet_right = wavelet_[index_w_r];
301 
302 #ifdef DEBUG_PEAK_PICKING
303  std::cout << " help+1 " << (help + 1)->getMZ() << " distance x, help+1" << distance << std::endl;
304  std::cout << "distance in wavelet_ " << index_w_r * spacing_ << std::endl;
305  std::cout << "wavelet_ at right " << wavelet_right << std::endl;
306 #endif
307 
308  v += fabs(help->getMZ() - (help + 1)->getMZ()) / 2. * (help->getIntensity() * wavelet_left + (help + 1)->getIntensity() * wavelet_right);
309  ++help;
310  }
311 
312 
313 #ifdef DEBUG_PEAK_PICKING
314  std::cout << "return" << (v / sqrt(scale_)) << std::endl;
315 #endif
316  return v / sqrt(scale_);
317  }
318 
320  double integrate_(const std::vector<double> & processed_input, double spacing_data, int index);
321 
323  inline double marr_(double x)
324  {
325  return (1 - x * x) * exp(-x * x / 2);
326  }
327 
328  };
329 } //namespace OpenMS
330 #endif
const double k
ContinuousWaveletTransform::PeakConstIterator PeakConstIterator
Raw data const iterator type.
Definition: ContinuousWaveletTransformNumIntegration.h:62
void transform(InputPeakIterator begin_input, InputPeakIterator end_input, float resolution, unsigned int zeros=0)
Computes the wavelet transform of a given raw data interval [begin_input,end_input) ...
Definition: ContinuousWaveletTransformNumIntegration.h:95
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:128
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:137
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:208
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:121
virtual ~ContinuousWaveletTransformNumIntegration()
Destructor.
Definition: ContinuousWaveletTransformNumIntegration.h:79
double marr_(double x)
Computes the marr wavelet at position x.
Definition: ContinuousWaveletTransformNumIntegration.h:323
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:96

OpenMS / TOPP release 2.0.0 Documentation generated on Fri May 29 2015 17:20:22 using doxygen 1.8.9.1