tracedetectionzivy.cpp

Go to the documentation of this file.

 
/*******************************************************************************
 * Copyright (c) 2015 Olivier Langella <Olivier.Langella@moulon.inra.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/>.
 *
 * Contributors:
 *     Olivier Langella <Olivier.Langella@moulon.inra.fr> - initial API and
 *implementation
 ******************************************************************************/
 
#include "tracedetectionzivy.h"
#include "tracepeak.h"
#include "pappsomspp/core/exception/exceptionoutofrange.h"
#include "pappsomspp/core/exception/exceptionnotpossible.h"
 
#include <QObject>
 
 
namespace pappso
{
TraceDetectionZivy::TraceDetectionZivy(unsigned int smoothing_half_window_length,
                                       unsigned int minmax_half_window_length,
                                       unsigned int maxmin_half_window_length,
                                       double detection_threshold_on_minmax,
                                       double detection_threshold_on_maxmin)
  : m_smooth(smoothing_half_window_length),
    m_minMax(minmax_half_window_length),
    m_maxMin(maxmin_half_window_length)
{
  m_detectionThresholdOnMaxMin = detection_threshold_on_maxmin;
  m_detectionThresholdOnMinMax = detection_threshold_on_minmax;
}
 
TraceDetectionZivy::TraceDetectionZivy(const TraceDetectionZivy &other)
  : m_smooth(other.m_smooth), m_minMax(other.m_minMax), m_maxMin(other.m_maxMin)
{
  m_detectionThresholdOnMaxMin = other.m_detectionThresholdOnMaxMin;
  m_detectionThresholdOnMinMax = other.m_detectionThresholdOnMinMax;
}
TraceDetectionZivy::~TraceDetectionZivy()
{
}
 
void
TraceDetectionZivy::setFilterMorphoMean(const FilterMorphoMean &smooth)
{
  m_smooth = smooth;
}
void
TraceDetectionZivy::setFilterMorphoMinMax(const FilterMorphoMinMax &minMax)
{
  m_minMax = minMax;
}
void
TraceDetectionZivy::setFilterMorphoMaxMin(const FilterMorphoMaxMin &maxMin)
{
  m_maxMin = maxMin;
}
 
void
TraceDetectionZivy::setDetectionThresholdOnMinmax(double detectionThresholdOnMinMax)
{
  m_detectionThresholdOnMinMax = detectionThresholdOnMinMax;
}
void
TraceDetectionZivy::setDetectionThresholdOnMaxmin(double detectionThresholdOnMaxMin)
{
  m_detectionThresholdOnMaxMin = detectionThresholdOnMaxMin;
}
unsigned int
TraceDetectionZivy::getSmoothingHalfEdgeWindows() const
{
  return m_smooth.getHalfWindowSize();
};
unsigned int
TraceDetectionZivy::getMaxMinHalfEdgeWindows() const
{
  return m_maxMin.getMaxMinHalfEdgeWindows();
};
 
unsigned int
TraceDetectionZivy::getMinMaxHalfEdgeWindows() const
{
  return m_minMax.getMinMaxHalfEdgeWindows();
};
pappso_double
TraceDetectionZivy::getDetectionThresholdOnMinmax() const
{
  return m_detectionThresholdOnMinMax;
};
pappso_double
TraceDetectionZivy::getDetectionThresholdOnMaxmin() const
{
  return m_detectionThresholdOnMaxMin;
};
 
 
void
TraceDetectionZivy::detect(const Trace &xic,
                           TraceDetectionSinkInterface &sink,
                           bool remove_peak_base) const
{
 
  // detect peak positions on close curve : a peak is an intensity value
  // strictly greater than the two surrounding values. In case of
  // equality (very rare, can happen with some old old spectrometers) we
  // take the last equal point to be the peak
  qDebug();
  std::size_t size = xic.size();
  if(size < 4)
    {
      qDebug() << QObject::tr("The original XIC is too small to detect peaks (%1)").arg(xic.size());
      return;
    }
  if(size <= m_maxMin.getMaxMinHalfEdgeWindows())
    return;
  if(size <= m_minMax.getMinMaxHalfEdgeWindows())
    return;
 
  Trace xic_minmax(xic); //"close" courbe du haut
  if(m_smooth.getHalfWindowSize() != 0)
    {
      qDebug() << "f";
      m_smooth.filter(xic_minmax);
    }
 
  qDebug();
  Trace xic_maxmin(xic_minmax); //"open" courbe du bas
 
  qDebug();
 
  try
    {
      qDebug() << "f1";
      m_minMax.filter(xic_minmax);
      qDebug() << "f2";
      m_maxMin.filter(xic_maxmin);
    }
  catch(const ExceptionOutOfRange &e)
    {
      qDebug() << QObject::tr("The original XIC is too small to detect peaks (%1) :\n%2")
                    .arg(xic.size())
                    .arg(e.qwhat());
      return;
    }
  qDebug() << "a";
 
  std::vector<DataPoint>::const_iterator previous_rt = xic_minmax.begin();
  std::vector<DataPoint>::const_iterator current_rt  = (previous_rt + 1);
  std::vector<DataPoint>::const_iterator last_rt     = (xic_minmax.end() - 1);
 
  std::vector<DataPoint>::const_iterator current_rt_on_maxmin = (xic_maxmin.begin() + 1);
 
  std::vector<DataPoint>::const_iterator xic_position = xic.begin();
  qDebug() << "b";
  while(current_rt != last_rt)
    // for (unsigned int i = 1, count = 0; i < xic_minmax.size() - 1; )
    {
      // conditions to have a peak
      if((previous_rt->y < current_rt->y) && (current_rt->y > m_detectionThresholdOnMinMax) &&
         (current_rt_on_maxmin->y > m_detectionThresholdOnMaxMin))
        {
          // here we test the last condition to have a peak
 
          // no peak case
          if(current_rt->y < (current_rt + 1)->y)
            {
              //++i;
              previous_rt = current_rt;
              current_rt++;
              current_rt_on_maxmin++;
            }
          // there is a peak here ! case
          else if(current_rt->y > (current_rt + 1)->y)
            {
              // peak.setMaxXicElement(*current_rt);
 
              // find left boundary
              std::vector<DataPoint>::const_iterator it_left =
                moveLowerYLeftDataPoint(xic_minmax, current_rt);
              // walk back
              it_left = moveLowerYRigthDataPoint(xic_minmax, it_left);
              // peak.setLeftBoundary(*it_left);
 
              // find right boundary
              std::vector<DataPoint>::const_iterator it_right =
                moveLowerYRigthDataPoint(xic_minmax, current_rt);
              // walk back
              it_right = moveLowerYLeftDataPoint(xic_minmax, it_right);
              // peak.setRightBoundary(*it_right);
 
              // integrate peak surface :
              auto it      = findFirstEqualOrGreaterX(xic_position, xic.end(), it_left->x);
              xic_position = findFirstEqualOrGreaterX(it, xic.end(), it_right->x) + 1;
              // peak.setArea(areaTrace(it, xic_position));
 
              // find the maximum :
              // peak.setMaxXicElement(*maxYDataPoint(it, xic_position));
 
              // areaTrace()
              TracePeak peak(it, xic_position, remove_peak_base);
              sink.setTracePeak(peak);
              // }
              //++i;
              previous_rt = current_rt;
              current_rt++;
              current_rt_on_maxmin++;
            }
          // equality case, skipping equal points
          else
            {
              // while (v_minmax[i] == v_minmax[i + 1]) {
              //++i;
              current_rt++;
              current_rt_on_maxmin++;
 
              //++count;
            }
        }
      // no chance to have a peak at all, continue looping
      else
        {
          //++i;
          previous_rt = current_rt;
          current_rt++;
          current_rt_on_maxmin++;
        }
    } // end loop for peaks
  qDebug();
}
} // namespace pappso
 
void
pappso::TraceDetectionZivy::guessParametersFromMsRunReader( [[maybe_unused]] const pappso::MsRunReader &reader)
{
  // reader.
}