KERNEL/PeakTypeEstimator.h

Go to the documentation of this file.

// Copyright (c) 2002-present, OpenMS Inc. -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Chris Bielow $
// $Authors: Chris Bielow $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/METADATA/SpectrumSettings.h>
 
#include <limits>
#include <cmath>
#include <numeric>
#include <vector>
 
namespace OpenMS
{
  class OPENMS_DLLAPI PeakTypeEstimator
  {
public:
    template <typename PeakConstIterator>
    static SpectrumSettings::SpectrumType estimateType(const PeakConstIterator& begin, const PeakConstIterator& end)
    {
      typedef typename PeakConstIterator::value_type PeakT;
      // abort if there are less than 5 peak in the iterator range
      if (end - begin < 5)
      {
        return SpectrumSettings::SpectrumType::UNKNOWN;
      }
 
      const int max_peaks = 5; // maximal number of peaks we are looking at
      int profile_evidence = 0; // number of peaks found to be profile
      int centroid_evidence = 0; // number of peaks found to be centroided
 
      // copy data, since we need to modify
      std::vector<PeakT> data(begin, end);
      // total intensity of spectrum
      double total_int = std::accumulate(begin, end, 0.0, [](double int_, const PeakT& p) { return int_ + p.getIntensity(); } );
      double explained_int = 0;
      // get the 5 highest peaks
      for (int i = 0; i < max_peaks; ++i)
      {
        // stop if we explained +50% of all intensity
        // (due to danger of interpreting noise - usually wrongly classified as centroided data)
        if (explained_int > 0.5 * total_int) break;
 
        double int_max = 0;
        Size idx = std::numeric_limits<Size>::max();
        // find highest peak position
        for (Size i = 0; i < data.size(); ++i)
        {
          if (data[i].getIntensity() > int_max)
          {
            int_max = data[i].getIntensity();
            idx = i;
          }
        }
        // no more peaks
        if (idx == std::numeric_limits<Size>::max()) break;
 
        // check left and right peak shoulders and count number of sample points
        typedef typename std::vector<PeakT>::iterator PeakIterator; // non-const version, since we need to modify the peaks
        PeakIterator it_max = data.begin() + idx;
        PeakIterator it = it_max;
        double int_last = int_max;
        while (it != data.begin()
               && it->getIntensity() <= int_last        // at most 100% of last sample point
               && it->getIntensity() > 0
               && (it->getIntensity() / int_last) > 0.1 // at least 10% of last sample point
               && it->getMZ() + 1 > it_max->getMZ())    // at most 1 Th away
        {
          int_last = it->getIntensity();
          explained_int += int_last;
          it->setIntensity(0); // remove peak from future consideration
          --it;
        }
        // if the current point is rising again, restore the intensity of the
        // previous 'sink' point (because it could belong to a neighbour peak)
        // e.g. imagine intensities: 1-2-3-2-1-2-4-2-1. We do not want to destroy the middle '1'
        if (it->getIntensity() > int_last) (it+1)->setIntensity(int_last);
 
        //std::cerr << "  Peak candidate: " << it_max->getMZ() << " ...";
        bool break_left = false;
        if (it_max - it < 2+1)  // 'it' does not fulfill the conditions, i.e. does not count
        { // fewer than two sampling points on left shoulder
          //std::cerr << " break left " << it_max - it << " points\n";
          break_left = true;
          // do not end loop here.. we still need to clean up the right side
        }
        it_max->setIntensity(int_max); // restore center intensity
        explained_int -= int_max;
        it = it_max;
        int_last = int_max;
        while (it != data.end()
               && it->getIntensity() <= int_last        // at most 100% of last sample point
               && it->getIntensity() > 0
               && (it->getIntensity() / int_last) > 0.1 // at least 10% of last sample point
               && it->getMZ() - 1 < it_max->getMZ())    // at most 1 Th away
        {
          int_last = it->getIntensity();
          explained_int += int_last;
          it->setIntensity(0); // remove peak from future consideration
          ++it;
        }
        // if the current point is rising again, restore the intensity of the
        // previous 'sink' point (because it could belong to a neighbour peak)
        // e.g. imagine intensities: 1-2-4-2-1-2-3-2-1. We do not want to destroy the middle '1'
        // (note: the sequence is not identical to the one of the left shoulder)
        if (it != data.end() && it->getIntensity() > int_last) (it-1)->setIntensity(int_last);
 
        if (break_left || it - it_max < 2+1)  // 'it' does not fulfill the conditions, i.e. does not count
        { // fewer than two sampling points on right shoulder
          //std::cerr << " break right " << it - it_max << " points\n";
          ++centroid_evidence;
          continue;
        }
        // peak has at least two sampling points on either side within 1 Th
        ++profile_evidence;
        //std::cerr << " PROFILE " << it - it_max << " points right\n";
 
      }
 
      float evidence_ratio = profile_evidence / float(profile_evidence + centroid_evidence);
      //std::cerr << "--> Evidence ratio: " << evidence_ratio;
 
      if (evidence_ratio > 0.75) // 80% are profile
      {
        //std::cerr << "  PROFILE\n";
        return SpectrumSettings::SpectrumType::PROFILE;
      }
      else
      {
        //std::cerr << "  CENTROID\n";
        return SpectrumSettings::SpectrumType::CENTROID;
      }
    }
  };
 
} // namespace OpenMS