localsignaltonoiseestimator.cpp

Go to the documentation of this file.

// Copyright 2026 Filippo Rusconi
// Inspired by work by Lars Nilse in OpenMS
 
 
/////////////////////// stdlib includes
 
 
/////////////////////// Qt includes
#include <QList>
#include <QDebug>
 
 
/////////////////////// pappsomspp includes
 
 
/////////////////////// Local includes
#include "pappsomspp/core/processing/detection/localsignaltonoiseestimator.h"
 
namespace pappso
{
LocalSignalToNoiseEstimator::LocalSignalToNoiseEstimator(const Parameters &parameters)
  : m_parameters(parameters)
{
}
 
LocalSignalToNoiseEstimator::~LocalSignalToNoiseEstimator()
{
}
 
// Calculate mean & stdev of intensities of a Trace
LocalSignalToNoiseEstimator::GaussianEstimateParams
LocalSignalToNoiseEstimator::estimateGaussian(
  const TraceIterator &iter_first_data_point,
  const TraceIterator &iter_last_data_point) const
{
  int size          = 0;
  // add up
  double v          = 0;
  double m          = 0;
  TraceIterator run = iter_first_data_point;
  while(run != iter_last_data_point)
    {
      m += (*run).y;
      ++size;
      ++run;
    }
  // average
  m = m / size;
 
  // determine variance
  run = iter_first_data_point;
  while(run != iter_last_data_point)
    {
      double tmp(m - (*run).y);
      v += tmp * tmp;
      ++run;
    }
  v = v / ((double)size); // divide by n
 
  GaussianEstimateParams estimate_params = {m, v};
  return estimate_params;
}
 
void
LocalSignalToNoiseEstimator::initialize(const Trace &trace)
{
  computeSignaToNoiseRatio(trace);
}
 
// Calculate signal-to-noise values for all data points given, by using a
//  sliding window approach
void
LocalSignalToNoiseEstimator::computeSignaToNoiseRatio(const Trace &trace)
{
  // First element in the Trace
  TraceIterator iter_first_data_point = trace.begin();
  // Last element in the Trace
  TraceIterator iter_last_data_point  = trace.end();
 
  // Reset counter for sparse windows
  m_sparseWindowPercentage      = 0;
  // Reset counter for histogram overflow
  m_histogramOverflowPercentage = 0;
 
  // Reset the results
  m_signalToNoiseEstimates.clear();
  m_signalToNoiseEstimates.resize(trace.size());
 
  // Maximal range of histogram needs to be calculated first
  if(m_parameters.maxIntMode == AUTO_MAX_BY_STDEV)
    {
      // Use MEAN+auto_m_parameters.maxIntensity*STDEV as threshold
      GaussianEstimateParams gaussian_estimate_global =
        estimateGaussian(iter_first_data_point, iter_last_data_point);
      m_parameters.maxIntensity =
        gaussian_estimate_global.mean +
        std::sqrt(gaussian_estimate_global.variance) * m_parameters.maxIntStdDevFactor;
 
      // qDebug() << "Computed max intensity:" << m_parameters.maxIntensity;
    }
  else if(m_parameters.maxIntMode == AUTO_MAX_BY_PERCENTILE)
    {
      // get value at "m_parameters.maxIntPercentile"th percentile
      // we use a histogram approach here as well.
      if((m_parameters.maxIntPercentile < 0) || (m_parameters.maxIntPercentile > 100))
        {
          qFatal() << "The m_parameters.maxIntPercentile value is not in the [0, 100] range:"
                   << m_parameters.maxIntPercentile;
        }
 
      // Histogram with 100 percent values all initialized to 0.
      QList<int> histogram_of_intensity_distribution(100, 0);
 
      // Find maximum of current scan
      auto max_intensity_data_point_iterator = std::max_element(
        trace.begin(), trace.end(), [](const DataPoint &a, const DataPoint &b) {
          return a.y > b.y;
        });
 
      double max_intensity = max_intensity_data_point_iterator->y;
 
      // qDebug() << "Max intensity value in the trace:" << max_intensity;
 
      double histogram_bin_size = max_intensity / 100;
 
      // Fill histogram
      for(const DataPoint &data_point : trace)
        {
          // Increment count at histogram bin that is the bin corresponding
          // to the iterated data point intensity in the percentage scale above.
          int bin = (int)((data_point.y - 1) / histogram_bin_size);
          ++histogram_of_intensity_distribution[bin];
        }
 
      // Add up element counts in histogram until ?th percentile is reached
      // For example if percentile is 95 percent, then the value elements_below_max_int_percentile
      // below is the number of data points that correspond to 95 % of the whole number of data
      // points in the spectrum.
      int elements_below_max_int_percentile =
        (int)(m_parameters.maxIntPercentile * trace.size() / 100);
 
      int elements_seen = 0;
      int i             = -1;
 
      TraceIterator the_iterator = iter_first_data_point;
 
      // This loop is equivalent to iterating in the bins of the histogram,
      // and for each bin sum the corresponding count to elements_seen.
      // When elements_seen reaches 95% of the number of DataPoint instances
      // in the mass spectrum, stop and record the bin number (i).
      while(the_iterator != iter_last_data_point &&
            elements_seen < elements_below_max_int_percentile)
        {
          ++i;
          elements_seen += histogram_of_intensity_distribution[i];
          ++the_iterator;
        }
 
      // Because each histogram_bin_size accounts for a hundredth of the max intensity of the
      // whole mass spectra data, returning the i * histogram_bin_size is equivalent to
      // returning the greatest intensity value encountered in the 95% of the mass spectrum.
      // Half a bin is added to center the returned value in the bin.
      m_parameters.maxIntensity = (((double)i) + 0.5) * histogram_bin_size;
 
      // qDebug() << "Computed max intensity:" << m_parameters.maxIntensity;
    }
  else // if (m_parameters.maxIntMode == MANUAL_BY_VALUE)
    {
      if(m_parameters.maxIntensity <= 0)
        {
          qFatal() << "The method to determine the maximal intensity is set to MANUAL"
                   << "However, the m_parameters.maxIntensity value is <= 0.";
        }
 
      // qDebug() << "The max intensity was set manually:" << m_parameters.maxIntensity;
    }
 
  if(m_parameters.maxIntensity < 0)
    {
      qFatal() << "The m_parameters.maxIntensity "
                  "value should be positive! ";
    }
 
  TraceIterator window_pos_center       = iter_first_data_point;
  TraceIterator window_pos_border_left  = iter_first_data_point;
  TraceIterator window_pos_border_right = iter_first_data_point;
 
  double window_half_size = m_parameters.windowSize / 2;
  // At least size of 1 for intensity bins
  double bin_size         = std::max(1.0, m_parameters.maxIntensity / m_parameters.binCount);
  int bin_count_minus_one = m_parameters.binCount - 1;
 
  // Allocate the histogram with the user-configured number of bins
  std::vector<int> histogram(m_parameters.binCount, 0);
 
  // This vector will store for each bin, the value in it.
  std::vector<double> bin_values(m_parameters.binCount, 0);
 
  // Calculate average intensity that is represented by a bin
  for(int bin = 0; bin < m_parameters.binCount; bin++)
    {
      histogram[bin] = 0;
      bin_values[bin] = (bin + 0.5) * bin_size;
    }
 
  // This will store the bin in which a datapoint would fall.
  int to_bin = 0;
 
  // index of bin where the median is located
  int median_bin        = 0;
  // additive number of elements from left to x in histogram
  int element_inc_count = 0;
 
  // tracks elements in current window, which may vary because of unevenly
  // spaced data
  int elements_in_window = 0;
  // number of windows
  int window_count       = 0;
 
  // number of elements where we find the median
  int element_in_window_half = 0;
 
  double noise; // noise value of a datapoint
 
  // Iterate in the mass spectrum and at each data point, determine the signal-to-noise ratio
  // and store it. The signal-to-noise ratio is computed over a rolling window of
  // a configured size that is centered on the iterated data point.
  // At each data point in the mass spectrum (that is, at
  // iterator window_pos_center), update the histogram data for the rolling window.
  // The data in the histogram will serve as the basis to determine the noise.
 
  while(window_pos_center != iter_last_data_point)
    {
      // Each time we iterate in a new mass spectrum DataPoint, that is, each
      // time we increment window_pos_center, we need to shift the rolling window
      // to center it at the new mass spectrum data point.
 
      // Erase all elements from histogram that will leave the rolling window on the
      // LEFT side
      while((*window_pos_border_left).x <
            (*window_pos_center).x - window_half_size)
        {
          to_bin = std::max(
            std::min<int>((int)((*window_pos_border_left).y / bin_size),
                          bin_count_minus_one),
            0);
 
          // Remove the current window's left half's DataPoint from
          // the histogram counts at bin corresponding to its intensity.
          --histogram[to_bin];
          // Remove the current window's left half's DataPoint from
          // the total count of DataPoint instances in the rolling window.
          --elements_in_window;
          // Now go on and shift the rolling window's left border iterator one position
          // to the right.
          ++window_pos_border_left;
 
          // Do this until the iterator of the left border of the rolling window points
          // to a mass spectrum DataPoint that has its m/z value matching the
          // current rolling window center iterator m/z value minus half of the
          // rolling window width.
        }
 
      // Add all elements to histogram that will enter the rolling window on the RIGHT
      // side
      while((window_pos_border_right != iter_last_data_point) &&
            ((*window_pos_border_right).x <=
             (*window_pos_center).x + window_half_size))
        {
          // std::cerr << (*window_pos_border_right).y << " " <<
          // bin_size << " " << m_parameters.binCountminus_1 << std::endl;
          to_bin = std::max(
            std::min<int>((int)((*window_pos_border_right).y / bin_size),
                          bin_count_minus_one),
            0);
          ++histogram[to_bin];
          ++elements_in_window;
          ++window_pos_border_right;
        }
 
      if(elements_in_window < m_parameters.minRequiredElements)
        {
          noise = m_parameters.noiseValueForEmptyWindow;
          ++m_sparseWindowPercentage;
        }
      else
        {
          // Find bin i where ceil[elements_in_window/2] <= sum_c(0..i){
          // histogram[c] }
          median_bin             = -1;
          element_inc_count      = 0;
          element_in_window_half = (elements_in_window + 1) / 2;
          while(median_bin < bin_count_minus_one &&
                element_inc_count < element_in_window_half)
            {
              ++median_bin;
              element_inc_count += histogram[median_bin];
            }
 
          // Increase the error count
          if(median_bin == bin_count_minus_one)
            {
              ++m_histogramOverflowPercentage;
            }
 
          // Just avoid division by 0
          noise = std::max(1.0, bin_values[median_bin]);
        }
 
      // Store result
      m_signalToNoiseEstimates[window_count] = (*window_pos_center).y / noise;
      // qDebug() << "Added new signal to noise estimate:"
      //          << m_signalToNoiseEstimates.at(window_count);
 
      // Advance the window center by one datapoint
      ++window_pos_center;
      ++window_count;
 
    } // end while
 
  // qDebug() << "Done filling in the signal to noise estimates, the are"
  //          << m_signalToNoiseEstimates.size() << "estimates.";
 
  m_sparseWindowPercentage = m_sparseWindowPercentage * 100 / window_count;
  m_histogramOverflowPercentage =
    m_histogramOverflowPercentage * 100 / window_count;
 
  // warn if percentage of sparse windows is above 20%
  if(m_sparseWindowPercentage > 20)
    {
      qWarning() << m_sparseWindowPercentage
                 << "% of all windows were sparse. You should consider increasing"
                 << "'m_parameters.windowSize' or decrease 'm_parameters.minRequiredElements'";
    }
 
  // warn if percentage of possibly wrong median estimates is above 1%
  if(m_histogramOverflowPercentage > 1)
    {
      qWarning() << m_histogramOverflowPercentage
                 << "% of all Signal-to-Noise estimates are too high, because the "
                    "median was found in the rightmost histogram-bin. "
                 << "You should consider increasing 'm_parameters.maxIntensityy' (and maybe "
                    "'m_parameters.binCount' with it, to keep bin width reasonable)";
    }
}
 
// Return to signal/noise estimate for date point @p index
// @note you will get a warning to stderr if more than 20% of the
//       noise estimates used sparse windows
double
LocalSignalToNoiseEstimator::getSignalToNoiseRatio(qsizetype index) const
{
  Q_ASSERT(index < m_signalToNoiseEstimates.size());
  return m_signalToNoiseEstimates[index];
}
 
} // namespace pappso