timsframetype1.cpp

Go to the documentation of this file.

/**
 * \file pappsomspp/vendors/tims/timsframetype1.cpp
 * \date 3/10/2021
 * \author Olivier Langella
 * \brief handle a single Bruker's TimsTof frame type 1 compression
 */
 
/*******************************************************************************
 * Copyright (c) 2021 Olivier Langella <Olivier.Langella@u-psud.fr>.
 *
 * This file is part of the PAPPSOms++ library.
 *
 *     PAPPSOms++ is free software: you can redistribute it and/or modify
 *     it under the terms of the GNU General Public License as published by
 *     the Free Software Foundation, either version 3 of the License, or
 *     (at your option) any later version.
 *
 *     PAPPSOms++ is distributed in the hope that it will be useful,
 *     but WITHOUT ANY WARRANTY; without even the implied warranty of
 *     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *     GNU General Public License for more details.
 *
 *     You should have received a copy of the GNU General Public License
 *     along with PAPPSOms++.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/
 
 
#include "pappsomspp/core/vendors/tims/timsframetype1.h"
#include "pappsomspp/core/pappsoexception.h"
#include "pappsomspp/core/exception/exceptionoutofrange.h"
#include "pappsomspp/core/exception/exceptionnotimplemented.h"
#include <QDebug>
#include <QObject>
#include <liblzf/lzf.h>
#include <cerrno>
 
 
namespace pappso
{
TimsFrameType1::TimsFrameType1(std::size_t timsId,
                               quint32 scanNum,
                               char *p_bytes,
                               std::size_t len)
  : TimsFrame(timsId, scanNum)
{
  qDebug() << timsId;
  m_binaryData.resize(len * 2);
 
  if(p_bytes != nullptr)
    {
      qDebug() << timsId;
      copyAndLzfDecompress(p_bytes, len);
      qDebug() << timsId;
    }
  else
    {
      if(m_scanCount == 0)
        {
 
          throw pappso::PappsoException(
            QObject::tr(
              "TimsFrameType1::TimsFrameType1(%1,%2,nullptr,%3) FAILED")
              .arg(m_frameId)
              .arg(m_scanCount)
              .arg(len));
        }
    }
}
 
TimsFrameType1::TimsFrameType1(const TimsFrameType1 &other) : TimsFrame(other)
{
}
 
TimsFrameType1::~TimsFrameType1()
{
}
 
 
void
TimsFrameType1::copyAndLzfDecompress(const char *src, std::size_t len)
{
 
  qDebug() << " m_scanNumber=" << m_scanCount << " len=" << len;
  // the start position offset for each scan and the length of the last scan
  // copy first m_scanNumber*4 bytes in qbyte array
  std::size_t count = (m_scanCount + 2) * 4;
 
  qDebug() << " count=" << count;
  if(m_binaryData.size() < (long)(count + count))
    {
      qDebug() << " m_timsDataFrame.size()=" << m_binaryData.size();
      m_binaryData.resize(count + count);
    }
 
  /*
    std::size_t decompressed_size =
      lzfDecompressScan(src + 3687 - 8,
                        9,
                        m_timsDataFrame.data() + 3660,
                        m_timsDataFrame.size() - 3660);
 
    qDebug() << "decompressed_size=" << decompressed_size;
    */
  // memcpy(m_timsDataFrame.data(), src, count);
 
  qDebug() << "offset begin at last :" << count + 4;
 
  // std::vector<std::size_t> compressed_len_list;
  std::size_t offset;
  std::size_t previous_offset = (*(quint32 *)(src));
  qDebug() << "first offset= " << previous_offset;
  std::size_t cumul_decompressed_size = 0;
 
 
  for(quint32 i = 1; i <= m_scanCount; i++)
    {
      offset = (*(quint32 *)(src + (i * 4)));
 
      std::size_t compressed_size = offset - previous_offset;
 
      qDebug() << "scan i=" << i << " previous_offset=" << previous_offset
               << " offset=" << offset << " length=" << compressed_size;
      // compressed_len_list.push_back(offset - previous_offset);
      std::size_t remaining_size = m_binaryData.size();
 
      if(cumul_decompressed_size < remaining_size)
        {
          remaining_size = remaining_size - cumul_decompressed_size;
        }
      else
        {
          remaining_size = 0;
        }
      qDebug() << " remaining_size=" << remaining_size;
      std::size_t decompressed_size =
        lzfDecompressScan(src + previous_offset - 8,
                          compressed_size,
                          m_binaryData.data() + cumul_decompressed_size,
                          remaining_size);
 
 
      m_scanOffsetList.push_back(cumul_decompressed_size);
      m_scanSizeList.push_back(decompressed_size / 4);
      cumul_decompressed_size += decompressed_size;
      qDebug() << " decompressed_size=" << decompressed_size;
 
 
      previous_offset = offset;
    }
  /*
    std::size_t last_offset = (*(quint32 *)(src + (m_scanNumber * 4)));
    qDebug() << "last scan length :" << last_offset;
 
    qDebug() << "last scan length bonus:"
             << (*(quint32 *)(src + (m_scanNumber + 1 * 4)));
 
    qDebug() << " m_scanOffsetList.size()=" << m_scanOffsetList.size()
             << " m_scanNumber=" << m_scanNumber;
             */
  /*
    throw PappsoException(
      QObject::tr("ERROR reading TimsFrameType1 ").arg(m_timsId));
      */
}
 
 
unsigned int
TimsFrameType1::lzfDecompressScan(const char *src,
                                  unsigned int src_len,
                                  char *dest,
                                  unsigned int dest_len)
{
  qDebug() << "src=" << src << " src_len=" << src_len
           << " dest_len=" << dest_len;
  if(src_len == 0)
    return 0;
  unsigned int decompressed_size;
  unsigned int more_space = src_len * 2;
  decompressed_size       = lzf_decompress(src, src_len, dest, dest_len);
  while(decompressed_size == 0)
    {
      qDebug() << "dest_len=" << dest_len;
      qDebug() << "decompressed_size=" << decompressed_size;
 
      if(errno == EINVAL)
        {
          throw PappsoException(
            QObject::tr("ERROR reading TimsFrameType1 %1 TIMS binary file %2: "
                        "LZF decompression error EINVAL")
              .arg(m_frameId));
        }
      else if(errno == E2BIG)
        {
          qDebug() << " m_timsDataFrame.size()=" << m_binaryData.size()
                   << " more_space=" << more_space;
          m_binaryData.resize(m_binaryData.size() + more_space);
          dest_len += more_space;
          qDebug();
          decompressed_size = lzf_decompress(src, src_len, dest, dest_len);
        }
      else
        {
          break;
        }
    }
  return decompressed_size;
}
 
