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

OpenMS / TOPP release 1.11.1 Documentation generated on Thu Nov 14 2013 11:19:17 using doxygen 1.8.5