MorphologicalFilter.h

Go to the documentation of this file.

// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2023.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Timo Sachsenberg $
// $Authors: $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/CONCEPT/ProgressLogger.h>
#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
#include <OpenMS/DATASTRUCTURES/ListUtils.h>
#include <OpenMS/KERNEL/MSExperiment.h>
#include <OpenMS/KERNEL/StandardTypes.h>
#include <OpenMS/MATH/MISC/MathFunctions.h>
#include <algorithm>
#include <iterator>
 
namespace OpenMS
{
 
  namespace Internal
  {
    template<typename IteratorT>
    class /* OPENMS_DLLAPI */ IntensityIteratorWrapper
    {
    public:
      typedef std::forward_iterator_tag iterator_category;
      typedef typename IteratorT::value_type::IntensityType value_type;
      typedef typename IteratorT::value_type::IntensityType& reference;
      typedef typename IteratorT::value_type::IntensityType* pointer;
      typedef typename IteratorT::difference_type difference_type;
 
      IntensityIteratorWrapper(const IteratorT& rhs) : base(rhs)
      {
      }
 
      value_type operator*()
      {
        return base->getIntensity();
      }
 
      template<typename IndexT>
      value_type operator[](const IndexT& index)
      {
        return base[index].getIntensity();
      }
 
      difference_type operator-(IntensityIteratorWrapper& rhs) const
      {
        return base - rhs.base;
      }
 
      IntensityIteratorWrapper& operator++()
      {
        ++base;
        return *this;
      }
 
      IntensityIteratorWrapper operator++(int)
      {
        IteratorT tmp = *this;
        ++(*this);
        return tmp;
      }
 
      bool operator==(const IntensityIteratorWrapper& rhs) const
      {
        return base == rhs.base;
      }
 
      bool operator!=(const IntensityIteratorWrapper& rhs) const
      {
        return base != rhs.base;
      }
 
    protected:
      IteratorT base;
    };
 
    template<typename IteratorT>
    IntensityIteratorWrapper<IteratorT> intensityIteratorWrapper(const IteratorT& rhs)
    {
      return IntensityIteratorWrapper<IteratorT>(rhs);
    }
 
  } // namespace Internal
 
  class OPENMS_DLLAPI MorphologicalFilter : public ProgressLogger, public DefaultParamHandler
  {
  public:
    MorphologicalFilter() : ProgressLogger(), DefaultParamHandler("MorphologicalFilter"), struct_size_in_datapoints_(0)
    {
      // structuring element
      defaults_.setValue("struc_elem_length", 3.0, "Length of the structuring element. This should be wider than the expected peak width.");
      defaults_.setValue("struc_elem_unit", "Thomson", "The unit of the 'struct_elem_length'.");
      defaults_.setValidStrings("struc_elem_unit", {"Thomson", "DataPoints"});
      // methods
      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.");
      defaults_.setValidStrings("method", {"identity", "erosion", "dilation", "opening", "closing", "gradient", "tophat", "bothat", "erosion_simple", "dilation_simple"});
 
      defaultsToParam_();
    }
 
    ~MorphologicalFilter() override
    {
    }
 
