SignalToNoiseEstimatorMedianRapid.h

Go to the documentation of this file.

// Copyright (c) 2002-2023, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Hannes Roest $
// $Authors: Hannes Roest $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/INTERFACES/DataStructures.h>
#include <OpenMS/INTERFACES/ISpectrumAccess.h>
#include <cassert>
#include <vector>
 
namespace OpenMS
{
 
  class OPENMS_DLLAPI SignalToNoiseEstimatorMedianRapid
  {
    double window_length_;
 
public:
 
    struct OPENMS_DLLAPI NoiseEstimator
    {
      int nr_windows;
      double mz_start;
      double window_length;
      std::vector<double> result_windows_even;
      std::vector<double> result_windows_odd;
 
      NoiseEstimator() {}
 
      NoiseEstimator(double nr_windows_, double mz_start_, double win_len_) :
        nr_windows(nr_windows_),
        mz_start(mz_start_),
        window_length(win_len_),
        result_windows_even(nr_windows_),
        result_windows_odd(nr_windows_+1)
      {}
 
      double get_noise_value (double mz)
      {
        // Take the average of the two stored values
        // Avoid division by 0 (since most clients will divide by the noise value)
        return std::max(1.0, (get_noise_even(mz)+get_noise_odd(mz))/2.0 );
      }
 
      double get_noise_even (double mz)
      {
        // PRECONDITION
        int window_nr = (int)((mz - mz_start)/window_length);
        assert(window_nr >= 0);
        assert(window_nr < (int)result_windows_even.size());
 
        double noise = result_windows_even[window_nr];
        return noise;
      }
 
      double get_noise_odd (double mz)
      {
        // PRECONDITION
        int window_nr = (int)((mz - mz_start + window_length/2.0)/window_length);
        assert(window_nr >= 0);
        assert(window_nr < (int)result_windows_odd.size());
 
        double noise = result_windows_odd[window_nr];
        return noise;
      }
    };
 
    SignalToNoiseEstimatorMedianRapid(double window_length) :
      window_length_(window_length)
    {
    }
 
    inline NoiseEstimator estimateNoise(OpenMS::Interfaces::SpectrumPtr spectrum)
    {
      return estimateNoise(spectrum->getMZArray()->data, spectrum->getIntensityArray()->data);
    }
 
    inline NoiseEstimator estimateNoise(OpenMS::Interfaces::ChromatogramPtr chrom)
    {
      return estimateNoise(chrom->getTimeArray()->data, chrom->getIntensityArray()->data);
    }
 
    NoiseEstimator estimateNoise(const std::vector<double>& mz_array, const std::vector<double>& int_array)
    {
      // PRECONDITION
      assert(mz_array.size() == int_array.size());
      assert(mz_array.size() > 2);
 
      int nr_windows = (int)((mz_array[mz_array.size()-1] - mz_array[0])/window_length_) + 1;
      NoiseEstimator eval(nr_windows, mz_array[0], window_length_);
 
      // Compute even windows
      computeNoiseInWindows_(mz_array, int_array, eval.result_windows_even, mz_array[0]);
      // Compute odd windows
      computeNoiseInWindows_(mz_array, int_array, eval.result_windows_odd, mz_array[0] - window_length_/2.0);
 
      return eval;
    }
 
private:
 
    void computeNoiseInWindows_(const std::vector<double>& mz_array, std::vector<double> int_array, std::vector<double> & result, double mz_start);
 
    double computeMedian_(std::vector<double>::iterator & first, std::vector<double>::iterator & last);
 
  };
 
} // namespace OpenMS