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-2015.
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: Clemens Groepl $
32 // $Authors: $
33 // --------------------------------------------------------------------------
34 
35 #ifndef OPENMS_FILTERING_BASELINE_MORPHOLOGICALFILTER_H
36 #define OPENMS_FILTERING_BASELINE_MORPHOLOGICALFILTER_H
37 
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 
73  value_type operator*()
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 
84  difference_type operator-(IntensityIteratorWrapper & rhs) const
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  template <typename PeakType>
286  void filter(MSSpectrum<PeakType> & 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  struct_size_in_datapoints_ =
298  UInt(
299  ceil(
300  (double)(param_.getValue("struc_elem_length"))
301  *
302  double(spectrum.size() - 1)
303  /
304  (spectrum.back().getMZ() - spectrum.begin()->getMZ())
305  )
306  );
307  }
308  else
309  {
310  struct_size_in_datapoints_ = (UInt)(double)param_.getValue("struc_elem_length");
311  }
312  //make it odd (needed for the algorithm)
313  if (!Math::isOdd(struct_size_in_datapoints_)) ++struct_size_in_datapoints_;
314 
315  //apply the filtering and overwrite the input data
316  std::vector<typename PeakType::IntensityType> output(spectrum.size());
317  filterRange(Internal::intensityIteratorWrapper(spectrum.begin()),
318  Internal::intensityIteratorWrapper(spectrum.end()),
319  output.begin()
320  );
321 
322  //overwrite output with data
323  for (Size i = 0; i < spectrum.size(); ++i)
324  {
325  spectrum[i].setIntensity(output[i]);
326  }
327  }
328 
335  template <typename PeakType>
337  {
338  startProgress(0, exp.size(), "filtering baseline");
339  for (UInt i = 0; i < exp.size(); ++i)
340  {
341  filter(exp[i]);
342  setProgress(i);
343  }
344  endProgress();
345  }
346 
347 protected:
348 
351 
356  template <typename InputIterator, typename OutputIterator>
357  void applyErosion_(Int struc_size, InputIterator input, InputIterator input_end, OutputIterator output)
358  {
359  typedef typename InputIterator::value_type ValueType;
360  const Int size = input_end - input;
361  const Int struc_size_half = struc_size / 2; // yes, integer division
362 
363  static std::vector<ValueType> buffer;
364  if (Int(buffer.size()) < struc_size) buffer.resize(struc_size);
365 
366  Int anchor; // anchoring position of the current block
367  Int i; // index relative to anchor, used for 'for' loops
368  Int ii = 0; // input index
369  Int oi = 0; // output index
370  ValueType current; // current value
371 
372  // we just can't get the case distinctions right in these cases, resorting to simple method.
373  if (size <= struc_size || size <= 5)
374  {
375  applyErosionSimple_(struc_size, input, input_end, output);
376  return;
377  }
378  {
379  // lower margin area
380  current = input[0];
381  for (++ii; ii < struc_size_half; ++ii) if (current > input[ii]) current = input[ii];
382  for (; ii < std::min(Int(struc_size), size); ++ii, ++oi)
383  {
384  if (current > input[ii]) current = input[ii];
385  output[oi] = current;
386  }
387  }
388  {
389  // middle (main) area
390  for (anchor = struc_size;
391  anchor <= size - struc_size;
392  anchor += struc_size
393  )
394  {
395  ii = anchor;
396  current = input[ii];
397  buffer[0] = current;
398  for (i = 1; i < struc_size; ++i, ++ii)
399  {
400  if (current > input[ii]) current = input[ii];
401  buffer[i] = current;
402  }
403  ii = anchor - 1;
404  oi = ii + struc_size_half;
405  current = input[ii];
406  for (i = 1; i < struc_size; ++i, --ii, --oi)
407  {
408  if (current > input[ii]) current = input[ii];
409  output[oi] = std::min(buffer[struc_size - i], current);
410  }
411  if (current > input[ii]) current = input[ii];
412  output[oi] = current;
413  }
414  }
415  {
416  // higher margin area
417  ii = size - 1;
418  oi = ii;
419  current = input[ii];
420  for (--ii; ii >= size - struc_size_half; --ii) if (current > input[ii]) current = input[ii];
421  for (; ii >= std::max(size - Int(struc_size), 0); --ii, --oi)
422  {
423  if (current > input[ii]) current = input[ii];
424  output[oi] = current;
425  }
426  anchor = size - struc_size;
427  ii = anchor;
428  current = input[ii];
429  buffer[0] = current;
430  for (i = 1; i < struc_size; ++i, ++ii)
431  {
432  if (current > input[ii]) current = input[ii];
433  buffer[i] = current;
434  }
435  ii = anchor - 1;
436  oi = ii + struc_size_half;
437  current = input[ii];
438  for (i = 1; (ii >= 0) && (i < struc_size); ++i, --ii, --oi)
439  {
440  if (current > input[ii]) current = input[ii];
441  output[oi] = std::min(buffer[struc_size - i], current);
442  }
443  if (ii >= 0)
444  {
445  if (current > input[ii]) current = input[ii];
446  output[oi] = current;
447  }
448  }
449  return;
450  }
451 
456  template <typename InputIterator, typename OutputIterator>
457  void applyDilation_(Int struc_size, InputIterator input, InputIterator input_end, OutputIterator output)
458  {
459  typedef typename InputIterator::value_type ValueType;
460  const Int size = input_end - input;
461  const Int struc_size_half = struc_size / 2; // yes, integer division
462 
463  static std::vector<ValueType> buffer;
464  if (Int(buffer.size()) < struc_size) buffer.resize(struc_size);
465 
466  Int anchor; // anchoring position of the current block
467  Int i; // index relative to anchor, used for 'for' loops
468  Int ii = 0; // input index
469  Int oi = 0; // output index
470  ValueType current; // current value
471 
472  // we just can't get the case distinctions right in these cases, resorting to simple method.
473  if (size <= struc_size || size <= 5)
474  {
475  applyDilationSimple_(struc_size, input, input_end, output);
476  return;
477  }
478  {
479  // lower margin area
480  current = input[0];
481  for (++ii; ii < struc_size_half; ++ii) if (current < input[ii]) current = input[ii];
482  for (; ii < std::min(Int(struc_size), size); ++ii, ++oi)
483  {
484  if (current < input[ii]) current = input[ii];
485  output[oi] = current;
486  }
487  }
488  {
489  // middle (main) area
490  for (anchor = struc_size;
491  anchor <= size - struc_size;
492  anchor += struc_size
493  )
494  {
495  ii = anchor;
496  current = input[ii];
497  buffer[0] = current;
498  for (i = 1; i < struc_size; ++i, ++ii)
499  {
500  if (current < input[ii]) current = input[ii];
501  buffer[i] = current;
502  }
503  ii = anchor - 1;
504  oi = ii + struc_size_half;
505  current = input[ii];
506  for (i = 1; i < struc_size; ++i, --ii, --oi)
507  {
508  if (current < input[ii]) current = input[ii];
509  output[oi] = std::max(buffer[struc_size - i], current);
510  }
511  if (current < input[ii]) current = input[ii];
512  output[oi] = current;
513  }
514  }
515  {
516  // higher margin area
517  ii = size - 1;
518  oi = ii;
519  current = input[ii];
520  for (--ii; ii >= size - struc_size_half; --ii) if (current < input[ii]) current = input[ii];
521  for (; ii >= std::max(size - Int(struc_size), 0); --ii, --oi)
522  {
523  if (current < input[ii]) current = input[ii];
524  output[oi] = current;
525  }
526  anchor = size - struc_size;
527  ii = anchor;
528  current = input[ii];
529  buffer[0] = current;
530  for (i = 1; i < struc_size; ++i, ++ii)
531  {
532  if (current < input[ii]) current = input[ii];
533  buffer[i] = current;
534  }
535  ii = anchor - 1;
536  oi = ii + struc_size_half;
537  current = input[ii];
538  for (i = 1; (ii >= 0) && (i < struc_size); ++i, --ii, --oi)
539  {
540  if (current < input[ii]) current = input[ii];
541  output[oi] = std::max(buffer[struc_size - i], current);
542  }
543  if (ii >= 0)
544  {
545  if (current < input[ii]) current = input[ii];
546  output[oi] = current;
547  }
548  }
549  return;
550  }
551 
553  template <typename InputIterator, typename OutputIterator>
554  void applyErosionSimple_(Int struc_size, InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
555  {
556  typedef typename InputIterator::value_type ValueType;
557  const int size = input_end - input_begin;
558  const Int struc_size_half = struc_size / 2; // yes integer division
559  for (Int index = 0; index < size; ++index)
560  {
561  Int start = std::max(0, index - struc_size_half);
562  Int stop = std::min(size - 1, index + struc_size_half);
563  ValueType value = input_begin[start];
564  for (Int i = start + 1; i <= stop; ++i) if (value > input_begin[i]) value = input_begin[i];
565  output_begin[index] = value;
566  }
567  return;
568  }
569 
571  template <typename InputIterator, typename OutputIterator>
572  void applyDilationSimple_(Int struc_size, InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
573  {
574  typedef typename InputIterator::value_type ValueType;
575  const int size = input_end - input_begin;
576  const Int struc_size_half = struc_size / 2; // yes integer division
577  for (Int index = 0; index < size; ++index)
578  {
579  Int start = std::max(0, index - struc_size_half);
580  Int stop = std::min(size - 1, index + struc_size_half);
581  ValueType value = input_begin[start];
582  for (Int i = start + 1; i <= stop; ++i) if (value < input_begin[i]) value = input_begin[i];
583  output_begin[index] = value;
584  }
585  return;
586  }
587 
588 private:
589 
592 
593  };
594 
595 } // namespace OpenMS
596 
597 #endif
UInt struct_size_in_datapoints_
Member for struct size in data points.
Definition: MorphologicalFilter.h:350
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:126
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's method. Only 3 min/max comparisons are required...
Definition: MorphologicalFilter.h:457
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
Size size() const
Definition: MSExperiment.h:117
bool operator==(const IntensityIteratorWrapper &rhs) const
Definition: MorphologicalFilter.h:102
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:88
virtual ~MorphologicalFilter()
Destructor.
Definition: MorphologicalFilter.h:184
bool operator!=(const IntensityIteratorWrapper &rhs) const
Definition: MorphologicalFilter.h:107
Raw data (also called profile data)
Definition: SpectrumSettings.h:75
Main OpenMS namespace.
Definition: FeatureDeconvolution.h:47
void filter(MSSpectrum< PeakType > &spectrum)
Applies the morphological filtering operation to an MSSpectrum.
Definition: MorphologicalFilter.h:286
void filterExperiment(MSExperiment< PeakType > &exp)
Applies the morphological filtering operation to an MSExperiment.
Definition: MorphologicalFilter.h:336
difference_type operator-(IntensityIteratorWrapper &rhs) const
Definition: MorphologicalFilter.h:84
In-Memory representation of a mass spectrometry experiment.
Definition: MSExperiment.h:69
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:572
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's method. Only 3 min/max comparisons are required ...
Definition: MorphologicalFilter.h:357
Base class for all classes that want to report their progress.
Definition: ProgressLogger.h:55
IteratorT::value_type::IntensityType & reference
Definition: MorphologicalFilter.h:64
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:92
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:554
int Int
Signed integer type.
Definition: Types.h:96
IteratorT::value_type::IntensityType value_type
Definition: MorphologicalFilter.h:63
IntensityIteratorWrapper operator++(int)
Definition: MorphologicalFilter.h:95
value_type operator*()
Definition: MorphologicalFilter.h:73

OpenMS / TOPP release 2.0.0 Documentation generated on Thu Aug 20 2015 01:44:25 using doxygen 1.8.9.1