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

OpenMS / TOPP release 2.3.0 Documentation generated on Tue Jan 9 2018 18:22:02 using doxygen 1.8.13