highrespeakpicker.cpp

Go to the documentation of this file.

// Copyright 2026 Filippo Rusconi
// Inspired by work by Lars Nilse in OpenMS
 
 
/////////////////////// stdlib includes
#include <limits>
 
 
/////////////////////// Qt includes
#include <QList>
#include <QDebug>
 
 
/////////////////////// pappsomspp includes
 
 
/////////////////////// Local includes
#include "pappsomspp/core/processing/detection/highrespeakpicker.h"
#include "pappsomspp/core/processing/detection/cubicsplinemodel.h"
#include "pappsomspp/core/processing/detection/localsignaltonoiseestimator.h"
 
namespace pappso
{
HighResPeakPicker::HighResPeakPicker(Parameters &parameters) : m_parameters(parameters)
{
}
 
HighResPeakPicker::~HighResPeakPicker()
{
}
 
void
HighResPeakPicker::pick(const Trace &trace,
                        Trace &picked_peaks,
                        QList<PeakRange> &peak_ranges)
{
  picked_peaks.clear();
 
  qsizetype trace_size = trace.size();
 
  if(trace_size < 5)
    return;
 
  // signal-to-noise estimation
  LocalSignalToNoiseEstimator::Parameters signal_to_noise_estimator_parameters;
  LocalSignalToNoiseEstimator local_signal_estimator(signal_to_noise_estimator_parameters);
 
  if(m_parameters.signalToNoise > 0.0)
    {
      local_signal_estimator.initialize(trace);
    }
 
  // Iterate in the trace and search for local maxima, that is, roughly,
  // DataPoint instances that are higher than previous and next point.
  for(qsizetype iter = 2; iter < trace_size - 2; ++iter)
    {
      double center_data_point_mz  = trace.at(iter).x;
      double center_data_point_int = trace.at(iter).y;
      double left_data_point_mz    = trace.at(iter - 1).x;
      double left_data_point_int   = trace.at(iter - 1).y;
      double right_data_point_mz   = trace.at(iter + 1).x;
      double right_data_point_int  = trace.at(iter + 1).y;
 
      // Do not interpolate when the left or right DataPoint is zero-intensity.
      if(std::fabs(left_data_point_int) <
         std::numeric_limits<double>::epsilon())
        {
          continue;
        }
      if(std::fabs(right_data_point_int) <
         std::numeric_limits<double>::epsilon())
        {
          continue;
        }
 
      double current_snt_estimate         = 0.0;
      double current_snt_estimate_left_1  = 0.0;
      double current_snt_estimate_right_1 = 0.0;
 
      if(m_parameters.signalToNoise > 0.0)
        {
          current_snt_estimate =
            local_signal_estimator.getSignalToNoiseRatio(iter);
 
          // qDebug() << "For current center point at iter:" << iter
          //          << "with mz:" << center_data_point_mz << "of intensity"
          //          << center_data_point_int
          //          << "the snt estimate is:" << current_snt_estimate;
 
          current_snt_estimate_left_1 =
            local_signal_estimator.getSignalToNoiseRatio(iter - 1);
          current_snt_estimate_right_1 =
            local_signal_estimator.getSignalToNoiseRatio(iter + 1);
        }
 
      // look for peak cores meeting MZ and intensity/SNT criteria
      if((center_data_point_int > left_data_point_int) &&
         (center_data_point_int > right_data_point_int) &&
         (current_snt_estimate >= m_parameters.signalToNoise) &&
         (current_snt_estimate_left_1 >= m_parameters.signalToNoise) &&
         (current_snt_estimate_right_1 >= m_parameters.signalToNoise))
        {
          // First check if the currently iterated center center point
          // is surrounded by more intense peaks on its left-left or on its
          // right-right, because that might indicate signal instability instead
          // of an actual maximum. In this case, we do not consider the data
          // point.
 
          double current_snt_estimate_left_2  = 0.0;
          double current_snt_estimate_right_2 = 0.0;
 
          if(m_parameters.signalToNoise > 0.0)
            {
              current_snt_estimate_left_2 =
                local_signal_estimator.getSignalToNoiseRatio(iter - 2);
              current_snt_estimate_right_2 =
                local_signal_estimator.getSignalToNoiseRatio(iter + 2);
            }
 
          // Checking signal-to-noise
          if((iter > 1) && (iter + 2 < static_cast<qsizetype>(trace.size())) &&
             (left_data_point_int < trace.at(iter - 2).y) &&
             (right_data_point_int < trace.at(iter + 2).y) &&
             (current_snt_estimate_left_2 >= m_parameters.signalToNoise) &&
             (current_snt_estimate_right_2 >= m_parameters.signalToNoise))
            {
              // If the point left of the left data point (with respect to
              // current center data point that was found to be a local maximum)
              // and if the point right of the right data point (again, with
              // respect to current center data point that was found to be a
              // local maximum) and both of these left_left and right_right
              // datapoints are above the signal to noise ratio, then we are
              // experiencing signal instability.
              ++iter;
              continue;
            }
 
          QMap<double, double> peak_raw_data;
 
          peak_raw_data[center_data_point_mz] = center_data_point_int;
          peak_raw_data[left_data_point_mz]   = left_data_point_int;
          peak_raw_data[right_data_point_mz]  = right_data_point_int;
 
          // peak core found, now extend it
          // to the left
          qsizetype k = 2;
 
          // No need to extend peak if previous intensity was zero
          bool previous_zero_left      = false;
          qsizetype missing_left_count = 0;
          // Index of the left boundary for the spline interpolation
          qsizetype left_boundary      = iter - 1;
 
          while((k <= iter) && // prevent underflow
                (iter - k + 1 > 0) && !previous_zero_left &&
                (trace.at(iter - k).y <= peak_raw_data.begin().value()))
            {
              double current_snt_estimate_left_k = 0.0;
 
              if(m_parameters.signalToNoise > 0.0)
                {
                  current_snt_estimate_left_k =
                    local_signal_estimator.getSignalToNoiseRatio(iter - k);
                }
 
              if(current_snt_estimate_left_k >= m_parameters.signalToNoise)
                {
                  peak_raw_data[trace.at(iter - k).x] = trace.at(iter - k).y;
                }
              {
                ++missing_left_count;
                if(missing_left_count <= m_parameters.missingPeakCount)
                  {
                    peak_raw_data[trace.at(iter - k).x] = trace.at(iter - k).y;
                  }
              }
 
              previous_zero_left = (trace.at(iter - k).y == 0);
              left_boundary      = iter - k;
              ++k;
            }
 
          // to the right
          k = 2;
 
          // No need to extend peak if previous intensity was zero
          bool previous_zero_right      = false;
          qsizetype missing_right_count = 0;
          // Index of the right boundary for the spline interpolation
          qsizetype right_boundary(iter + 1);
 
          while((iter + k < static_cast<qsizetype>(trace.size())) &&
                !previous_zero_right &&
                (trace.at(iter + k).y <= peak_raw_data.last()))
            {
              double current_snt_estimate_right_k = 0.0;
 
              if(m_parameters.signalToNoise > 0.0)
                {
                  current_snt_estimate_right_k =
                    local_signal_estimator.getSignalToNoiseRatio(iter + k);
                }
 
              if(current_snt_estimate_right_k >= m_parameters.signalToNoise)
                {
                  peak_raw_data[trace.at(iter + k).x] = trace.at(iter + k).y;
                }
              else
                {
                  ++missing_right_count;
                  if(missing_right_count <= m_parameters.missingPeakCount)
                    {
                      peak_raw_data[trace.at(iter + k).x] =
                        trace.at(iter + k).y;
                    }
                }
 
              previous_zero_right = (trace.at(iter + k).y == 0);
              right_boundary      = iter + k;
              ++k;
            }
 
          // skip if the minimal number of 3 points for fitting is not reached
          if(peak_raw_data.size() < 3)
            {
              continue;
            }
 
          // qDebug() << "peak_raw_data size:" << peak_raw_data.size();
 
          CubicSplineModel cubic_spline_model(peak_raw_data);
          // qDebug() << "Stack-allocated:" << &cubic_spline_model;
 
 
          // qDebug() << "The cubic spline model has "
          //          << cubic_spline_model.getKnots().size() << "knots.";
 
          // calculate maximum by evaluating the spline's 1st derivative
          // (bisection method)
          double max_peak_mz  = center_data_point_mz;
          double max_peak_int = center_data_point_int;
          double threshold    = 1e-6;
          spline_bisection(cubic_spline_model,
                           left_data_point_mz,
                           right_data_point_mz,
                           max_peak_mz,
                           max_peak_int,
                           threshold);
 
          // save picked peak into output spectrum
          DataPoint data_point;
 
          PeakRange peak_range;
 
          data_point.x        = max_peak_mz;
          data_point.y        = max_peak_int;
          peak_range.mz_start = trace.at(left_boundary).x;
          peak_range.mz_stop  = trace.at(right_boundary).x;
 
          // qDebug() << "Now appending new data point:" <<
          // data_point.toString();
 
          picked_peaks.append(data_point);
 
          peak_ranges.push_back(peak_range);
 
          // jump over profile data points that have been considered already
          iter += k - 1;
        }
    }
 
  // qDebug() << "Now returning" << picked_peaks.size()
  //          << "picked centroids:" << picked_peaks.toString();
}
 
 
} // namespace pappso