std::size_t
TimsFrameType1::getScanPeakCount(std::size_t scanNum) const
{
  pappso::MassSpectrumSPtr mass_spectrum_sptr = getMassSpectrumSPtr(scanNum);
  return mass_spectrum_sptr.get()->size();
}
 
 
void
TimsFrameType1::cumulateScan(std::size_t scanNum,
                             TimsDataFastMap &accumulate_into) const
{
  if(m_binaryData.size() == 0)
    return;
  // checkScanNum(scanNum);
 
 
  std::size_t size = m_scanSizeList[scanNum];
 
  std::size_t offset = m_scanOffsetList[scanNum];
 
  // qDebug() << "begin offset=" << offset << " size=" << size;
  qint32 value     = 0;
  qint32 tof_index = 0;
  for(std::size_t i = 0; i < size; i++)
    {
      value = (*(qint32 *)(m_binaryData.constData() + offset + (i * 4)));
      // qDebug() << " i=" << i << " value=" << value;
 
      if(value < 0)
        {
          tof_index += -1 * value;
        }
      else
        {
 
          quint32 x = tof_index;
          quint32 y = value;
 
          accumulate_into.accumulateIntensity(x, y);
 
          tof_index++;
        }
    }
  qDebug() << "end";
}
 
 
void
TimsFrameType1::cumulateScan2(std::size_t scanNum,
                              TimsDataFastMap &accumulate_into,
                              quint32 accepted_tof_index_range_begin,
                              quint32 accepted_tof_index_range_end) const
{
  if(m_binaryData.size() == 0)
    return;
  // checkScanNum(scanNum);
 
 
  std::size_t size = m_scanSizeList[scanNum];
 
  std::size_t offset = m_scanOffsetList[scanNum];
 
  // qDebug() << "begin offset=" << offset << " size=" << size;
  qint32 value     = 0;
  qint32 tof_index = 0;
  for(std::size_t i = 0; i < size; i++)
    {
      value = (*(qint32 *)(m_binaryData.constData() + offset + (i * 4)));
      // qDebug() << " i=" << i << " value=" << value;
 
      if(value < 0)
        {
          tof_index += -1 * value;
        }
      else
        {
 
          quint32 x = tof_index;
          quint32 y = value;
 
 
          if(x < accepted_tof_index_range_begin)
            {
              tof_index++;
              continue;
            }
          if(x > accepted_tof_index_range_end)
            {
              break;
            }
          accumulate_into.accumulateIntensity(x, y);
 
          tof_index++;
        }
    }
  qDebug() << "end";
}
 
