PeakIntegrator.h

Go to the documentation of this file.

// Copyright (c) 2002-2023, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Douglas McCloskey, Pasquale Domenico Colaianni $
// $Authors: Douglas McCloskey, Pasquale Domenico Colaianni $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/config.h> // OPENMS_DLLAPI
#include <OpenMS/CONCEPT/LogStream.h>
#include <OpenMS/DATASTRUCTURES/ConvexHull2D.h>
#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
#include <OpenMS/DATASTRUCTURES/String.h>
#include <OpenMS/KERNEL/MSChromatogram.h>
#include <OpenMS/KERNEL/MSSpectrum.h>
#include <OpenMS/MATH/MISC/EmgGradientDescent.h>
 
namespace OpenMS
{
  class OPENMS_DLLAPI PeakIntegrator :
    public DefaultParamHandler
  {
public:
    PeakIntegrator();
    ~PeakIntegrator() override;
 
    struct PeakArea
    {
      double area = 0.0;
      double height = 0.0;
      double apex_pos = 0.0;
      ConvexHull2D::PointArrayType hull_points;
    };
 
    struct PeakBackground
    {
      double area = 0.0;
      double height = 0.0;
    };
 
    struct PeakShapeMetrics
    {
      double width_at_5 = 0.0;
      double width_at_10 = 0.0;
      double width_at_50 = 0.0;
      double start_position_at_5 = 0.0;
      double start_position_at_10 = 0.0;
      double start_position_at_50 = 0.0;
      double end_position_at_5 = 0.0;
      double end_position_at_10 = 0.0;
      double end_position_at_50 = 0.0;
      double total_width = 0.0;
      double tailing_factor = 0.0;
      double asymmetry_factor = 0.0;
      double slope_of_baseline = 0.0;
      double baseline_delta_2_height = 0.0;
      Int points_across_baseline = 0;
      Int points_across_half_height = 0;
    };
 
    static constexpr const char* INTEGRATION_TYPE_INTENSITYSUM = "intensity_sum";
    static constexpr const char* INTEGRATION_TYPE_TRAPEZOID = "trapezoid";
    static constexpr const char* INTEGRATION_TYPE_SIMPSON = "simpson";
    static constexpr const char* BASELINE_TYPE_BASETOBASE = "base_to_base";
    static constexpr const char* BASELINE_TYPE_VERTICALDIVISION = "vertical_division";
    static constexpr const char* BASELINE_TYPE_VERTICALDIVISION_MIN = "vertical_division_min";
    static constexpr const char* BASELINE_TYPE_VERTICALDIVISION_MAX = "vertical_division_max";
 
    PeakArea integratePeak(
      const MSChromatogram& chromatogram, const double left, const double right
    ) const;
 
    PeakArea integratePeak(
      const MSChromatogram& chromatogram, MSChromatogram::ConstIterator& left, MSChromatogram::ConstIterator& right
    ) const;
 
    PeakArea integratePeak(
      const MSSpectrum& spectrum, const double left, const double right
    ) const;
 
    PeakArea integratePeak(
      const MSSpectrum& spectrum, MSSpectrum::ConstIterator& left, MSSpectrum::ConstIterator& right
    ) const;
 
    PeakBackground estimateBackground(
      const MSChromatogram& chromatogram, const double left, const double right,
      const double peak_apex_pos
    ) const;
 
    PeakBackground estimateBackground(
      const MSChromatogram& chromatogram, MSChromatogram::ConstIterator& left, MSChromatogram::ConstIterator& right,
      const double peak_apex_pos
    ) const;
 
    PeakBackground estimateBackground(
      const MSSpectrum& spectrum, const double left, const double right,
      const double peak_apex_pos
    ) const;
 
    PeakBackground estimateBackground(
      const MSSpectrum& spectrum, MSSpectrum::ConstIterator& left, MSSpectrum::ConstIterator& right,
      const double peak_apex_pos
    ) const;
 
    PeakShapeMetrics calculatePeakShapeMetrics(
      const MSChromatogram& chromatogram, const double left, const double right,
      const double peak_height, const double peak_apex_pos
    ) const;
 
    PeakShapeMetrics calculatePeakShapeMetrics(
      const MSChromatogram& chromatogram, MSChromatogram::ConstIterator& left, MSChromatogram::ConstIterator& right,
      const double peak_height, const double peak_apex_pos
    ) const;
 
    PeakShapeMetrics calculatePeakShapeMetrics(
      const MSSpectrum& spectrum, const double left, const double right,
      const double peak_height, const double peak_apex_pos
    ) const;
 
    PeakShapeMetrics calculatePeakShapeMetrics(
      const MSSpectrum& spectrum, MSSpectrum::ConstIterator& left, MSSpectrum::ConstIterator& right,
      const double peak_height, const double peak_apex_pos
    ) const;
 
    void getDefaultParameters(Param& params);
 
protected:
    void updateMembers_() override;
 
    template <typename PeakContainerT>
    PeakArea integratePeak_(const PeakContainerT& pc, double left, double right) const
    {
      OPENMS_PRECONDITION(left <= right, "Left peak boundary must be smaller than right boundary!") // otherwise the code below will segfault (due to PosBegin/PosEnd)
      PeakContainerT emg_pc;
      const PeakContainerT& p = EMGPreProcess_(pc, emg_pc, left, right);
 
      std::function<double(const double, const double)>
      compute_peak_area_trapezoid = [&p](const double left, const double right)
      {
        double peak_area { 0.0 };
        for (typename PeakContainerT::ConstIterator it = p.PosBegin(left); it != p.PosEnd(right) - 1; ++it)
        {
          peak_area += ((it + 1)->getPos() - it->getPos()) * ((it->getIntensity() + (it + 1)->getIntensity()) / 2.0);
        }
        return peak_area;
      };
 
      std::function<double(const double, const double)>
      compute_peak_area_intensity_sum = [&p](const double left, const double right)
      {
        // OPENMS_LOG_WARN << "WARNING: intensity_sum method is being used." << std::endl;
        double peak_area { 0.0 };
        for (typename PeakContainerT::ConstIterator it = p.PosBegin(left); it != p.PosEnd(right); ++it)
        {
          peak_area += it->getIntensity();
        }
        return peak_area;
      };
 
      PeakArea pa;
      pa.apex_pos = (left + right) / 2; // initial estimate, to avoid apex being outside of [left,right]
      UInt n_points = std::distance(p.PosBegin(left), p.PosEnd(right));
      for (auto it = p.PosBegin(left); it != p.PosEnd(right); ++it) //OMS_CODING_TEST_EXCLUDE
      {
        pa.hull_points.push_back(DPosition<2>(it->getPos(), it->getIntensity()));
        if (pa.height < it->getIntensity())
        {
          pa.height = it->getIntensity();
          pa.apex_pos = it->getPos();
        }
      }
 
      if (integration_type_ == INTEGRATION_TYPE_TRAPEZOID)
      {
        if (n_points >= 2)
        {
          pa.area = compute_peak_area_trapezoid(left, right);
        }
      }
      else if (integration_type_ == INTEGRATION_TYPE_SIMPSON)
      {
        if (n_points == 2)
        {
          OPENMS_LOG_WARN << std::endl << "PeakIntegrator::integratePeak:"
            "number of points is 2, falling back to `trapezoid`." << std::endl;
          pa.area = compute_peak_area_trapezoid(left, right);
        }
        else if (n_points > 2)
        {
          if (n_points % 2)
          {
            pa.area = simpson_(p.PosBegin(left), p.PosEnd(right));
          }
          else
          {
            double areas[4] = {-1.0, -1.0, -1.0, -1.0};
            areas[0] = simpson_(p.PosBegin(left), p.PosEnd(right) - 1);   // without last point
            areas[1] = simpson_(p.PosBegin(left) + 1, p.PosEnd(right));   // without first point
            if (p.begin() <= p.PosBegin(left) - 1)
            {
              areas[2] = simpson_(p.PosBegin(left) - 1, p.PosEnd(right)); // with one more point on the left
            }
            if (p.PosEnd(right) < p.end())
            {
              areas[3] = simpson_(p.PosBegin(left), p.PosEnd(right) + 1); // with one more point on the right
            }
            UInt valids = 0;
            for (const auto& area : areas)
            {
              if (area != -1.0)
              {
                pa.area += area;
                ++valids;
              }
            }
            pa.area /= valids;
          }
        }
      }
      else if (integration_type_ == INTEGRATION_TYPE_INTENSITYSUM)
      {
        pa.area = compute_peak_area_intensity_sum(left, right);
      }
      else
      {
        throw Exception::InvalidParameter(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Please set a valid value for the parameter \"integration_type\".");
      }
    
      return pa;
    }
 
 
 
    template <typename PeakContainerT>
    PeakBackground estimateBackground_(
      const PeakContainerT& pc, double left, double right,
      const double peak_apex_pos
    ) const
    {
      PeakContainerT emg_pc;
      const PeakContainerT& p = EMGPreProcess_(pc, emg_pc, left, right);
 
      const double int_l = p.PosBegin(left)->getIntensity();
      const double int_r = (p.PosEnd(right) - 1)->getIntensity();
      const double delta_int = int_r - int_l;
      const double delta_pos = (p.PosEnd(right) - 1)->getPos() - p.PosBegin(left)->getPos();
      const double min_int_pos = int_r <= int_l ? (p.PosEnd(right) - 1)->getPos() : p.PosBegin(left)->getPos();
      const double delta_int_apex = std::fabs(delta_int) * std::fabs(min_int_pos - peak_apex_pos) / delta_pos;
      double area {0.0};
      double height {0.0};
      if (baseline_type_ == BASELINE_TYPE_BASETOBASE)
      {
        height = std::min(int_r, int_l) + delta_int_apex;
        if (integration_type_ == INTEGRATION_TYPE_TRAPEZOID || integration_type_ == INTEGRATION_TYPE_SIMPSON)
        {
          // formula for calculating the background using the trapezoidal rule
          // area = intensity_min*delta_pos + 0.5*delta_int*delta_pos;
          area = delta_pos * (std::min(int_r, int_l) + 0.5 * std::fabs(delta_int));
        }
        else if (integration_type_ == INTEGRATION_TYPE_INTENSITYSUM)
        {
          // calculate the background using an estimator of the form
          //    y = mx + b
          //    where x = rt or mz, m = slope, b = left intensity
          // sign of delta_int will determine line direction
          // area += delta_int / delta_pos * (it->getPos() - left) + int_l;
          double pos_sum = 0.0; // rt or mz
          for (auto it = p.PosBegin(left); it != p.PosEnd(right); ++it) //OMS_CODING_TEST_EXCLUDE
          {
            pos_sum += it->getPos();
          }
          UInt n_points = std::distance(p.PosBegin(left), p.PosEnd(right));
 
          // We construct the background area as the sum of a rectangular part
          // and a triangle on top. The triangle is constructed as the sum of the
          // line's y value at each sampled point: \sum_{i=0}^{n} (x_i - x_0)  * m
          const double rectangle_area = n_points * int_l;
          const double slope = delta_int / delta_pos;
          const double triangle_area = (pos_sum - n_points * p.PosBegin(left)->getPos()) * slope;
          area = triangle_area + rectangle_area;
        }
      }
      else if (baseline_type_ == BASELINE_TYPE_VERTICALDIVISION || baseline_type_ == BASELINE_TYPE_VERTICALDIVISION_MIN)
      {
        height = std::min(int_r, int_l);
        if (integration_type_ == INTEGRATION_TYPE_TRAPEZOID || integration_type_ == INTEGRATION_TYPE_SIMPSON)
        {
          area = delta_pos * std::min(int_r, int_l);
        }
        else if (integration_type_ == INTEGRATION_TYPE_INTENSITYSUM)
        {
          area = std::min(int_r, int_l) * std::distance(p.PosBegin(left), p.PosEnd(right));;
        }
      }
      else if (baseline_type_ == BASELINE_TYPE_VERTICALDIVISION_MAX)
      {
        height = std::max(int_r, int_l);
        if (integration_type_ == INTEGRATION_TYPE_TRAPEZOID || integration_type_ == INTEGRATION_TYPE_SIMPSON)
        {
          area = delta_pos * std::max(int_r, int_l);
        }
        else if (integration_type_ == INTEGRATION_TYPE_INTENSITYSUM)
        {
          area = std::max(int_r, int_l) * std::distance(p.PosBegin(left), p.PosEnd(right));
        }
      }
      else
      {
        throw Exception::InvalidParameter(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION, "Please set a valid value for the parameter \"baseline_type\".");
      }
      PeakBackground pb;
      pb.area = area;
      pb.height = height;
      return pb;
    }
 
    template <typename PeakContainerConstIteratorT>
    double simpson_(PeakContainerConstIteratorT it_begin, PeakContainerConstIteratorT it_end) const
    {
      double integral = 0.0;
      for (auto it = it_begin + 1; it < it_end - 1; it = it + 2) //OMS_CODING_TEST_EXCLUDE
      {
        const double h = it->getPos() - (it - 1)->getPos();
        const double k = (it + 1)->getPos() - it->getPos();
        const double y_h = (it - 1)->getIntensity();
        const double y_0 = it->getIntensity();
        const double y_k = (it + 1)->getIntensity();
        integral += (1.0 / 6.0) * (h + k) * ((2.0 - k / h) * y_h + (pow(h + k, 2) / (h * k)) * y_0 + (2.0 - h / k) * y_k);
      }
      return integral;
    }
 
 
    template <typename PeakContainerT>
    PeakShapeMetrics calculatePeakShapeMetrics_(
      const PeakContainerT& pc, double left, double right,
      const double peak_height, const double peak_apex_pos
    ) const
    {
      PeakShapeMetrics psm;
 
      if (pc.empty()) return psm; // return all '0'
 
      // enforce order: left <= peakapex <= right
      if (!(left <= peak_apex_pos && peak_apex_pos <= right)) throw Exception::InvalidRange(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION);
 
      PeakContainerT emg_pc;
      const PeakContainerT& p = EMGPreProcess_(pc, emg_pc, left, right);
    
      typename PeakContainerT::ConstIterator it_PosBegin_l = p.PosBegin(left);
      typename PeakContainerT::ConstIterator it_PosEnd_apex = p.PosBegin(peak_apex_pos); // if peak_apex_pos is correct, this will get the underlying iterator
      typename PeakContainerT::ConstIterator it_PosEnd_r = p.PosEnd(right); // past the end. Do not dereference (might be the true .end())
      for (auto it = it_PosBegin_l; it != it_PosEnd_r; ++it) //OMS_CODING_TEST_EXCLUDE
      {
        // points across the peak
        ++(psm.points_across_baseline);
        if (it->getIntensity() >= 0.5 * peak_height)
        {
          ++(psm.points_across_half_height);
        }
      }
      // positions at peak heights
      psm.start_position_at_5 = findPosAtPeakHeightPercent_(it_PosBegin_l, it_PosEnd_apex, p.end(), peak_height, 0.05, true);
      psm.start_position_at_10 = findPosAtPeakHeightPercent_(it_PosBegin_l, it_PosEnd_apex, p.end(), peak_height, 0.1, true);
      psm.start_position_at_50 = findPosAtPeakHeightPercent_(it_PosBegin_l, it_PosEnd_apex, p.end(), peak_height, 0.5, true);
      psm.end_position_at_5 = findPosAtPeakHeightPercent_(it_PosEnd_apex, it_PosEnd_r, p.end(), peak_height, 0.05, false);
      psm.end_position_at_10 = findPosAtPeakHeightPercent_(it_PosEnd_apex, it_PosEnd_r, p.end(), peak_height, 0.1, false);
      psm.end_position_at_50 = findPosAtPeakHeightPercent_(it_PosEnd_apex, it_PosEnd_r, p.end(), peak_height, 0.5, false);
      // peak widths
      psm.width_at_5 = psm.end_position_at_5 - psm.start_position_at_5;
      psm.width_at_10 = psm.end_position_at_10 - psm.start_position_at_10;
      psm.width_at_50 = psm.end_position_at_50 - psm.start_position_at_50;
      psm.total_width = (p.PosEnd(right) - 1)->getPos() - p.PosBegin(left)->getPos();
      psm.slope_of_baseline = (p.PosEnd(right) - 1)->getIntensity() - p.PosBegin(left)->getIntensity();
      psm.baseline_delta_2_height = psm.slope_of_baseline / peak_height;
      // Source of tailing_factor and asymmetry_factor formulas:
      // USP 40 - NF 35 The United States Pharmacopeia and National Formulary - Supplementary
      psm.tailing_factor = psm.width_at_5 / (2*(peak_apex_pos - psm.start_position_at_5));
      psm.asymmetry_factor = (psm.end_position_at_10 - peak_apex_pos) / (peak_apex_pos - psm.start_position_at_10);
      return psm;
    }
 
 
 
    template <typename PeakContainerConstIteratorT>
    double findPosAtPeakHeightPercent_(
      PeakContainerConstIteratorT it_left,  // must not be past the end
      PeakContainerConstIteratorT it_right, // might be past the end
      PeakContainerConstIteratorT it_end,   // definitely past-the-end
      const double peak_height,
      const double percent,
      const bool is_left_half) const
    {
      // no points in range
      if (it_left == it_end) throw Exception::InvalidRange(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION);
 
      // only one point in range
      if (it_left == it_right) return it_left->getPos();
 
      const double percent_intensity = peak_height * percent;
      PeakContainerConstIteratorT closest;
      if (is_left_half)
      {
        closest = it_left;
        for (
          PeakContainerConstIteratorT it = it_left;
          it < it_right && it->getIntensity() <= percent_intensity;
          closest = it++
        ) {}
      }
      else // right half; search from right to left
      {
        closest = it_right - 1; // make sure we can deference it
        for (
          PeakContainerConstIteratorT it = it_right - 1;
          it >= it_left && it->getIntensity() <= percent_intensity;
          closest = it--
        ) {}
      }
      return closest->getPos();
    }
 
private:
 
    String integration_type_ = INTEGRATION_TYPE_INTENSITYSUM;
    String baseline_type_ = BASELINE_TYPE_BASETOBASE;
 
    bool fit_EMG_;
    EmgGradientDescent emg_;
 
 
    template <typename PeakContainerT>
    const PeakContainerT& EMGPreProcess_(
      const PeakContainerT& pc,
      PeakContainerT& emg_pc,
      double& left,
      double& right
    ) const
    {
      if (fit_EMG_)
      {
        emg_.fitEMGPeakModel(pc, emg_pc, left, right);
        left = emg_pc.front().getPos();
        right = emg_pc.back().getPos();
        return emg_pc;
      }
      return pc;
    }
  };
}