LinearResamplerAlign.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, Luis Jacob Keller, Alen Saric$
// --------------------------------------------------------------------------
 
#pragma once
 
#include <limits>
#include <atomic>
#include <OpenMS/CONCEPT/LogStream.h>
#include <OpenMS/CONCEPT/Macros.h>
#include <OpenMS/CONCEPT/ProgressLogger.h>
#include <OpenMS/DATASTRUCTURES/DefaultParamHandler.h>
#include <OpenMS/KERNEL/MSExperiment.h>
#include <OpenMS/KERNEL/StandardTypes.h>
 
namespace OpenMS
{
  extern OPENMS_DLLAPI std::atomic<bool> suppress_resampling_spacing_warning;
 
  namespace Internal
  {
    class ScopedResamplingWarningSuppression
    {
    public:
      explicit ScopedResamplingWarningSuppression(bool suppress = true) :
        previous_(suppress_resampling_spacing_warning.exchange(suppress))
      {
      }
 
      ~ScopedResamplingWarningSuppression()
      {
        suppress_resampling_spacing_warning.store(previous_);
      }
 
      ScopedResamplingWarningSuppression(const ScopedResamplingWarningSuppression&) = delete;
      ScopedResamplingWarningSuppression& operator=(const ScopedResamplingWarningSuppression&) = delete;
 
    private:
      bool previous_;
    };
  } // namespace Internal
 