    template<typename InputIterator, typename OutputIterator>
    void filterRange(InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
    {
      // the buffer is static only to avoid reallocation
      static std::vector<typename InputIterator::value_type> buffer;
      const UInt size = input_end - input_begin;
 
      // determine the struct size in data points if not already set
      if (struct_size_in_datapoints_ == 0)
      {
        struct_size_in_datapoints_ = (UInt)(double)param_.getValue("struc_elem_length");
      }
 
      // apply the filtering
      std::string method = param_.getValue("method");
      if (method == "identity")
      {
        std::copy(input_begin, input_end, output_begin);
      }
      else if (method == "erosion")
      {
        applyErosion_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
      }
      else if (method == "dilation")
      {
        applyDilation_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
      }
      else if (method == "opening")
      {
        if (buffer.size() < size)
          buffer.resize(size);
        applyErosion_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
        applyDilation_(struct_size_in_datapoints_, buffer.begin(), buffer.begin() + size, output_begin);
      }
      else if (method == "closing")
      {
        if (buffer.size() < size)
          buffer.resize(size);
        applyDilation_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
        applyErosion_(struct_size_in_datapoints_, buffer.begin(), buffer.begin() + size, output_begin);
      }
      else if (method == "gradient")
      {
        if (buffer.size() < size)
          buffer.resize(size);
        applyErosion_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
        applyDilation_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
        for (UInt i = 0; i < size; ++i)
          output_begin[i] -= buffer[i];
      }
      else if (method == "tophat")
      {
        if (buffer.size() < size)
          buffer.resize(size);
        applyErosion_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
        applyDilation_(struct_size_in_datapoints_, buffer.begin(), buffer.begin() + size, output_begin);
        for (UInt i = 0; i < size; ++i)
          output_begin[i] = input_begin[i] - output_begin[i];
      }
      else if (method == "bothat")
      {
        if (buffer.size() < size)
          buffer.resize(size);
        applyDilation_(struct_size_in_datapoints_, input_begin, input_end, buffer.begin());
        applyErosion_(struct_size_in_datapoints_, buffer.begin(), buffer.begin() + size, output_begin);
        for (UInt i = 0; i < size; ++i)
          output_begin[i] = input_begin[i] - output_begin[i];
      }
      else if (method == "erosion_simple")
      {
        applyErosionSimple_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
      }
      else if (method == "dilation_simple")
      {
        applyDilationSimple_(struct_size_in_datapoints_, input_begin, input_end, output_begin);
      }
 
      struct_size_in_datapoints_ = 0;
    }
 
    void filter(MSSpectrum& spectrum)
    {
      // make sure the right peak type is set
      spectrum.setType(SpectrumSettings::PROFILE);
 
      // Abort if there is nothing to do
      if (spectrum.size() <= 1)
      {
        return;
      }
 
      // Determine structuring element size in datapoints (depending on the unit)
      if (param_.getValue("struc_elem_unit") == "Thomson")
      {
        const double struc_elem_length = (double)param_.getValue("struc_elem_length");
        const double mz_diff = spectrum.back().getMZ() - spectrum.begin()->getMZ();
        struct_size_in_datapoints_ = (UInt)(ceil(struc_elem_length * (double)(spectrum.size() - 1) / mz_diff));
      }
      else
      {
        struct_size_in_datapoints_ = (UInt)(double)param_.getValue("struc_elem_length");
      }
      // make it odd (needed for the algorithm)
      if (!Math::isOdd(struct_size_in_datapoints_))
        ++struct_size_in_datapoints_;
 
      // apply the filtering and overwrite the input data
      std::vector<Peak1D::IntensityType> output(spectrum.size());
      filterRange(Internal::intensityIteratorWrapper(spectrum.begin()), Internal::intensityIteratorWrapper(spectrum.end()), output.begin());
 
      // overwrite output with data
      for (Size i = 0; i < spectrum.size(); ++i)
      {
        spectrum[i].setIntensity(output[i]);
      }
    }
 
    void filterExperiment(PeakMap& exp)
    {
      startProgress(0, exp.size(), "filtering baseline");
      for (UInt i = 0; i < exp.size(); ++i)
      {
        filter(exp[i]);
        setProgress(i);
      }
      endProgress();
    }
 
  protected:
    UInt struct_size_in_datapoints_;
 
    template<typename InputIterator, typename OutputIterator>
    void applyErosion_(Int struc_size, InputIterator input, InputIterator input_end, OutputIterator output)
    {
      typedef typename InputIterator::value_type ValueType;
      const Int size = input_end - input;
      const Int struc_size_half = struc_size / 2; // yes, integer division
 
      static std::vector<ValueType> buffer;
      if (Int(buffer.size()) < struc_size)
        buffer.resize(struc_size);
 
      Int anchor;        // anchoring position of the current block
      Int i;             // index relative to anchor, used for 'for' loops
      Int ii = 0;        // input index
      Int oi = 0;        // output index
      ValueType current; // current value
 
      // we just can't get the case distinctions right in these cases, resorting to simple method.
      if (size <= struc_size || size <= 5)
      {
        applyErosionSimple_(struc_size, input, input_end, output);
        return;
      }
      {
        // lower margin area
        current = input[0];
        for (++ii; ii < struc_size_half; ++ii)
          if (current > input[ii])
            current = input[ii];
        for (; ii < std::min(Int(struc_size), size); ++ii, ++oi)
        {
          if (current > input[ii])
            current = input[ii];
          output[oi] = current;
        }
      }
      {
        // middle (main) area
        for (anchor = struc_size; anchor <= size - struc_size; anchor += struc_size)
        {
          ii = anchor;
          current = input[ii];
          buffer[0] = current;
          for (i = 1; i < struc_size; ++i, ++ii)
          {
            if (current > input[ii])
              current = input[ii];
            buffer[i] = current;
          }
          ii = anchor - 1;
          oi = ii + struc_size_half;
          current = input[ii];
          for (i = 1; i < struc_size; ++i, --ii, --oi)
          {
            if (current > input[ii])
              current = input[ii];
            output[oi] = std::min(buffer[struc_size - i], current);
          }
          if (current > input[ii])
            current = input[ii];
          output[oi] = current;
        }
      }
      {
        // higher margin area
        ii = size - 1;
        oi = ii;
        current = input[ii];
        for (--ii; ii >= size - struc_size_half; --ii)
          if (current > input[ii])
            current = input[ii];
        for (; ii >= std::max(size - Int(struc_size), 0); --ii, --oi)
        {
          if (current > input[ii])
            current = input[ii];
          output[oi] = current;
        }
        anchor = size - struc_size;
        ii = anchor;
        current = input[ii];
        buffer[0] = current;
        for (i = 1; i < struc_size; ++i, ++ii)
        {
          if (current > input[ii])
            current = input[ii];
          buffer[i] = current;
        }
        ii = anchor - 1;
        oi = ii + struc_size_half;
        current = input[ii];
        for (i = 1; (ii >= 0) && (i < struc_size); ++i, --ii, --oi)
        {
          if (current > input[ii])
            current = input[ii];
          output[oi] = std::min(buffer[struc_size - i], current);
        }
        if (ii >= 0)
        {
          if (current > input[ii])
            current = input[ii];
          output[oi] = current;
        }
      }
      return;
    }
 
    template<typename InputIterator, typename OutputIterator>
    void applyDilation_(Int struc_size, InputIterator input, InputIterator input_end, OutputIterator output)
    {
      typedef typename InputIterator::value_type ValueType;
      const Int size = input_end - input;
      const Int struc_size_half = struc_size / 2; // yes, integer division
 
      static std::vector<ValueType> buffer;
      if (Int(buffer.size()) < struc_size)
        buffer.resize(struc_size);
 
      Int anchor;        // anchoring position of the current block
      Int i;             // index relative to anchor, used for 'for' loops
      Int ii = 0;        // input index
      Int oi = 0;        // output index
      ValueType current; // current value
 
      // we just can't get the case distinctions right in these cases, resorting to simple method.
      if (size <= struc_size || size <= 5)
      {
        applyDilationSimple_(struc_size, input, input_end, output);
        return;
      }
      {
        // lower margin area
        current = input[0];
        for (++ii; ii < struc_size_half; ++ii)
          if (current < input[ii])
            current = input[ii];
        for (; ii < std::min(Int(struc_size), size); ++ii, ++oi)
        {
          if (current < input[ii])
            current = input[ii];
          output[oi] = current;
        }
      }
      {
        // middle (main) area
        for (anchor = struc_size; anchor <= size - struc_size; anchor += struc_size)
        {
          ii = anchor;
          current = input[ii];
          buffer[0] = current;
          for (i = 1; i < struc_size; ++i, ++ii)
          {
            if (current < input[ii])
              current = input[ii];
            buffer[i] = current;
          }
          ii = anchor - 1;
          oi = ii + struc_size_half;
          current = input[ii];
          for (i = 1; i < struc_size; ++i, --ii, --oi)
          {
            if (current < input[ii])
              current = input[ii];
            output[oi] = std::max(buffer[struc_size - i], current);
          }
          if (current < input[ii])
            current = input[ii];
          output[oi] = current;
        }
      }
      {
        // higher margin area
        ii = size - 1;
        oi = ii;
        current = input[ii];
        for (--ii; ii >= size - struc_size_half; --ii)
          if (current < input[ii])
            current = input[ii];
        for (; ii >= std::max(size - Int(struc_size), 0); --ii, --oi)
        {
          if (current < input[ii])
            current = input[ii];
          output[oi] = current;
        }
        anchor = size - struc_size;
        ii = anchor;
        current = input[ii];
        buffer[0] = current;
        for (i = 1; i < struc_size; ++i, ++ii)
        {
          if (current < input[ii])
            current = input[ii];
          buffer[i] = current;
        }
        ii = anchor - 1;
        oi = ii + struc_size_half;
        current = input[ii];
        for (i = 1; (ii >= 0) && (i < struc_size); ++i, --ii, --oi)
        {
          if (current < input[ii])
            current = input[ii];
          output[oi] = std::max(buffer[struc_size - i], current);
        }
        if (ii >= 0)
        {
          if (current < input[ii])
            current = input[ii];
          output[oi] = current;
        }
      }
      return;
    }
 
    template<typename InputIterator, typename OutputIterator>
    void applyErosionSimple_(Int struc_size, InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
    {
      typedef typename InputIterator::value_type ValueType;
      const int size = input_end - input_begin;
      const Int struc_size_half = struc_size / 2; // yes integer division
      for (Int index = 0; index < size; ++index)
      {
        Int start = std::max(0, index - struc_size_half);
        Int stop = std::min(size - 1, index + struc_size_half);
        ValueType value = input_begin[start];
        for (Int i = start + 1; i <= stop; ++i)
          if (value > input_begin[i])
            value = input_begin[i];
        output_begin[index] = value;
      }
      return;
    }
 
    template<typename InputIterator, typename OutputIterator>
    void applyDilationSimple_(Int struc_size, InputIterator input_begin, InputIterator input_end, OutputIterator output_begin)
    {
      typedef typename InputIterator::value_type ValueType;
      const int size = input_end - input_begin;
      const Int struc_size_half = struc_size / 2; // yes integer division
      for (Int index = 0; index < size; ++index)
      {
        Int start = std::max(0, index - struc_size_half);
        Int stop = std::min(size - 1, index + struc_size_half);
        ValueType value = input_begin[start];
        for (Int i = start + 1; i <= stop; ++i)
          if (value < input_begin[i])
            value = input_begin[i];
        output_begin[index] = value;
      }
      return;
    }
 
  private:
    MorphologicalFilter(const MorphologicalFilter& source);
  };
 
} // namespace OpenMS