SignalToNoiseEstimatorMedianRapid.h

Go to the documentation of this file.

// Copyright (c) 2002-present, OpenMS Inc. -- 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/CONCEPT/Types.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)
      {
        return result_windows_even[windowIndex_((mz - mz_start) / window_length, result_windows_even.size())];
      }
 
      double get_noise_odd (double mz)
      {
        return result_windows_odd[windowIndex_((mz - mz_start + window_length / 2.0) / window_length, result_windows_odd.size())];
      }
 
    private:
      static Size windowIndex_(double window, Size n_windows)
      {
        if (!(window > -1.0 && window < static_cast<double>(n_windows)))
        {
          throwWindowOutOfRange_(window, n_windows);
        }
        // convert via a signed integer: a single instruction on x86, unlike double -> unsigned
        return static_cast<Size>(static_cast<SignedSize>(window));
      }
 
      [[noreturn]] static void throwWindowOutOfRange_(double window, Size n_windows);
    };
 
    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
      checkInputArrays_(mz_array, int_array);
      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:
 
    static void checkInputArrays_(const std::vector<double>& mz_array, const std::vector<double>& int_array);
 
    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