std::vector<quint32>
TimsFrameType1::getScanTofIndexList(std::size_t scanNum) const
{
  qDebug();
  checkScanNum(scanNum);
 
  std::vector<quint32> mzindex_values;
 
  try
    {
      qDebug();
 
 
      if(m_binaryData.size() == 0)
        return mzindex_values;
      qDebug();
 
      std::size_t size = m_scanSizeList[scanNum];
 
      std::size_t offset = m_scanOffsetList[scanNum];
 
      qDebug() << " offset=" << offset << " size=" << size;
      if(size == 0)
        return mzindex_values;
 
      qint32 value     = 0;
      qint32 tof_index = 0;
      // std::vector<quint32> index_list;
      for(std::size_t i = 0; i < size; i++)
        {
          value = (*(qint32 *)(m_binaryData.constData() + offset + (i * 4)));
 
          if(value < 0)
            {
              tof_index += -1 * value;
            }
          else
            {
              mzindex_values.push_back(tof_index);
              tof_index++;
            }
        }
 
 
      qDebug();
      return mzindex_values;
    }
  catch(PappsoException &error)
    {
      throw pappso::PappsoException(QObject::tr("Error %1 frameId=%2 "
                                                "scanNum=%3 :\n%4")
                                      .arg(__FUNCTION__)
                                      .arg(getId())
                                      .arg(scanNum)
                                      .arg(error.qwhat()));
    }
  qDebug();
}
 
std::vector<quint32>
TimsFrameType1::getScanIntensityList(std::size_t scanNum) const
{
 
 
  qDebug() << " scanNum=" << scanNum;
 
  checkScanNum(scanNum);
 
  std::vector<quint32> int_values;
 
  try
    {
      qDebug();
 
 
      if(m_binaryData.size() == 0)
        return int_values;
      qDebug();
 
      std::size_t size = m_scanSizeList[scanNum];
 
      std::size_t offset = m_scanOffsetList[scanNum];
 
      qDebug() << " offset=" << offset << " size=" << size;
      if(size == 0)
        return int_values;
 
      qint32 value     = 0;
      qint32 tof_index = 0;
      // std::vector<quint32> index_list;
      for(std::size_t i = 0; i < size; i++)
        {
          value = (*(qint32 *)(m_binaryData.constData() + offset + (i * 4)));
 
          if(value < 0)
            {
              tof_index += -1 * value;
            }
          else
            {
              int_values.push_back(value);
              tof_index++;
            }
        }
 
 
      qDebug();
      return int_values;
    }
  catch(PappsoException &error)
    {
      throw pappso::PappsoException(QObject::tr("Error %1 frameId=%2 "
                                                "scanNum=%3 :\n%4")
                                      .arg(__FUNCTION__)
                                      .arg(getId())
                                      .arg(scanNum)
                                      .arg(error.qwhat()));
    }
}
 
