BernNorm.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: Mathias Walzer $
// $Authors: $
// --------------------------------------------------------------------------
//
#pragma once
 
#include <OpenMS/KERNEL/StandardTypes.h>
 
#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
 
#include <map>
 
namespace OpenMS
{
  class OPENMS_DLLAPI BernNorm :
    public DefaultParamHandler
  {
public:
 
    // @name Constructors and Destructors
    BernNorm();
 
    BernNorm(const BernNorm & source);
 
    ~BernNorm() override;
 
    // @name Operators
    // @{
    BernNorm & operator=(const BernNorm & source);
 
    // @name Accessors
    // @{
 
    template <typename SpectrumType>
    void filterSpectrum(SpectrumType & spectrum)
    {
      typedef typename SpectrumType::Iterator Iterator;
      typedef typename SpectrumType::ConstIterator ConstIterator;
 
      c1_ = (double)param_.getValue("C1");
      c2_ = (double)param_.getValue("C2");
      th_ = (double)param_.getValue("threshold");
 
      spectrum.sortByPosition();
 
      // find highest peak and ranking
      double maxint = 0;
      std::map<double, Size> peakranks;
      for (ConstIterator it = spectrum.begin(); it != spectrum.end(); ++it)
      {
        peakranks[it->getIntensity()] = 0;
        if (it->getIntensity() > maxint)
        {
          maxint = it->getIntensity();
        }
      }
      UInt rank = 0;
      for (std::map<double, Size>::reverse_iterator mit = peakranks.rbegin(); mit != peakranks.rend(); ++mit)
      {
        mit->second = ++rank;
      }
 
      // find maxmz i.e. significant (> threshold * maxpeak) peak with highest m/z
      double maxmz = 0;
      for (SignedSize i = spectrum.size() - 1; i >= 0; --i)
      {
        if (spectrum[i].getIntensity() > maxint * th_)
        {
          maxmz = spectrum[i].getMZ();
          break;
        }
      }
 
      // rank
      for (Iterator it = spectrum.begin(); it != spectrum.end(); )
      {
        double newint = c1_ - (c2_ / maxmz) * peakranks[it->getIntensity()];
        if (newint < 0)
        {
          it = spectrum.erase(it);
        }
        else
        {
          it->setIntensity(newint);
          ++it;
        }
      }
      return;
    }
 
    void filterPeakSpectrum(PeakSpectrum & spectrum);
 
    void filterPeakMap(PeakMap & exp);
    //TODO reimplement DefaultParamHandler::updateMembers_()
 
private:
    double c1_;
    double c2_;
    double th_;
 
    // @}
 
  };
 
} // namespace OpenMS