OpenMS  2.4.0
MorphologicalFilter.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-2018.
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: $
33 // --------------------------------------------------------------------------
34 
35 #pragma once
36 
43 
44 #include <algorithm>
45 #include <iterator>
46 
47 namespace OpenMS
48 {
49 
50  namespace Internal
51  {
58  template <typename IteratorT>
59  class /* OPENMS_DLLAPI */ IntensityIteratorWrapper :
60  public std::iterator<std::forward_iterator_tag, typename IteratorT::value_type::IntensityType>
61  {
62 public:
63  typedef typename IteratorT::value_type::IntensityType value_type;
64  typedef typename IteratorT::value_type::IntensityType & reference;
65  typedef typename IteratorT::value_type::IntensityType * pointer;
66  typedef typename IteratorT::difference_type difference_type;
67 
68  IntensityIteratorWrapper(const IteratorT & rhs) :
69  base(rhs)
70  {
71  }
72 
74  {
75  return base->getIntensity();
76  }
77 
78  template <typename IndexT>
79  value_type operator[](const IndexT & index)
80  {
81  return base[index].getIntensity();
82  }
83 
85  {
86  return base - rhs.base;
87  }
88 
90  {
91  ++base;
92  return *this;
93  }
94 
96  {
97  IteratorT tmp = *this;
98  ++(*this);
99  return tmp;
100  }
101 
102  bool operator==(const IntensityIteratorWrapper & rhs) const
103  {
104  return base == rhs.base;
105  }
106 
107  bool operator!=(const IntensityIteratorWrapper & rhs) const
108  {
109  return base != rhs.base;
110  }
111 
112 protected:
113  IteratorT base;
114  };
115 
117  template <typename IteratorT>
119  {
121  }
122 
123  }
124 
160  class OPENMS_DLLAPI MorphologicalFilter :
161  public ProgressLogger,
162  public DefaultParamHandler
163  {
164 public:
165 
168  ProgressLogger(),
169  DefaultParamHandler("MorphologicalFilter"),
170  struct_size_in_datapoints_(0)
171  {
172  //structuring element
173  defaults_.setValue("struc_elem_length", 3.0, "Length of the structuring element. This should be wider than the expected peak width.");
174  defaults_.setValue("struc_elem_unit", "Thomson", "The unit of the 'struct_elem_length'.");
175  defaults_.setValidStrings("struc_elem_unit", ListUtils::create<String>("Thomson,DataPoints"));
176  //methods
177  defaults_.setValue("method", "tophat", "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 documentation of MorpthologicalFilter.");
178  defaults_.setValidStrings("method", ListUtils::create<String>("identity,erosion,dilation,opening,closing,gradient,tophat,bothat,erosion_simple,dilation_simple"));
179 
180  defaultsToParam_();
181  }
182 
185  {
186  }
187 
199  template <typename InputIterator, typename OutputIterator>
200  void filterRange(InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
201  {
202  // the buffer is static only to avoid reallocation
203  static std::vector<typename InputIterator::value_type> buffer;
204  const UInt size = input_end - input_begin;
205 
206  //determine the struct size in data points if not already set
207  if (struct_size_in_datapoints_ == 0)
208  {
209  struct_size_in_datapoints_ = (UInt)(double)param_.getValue("struc_elem_length");
210  }
211 
212  //apply the filtering
213  String method = param_.getValue("method");
214  if (method == "identity")
215  {
216  std::copy(input_begin, input_end, output_begin);
217  }
218  else if (method == "erosion")
219  {
220  applyErosion_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
221  }
222  else if (method == "dilation")
223  {
224  applyDilation_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
225  }
226  else if (method == "opening")
227  {
228  if (buffer.size() < size) buffer.resize(size);
229  applyErosion_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
230  applyDilation_(struct_size_in_datapoints_, buffer.begin(), buffer.begin() + size, output_begin);
231  }
232  else if (method == "closing")
233  {
234  if (buffer.size() < size) buffer.resize(size);
235  applyDilation_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
236  applyErosion_(struct_size_in_datapoints_, buffer.begin(), buffer.begin() + size, output_begin);
237  }
238  else if (method == "gradient")
239  {
240  if (buffer.size() < size) buffer.resize(size);
241  applyErosion_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
242  applyDilation_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
243  for (UInt i = 0; i < size; ++i) output_begin[i] -= buffer[i];
244  }
245  else if (method == "tophat")
246  {
247  if (buffer.size() < size) buffer.resize(size);
248  applyErosion_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
249  applyDilation_(struct_size_in_datapoints_, buffer.begin(), buffer.begin() + size, output_begin);
250  for (UInt i = 0; i < size; ++i) output_begin[i] = input_begin[i] - output_begin[i];
251  }
252  else if (method == "bothat")
253  {
254  if (buffer.size() < size) buffer.resize(size);
255  applyDilation_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
256  applyErosion_(struct_size_in_datapoints_, buffer.begin(), buffer.begin() + size, output_begin);
257  for (UInt i = 0; i < size; ++i) output_begin[i] = input_begin[i] - output_begin[i];
258  }
259  else if (method == "erosion_simple")
260  {
261  applyErosionSimple_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
262  }
263  else if (method == "dilation_simple")
264  {
265  applyDilationSimple_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
266  }
267 
268  struct_size_in_datapoints_ = 0;
269  }
270 
285  void filter(MSSpectrum & spectrum)
286  {
287  //make sure the right peak type is set
289 
290  //Abort if there is nothing to do
291  if (spectrum.size() <= 1) { return; }
292 
293  //Determine structuring element size in datapoints (depending on the unit)
294  if ((String)(param_.getValue("struc_elem_unit")) == "Thomson")
295  {
296  const double struc_elem_length = (double)param_.getValue("struc_elem_length");
297  const double mz_diff = spectrum.back().getMZ() - spectrum.begin()->getMZ();
298  struct_size_in_datapoints_ = (UInt)(ceil(struc_elem_length*(double)(spectrum.size() - 1)/mz_diff));
299  }
300  else
301  {
302  struct_size_in_datapoints_ = (UInt)(double)param_.getValue("struc_elem_length");
303  }
304  //make it odd (needed for the algorithm)
305  if (!Math::isOdd(struct_size_in_datapoints_)) ++struct_size_in_datapoints_;
306 
307  //apply the filtering and overwrite the input data
308  std::vector<Peak1D::IntensityType> output(spectrum.size());
309  filterRange(Internal::intensityIteratorWrapper(spectrum.begin()),
310  Internal::intensityIteratorWrapper(spectrum.end()),
311  output.begin()
312  );
313 
314  //overwrite output with data
315  for (Size i = 0; i < spectrum.size(); ++i)
316  {
317  spectrum[i].setIntensity(output[i]);
318  }
319  }
320 
328  {
329  startProgress(0, exp.size(), "filtering baseline");
330  for (UInt i = 0; i < exp.size(); ++i)
331  {
332  filter(exp[i]);
333  setProgress(i);
334  }
335  endProgress();
336  }
337 
338 protected:
339 
342 
347  template <typename InputIterator, typename OutputIterator>
348  void applyErosion_(Int struc_size, InputIterator input, InputIterator input_end, OutputIterator output)
349  {
350  typedef typename InputIterator::value_type ValueType;
351  const Int size = input_end - input;
352  const Int struc_size_half = struc_size / 2; // yes, integer division
353 
354  static std::vector<ValueType> buffer;
355  if (Int(buffer.size()) < struc_size) buffer.resize(struc_size);
356 
357  Int anchor; // anchoring position of the current block
358  Int i; // index relative to anchor, used for 'for' loops
359  Int ii = 0; // input index
360  Int oi = 0; // output index
361  ValueType current; // current value
362 
363  // we just can't get the case distinctions right in these cases, resorting to simple method.
364  if (size <= struc_size || size <= 5)
365  {
366  applyErosionSimple_(struc_size, input, input_end, output);
367  return;
368  }
369  {
370  // lower margin area
371  current = input[0];
372  for (++ii; ii < struc_size_half; ++ii) if (current > input[ii]) current = input[ii];
373  for (; ii < std::min(Int(struc_size), size); ++ii, ++oi)
374  {
375  if (current > input[ii]) current = input[ii];
376  output[oi] = current;
377  }
378  }
379  {
380  // middle (main) area
381  for (anchor = struc_size;
382  anchor <= size - struc_size;
383  anchor += struc_size
384  )
385  {
386  ii = anchor;
387  current = input[ii];
388  buffer[0] = current;
389  for (i = 1; i < struc_size; ++i, ++ii)
390  {
391  if (current > input[ii]) current = input[ii];
392  buffer[i] = current;
393  }
394  ii = anchor - 1;
395  oi = ii + struc_size_half;
396  current = input[ii];
397  for (i = 1; i < struc_size; ++i, --ii, --oi)
398  {
399  if (current > input[ii]) current = input[ii];
400  output[oi] = std::min(buffer[struc_size - i], current);
401  }
402  if (current > input[ii]) current = input[ii];
403  output[oi] = current;
404  }
405  }
406  {
407  // higher margin area
408  ii = size - 1;
409  oi = ii;
410  current = input[ii];
411  for (--ii; ii >= size - struc_size_half; --ii) if (current > input[ii]) current = input[ii];
412  for (; ii >= std::max(size - Int(struc_size), 0); --ii, --oi)
413  {
414  if (current > input[ii]) current = input[ii];
415  output[oi] = current;
416  }
417  anchor = size - struc_size;
418  ii = anchor;
419  current = input[ii];
420  buffer[0] = current;
421  for (i = 1; i < struc_size; ++i, ++ii)
422  {
423  if (current > input[ii]) current = input[ii];
424  buffer[i] = current;
425  }
426  ii = anchor - 1;
427  oi = ii + struc_size_half;
428  current = input[ii];
429  for (i = 1; (ii >= 0) && (i < struc_size); ++i, --ii, --oi)
430  {
431  if (current > input[ii]) current = input[ii];
432  output[oi] = std::min(buffer[struc_size - i], current);
433  }
434  if (ii >= 0)
435  {
436  if (current > input[ii]) current = input[ii];
437  output[oi] = current;
438  }
439  }
440  return;
441  }
442 
447  template <typename InputIterator, typename OutputIterator>
448  void applyDilation_(Int struc_size, InputIterator input, InputIterator input_end, OutputIterator output)
449  {
450  typedef typename InputIterator::value_type ValueType;
451  const Int size = input_end - input;
452  const Int struc_size_half = struc_size / 2; // yes, integer division
453 
454  static std::vector<ValueType> buffer;
455  if (Int(buffer.size()) < struc_size) buffer.resize(struc_size);
456 
457  Int anchor; // anchoring position of the current block
458  Int i; // index relative to anchor, used for 'for' loops
459  Int ii = 0; // input index
460  Int oi = 0; // output index
461  ValueType current; // current value
462 
463  // we just can't get the case distinctions right in these cases, resorting to simple method.
464  if (size <= struc_size || size <= 5)
465  {
466  applyDilationSimple_(struc_size, input, input_end, output);
467  return;
468  }
469  {
470  // lower margin area
471  current = input[0];
472  for (++ii; ii < struc_size_half; ++ii) if (current < input[ii]) current = input[ii];
473  for (; ii < std::min(Int(struc_size), size); ++ii, ++oi)
474  {
475  if (current < input[ii]) current = input[ii];
476  output[oi] = current;
477  }
478  }
479  {
480  // middle (main) area
481  for (anchor = struc_size;
482  anchor <= size - struc_size;
483  anchor += struc_size
484  )
485  {
486  ii = anchor;
487  current = input[ii];
488  buffer[0] = current;
489  for (i = 1; i < struc_size; ++i, ++ii)
490  {
491  if (current < input[ii]) current = input[ii];
492  buffer[i] = current;
493  }
494  ii = anchor - 1;
495  oi = ii + struc_size_half;
496  current = input[ii];
497  for (i = 1; i < struc_size; ++i, --ii, --oi)
498  {
499  if (current < input[ii]) current = input[ii];
500  output[oi] = std::max(buffer[struc_size - i], current);
501  }
502  if (current < input[ii]) current = input[ii];
503  output[oi] = current;
504  }
505  }
506  {
507  // higher margin area
508  ii = size - 1;
509  oi = ii;
510  current = input[ii];
511  for (--ii; ii >= size - struc_size_half; --ii) if (current < input[ii]) current = input[ii];
512  for (; ii >= std::max(size - Int(struc_size), 0); --ii, --oi)
513  {
514  if (current < input[ii]) current = input[ii];
515  output[oi] = current;
516  }
517  anchor = size - struc_size;
518  ii = anchor;
519  current = input[ii];
520  buffer[0] = current;
521  for (i = 1; i < struc_size; ++i, ++ii)
522  {
523  if (current < input[ii]) current = input[ii];
524  buffer[i] = current;
525  }
526  ii = anchor - 1;
527  oi = ii + struc_size_half;
528  current = input[ii];
529  for (i = 1; (ii >= 0) && (i < struc_size); ++i, --ii, --oi)
530  {
531  if (current < input[ii]) current = input[ii];
532  output[oi] = std::max(buffer[struc_size - i], current);
533  }
534  if (ii >= 0)
535  {
536  if (current < input[ii]) current = input[ii];
537  output[oi] = current;
538  }
539  }
540  return;
541  }
542 
544  template <typename InputIterator, typename OutputIterator>
545  void applyErosionSimple_(Int struc_size, InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
546  {
547  typedef typename InputIterator::value_type ValueType;
548  const int size = input_end - input_begin;
549  const Int struc_size_half = struc_size / 2; // yes integer division
550  for (Int index = 0; index < size; ++index)
551  {
552  Int start = std::max(0, index - struc_size_half);
553  Int stop = std::min(size - 1, index + struc_size_half);
554  ValueType value = input_begin[start];
555  for (Int i = start + 1; i <= stop; ++i) if (value > input_begin[i]) value = input_begin[i];
556  output_begin[index] = value;
557  }
558  return;
559  }
560 
562  template <typename InputIterator, typename OutputIterator>
563  void applyDilationSimple_(Int struc_size, InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
564  {
565  typedef typename InputIterator::value_type ValueType;
566  const int size = input_end - input_begin;
567  const Int struc_size_half = struc_size / 2; // yes integer division
568  for (Int index = 0; index < size; ++index)
569  {
570  Int start = std::max(0, index - struc_size_half);
571  Int stop = std::min(size - 1, index + struc_size_half);
572  ValueType value = input_begin[start];
573  for (Int i = start + 1; i <= stop; ++i) if (value < input_begin[i]) value = input_begin[i];
574  output_begin[index] = value;
575  }
576  return;
577  }
578 
579 private:
580 
583 
584  };
585 
586 } // namespace OpenMS
587 
UInt struct_size_in_datapoints_
Member for struct size in data points.
Definition: MorphologicalFilter.h:341
This class implements baseline filtering operations using methods from mathematical morphology...
Definition: MorphologicalFilter.h:160
bool isOdd(UInt x)
Returns true if the given integer is odd.
Definition: MathFunctions.h:127
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:118
void applyDilation_(Int struc_size, InputIterator input, InputIterator input_end, OutputIterator output)
Applies dilation. This implementation uses van Herk&#39;s method. Only 3 min/max comparisons are required...
Definition: MorphologicalFilter.h:448
IntensityIteratorWrapper(const IteratorT &rhs)
Definition: MorphologicalFilter.h:68
An iterator wrapper to access peak intensities instead of the peak itself.
Definition: MorphologicalFilter.h:59
value_type operator[](const IndexT &index)
Definition: MorphologicalFilter.h:79
A more convenient string class.
Definition: String.h:57
IteratorT::difference_type difference_type
Definition: MorphologicalFilter.h:66
MorphologicalFilter()
Constructor.
Definition: MorphologicalFilter.h:167
void filterRange(InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
Applies the morphological filtering operation to an iterator range.
Definition: MorphologicalFilter.h:200
IntensityIteratorWrapper & operator++()
Definition: MorphologicalFilter.h:89
IteratorT base
Definition: MorphologicalFilter.h:113
unsigned int UInt
Unsigned integer type.
Definition: Types.h:94
profile data
Definition: SpectrumSettings.h:74
~MorphologicalFilter() override
Destructor.
Definition: MorphologicalFilter.h:184
bool operator==(const IntensityIteratorWrapper &rhs) const
Definition: MorphologicalFilter.h:102
Size size() const
Definition: MSExperiment.h:127
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:46
The representation of a 1D spectrum.
Definition: MSSpectrum.h:66
difference_type operator-(IntensityIteratorWrapper &rhs) const
Definition: MorphologicalFilter.h:84
void filter(MSSpectrum &spectrum)
Applies the morphological filtering operation to an MSSpectrum.
Definition: MorphologicalFilter.h:285
In-Memory representation of a mass spectrometry experiment.
Definition: MSExperiment.h:77
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, and used in some special cases.
Definition: MorphologicalFilter.h:563
void setType(SpectrumType type)
sets the spectrum type
void applyErosion_(Int struc_size, InputIterator input, InputIterator input_end, OutputIterator output)
Applies erosion. This implementation uses van Herk&#39;s method. Only 3 min/max comparisons are required ...
Definition: MorphologicalFilter.h:348
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:127
Base class for all classes that want to report their progress.
Definition: ProgressLogger.h:54
IteratorT::value_type::IntensityType & reference
Definition: MorphologicalFilter.h:64
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:91
IteratorT::value_type::IntensityType * pointer
Definition: MorphologicalFilter.h:65
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, and used in some special cases.
Definition: MorphologicalFilter.h:545
void filterExperiment(PeakMap &exp)
Applies the morphological filtering operation to an MSExperiment.
Definition: MorphologicalFilter.h:327
bool operator!=(const IntensityIteratorWrapper &rhs) const
Definition: MorphologicalFilter.h:107
int Int
Signed integer type.
Definition: Types.h:102
IteratorT::value_type::IntensityType value_type
Definition: MorphologicalFilter.h:63
IntensityIteratorWrapper operator++(int)
Definition: MorphologicalFilter.h:95
value_type operator*()
Definition: MorphologicalFilter.h:73