Loading...
Searching...
No Matches
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
pappso::CubicSplineModel
Definition cubicsplinemodel.h:24
pappso::Trace
A simple container of DataPoint instances.
Definition trace.h:152
cubicsplinemodel.h
highrespeakpicker.h
pappso
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition aa.cpp:39
pappso::spline_bisection
void spline_bisection(const CubicSplineModel &spline_model, double const mz_at_left, double const mz_at_right, double &center_peak_mz, double &center_peak_intensity, double const threshold)
Definition cubicsplinemodel.cpp:241
pappso::PeptideIonCter::y
@ y
Definition peptide.h:87
pappso::DataPoint
Definition datapoint.h:22
pappso::DataPoint::x
pappso_double x
Definition datapoint.h:24