| OpenMS 2.5.0

LinearResamplerAlign.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: Hannes Roest $
32 // $Authors: Hannes Roest $
33 // --------------------------------------------------------------------------
34
35 #pragma once
36
37 #include <OpenMS/FILTERING/TRANSFORMERS/LinearResampler.h>
38 #include <OpenMS/CONCEPT/Macros.h>
39
40 namespace OpenMS
41 {
42
57 class OPENMS_DLLAPI LinearResamplerAlign :
58 public LinearResampler
59 {
60
61 public:
62
63 LinearResamplerAlign()
64 {
65 defaults_.setValue("spacing", 0.05, "Spacing of the resampled output peaks.");
66 defaults_.setValue("ppm", "false", "Whether spacing is in ppm or Th");
67 defaultsToParam_();
68 }
69
79 template <class SpecT>
80 void raster(SpecT& container)
81 {
82 //return if nothing to do
83 if (container.empty()) return;
84
85 typename SpecT::iterator first = container.begin();
86 typename SpecT::iterator last = container.end();
87
88 double end_pos = (last - 1)->getMZ();
89 double start_pos = first->getMZ();
90 int number_resampled_points = (int)(ceil((end_pos - start_pos) / spacing_ + 1));
91
92 std::vector<typename SpecT::PeakType> resampled_peak_container;
93 populate_raster_(resampled_peak_container, start_pos, end_pos, number_resampled_points);
94
95 raster(container.begin(), container.end(), resampled_peak_container.begin(), resampled_peak_container.end());
96
97 container.swap(resampled_peak_container);
98 }
99
115 template <typename SpecT>
116 void raster_align(SpecT& container, double start_pos, double end_pos)
117 {
118 //return if nothing to do
119 if (container.empty()) return;
120
121 if (end_pos < start_pos)
122 {
123 std::vector<typename SpecT::PeakType> empty;
124 container.swap(empty);
125 return;
126 }
127
128 typename SpecT::iterator first = container.begin();
129 typename SpecT::iterator last = container.end();
130
131 // get the iterators just before / after the two points start_pos / end_pos
132 while (first != container.end() && (first)->getMZ() < start_pos) {++first;}
133 while (last != first && (last - 1)->getMZ() > end_pos) {--last;}
134
135 int number_resampled_points = (int)(ceil((end_pos - start_pos) / spacing_ + 1));
136
137 std::vector<typename SpecT::PeakType> resampled_peak_container;
138 populate_raster_(resampled_peak_container, start_pos, end_pos, number_resampled_points);
139
140 raster(first, last, resampled_peak_container.begin(), resampled_peak_container.end());
141
142 container.swap(resampled_peak_container);
143 }
144
166 template <typename PeakTypeIterator, typename ConstPeakTypeIterator>
167 void raster(ConstPeakTypeIterator raw_it, ConstPeakTypeIterator raw_end, PeakTypeIterator resample_it, PeakTypeIterator resample_end)
168 {
169 OPENMS_PRECONDITION(resample_it != resample_end, "Output iterators cannot be identical") // as we use +1
170 // OPENMS_PRECONDITION(raw_it != raw_end, "Input iterators cannot be identical")
171
172 PeakTypeIterator resample_start = resample_it;
173
174 // need to get the raw iterator between two resampled iterators of the raw data
175 while (raw_it != raw_end && raw_it->getMZ() < resample_it->getMZ())
176 {
177 resample_it->setIntensity(resample_it->getIntensity() + raw_it->getIntensity());
178 raw_it++;
179 }
180
181 while (raw_it != raw_end)
182 {
183 //advance the resample iterator until our raw point is between two resampled iterators
184 while (resample_it != resample_end && resample_it->getMZ() < raw_it->getMZ()) {resample_it++;}
185 if (resample_it != resample_start) {resample_it--;}
186
187 // if we have the last datapoint we break
188 if ((resample_it + 1) == resample_end) {break;}
189
190 double dist_left = fabs(raw_it->getMZ() - resample_it->getMZ());
191 double dist_right = fabs(raw_it->getMZ() - (resample_it + 1)->getMZ());
192
193 // distribute the intensity of the raw point according to the distance to resample_it and resample_it+1
194 resample_it->setIntensity(resample_it->getIntensity() + raw_it->getIntensity() * dist_right / (dist_left + dist_right));
195 (resample_it + 1)->setIntensity((resample_it + 1)->getIntensity() + raw_it->getIntensity() * dist_left / (dist_left + dist_right));
196
197 raw_it++;
198 }
199
200 // add the final intensity to the right
201 while (raw_it != raw_end)
202 {
203 resample_it->setIntensity(resample_it->getIntensity() + raw_it->getIntensity());
204 raw_it++;
205 }
206 }
207
233 template <typename PeakTypeIterator, typename ConstPeakTypeIterator>
234 #ifdef OPENMS_ASSERTIONS
235 void raster(ConstPeakTypeIterator mz_raw_it, ConstPeakTypeIterator mz_raw_end,
236 ConstPeakTypeIterator int_raw_it, ConstPeakTypeIterator int_raw_end,
237 PeakTypeIterator mz_resample_it, PeakTypeIterator mz_resample_end,
238 PeakTypeIterator int_resample_it, PeakTypeIterator int_resample_end)
239 #else
240 void raster(ConstPeakTypeIterator mz_raw_it, ConstPeakTypeIterator mz_raw_end,
241 ConstPeakTypeIterator int_raw_it, ConstPeakTypeIterator /* int_raw_end */,
242 PeakTypeIterator mz_resample_it, PeakTypeIterator mz_resample_end,
243 PeakTypeIterator int_resample_it, PeakTypeIterator /* int_resample_end */)
244 #endif
245 {
246 OPENMS_PRECONDITION(mz_resample_it != mz_resample_end, "Output iterators cannot be identical") // as we use +1
247 OPENMS_PRECONDITION(std::distance(mz_resample_it, mz_resample_end) == std::distance(int_resample_it, int_resample_end),
248 "Resample m/z and intensity iterators need to cover the same distance")
249 OPENMS_PRECONDITION(std::distance(mz_raw_it, mz_raw_end) == std::distance(int_raw_it, int_raw_end),
250 "Raw m/z and intensity iterators need to cover the same distance")
251 // OPENMS_PRECONDITION(raw_it != raw_end, "Input iterators cannot be identical")
252
253 PeakTypeIterator mz_resample_start = mz_resample_it;
254
255 // need to get the raw iterator between two resampled iterators of the raw data
256 while (mz_raw_it != mz_raw_end && (*mz_raw_it) < (*mz_resample_it) )
257 {
258 (*int_resample_it) = *int_resample_it + *int_raw_it;
259 ++mz_raw_it;
260 ++int_raw_it;
261 }
262
263 while (mz_raw_it != mz_raw_end)
264 {
265 //advance the resample iterator until our raw point is between two resampled iterators
266 while (mz_resample_it != mz_resample_end && *mz_resample_it < *mz_raw_it)
267 {
268 ++mz_resample_it; ++int_resample_it;
269 }
270 if (mz_resample_it != mz_resample_start)
271 {
272 --mz_resample_it; --int_resample_it;
273 }
274
275 // if we have the last datapoint we break
276 if ((mz_resample_it + 1) == mz_resample_end) {break;}
277
278 double dist_left = fabs(*mz_raw_it - *mz_resample_it);
279 double dist_right = fabs(*mz_raw_it - *(mz_resample_it + 1));
280
281 // distribute the intensity of the raw point according to the distance to resample_it and resample_it+1
282 *(int_resample_it) = *int_resample_it + (*int_raw_it) * dist_right / (dist_left + dist_right);
283 *(int_resample_it + 1) = *(int_resample_it + 1) + (*int_raw_it) * dist_left / (dist_left + dist_right);
284
285 ++mz_raw_it;
286 ++int_raw_it;
287 }
288
289 // add the final intensity to the right
290 while (mz_raw_it != mz_raw_end)
291 {
292 *int_resample_it = *int_resample_it + (*int_raw_it);
293 ++mz_raw_it;
294 ++int_raw_it;
295 }
296 }
297
314 template <typename PeakTypeIterator>
315 void raster_interpolate(PeakTypeIterator raw_it, PeakTypeIterator raw_end, PeakTypeIterator resample_it, PeakTypeIterator resampled_end)
316 {
317 // OPENMS_PRECONDITION(resample_it != resampled_end, "Output iterators cannot be identical")
318 OPENMS_PRECONDITION(raw_it != raw_end, "Input iterators cannot be identical") // as we use +1
319
320 PeakTypeIterator raw_start = raw_it;
321
322 // need to get the resampled iterator between two iterators of the raw data
323 while (resample_it != resampled_end && resample_it->getMZ() < raw_it->getMZ()) {resample_it++;}
324
325 while (resample_it != resampled_end)
326 {
327 //advance the raw_iterator until our current point we want to interpolate is between them
328 while (raw_it != raw_end && raw_it->getMZ() < resample_it->getMZ()) {raw_it++;}
329 if (raw_it != raw_start) {raw_it--;}
330
331 // if we have the last datapoint we break
332 if ((raw_it + 1) == raw_end) {break;}
333
334 // use a linear interpolation between raw_it and raw_it+1
335 double m = ((raw_it + 1)->getIntensity() - raw_it->getIntensity()) / ((raw_it + 1)->getMZ() - raw_it->getMZ());
336 resample_it->setIntensity(raw_it->getIntensity() + (resample_it->getMZ() - raw_it->getMZ()) * m);
337 resample_it++;
338 }
339
340 }
341
342 protected:
343
345 bool ppm_;
346
347 void updateMembers_() override
348 {
349 spacing_ = param_.getValue("spacing");
350 ppm_ = (bool)param_.getValue("ppm").toBool();
351 }
352
354 template <typename PeakType>
355 void populate_raster_(std::vector<PeakType>& resampled_peak_container,
356 double start_pos, double end_pos, int number_resampled_points)
357 {
358 if (!ppm_)
359 {
360 // generate the resampled peaks at positions origin+i*spacing_
361 resampled_peak_container.resize(number_resampled_points);
362 typename std::vector<PeakType>::iterator it = resampled_peak_container.begin();
363 for (int i = 0; i < number_resampled_points; ++i)
364 {
365 it->setMZ(start_pos + i * spacing_);
366 ++it;
367 }
368 }
369 else
370 {
371 // generate resampled peaks with ppm distance (not fixed)
372 double current_mz = start_pos;
373 while (current_mz < end_pos)
374 {
375 PeakType p;
376 p.setIntensity(0);
377 p.setMZ(current_mz);
378 resampled_peak_container.push_back(p);
379
380 // increment current_mz
381 current_mz += current_mz * (spacing_ / 1e6);
382 }
383 }
384 }
385 };
386
387 }
388
389
OpenMS::LinearResamplerAlign::raster_align
void raster_align(SpecT &container, double start_pos, double end_pos)
Applies the resampling algorithm to a container (MSSpectrum or MSChromatogram) with fixed coordinates...
Definition: LinearResamplerAlign.h:116
OpenMS::Peak2D
A 2-dimensional raw data point or peak.
Definition: Peak2D.h:54
OpenMS::Peak2D::setMZ
void setMZ(CoordinateType coordinate)
Mutable access to the m/z coordinate (index 1)
Definition: Peak2D.h:202
OpenMS::LinearResamplerAlign::populate_raster_
void populate_raster_(std::vector< PeakType > &resampled_peak_container, double start_pos, double end_pos, int number_resampled_points)
Generate raster for resampled peak container.
Definition: LinearResamplerAlign.h:355
OpenMS::LinearResamplerAlign::ppm_
bool ppm_
Spacing of the resampled data.
Definition: LinearResamplerAlign.h:345
OpenMS::LinearResamplerAlign::raster
void raster(ConstPeakTypeIterator raw_it, ConstPeakTypeIterator raw_end, PeakTypeIterator resample_it, PeakTypeIterator resample_end)
Applies the resampling algorithm.
Definition: LinearResamplerAlign.h:167
OpenMS::LinearResamplerAlign
Linear Resampling of raw data with alignment.
Definition: LinearResamplerAlign.h:57
OpenMS::LinearResamplerAlign::updateMembers_
void updateMembers_() override
This method is used to update extra member variables at the end of the setParameters() method...
Definition: LinearResamplerAlign.h:347
OpenMS::LinearResamplerAlign::raster
void raster(ConstPeakTypeIterator mz_raw_it, ConstPeakTypeIterator mz_raw_end, ConstPeakTypeIterator int_raw_it, ConstPeakTypeIterator, PeakTypeIterator mz_resample_it, PeakTypeIterator mz_resample_end, PeakTypeIterator int_resample_it, PeakTypeIterator)
Applies the resampling algorithm.
Definition: LinearResamplerAlign.h:240
OpenMS::LinearResamplerAlign::raster_interpolate
void raster_interpolate(PeakTypeIterator raw_it, PeakTypeIterator raw_end, PeakTypeIterator resample_it, PeakTypeIterator resampled_end)
Applies the resampling algorithm using a linear interpolation.
Definition: LinearResamplerAlign.h:315
OPENMS_PRECONDITION
#define OPENMS_PRECONDITION(condition, message)
Precondition macro.
Definition: openms/include/OpenMS/CONCEPT/Macros.h:136
OpenMS::LinearResamplerAlign::raster
void raster(SpecT &container)
Applies the resampling algorithm to a container (MSSpectrum or MSChromatogram).
Definition: LinearResamplerAlign.h:80
OpenMS
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:46
OpenMS::LinearResamplerAlign::LinearResamplerAlign
LinearResamplerAlign()
Definition: LinearResamplerAlign.h:63
LinearResampler.h
OpenMS::LinearResampler
Linear Resampling of raw data.
Definition: LinearResampler.h:61
int
Macros.h
OpenMS::Peak2D::setIntensity
void setIntensity(IntensityType intensity)
Non-mutable access to the data point intensity (height)
Definition: Peak2D.h:172