pappso::MassSpectrumSPtr
TimsFrameType1::getMassSpectrumSPtr(std::size_t scanNum) const
{
 
  qDebug() << " scanNum=" << scanNum;
 
  checkScanNum(scanNum);
 
  try
    {
      qDebug();
 
      pappso::MassSpectrumSPtr mass_spectrum_sptr =
        std::make_shared<pappso::MassSpectrum>();
      // std::vector<DataPoint>
 
      if(m_binaryData.size() == 0)
        return mass_spectrum_sptr;
      qDebug();
 
      std::size_t size = m_scanSizeList[scanNum];
 
      std::size_t offset = m_scanOffsetList[scanNum];
 
      qDebug() << " offset=" << offset << " size=" << size;
      if(size == 0)
        return mass_spectrum_sptr;
 
 
      MzCalibrationInterface *mz_calibration_p =
        getMzCalibrationInterfaceSPtr().get();
 
 
      qint32 value     = 0;
      qint32 tof_index = 0;
      // std::vector<quint32> index_list;
      DataPoint data_point;
      for(std::size_t i = 0; i < size; i++)
        {
          value = (*(qint32 *)(m_binaryData.constData() + offset + (i * 4)));
 
          if(value < 0)
            {
              tof_index += -1 * value;
            }
          else
            {
              data_point.y = value;
 
              // intensity normalization
              data_point.y *= 100.0 / m_acqDurationInMilliseconds;
 
 
              // mz calibration
              data_point.x = mz_calibration_p->getMzFromTofIndex(tof_index);
              mass_spectrum_sptr.get()->push_back(data_point);
              tof_index++;
            }
        }
 
 
      qDebug() << mass_spectrum_sptr.get()->toString();
      return mass_spectrum_sptr;
    }
  catch(PappsoException &error)
    {
      throw pappso::PappsoException(
        QObject::tr("Error TimsFrameType1::getMassSpectrumSPtr frameId=%1 "
                    "scanNum=%2 :\n%3")
          .arg(getId())
          .arg(scanNum)
          .arg(error.qwhat()));
    }
}
 
 
pappso::TraceSPtr
TimsFrameType1::getRawTraceSPtr(std::size_t scanNum) const
{
 
  // qDebug();
 
  pappso::TraceSPtr trace_sptr = std::make_shared<pappso::Trace>();
  // std::vector<DataPoint>
 
  if(m_binaryData.size() == 0)
    return trace_sptr;
  qDebug();
 
  std::size_t size = m_scanSizeList[scanNum];
 
  std::size_t offset = m_scanOffsetList[scanNum];
 
  qDebug() << " offset=" << offset << " size=" << size;
  if(size == 0)
    return trace_sptr;
 
  // qDebug();
  qint32 value     = 0;
  qint32 tof_index = 0;
 
  // std::vector<quint32> index_list;
  DataPoint data_point;
  for(std::size_t i = 0; i < size; i++)
    {
      value = (*(qint32 *)(m_binaryData.constData() + offset + (i * 4)));
 
      if(value < 0)
        {
          tof_index += -1 * value;
        }
      else
        {
          data_point.y = value;
 
          // intensity normalization
          data_point.y *= 100.0 / m_acqDurationInMilliseconds;
 
 
          // mz calibration
          data_point.x = tof_index;
          trace_sptr.get()->push_back(data_point);
          tof_index++;
        }
    }
 
 
  // qDebug();
  return trace_sptr;
}
 
 
std::vector<TimsFrameType1::TofIndexIntensityPair>
TimsFrameType1::getRawValuePairList(std::size_t scanNum,
                                    quint32 accepted_tof_index_range_begin,
                                    quint32 accepted_tof_index_range_end) const
{
 
  // qDebug();
 
  std::vector<TimsFrame::TofIndexIntensityPair> trace_sptr;
  // std::vector<DataPoint>
 
  if(m_binaryData.size() == 0)
    return trace_sptr;
  qDebug();
 
  std::size_t size = m_scanSizeList[scanNum];
 
  std::size_t offset = m_scanOffsetList[scanNum];
 
  qDebug() << " offset=" << offset << " size=" << size;
  if(size == 0)
    return trace_sptr;
 
  // qDebug();
  qint32 value     = 0;
  qint32 tof_index = 0;
 
  // std::vector<quint32> index_list;
  TimsFrame::TofIndexIntensityPair data_point;
  for(std::size_t i = 0; i < size; i++)
    {
      value = (*(qint32 *)(m_binaryData.constData() + offset + (i * 4)));
 
      if(value < 0)
        {
          tof_index += -1 * value;
        }
      else
        {
          data_point.intensity_index = value;
 
          if(data_point.intensity_index < accepted_tof_index_range_begin)
            {
              tof_index++;
              continue;
            }
          if(data_point.intensity_index > accepted_tof_index_range_end)
            {
              break;
            }
          // intensity normalization
          data_point.intensity_index *= 100.0 / m_acqDurationInMilliseconds;
 
 
          // mz calibration
          data_point.tof_index = tof_index;
          trace_sptr.push_back(data_point);
          tof_index++;
        }
    }
 
 
  // qDebug();
  return trace_sptr;
}
 
} // namespace pappso