IonImageExtraction.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: Timo Sachsenberg $
// $Authors: Timo Sachsenberg, Patrick Boschmann, Aditya Sarna $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/CONCEPT/Exception.h>
#include <OpenMS/CONCEPT/Types.h>
#include <OpenMS/DATASTRUCTURES/StringUtils.h>
#include <OpenMS/IMAGING/IonImage.h>
#include <OpenMS/IMAGING/MSImagingGeometry.h>
#include <OpenMS/KERNEL/MSSpectrum.h>
 
#include <cmath>
#include <vector>
 
namespace OpenMS
{
namespace Internal
{
 
template <typename GetSpectrum>
IonImage extractIonImage(const MSImagingGeometry& geom,
                         double mz, double tolerance_ppm,
                         const std::vector<Size>& pixel_indices,
                         Size n_spectra,
                         GetSpectrum&& get_spectrum)
{
  if (!std::isfinite(mz) || !std::isfinite(tolerance_ppm) || mz < 0.0 || tolerance_ppm < 0.0)
  {
    throw Exception::InvalidValue(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
                                  "mz and tolerance_ppm must be finite and non-negative",
                                  "mz=" + StringUtils::toStr(mz) + ", tolerance_ppm=" + StringUtils::toStr(tolerance_ppm));
  }
 
  const double dm    = mz * tolerance_ppm * 1e-6;
  const double mz_lo = mz - dm;
  const double mz_hi = mz + dm;
 
  IonImage image(geom.getWidth(), geom.getHeight());
  image.setMzRange(RangeMZ(mz_lo, mz_hi));
 
  const auto& pixels = geom.getPixels();
  for (Size i : pixel_indices)
  {
    const auto& p = pixels[i];
    if (p.spectrum_index >= n_spectra)
    {
      throw Exception::InvalidValue(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
                                    "Pixel references missing spectrum",
                                    StringUtils::toStr(p.spectrum_index));
    }
    // const auto& spec: for a lambda returning const MSSpectrum& (in-memory) this
    // binds the reference with no copy; for a lambda returning MSSpectrum by value
    // (on-disc) this lifetime-extends the temporary for the duration of the loop body.
    const auto& spec = get_spectrum(p.spectrum_index);
    double sum = 0.0;
    for (auto it = spec.MZBegin(mz_lo); it != spec.MZEnd(mz_hi); ++it)
    {
      sum += static_cast<double>(it->getIntensity());
    }
    image.setIntensity(p.x, p.y, sum);
  }
  return image;
}
 
} // namespace Internal
} // namespace OpenMS