  class OPENMS_DLLAPI LinearResamplerAlign :
    public DefaultParamHandler,
    public ProgressLogger
  {
 
  public:
    LinearResamplerAlign() : DefaultParamHandler("LinearResamplerAlign")
    {
      defaults_.setValue("spacing", 0.05, "Spacing of the resampled output peaks.");
      defaults_.setValue("ppm", "false", "Whether spacing is in ppm or Th");
      defaultsToParam_();
    }
 
    template <class PeakContainerT>
    void raster(PeakContainerT& container)
    {
      //return if nothing to do
      if (container.empty()) return;
 
      auto first = container.begin();
      auto last = container.end();
 
      double end_pos = (last - 1)->getPos();
      double start_pos = first->getPos();
      int number_resampled_points = (int)(ceil((end_pos - start_pos) / spacing_ + 1));
 
      std::vector<typename PeakContainerT::PeakType> resampled_peak_container;
      populate_raster_(resampled_peak_container, start_pos, end_pos, number_resampled_points);
 
      raster(container.begin(), container.end(), resampled_peak_container.begin(), resampled_peak_container.end());
 
      container.swap(resampled_peak_container);
    }
 
    template <typename PeakContainerT>
    void raster_align(PeakContainerT& container, double start_pos, double end_pos)
    {
      //return if nothing to do
      if (container.empty()) return;
 
      if (end_pos < start_pos)
      {
        std::vector<typename PeakContainerT::PeakType> empty;
        container.swap(empty);
        return;
      }
 
      auto first = container.begin();
      auto last = container.end();
 
      verifySpacing_(first, last, [](auto x){return x->getPos();});
 
      // get the iterators just before / after the two points start_pos / end_pos
      while (first != container.end() && (first)->getPos() < start_pos) {++first;}
      while (last != first && (last - 1)->getPos() > end_pos) {--last;}
 
      int number_resampled_points = (int)(ceil((end_pos - start_pos) / spacing_ + 1));
 
      std::vector<typename PeakContainerT::PeakType> resampled_peak_container;
      populate_raster_(resampled_peak_container, start_pos, end_pos, number_resampled_points);
 
      raster(first, last, resampled_peak_container.begin(), resampled_peak_container.end());
 
      container.swap(resampled_peak_container);
    }
 
    template <typename PeakTypeIterator, typename ConstPeakTypeIterator>
    void raster(ConstPeakTypeIterator raw_it, ConstPeakTypeIterator raw_end, PeakTypeIterator resampled_begin, PeakTypeIterator resampled_end)
    {
      OPENMS_PRECONDITION(resampled_begin != resampled_end, "Output iterators cannot be identical") // as we use +1
      // OPENMS_PRECONDITION(raw_it != raw_end, "Input iterators cannot be identical")
 
      verifySpacing_(raw_it, raw_end, [](auto x){return x->getPos();});
 
      PeakTypeIterator resample_start = resampled_begin;
 
      // need to get the raw iterator between two resampled iterators of the raw data
      while (raw_it != raw_end && raw_it->getPos() < resampled_begin->getPos())
      {
        resampled_begin->setIntensity(resampled_begin->getIntensity() + raw_it->getIntensity());
        raw_it++;
      }
 
      while (raw_it != raw_end)
      {
        //advance the resample iterator until our raw point is between two resampled iterators
        while (resampled_begin != resampled_end && resampled_begin->getPos() < raw_it->getPos()) {resampled_begin++;}
        if (resampled_begin != resample_start) {resampled_begin--;}
 
        // if we have the last datapoint we break
        if ((resampled_begin + 1) == resampled_end) {break;}
 
        double dist_left =  fabs(raw_it->getPos() - resampled_begin->getPos());
        double dist_right = fabs(raw_it->getPos() - (resampled_begin + 1)->getPos());
 
        // distribute the intensity of the raw point according to the distance to resample_it and resample_it+1
        resampled_begin->setIntensity(resampled_begin->getIntensity() + raw_it->getIntensity() * dist_right / (dist_left + dist_right));
        (resampled_begin + 1)->setIntensity((resampled_begin + 1)->getIntensity() + raw_it->getIntensity() * dist_left / (dist_left + dist_right));
 
        raw_it++;
      }
 
      // add the final intensity to the right
      while (raw_it != raw_end)
      {
        resampled_begin->setIntensity(resampled_begin->getIntensity() + raw_it->getIntensity());
        raw_it++;
      }
    }
 
    template <typename PeakTypeIterator, typename ConstPeakTypeIterator>
    void raster(ConstPeakTypeIterator mz_raw_it, ConstPeakTypeIterator mz_raw_end,
        ConstPeakTypeIterator int_raw_it, ConstPeakTypeIterator int_raw_end,
        ConstPeakTypeIterator mz_resample_it, ConstPeakTypeIterator mz_resample_end,
        PeakTypeIterator int_resample_it, PeakTypeIterator int_resample_end)
    {
      (void)int_raw_end;      // avoid 'unused parameter' compile error
      (void)int_resample_end; // avoid 'unused parameter' compile error
      OPENMS_PRECONDITION(mz_resample_it != mz_resample_end, "Output iterators cannot be identical") // as we use +1
      OPENMS_PRECONDITION(std::distance(mz_resample_it, mz_resample_end) == std::distance(int_resample_it, int_resample_end),
          "Resample m/z and intensity iterators need to cover the same distance")
      OPENMS_PRECONDITION(std::distance(mz_raw_it, mz_raw_end) == std::distance(int_raw_it, int_raw_end),
          "Raw m/z and intensity iterators need to cover the same distance")
      // OPENMS_PRECONDITION(raw_it != raw_end, "Input iterators cannot be identical")
 
      verifySpacing_(mz_raw_it, mz_raw_end, [](auto x){return *x;});
 
      PeakTypeIterator mz_resample_start = mz_resample_it;
 
      // need to get the raw iterator between two resampled iterators of the raw data
      while (mz_raw_it != mz_raw_end && (*mz_raw_it) < (*mz_resample_it) )
      {
        (*int_resample_it) = *int_resample_it + *int_raw_it;
        ++mz_raw_it;
        ++int_raw_it;
      }
 
      while (mz_raw_it != mz_raw_end)
      {
        //advance the resample iterator until our raw point is between two resampled iterators
        while (mz_resample_it != mz_resample_end && *mz_resample_it < *mz_raw_it)
        {
          ++mz_resample_it; ++int_resample_it;
        }
        if (mz_resample_it != mz_resample_start)
        {
          --mz_resample_it; --int_resample_it;
        }
 
        // if we have the last datapoint we break
        if ((mz_resample_it + 1) == mz_resample_end) {break;}
 
        double dist_left =  fabs(*mz_raw_it - *mz_resample_it);
        double dist_right = fabs(*mz_raw_it - *(mz_resample_it + 1));
 
        // distribute the intensity of the raw point according to the distance to resample_it and resample_it+1
        *(int_resample_it) = *int_resample_it + (*int_raw_it) * dist_right / (dist_left + dist_right);
        *(int_resample_it + 1) = *(int_resample_it + 1) + (*int_raw_it) * dist_left / (dist_left + dist_right);
 
        ++mz_raw_it;
        ++int_raw_it;
      }
 
      // add the final intensity to the right
      while (mz_raw_it != mz_raw_end)
      {
        *int_resample_it = *int_resample_it + (*int_raw_it);
        ++mz_raw_it;
        ++int_raw_it;
      }
    }
 
    template <typename PeakTypeIterator>
    void raster_interpolate(PeakTypeIterator raw_it, PeakTypeIterator raw_end, PeakTypeIterator resampled_start, PeakTypeIterator resampled_end)
    {
      // OPENMS_PRECONDITION(resampled_start != resampled_end, "Output iterators cannot be identical")
      OPENMS_PRECONDITION(raw_it != raw_end, "Input iterators cannot be identical") // as we use +1
 
      verifySpacing_(raw_it, raw_end, [](PeakTypeIterator x){return x->getPos();});
 
      PeakTypeIterator raw_start = raw_it;
 
      // need to get the resampled iterator between two iterators of the raw data
      while (resampled_start != resampled_end && resampled_start->getPos() < raw_it->getPos()) {resampled_start++;}
 
      while (resampled_start != resampled_end)
      {
        //advance the raw_iterator until our current point we want to interpolate is between them
        while (raw_it != raw_end && raw_it->getPos() < resampled_start->getPos()) {raw_it++;}
        if (raw_it != raw_start) {raw_it--;}
 
        // if we have the last datapoint we break
        if ((raw_it + 1) == raw_end) {break;}
 
        // use a linear interpolation between raw_it and raw_it+1
        double m = ((raw_it + 1)->getIntensity() - raw_it->getIntensity()) / ((raw_it + 1)->getPos() - raw_it->getPos());
        resampled_start->setIntensity(raw_it->getIntensity() + (resampled_start->getPos() - raw_it->getPos()) * m);
        resampled_start++;
      }
 
    }
 
    void rasterExperiment(PeakMap& exp)
    {
      startProgress(0, exp.size(), "resampling of data");
      for (Size i = 0; i < exp.size(); ++i)
      {
        raster(exp[i]);
        setProgress(i);
      }
      endProgress();
    }
 
 
protected:
 
    double spacing_;
 
    bool ppm_;
 
    void updateMembers_() override
    {
      spacing_ =  param_.getValue("spacing");
      ppm_ =  (bool)param_.getValue("ppm").toBool();
    }
 
    template <typename PeakType>
    void populate_raster_(std::vector<PeakType>& resampled_peak_container,
        double start_pos, double end_pos, int number_resampled_points)
    {
      if (!ppm_)
      {
        // generate the resampled peaks at positions origin+i*spacing_
        resampled_peak_container.resize(number_resampled_points);
        typename std::vector<PeakType>::iterator it = resampled_peak_container.begin();
        for (int i = 0; i < number_resampled_points; ++i)
        {
          it->setPos(start_pos + i * spacing_);
          ++it;
        }
      }
      else
      {
        // generate resampled peaks with ppm distance (not fixed)
        double current_mz = start_pos;
        while (current_mz < end_pos)
        {
          PeakType p;
          p.setIntensity(0);
          p.setPos(current_mz);
          resampled_peak_container.push_back(p);
 
          // increment current_mz
          current_mz += current_mz * (spacing_ / 1e6);
        }
      }
    }
 
 
    template <typename PeakTypeIterator>
    void verifySpacing_(PeakTypeIterator it, PeakTypeIterator end, auto access)
    {
      // ppm_ spacing is relative (parts-per-million) and cannot be compared
      // directly against the absolute neighbour distance computed below.
      if (ppm_) return;
      if (it == end || std::next(it) == end) return;
      double min_dist = std::numeric_limits<double>::infinity();
      double current_dist{};
 
      while (std::next(it) != end)
      {
        current_dist = (access(std::next(it)) - access(it));
        if (min_dist > current_dist) min_dist = current_dist;
        ++it;
      }
 
      if (spacing_ < min_dist && !suppress_resampling_spacing_warning.load())
      {
        OPENMS_LOG_WARN << "Resampling spacing (" << spacing_
                        << ") is smaller than the smallest distance between data points ("
                        << min_dist << "). This approximates the detector dead time and may produce spurious peaks.\n";
      }
    }
 
  };
}