OpenMS
MorphologicalFilter.h
Go to the documentation of this file.
1 // Copyright (c) 2002-2023, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin
2 // SPDX-License-Identifier: BSD-3-Clause
3 //
4 // --------------------------------------------------------------------------
5 // $Maintainer: Timo Sachsenberg $
6 // $Authors: $
7 // --------------------------------------------------------------------------
8 
9 #pragma once
10 
17 #include <algorithm>
18 #include <iterator>
19 
20 namespace OpenMS
21 {
22 
23  namespace Internal
24  {
31  template<typename IteratorT>
32  class /* OPENMS_DLLAPI */ IntensityIteratorWrapper
33  {
34  public:
35  typedef std::forward_iterator_tag iterator_category;
36  typedef typename IteratorT::value_type::IntensityType value_type;
37  typedef typename IteratorT::value_type::IntensityType& reference;
38  typedef typename IteratorT::value_type::IntensityType* pointer;
39  typedef typename IteratorT::difference_type difference_type;
40 
41  IntensityIteratorWrapper(const IteratorT& rhs) : base(rhs)
42  {
43  }
44 
46  {
47  return base->getIntensity();
48  }
49 
50  template<typename IndexT>
51  value_type operator[](const IndexT& index)
52  {
53  return base[index].getIntensity();
54  }
55 
57  {
58  return base - rhs.base;
59  }
60 
62  {
63  ++base;
64  return *this;
65  }
66 
68  {
69  IteratorT tmp = *this;
70  ++(*this);
71  return tmp;
72  }
73 
74  bool operator==(const IntensityIteratorWrapper& rhs) const
75  {
76  return base == rhs.base;
77  }
78 
79  bool operator!=(const IntensityIteratorWrapper& rhs) const
80  {
81  return base != rhs.base;
82  }
83 
84  protected:
85  IteratorT base;
86  };
87 
89  template<typename IteratorT>
91  {
93  }
94 
95  } // namespace Internal
96 
132  class OPENMS_DLLAPI MorphologicalFilter : public ProgressLogger, public DefaultParamHandler
133  {
134  public:
136  MorphologicalFilter() : ProgressLogger(), DefaultParamHandler("MorphologicalFilter"), struct_size_in_datapoints_(0)
137  {
138  // structuring element
139  defaults_.setValue("struc_elem_length", 3.0, "Length of the structuring element. This should be wider than the expected peak width.");
140  defaults_.setValue("struc_elem_unit", "Thomson", "The unit of the 'struct_elem_length'.");
141  defaults_.setValidStrings("struc_elem_unit", {"Thomson", "DataPoints"});
142  // methods
143  defaults_.setValue("method", "tophat",
144  "Method to use, the default is 'tophat'. Do not change this unless you know what you are doing. The other methods may be useful for tuning the parameters, see the class "
145  "documentation of MorpthologicalFilter.");
146  defaults_.setValidStrings("method", {"identity", "erosion", "dilation", "opening", "closing", "gradient", "tophat", "bothat", "erosion_simple", "dilation_simple"});
147 
148  defaultsToParam_();
149  }
150 
153  {
154  }
155 
167  template<typename InputIterator, typename OutputIterator>
168  void filterRange(InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
169  {
170  // the buffer is static only to avoid reallocation
171  static std::vector<typename InputIterator::value_type> buffer;
172  const UInt size = input_end - input_begin;
173 
174  // determine the struct size in data points if not already set
175  if (struct_size_in_datapoints_ == 0)
176  {
177  struct_size_in_datapoints_ = (UInt)(double)param_.getValue("struc_elem_length");
178  }
179 
180  // apply the filtering
181  std::string method = param_.getValue("method");
182  if (method == "identity")
183  {
184  std::copy(input_begin, input_end, output_begin);
185  }
186  else if (method == "erosion")
187  {
188  applyErosion_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
189  }
190  else if (method == "dilation")
191  {
192  applyDilation_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
193  }
194  else if (method == "opening")
195  {
196  if (buffer.size() < size)
197  buffer.resize(size);
198  applyErosion_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
199  applyDilation_(struct_size_in_datapoints_, buffer.begin(), buffer.begin() + size, output_begin);
200  }
201  else if (method == "closing")
202  {
203  if (buffer.size() < size)
204  buffer.resize(size);
205  applyDilation_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
206  applyErosion_(struct_size_in_datapoints_, buffer.begin(), buffer.begin() + size, output_begin);
207  }
208  else if (method == "gradient")
209  {
210  if (buffer.size() < size)
211  buffer.resize(size);
212  applyErosion_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
213  applyDilation_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
214  for (UInt i = 0; i < size; ++i)
215  output_begin[i] -= buffer[i];
216  }
217  else if (method == "tophat")
218  {
219  if (buffer.size() < size)
220  buffer.resize(size);
221  applyErosion_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
222  applyDilation_(struct_size_in_datapoints_, buffer.begin(), buffer.begin() + size, output_begin);
223  for (UInt i = 0; i < size; ++i)
224  output_begin[i] = input_begin[i] - output_begin[i];
225  }
226  else if (method == "bothat")
227  {
228  if (buffer.size() < size)
229  buffer.resize(size);
230  applyDilation_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
231  applyErosion_(struct_size_in_datapoints_, buffer.begin(), buffer.begin() + size, output_begin);
232  for (UInt i = 0; i < size; ++i)
233  output_begin[i] = input_begin[i] - output_begin[i];
234  }
235  else if (method == "erosion_simple")
236  {
237  applyErosionSimple_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
238  }
239  else if (method == "dilation_simple")
240  {
241  applyDilationSimple_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
242  }
243 
244  struct_size_in_datapoints_ = 0;
245  }
246 
261  void filter(MSSpectrum& spectrum)
262  {
263  // make sure the right peak type is set
265 
266  // Abort if there is nothing to do
267  if (spectrum.size() <= 1)
268  {
269  return;
270  }
271 
272  // Determine structuring element size in datapoints (depending on the unit)
273  if (param_.getValue("struc_elem_unit") == "Thomson")
274  {
275  const double struc_elem_length = (double)param_.getValue("struc_elem_length");
276  const double mz_diff = spectrum.back().getMZ() - spectrum.begin()->getMZ();
277  struct_size_in_datapoints_ = (UInt)(ceil(struc_elem_length * (double)(spectrum.size() - 1) / mz_diff));
278  }
279  else
280  {
281  struct_size_in_datapoints_ = (UInt)(double)param_.getValue("struc_elem_length");
282  }
283  // make it odd (needed for the algorithm)
284  if (!Math::isOdd(struct_size_in_datapoints_))
285  ++struct_size_in_datapoints_;
286 
287  // apply the filtering and overwrite the input data
288  std::vector<Peak1D::IntensityType> output(spectrum.size());
289  filterRange(Internal::intensityIteratorWrapper(spectrum.begin()), Internal::intensityIteratorWrapper(spectrum.end()), output.begin());
290 
291  // overwrite output with data
292  for (Size i = 0; i < spectrum.size(); ++i)
293  {
294  spectrum[i].setIntensity(output[i]);
295  }
296  }
297 
305  {
306  startProgress(0, exp.size(), "filtering baseline");
307  for (UInt i = 0; i < exp.size(); ++i)
308  {
309  filter(exp[i]);
310  setProgress(i);
311  }
312  endProgress();
313  }
314 
315  protected:
318 
323  template<typename InputIterator, typename OutputIterator>
324  void applyErosion_(Int struc_size, InputIterator input, InputIterator input_end, OutputIterator output)
325  {
326  typedef typename InputIterator::value_type ValueType;
327  const Int size = input_end - input;
328  const Int struc_size_half = struc_size / 2; // yes, integer division
329 
330  static std::vector<ValueType> buffer;
331  if (Int(buffer.size()) < struc_size)
332  buffer.resize(struc_size);
333 
334  Int anchor; // anchoring position of the current block
335  Int i; // index relative to anchor, used for 'for' loops
336  Int ii = 0; // input index
337  Int oi = 0; // output index
338  ValueType current; // current value
339 
340  // we just can't get the case distinctions right in these cases, resorting to simple method.
341  if (size <= struc_size || size <= 5)
342  {
343  applyErosionSimple_(struc_size, input, input_end, output);
344  return;
345  }
346  {
347  // lower margin area
348  current = input[0];
349  for (++ii; ii < struc_size_half; ++ii)
350  if (current > input[ii])
351  current = input[ii];
352  for (; ii < std::min(Int(struc_size), size); ++ii, ++oi)
353  {
354  if (current > input[ii])
355  current = input[ii];
356  output[oi] = current;
357  }
358  }
359  {
360  // middle (main) area
361  for (anchor = struc_size; anchor <= size - struc_size; anchor += struc_size)
362  {
363  ii = anchor;
364  current = input[ii];
365  buffer[0] = current;
366  for (i = 1; i < struc_size; ++i, ++ii)
367  {
368  if (current > input[ii])
369  current = input[ii];
370  buffer[i] = current;
371  }
372  ii = anchor - 1;
373  oi = ii + struc_size_half;
374  current = input[ii];
375  for (i = 1; i < struc_size; ++i, --ii, --oi)
376  {
377  if (current > input[ii])
378  current = input[ii];
379  output[oi] = std::min(buffer[struc_size - i], current);
380  }
381  if (current > input[ii])
382  current = input[ii];
383  output[oi] = current;
384  }
385  }
386  {
387  // higher margin area
388  ii = size - 1;
389  oi = ii;
390  current = input[ii];
391  for (--ii; ii >= size - struc_size_half; --ii)
392  if (current > input[ii])
393  current = input[ii];
394  for (; ii >= std::max(size - Int(struc_size), 0); --ii, --oi)
395  {
396  if (current > input[ii])
397  current = input[ii];
398  output[oi] = current;
399  }
400  anchor = size - struc_size;
401  ii = anchor;
402  current = input[ii];
403  buffer[0] = current;
404  for (i = 1; i < struc_size; ++i, ++ii)
405  {
406  if (current > input[ii])
407  current = input[ii];
408  buffer[i] = current;
409  }
410  ii = anchor - 1;
411  oi = ii + struc_size_half;
412  current = input[ii];
413  for (i = 1; (ii >= 0) && (i < struc_size); ++i, --ii, --oi)
414  {
415  if (current > input[ii])
416  current = input[ii];
417  output[oi] = std::min(buffer[struc_size - i], current);
418  }
419  if (ii >= 0)
420  {
421  if (current > input[ii])
422  current = input[ii];
423  output[oi] = current;
424  }
425  }
426  return;
427  }
428 
433  template<typename InputIterator, typename OutputIterator>
434  void applyDilation_(Int struc_size, InputIterator input, InputIterator input_end, OutputIterator output)
435  {
436  typedef typename InputIterator::value_type ValueType;
437  const Int size = input_end - input;
438  const Int struc_size_half = struc_size / 2; // yes, integer division
439 
440  static std::vector<ValueType> buffer;
441  if (Int(buffer.size()) < struc_size)
442  buffer.resize(struc_size);
443 
444  Int anchor; // anchoring position of the current block
445  Int i; // index relative to anchor, used for 'for' loops
446  Int ii = 0; // input index
447  Int oi = 0; // output index
448  ValueType current; // current value
449 
450  // we just can't get the case distinctions right in these cases, resorting to simple method.
451  if (size <= struc_size || size <= 5)
452  {
453  applyDilationSimple_(struc_size, input, input_end, output);
454  return;
455  }
456  {
457  // lower margin area
458  current = input[0];
459  for (++ii; ii < struc_size_half; ++ii)
460  if (current < input[ii])
461  current = input[ii];
462  for (; ii < std::min(Int(struc_size), size); ++ii, ++oi)
463  {
464  if (current < input[ii])
465  current = input[ii];
466  output[oi] = current;
467  }
468  }
469  {
470  // middle (main) area
471  for (anchor = struc_size; anchor <= size - struc_size; anchor += struc_size)
472  {
473  ii = anchor;
474  current = input[ii];
475  buffer[0] = current;
476  for (i = 1; i < struc_size; ++i, ++ii)
477  {
478  if (current < input[ii])
479  current = input[ii];
480  buffer[i] = current;
481  }
482  ii = anchor - 1;
483  oi = ii + struc_size_half;
484  current = input[ii];
485  for (i = 1; i < struc_size; ++i, --ii, --oi)
486  {
487  if (current < input[ii])
488  current = input[ii];
489  output[oi] = std::max(buffer[struc_size - i], current);
490  }
491  if (current < input[ii])
492  current = input[ii];
493  output[oi] = current;
494  }
495  }
496  {
497  // higher margin area
498  ii = size - 1;
499  oi = ii;
500  current = input[ii];
501  for (--ii; ii >= size - struc_size_half; --ii)
502  if (current < input[ii])
503  current = input[ii];
504  for (; ii >= std::max(size - Int(struc_size), 0); --ii, --oi)
505  {
506  if (current < input[ii])
507  current = input[ii];
508  output[oi] = current;
509  }
510  anchor = size - struc_size;
511  ii = anchor;
512  current = input[ii];
513  buffer[0] = current;
514  for (i = 1; i < struc_size; ++i, ++ii)
515  {
516  if (current < input[ii])
517  current = input[ii];
518  buffer[i] = current;
519  }
520  ii = anchor - 1;
521  oi = ii + struc_size_half;
522  current = input[ii];
523  for (i = 1; (ii >= 0) && (i < struc_size); ++i, --ii, --oi)
524  {
525  if (current < input[ii])
526  current = input[ii];
527  output[oi] = std::max(buffer[struc_size - i], current);
528  }
529  if (ii >= 0)
530  {
531  if (current < input[ii])
532  current = input[ii];
533  output[oi] = current;
534  }
535  }
536  return;
537  }
538 
540  template<typename InputIterator, typename OutputIterator>
541  void applyErosionSimple_(Int struc_size, InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
542  {
543  typedef typename InputIterator::value_type ValueType;
544  const int size = input_end - input_begin;
545  const Int struc_size_half = struc_size / 2; // yes integer division
546  for (Int index = 0; index < size; ++index)
547  {
548  Int start = std::max(0, index - struc_size_half);
549  Int stop = std::min(size - 1, index + struc_size_half);
550  ValueType value = input_begin[start];
551  for (Int i = start + 1; i <= stop; ++i)
552  if (value > input_begin[i])
553  value = input_begin[i];
554  output_begin[index] = value;
555  }
556  return;
557  }
558 
560  template<typename InputIterator, typename OutputIterator>
561  void applyDilationSimple_(Int struc_size, InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
562  {
563  typedef typename InputIterator::value_type ValueType;
564  const int size = input_end - input_begin;
565  const Int struc_size_half = struc_size / 2; // yes integer division
566  for (Int index = 0; index < size; ++index)
567  {
568  Int start = std::max(0, index - struc_size_half);
569  Int stop = std::min(size - 1, index + struc_size_half);
570  ValueType value = input_begin[start];
571  for (Int i = start + 1; i <= stop; ++i)
572  if (value < input_begin[i])
573  value = input_begin[i];
574  output_begin[index] = value;
575  }
576  return;
577  }
578 
579  private:
582  };
583 
584 } // namespace OpenMS
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:66
An iterator wrapper to access peak intensities instead of the peak itself.
Definition: MorphologicalFilter.h:33
bool operator==(const IntensityIteratorWrapper &rhs) const
Definition: MorphologicalFilter.h:74
IteratorT base
Definition: MorphologicalFilter.h:85
std::forward_iterator_tag iterator_category
Definition: MorphologicalFilter.h:35
bool operator!=(const IntensityIteratorWrapper &rhs) const
Definition: MorphologicalFilter.h:79
IteratorT::difference_type difference_type
Definition: MorphologicalFilter.h:39
IteratorT::value_type::IntensityType value_type
Definition: MorphologicalFilter.h:36
IteratorT::value_type::IntensityType * pointer
Definition: MorphologicalFilter.h:38
difference_type operator-(IntensityIteratorWrapper &rhs) const
Definition: MorphologicalFilter.h:56
value_type operator*()
Definition: MorphologicalFilter.h:45
IntensityIteratorWrapper operator++(int)
Definition: MorphologicalFilter.h:67
IntensityIteratorWrapper & operator++()
Definition: MorphologicalFilter.h:61
IteratorT::value_type::IntensityType & reference
Definition: MorphologicalFilter.h:37
value_type operator[](const IndexT &index)
Definition: MorphologicalFilter.h:51
IntensityIteratorWrapper(const IteratorT &rhs)
Definition: MorphologicalFilter.h:41
In-Memory representation of a mass spectrometry run.
Definition: MSExperiment.h:46
Size size() const
The number of spectra.
Definition: MSExperiment.h:121
The representation of a 1D spectrum.
Definition: MSSpectrum.h:44
This class implements baseline filtering operations using methods from mathematical morphology.
Definition: MorphologicalFilter.h:133
void applyDilation_(Int struc_size, InputIterator input, InputIterator input_end, OutputIterator output)
Applies dilation. This implementation uses van Herk's method. Only 3 min/max comparisons are required...
Definition: MorphologicalFilter.h:434
MorphologicalFilter(const MorphologicalFilter &source)
copy constructor not implemented
void filterRange(InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
Applies the morphological filtering operation to an iterator range.
Definition: MorphologicalFilter.h:168
void filter(MSSpectrum &spectrum)
Applies the morphological filtering operation to an MSSpectrum.
Definition: MorphologicalFilter.h:261
void filterExperiment(PeakMap &exp)
Applies the morphological filtering operation to an MSExperiment.
Definition: MorphologicalFilter.h:304
UInt struct_size_in_datapoints_
Member for struct size in data points.
Definition: MorphologicalFilter.h:317
MorphologicalFilter()
Constructor.
Definition: MorphologicalFilter.h:136
void applyDilationSimple_(Int struc_size, InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
Applies dilation. Simple implementation, possibly faster if struc_size is very small,...
Definition: MorphologicalFilter.h:561
void applyErosion_(Int struc_size, InputIterator input, InputIterator input_end, OutputIterator output)
Applies erosion. This implementation uses van Herk's method. Only 3 min/max comparisons are required ...
Definition: MorphologicalFilter.h:324
void applyErosionSimple_(Int struc_size, InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
Applies erosion. Simple implementation, possibly faster if struc_size is very small,...
Definition: MorphologicalFilter.h:541
~MorphologicalFilter() override
Destructor.
Definition: MorphologicalFilter.h:152
Base class for all classes that want to report their progress.
Definition: ProgressLogger.h:27
void setType(SpectrumType type)
sets the spectrum type
@ PROFILE
profile data
Definition: SpectrumSettings.h:48
int Int
Signed integer type.
Definition: Types.h:76
unsigned int UInt
Unsigned integer type.
Definition: Types.h:68
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:101
bool isOdd(UInt x)
Returns true if the given integer is odd.
Definition: MathFunctions.h:173
IntensityIteratorWrapper< IteratorT > intensityIteratorWrapper(const IteratorT &rhs)
make-function so that we need no write out all those type names to get the wrapped iterator.
Definition: MorphologicalFilter.h:90
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:22