localsignaltonoiseestimator.h

Go to the documentation of this file.

// Copyright 2026 Filippo Rusconi
// Largely inspired by code in OpenMS, by Chris Bielow
 
#pragma once
 
/////////////////////// stdlib includes
 
 
/////////////////////// Qt includes
#include <QList>
#include <QMap>
 
/////////////////////// pappsomspp includes
#include "pappsomspp/core/trace/trace.h"
 
 
/////////////////////// Local includes
#include "pappsomspp/export-import-config.h"
 
namespace pappso
{
 
using TraceIterator = std::vector<DataPoint>::const_iterator;
 
class PMSPP_LIB_DECL LocalSignalToNoiseEstimator
{
 
  public:
  // We will need that in the Parameters below.
  enum IntensityThresholdCalculation
  {
    MANUAL_BY_VALUE = 0x0,
    AUTO_MAX_BY_STDEV,
    AUTO_MAX_BY_PERCENTILE
  };
 
  struct Parameters
  {
    // Maximal intensity considered during binning (values above get discarded)
    // This ensures that absurdly-intense peaks are not used to compute the
    // stats about the maximum intensity of the trace.
    double maxIntensity = -1;
 
    // parameter for initial automatic estimation of "max_intensity_": a stdev
    // multiplier
    double maxIntStdDevFactor = 3;
 
    // parameter for initial automatic estimation of "max_intensity_" percentile
    // or a stde
    double maxIntPercentile = 95;
 
    // Which method is used for estimating m_maxIntensity "max intensity of
    // trace". valid are MANUAL, AUTOMABYSTDEV or AUTOMAXBYPERCENTILE
    int maxIntMode = IntensityThresholdCalculation::AUTO_MAX_BY_STDEV;
 
    // range of data points which belong to a window in Thomson
    double windowSize = 200;
 
    // number of bins in the histogram
    int binCount = 30;
 
    // minimal number of elements a window needs to cover to be used
    int minRequiredElements = 10;
 
    // used as noise value for windows which cover less than
    // "min_required_elements_" use a very high value if you want to get a low
    // S/N result
    double noiseValueForEmptyWindow = std::pow(10.0, 20);
 
    Parameters() {};
    Parameters(double max_intensity,
               double max_int_std_dev_factor,
               double max_int_percentile,
               int max_int_mode,
               double window_size,
               int bin_count,
               int min_required_elements,
               double noise_value_for_empty_window)
    {
      maxIntensity             = max_intensity;
      maxIntStdDevFactor       = max_int_std_dev_factor;
      maxIntPercentile         = max_int_percentile;
      maxIntMode               = max_int_mode;
      windowSize               = window_size;
      binCount                 = bin_count;
      minRequiredElements      = min_required_elements;
      noiseValueForEmptyWindow = noise_value_for_empty_window;
    }
 
    Parameters(const Parameters &other)
    {
      maxIntensity             = other.maxIntensity;
      maxIntStdDevFactor       = other.maxIntStdDevFactor;
      maxIntPercentile         = other.maxIntPercentile;
      maxIntMode               = other.maxIntMode;
      windowSize               = other.windowSize;
      binCount                 = other.binCount;
      minRequiredElements      = other.minRequiredElements;
      noiseValueForEmptyWindow = other.noiseValueForEmptyWindow;
    }
  };
 
  explicit LocalSignalToNoiseEstimator(const Parameters &parameters);
  virtual ~LocalSignalToNoiseEstimator();
 
  struct GaussianEstimateParams
  {
    // mean of estimated Gaussian
    double mean;
    // variance of estimated Gaussian
    double variance;
  };
 
  GaussianEstimateParams
  estimateGaussian(const TraceIterator &iter_first_data_point,
                   const TraceIterator &iter_last_data_point) const;
 
  void initialize(const Trace &trace);
  void computeSignaToNoiseRatio(const Trace &trace);
  double getSignalToNoiseRatio(qsizetype index) const;
 
  protected:
  Parameters m_parameters;
 
  // counter for sparse windows
  double m_sparseWindowPercentage;
  // counter for histogram overflow
  double m_histogramOverflowPercentage;
 
  QList<double> m_signalToNoiseEstimates;
};
 
 
} // namespace pappso