36 #ifndef OPENMS_TRANSFORMATIONS_RAW2PEAK_CONTINUOUSWAVELETTRANSFORMNUMINTEGRATION_H
37 #define OPENMS_TRANSFORMATIONS_RAW2PEAK_CONTINUOUSWAVELETTRANSFORMNUMINTEGRATION_H
45 #ifdef DEBUG_PEAK_PICKING
94 template <
typename InputPeakIterator>
96 InputPeakIterator end_input,
98 unsigned int zeros = 0)
101 #ifdef DEBUG_PEAK_PICKING
102 std::cout <<
"ContinuousWaveletTransformNumIntegration::transform: start " << begin_input->getMZ() <<
" until " << (end_input - 1)->getMZ() << std::endl;
104 if (fabs(resolution - 1) < 0.0001)
107 SignedSize n = distance(begin_input, end_input);
114 #ifdef DEBUG_PEAK_PICKING
115 std::cout <<
"---------START TRANSFORM---------- \n";
117 InputPeakIterator help = begin_input;
118 for (
int i = 0; i < n; ++i)
120 signal_[i].setMZ(help->getMZ());
124 #ifdef DEBUG_PEAK_PICKING
125 std::cout <<
"---------END TRANSFORM----------" << std::endl;
128 begin_right_padding_ = n;
129 end_left_padding_ = -1;
134 double origin = begin_input->getMZ();
135 double spacing = ((end_input - 1)->getMZ() - origin) / (n - 1);
144 std::vector<double> processed_input(n);
148 InputPeakIterator it_help = begin_input;
151 processed_input[0] = it_help->getMZ() - zeros * spacing;
152 for (
unsigned int i = 0; i < zeros; ++i) processed_input[i] = 0;
154 else processed_input[0] = it_help->getIntensity();
159 x = origin +
k * spacing;
161 while (((it_help + 1) < end_input) && ((it_help + 1)->getMZ() < x))
165 processed_input[
k] = getInterpolatedValue_(x, it_help);
169 for (
unsigned int i = 0; i < zeros; ++i) processed_input[n - zeros + i] = 0;
173 for (
Int i = 0; i < n; ++i)
175 signal_[i].setMZ(origin + i * spacing);
181 begin_right_padding_ = n;
182 end_left_padding_ = -1;
186 begin_right_padding_ = n - zeros;
187 end_left_padding_ = zeros - 1;
202 virtual void init(
double scale,
double spacing);
207 template <
typename InputPeakIterator>
208 double integrate_(InputPeakIterator x, InputPeakIterator first, InputPeakIterator last)
210 #ifdef DEBUG_PEAK_PICKING
211 std::cout <<
"integrate_" << std::endl;
215 Size middle = wavelet_.size();
217 double start_pos = ((x->getMZ() - (middle * spacing_)) > first->getMZ()) ? (x->getMZ() - (middle * spacing_))
219 double end_pos = ((x->getMZ() + (middle * spacing_)) < (last - 1)->getMZ()) ? (x->getMZ() + (middle * spacing_))
220 : (last - 1)->getMZ();
222 InputPeakIterator help = x;
224 #ifdef DEBUG_PEAK_PICKING
225 std::cout <<
"integrate from middle to start_pos " << help->getMZ() <<
" until " << start_pos << std::endl;
229 while ((help != first) && ((help - 1)->getMZ() > start_pos))
232 double distance = fabs(x->getMZ() - help->getMZ());
234 if (index_w_r >= wavelet_.size())
236 index_w_r = wavelet_.size() - 1;
238 double wavelet_right = wavelet_[index_w_r];
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;
247 distance = fabs(x->getMZ() - (help - 1)->getMZ());
249 if (index_w_l >= wavelet_.size())
251 index_w_l = wavelet_.size() - 1;
253 double wavelet_left = wavelet_[index_w_l];
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;
262 std::cout <<
" intensity " << fabs((help - 1)->getMZ() - help->getMZ()) / 2. <<
" * " << (help - 1)->getIntensity() <<
" * " << wavelet_left <<
" + " << (help)->getIntensity() <<
"* " << wavelet_right
266 v += fabs((help - 1)->getMZ() - help->getMZ()) / 2. * ((help - 1)->getIntensity() * wavelet_left + help->getIntensity() * wavelet_right);
273 #ifdef DEBUG_PEAK_PICKING
274 std::cout <<
"integrate from middle to endpos " << (help)->getMZ() <<
" until " << end_pos << std::endl;
276 while ((help != (last - 1)) && ((help + 1)->getMZ() < end_pos))
279 double distance = fabs(x->getMZ() - help->getMZ());
281 if (index_w_l >= wavelet_.size())
283 index_w_l = wavelet_.size() - 1;
285 double wavelet_left = wavelet_[index_w_l];
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;
294 distance = fabs(x->getMZ() - (help + 1)->getMZ());
296 if (index_w_r >= wavelet_.size())
298 index_w_r = wavelet_.size() - 1;
300 double wavelet_right = wavelet_[index_w_r];
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;
308 v += fabs(help->getMZ() - (help + 1)->getMZ()) / 2. * (help->getIntensity() * wavelet_left + (help + 1)->getIntensity() * wavelet_right);
313 #ifdef DEBUG_PEAK_PICKING
314 std::cout <<
"return" << (v / sqrt(scale_)) << std::endl;
316 return v / sqrt(scale_);
320 double integrate_(
const std::vector<double> & processed_input,
double spacing_data,
int index);
325 return (1 - x * x) * exp(-x * x / 2);
ptrdiff_t SignedSize
Signed Size type e.g. used as pointer difference.
Definition: Types.h:128
T round(T x)
Rounds the value.
Definition: MathFunctions.h:137
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:121
int Int
Signed integer type.
Definition: Types.h:96