SignalToNoiseEstimatorMedianRapid.h

Go to the documentation of this file.

// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2017.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Hannes Roest $
// $Authors: Hannes Roest $
// --------------------------------------------------------------------------

#ifndef OPENMS_FILTERING_NOISEESTIMATION_SIGNALTONOISEESTIMATORMEDIANRAPID_H
#define OPENMS_FILTERING_NOISEESTIMATION_SIGNALTONOISEESTIMATORMEDIANRAPID_H

#include <OpenMS/INTERFACES/DataStructures.h>
#include <OpenMS/INTERFACES/ISpectrumAccess.h>
#include <assert.h>
#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

#endif // OPENMS_FILTERING_NOISEESTIMATION_SIGNALTONOISEESTIMATORMEDIANRAPID_H