signaltonoiseestimatormedian.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"
 
namespace pappso
{
SignalToNoiseEstimatorMedian::SignalToNoiseEstimatorMedian()
{
}
 
SignalToNoiseEstimatorMedian::~SignalToNoiseEstimatorMedian()
{
}
 
void
SignalToNoiseEstimatorMedian::pick(const Trace &trace,
                        Trace &picked_peaks,
                        QList<PeakRange> &peak_ranges)
{
  picked_peaks.clear();
 
  qsizetype trace_size = trace.size();
 
  if(trace_size < 5)
    return;
 
  for(qsizetype iter = 2; iter < trace_size - 2; ++iter)
    {
      double center_data_point_mz  = trace.at(iter).x,
             center_data_point_int = trace.at(iter).y;
      double left_data_point_mz    = trace.at(iter - 1).x,
             left_data_point_int   = trace.at(iter - 1).y;
      double right_data_point_mz   = trace.at(iter + 1).x,
             right_data_point_int  = trace.at(iter + 1).y;
 
      // do not interpolate when the left or right support is a zero-data-point
      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 act_snt = 0.0, act_snt_l1 = 0.0, act_snt_r1 = 0.0;
 
      if(m_signalToNoise > 0.0)
        {
          // FIXME
        }
 
      // 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) &&
         (act_snt >= m_signalToNoise) && (act_snt_l1 >= m_signalToNoise) &&
         (act_snt_r1 >= m_signalToNoise))
        {
          // special case: if a peak core is surrounded by more intense
          // satellite peaks (indicates oscillation rather than
          // real peaks) -> remove
 
          double act_snt_l2 = 0.0, act_snt_r2 = 0.0;
 
          if(m_signalToNoise > 0.0)
            {
              // FIXME
            }
 
          // 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) &&
             (act_snt_l2 >= m_signalToNoise) && (act_snt_r2 >= m_signalToNoise))
            {
              ++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);
 
          // 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 act_snt_lk = 0.0;
 
              if(m_signalToNoise > 0.0)
                {
                  // FIXME
                }
 
              if(act_snt_lk >= m_signalToNoise)
                {
                  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);
 
          // 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 act_snt_rk = 0.0;
 
              if(m_signalToNoise > 0.0)
                {
                  // FIXME
                }
 
              if(act_snt_rk >= m_signalToNoise)
                {
                  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;
        }
    }
}
 
 
} // namespace pappso