| OpenMS 2.5.0

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-2020.
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 #pragma once
37
38 #include <cmath>
39
40 #include <OpenMS/MATH/MISC/MathFunctions.h>
41 #include <OpenMS/TRANSFORMATIONS/RAW2PEAK/ContinuousWaveletTransform.h>
42
43
44 #ifdef DEBUG_PEAK_PICKING
45 #include <iostream>
46 #include <fstream>
47 #endif
48
49 namespace OpenMS
50 {
56 class OPENMS_DLLAPI ContinuousWaveletTransformNumIntegration :
57 public ContinuousWaveletTransform
58 {
59 public:
61 typedef ContinuousWaveletTransform::PeakConstIterator PeakConstIterator;
62
63 using ContinuousWaveletTransform::signal_;
64 using ContinuousWaveletTransform::wavelet_;
65 using ContinuousWaveletTransform::scale_;
66 using ContinuousWaveletTransform::spacing_;
67 using ContinuousWaveletTransform::end_left_padding_;
68 using ContinuousWaveletTransform::begin_right_padding_;
69 using ContinuousWaveletTransform::signal_length_;
70
71
73 ContinuousWaveletTransformNumIntegration() :
74 ContinuousWaveletTransform()
75 {}
76
78 ~ContinuousWaveletTransformNumIntegration() override {}
79
93 template <typename InputPeakIterator>
94 void transform(InputPeakIterator begin_input,
95 InputPeakIterator end_input,
96 float resolution)
97 {
98 #ifdef DEBUG_PEAK_PICKING
99 std::cout << "ContinuousWaveletTransformNumIntegration::transform: start " << begin_input->getMZ() << " until " << (end_input - 1)->getMZ() << std::endl;
100 #endif
101 if (fabs(resolution - 1) < 0.0001)
102 {
103 // resolution = 1 corresponds to the cwt at supporting points which have a distance corresponding to the minimal spacing in [begin_input,end_input)
104 SignedSize n = distance(begin_input, end_input);
105 signal_length_ = n;
106
107 signal_.clear();
108 signal_.resize(n);
109
110 // TODO avoid to compute the cwt for the zeros in signal
111 #ifdef DEBUG_PEAK_PICKING
112 std::cout << "---------START TRANSFORM---------- \n";
113 #endif
114 InputPeakIterator help = begin_input;
115 for (int i = 0; i < n; ++i)
116 {
117 signal_[i].setMZ(help->getMZ());
118 signal_[i].setIntensity((Peak1D::IntensityType)integrate_(help, begin_input, end_input));
119 ++help;
120 }
121 #ifdef DEBUG_PEAK_PICKING
122 std::cout << "---------END TRANSFORM----------" << std::endl;
123 #endif
124 // no zeropadding
125 begin_right_padding_ = n;
126 end_left_padding_ = -1;
127 }
128 else
129 {
130 SignedSize n = SignedSize(resolution * distance(begin_input, end_input));
131 double origin = begin_input->getMZ();
132 double spacing = ((end_input - 1)->getMZ() - origin) / (n - 1);
133
134 std::vector<double> processed_input(n);
135 signal_.clear();
136 signal_.resize(n);
137
138 InputPeakIterator it_help = begin_input;
139 processed_input[0] = it_help->getIntensity();
140
141 double x;
142 for (SignedSize k = 1; k < n; ++k)
143 {
144 x = origin + k * spacing;
145 // go to the real data point next to x
146 while (((it_help + 1) < end_input) && ((it_help + 1)->getMZ() < x))
147 {
148 ++it_help;
149 }
150 processed_input[k] = getInterpolatedValue_(x, it_help);
151 }
152
153 // TODO avoid to compute the cwt for the zeros in signal
154 for (Int i = 0; i < n; ++i)
155 {
156 signal_[i].setMZ(origin + i * spacing);
157 signal_[i].setIntensity((Peak1D::IntensityType)integrate_(processed_input, spacing, i));
158 }
159
160 begin_right_padding_ = n;
161 end_left_padding_ = -1;
162 }
163 }
164
175 void init(double scale, double spacing) override;
176
177 protected:
178
180 template <typename InputPeakIterator>
181 double integrate_(InputPeakIterator x, InputPeakIterator first, InputPeakIterator last)
182 {
183 #ifdef DEBUG_PEAK_PICKING
184 std::cout << "integrate_" << std::endl;
185 #endif
186
187 double v = 0.;
188 double middle_spacing = wavelet_.size() * spacing_;
189
190 double start_pos = ((x->getMZ() - middle_spacing) > first->getMZ()) ? (x->getMZ() - middle_spacing)
191 : first->getMZ();
192 double end_pos = ((x->getMZ() + middle_spacing) < (last - 1)->getMZ()) ? (x->getMZ() + middle_spacing)
193 : (last - 1)->getMZ();
194
195 InputPeakIterator help = x;
196
197 #ifdef DEBUG_PEAK_PICKING
198 std::cout << "integrate from middle to start_pos " << help->getMZ() << " until " << start_pos << std::endl;
199 #endif
200
201 //integrate from middle to start_pos
202 while ((help != first) && ((help - 1)->getMZ() > start_pos))
203 {
204 // search for the corresponding data point of help in the wavelet (take the left most adjacent point)
205 double distance = fabs(x->getMZ() - help->getMZ());
206 Size index_w_r = (Size) Math::round(distance / spacing_);
207 if (index_w_r >= wavelet_.size())
208 {
209 index_w_r = wavelet_.size() - 1;
210 }
211 double wavelet_right = wavelet_[index_w_r];
212
213 #ifdef DEBUG_PEAK_PICKING
214 std::cout << "distance x help " << distance << std::endl;
215 std::cout << "distance in wavelet_ " << index_w_r * spacing_ << std::endl;
216 std::cout << "wavelet_right " << wavelet_right << std::endl;
217 #endif
218
219 // search for the corresponding datapoint for (help-1) in the wavelet (take the left most adjacent point)
220 distance = fabs(x->getMZ() - (help - 1)->getMZ());
221 Size index_w_l = (Size) Math::round(distance / spacing_);
222 if (index_w_l >= wavelet_.size())
223 {
224 index_w_l = wavelet_.size() - 1;
225 }
226 double wavelet_left = wavelet_[index_w_l];
227
228 // start the interpolation for the true value in the wavelet
229
230 #ifdef DEBUG_PEAK_PICKING
231 std::cout << " help-1 " << (help - 1)->getMZ() << " distance x, help-1" << distance << std::endl;
232 std::cout << "distance in wavelet_ " << index_w_l * spacing_ << std::endl;
233 std::cout << "wavelet_ at left " << wavelet_left << std::endl;
234
235 std::cout << " intensity " << fabs((help - 1)->getMZ() - help->getMZ()) / 2. << " * " << (help - 1)->getIntensity() << " * " << wavelet_left << " + " << (help)->getIntensity() << "* " << wavelet_right
236 << std::endl;
237 #endif
238
239 v += fabs((help - 1)->getMZ() - help->getMZ()) / 2. * ((help - 1)->getIntensity() * wavelet_left + help->getIntensity() * wavelet_right);
240 --help;
241 }
242
243
244 //integrate from middle to end_pos
245 help = x;
246 #ifdef DEBUG_PEAK_PICKING
247 std::cout << "integrate from middle to endpos " << (help)->getMZ() << " until " << end_pos << std::endl;
248 #endif
249 while ((help != (last - 1)) && ((help + 1)->getMZ() < end_pos))
250 {
251 // search for the corresponding datapoint for help in the wavelet (take the left most adjacent point)
252 double distance = fabs(x->getMZ() - help->getMZ());
253 Size index_w_l = (Size) Math::round(distance / spacing_);
254 if (index_w_l >= wavelet_.size())
255 {
256 index_w_l = wavelet_.size() - 1;
257 }
258 double wavelet_left = wavelet_[index_w_l];
259
260 #ifdef DEBUG_PEAK_PICKING
261 std::cout << " help " << (help)->getMZ() << " distance x, help" << distance << std::endl;
262 std::cout << "distance in wavelet_ " << index_w_l * spacing_ << std::endl;
263 std::cout << "wavelet_ at left " << wavelet_left << std::endl;
264 #endif
265
266 // search for the corresponding datapoint for (help+1) in the wavelet (take the left most adjacent point)
267 distance = fabs(x->getMZ() - (help + 1)->getMZ());
268 Size index_w_r = (Size) Math::round(distance / spacing_);
269 if (index_w_r >= wavelet_.size())
270 {
271 index_w_r = wavelet_.size() - 1;
272 }
273 double wavelet_right = wavelet_[index_w_r];
274
275 #ifdef DEBUG_PEAK_PICKING
276 std::cout << " help+1 " << (help + 1)->getMZ() << " distance x, help+1" << distance << std::endl;
277 std::cout << "distance in wavelet_ " << index_w_r * spacing_ << std::endl;
278 std::cout << "wavelet_ at right " << wavelet_right << std::endl;
279 #endif
280
281 v += fabs(help->getMZ() - (help + 1)->getMZ()) / 2. * (help->getIntensity() * wavelet_left + (help + 1)->getIntensity() * wavelet_right);
282 ++help;
283 }
284
285
286 #ifdef DEBUG_PEAK_PICKING
287 std::cout << "return" << (v / sqrt(scale_)) << std::endl;
288 #endif
289 return v / sqrt(scale_);
290 }
291
293 double integrate_(const std::vector<double> & processed_input, double spacing_data, int index);
294
296 inline double marr_(const double x) const
297 {
298 return (1 - x * x) * exp(-x * x / 2);
299 }
300
301 };
302 } //namespace OpenMS
OpenMS::SignedSize
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:134
OpenMS::ContinuousWaveletTransform::begin_right_padding_
SignedSize begin_right_padding_
Definition: ContinuousWaveletTransform.h:231
OpenMS::ContinuousWaveletTransformNumIntegration::marr_
double marr_(const double x) const
Computes the Marr wavelet at position x.
Definition: ContinuousWaveletTransformNumIntegration.h:296
OpenMS::ContinuousWaveletTransform
This class is the base class of the continuous wavelet transformation.
Definition: ContinuousWaveletTransform.h:46
OpenMS::Constants::k
const double k
MathFunctions.h
OpenMS::ContinuousWaveletTransformNumIntegration::ContinuousWaveletTransformNumIntegration
ContinuousWaveletTransformNumIntegration()
Constructor.
Definition: ContinuousWaveletTransformNumIntegration.h:73
OpenMS::ContinuousWaveletTransform::PeakConstIterator
std::vector< Peak1D >::const_iterator PeakConstIterator
Raw data const iterator type.
Definition: ContinuousWaveletTransform.h:50
OpenMS
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:46
OpenMS::ContinuousWaveletTransformNumIntegration::~ContinuousWaveletTransformNumIntegration
~ContinuousWaveletTransformNumIntegration() override
Destructor.
Definition: ContinuousWaveletTransformNumIntegration.h:78
OpenMS::ContinuousWaveletTransformNumIntegration
This class computes the continuous wavelet transformation using a marr wavelet.
Definition: ContinuousWaveletTransformNumIntegration.h:56
OpenMS::ContinuousWaveletTransform::signal_length_
SignedSize signal_length_
Definition: ContinuousWaveletTransform.h:225
OpenMS::ContinuousWaveletTransform::scale_
double scale_
Spacing and scale of the wavelet and length of the signal.
Definition: ContinuousWaveletTransform.h:223
OpenMS::ContinuousWaveletTransform::end_left_padding_
SignedSize end_left_padding_
Definition: ContinuousWaveletTransform.h:230
ContinuousWaveletTransform.h
float
OpenMS::ContinuousWaveletTransformNumIntegration::PeakConstIterator
ContinuousWaveletTransform::PeakConstIterator PeakConstIterator
Profile data const iterator type.
Definition: ContinuousWaveletTransformNumIntegration.h:61
int
OpenMS::ContinuousWaveletTransform::spacing_
double spacing_
Definition: ContinuousWaveletTransform.h:224
OpenMS::ContinuousWaveletTransformNumIntegration::integrate_
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:181
OpenMS::ContinuousWaveletTransform::wavelet_
std::vector< double > wavelet_
The pre-tabulated wavelet used for the transform.
Definition: ContinuousWaveletTransform.h:220
OpenMS::ContinuousWaveletTransform::signal_
std::vector< Peak1D > signal_
The transformed signal.
Definition: ContinuousWaveletTransform.h:217
OpenMS::Size
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:127
OpenMS::Math::round
T round(T x)
Rounds the value.
Definition: MathFunctions.h:138
OpenMS::ContinuousWaveletTransformNumIntegration::transform
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:94