libpappsomspp
Library for mass spectrometry
tracedetectionzivy.cpp
Go to the documentation of this file.
1 
2 /*******************************************************************************
3  * Copyright (c) 2015 Olivier Langella <Olivier.Langella@moulon.inra.fr>.
4  *
5  * This file is part of the PAPPSOms++ library.
6  *
7  * PAPPSOms++ is free software: you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation, either version 3 of the License, or
10  * (at your option) any later version.
11  *
12  * PAPPSOms++ is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with PAPPSOms++. If not, see <http://www.gnu.org/licenses/>.
19  *
20  * Contributors:
21  * Olivier Langella <Olivier.Langella@moulon.inra.fr> - initial API and
22  *implementation
23  ******************************************************************************/
24 
25 #include "tracedetectionzivy.h"
26 #include "tracepeak.h"
27 #include "../../exception/exceptionoutofrange.h"
28 #include "../../exception/exceptionnotpossible.h"
29 
30 #include <QObject>
31 
32 
33 namespace pappso
34 {
36  unsigned int smoothing_half_window_length,
37  unsigned int minmax_half_window_length,
38  unsigned int maxmin_half_window_length,
39  pappso_double detection_threshold_on_minmax,
40  pappso_double detection_threshold_on_maxmin)
41  : m_smooth(smoothing_half_window_length),
42  m_minMax(minmax_half_window_length),
43  m_maxMin(maxmin_half_window_length)
44 {
45  m_detectionThresholdOnMaxMin = detection_threshold_on_maxmin;
46  m_detectionThresholdOnMinMax = detection_threshold_on_minmax;
47 }
49 {
50 }
51 
52 void
54 {
55  m_smooth = smooth;
56 }
57 void
59 {
60  m_minMax = minMax;
61 }
62 void
64 {
65  m_maxMin = maxMin;
66 }
67 
68 void
70  double detectionThresholdOnMinMax)
71 {
72  m_detectionThresholdOnMinMax = detectionThresholdOnMinMax;
73 }
74 void
76  double detectionThresholdOnMaxMin)
77 {
78  m_detectionThresholdOnMaxMin = detectionThresholdOnMaxMin;
79 }
80 unsigned int
82 {
84 };
85 unsigned int
87 {
89 };
90 
91 unsigned int
93 {
95 };
98 {
100 };
103 {
105 };
106 
107 
108 void
111  bool remove_peak_base) const
112 {
113 
114  // detect peak positions on close curve : a peak is an intensity value
115  // strictly greater than the two surrounding values. In case of
116  // equality (very rare, can happen with some old old spectrometers) we
117  // take the last equal point to be the peak
118  qDebug();
119  std::size_t size = xic.size();
120  if(size < 4)
121  {
122  qDebug() << QObject::tr(
123  "The original XIC is too small to detect peaks (%1)")
124  .arg(xic.size());
125  return;
126  }
127  if(size <= m_maxMin.getMaxMinHalfEdgeWindows())
128  return;
129  if(size <= m_minMax.getMinMaxHalfEdgeWindows())
130  return;
131 
132  Trace xic_minmax(xic); //"close" courbe du haut
134  {
135  qDebug() << "f";
136  m_smooth.filter(xic_minmax);
137  }
138 
139  qDebug();
140  Trace xic_maxmin(xic_minmax); //"open" courbe du bas
141 
142  qDebug();
143 
144  try
145  {
146  qDebug() << "f1";
147  m_minMax.filter(xic_minmax);
148  qDebug() << "f2";
149  m_maxMin.filter(xic_maxmin);
150  }
151  catch(const ExceptionOutOfRange &e)
152  {
153  qDebug() << QObject::tr(
154  "The original XIC is too small to detect peaks (%1) :\n%2")
155  .arg(xic.size())
156  .arg(e.qwhat());
157  return;
158  }
159  qDebug() << "a";
160 
161  std::vector<DataPoint>::const_iterator previous_rt = xic_minmax.begin();
162  std::vector<DataPoint>::const_iterator current_rt = (previous_rt + 1);
163  std::vector<DataPoint>::const_iterator last_rt = (xic_minmax.end() - 1);
164 
165  std::vector<DataPoint>::const_iterator current_rt_on_maxmin =
166  (xic_maxmin.begin() + 1);
167 
168  std::vector<DataPoint>::const_iterator xic_position = xic.begin();
169  qDebug() << "b";
170  while(current_rt != last_rt)
171  // for (unsigned int i = 1, count = 0; i < xic_minmax.size() - 1; )
172  {
173  // conditions to have a peak
174  if((previous_rt->y < current_rt->y) &&
175  (current_rt->y > m_detectionThresholdOnMinMax) &&
176  (current_rt_on_maxmin->y > m_detectionThresholdOnMaxMin))
177  {
178  // here we test the last condition to have a peak
179 
180  // no peak case
181  if(current_rt->y < (current_rt + 1)->y)
182  {
183  //++i;
184  previous_rt = current_rt;
185  current_rt++;
186  current_rt_on_maxmin++;
187  }
188  // there is a peak here ! case
189  else if(current_rt->y > (current_rt + 1)->y)
190  {
191  // peak.setMaxXicElement(*current_rt);
192 
193  // find left boundary
194  std::vector<DataPoint>::const_iterator it_left =
195  moveLowerYLeftDataPoint(xic_minmax, current_rt);
196  // walk back
197  it_left = moveLowerYRigthDataPoint(xic_minmax, it_left);
198  // peak.setLeftBoundary(*it_left);
199 
200  // find right boundary
201  std::vector<DataPoint>::const_iterator it_right =
202  moveLowerYRigthDataPoint(xic_minmax, current_rt);
203  // walk back
204  it_right = moveLowerYLeftDataPoint(xic_minmax, it_right);
205  // peak.setRightBoundary(*it_right);
206 
207  // integrate peak surface :
208  auto it =
209  findFirstEqualOrGreaterX(xic_position, xic.end(), it_left->x);
210  xic_position =
211  findFirstEqualOrGreaterX(it, xic.end(), it_right->x) + 1;
212  // peak.setArea(areaTrace(it, xic_position));
213 
214  // find the maximum :
215  // peak.setMaxXicElement(*maxYDataPoint(it, xic_position));
216 
217  // areaTrace()
218  TracePeak peak(it, xic_position, remove_peak_base);
219  sink.setTracePeak(peak);
220  // }
221  //++i;
222  previous_rt = current_rt;
223  current_rt++;
224  current_rt_on_maxmin++;
225  }
226  // equality case, skipping equal points
227  else
228  {
229  // while (v_minmax[i] == v_minmax[i + 1]) {
230  //++i;
231  current_rt++;
232  current_rt_on_maxmin++;
233 
234  //++count;
235  }
236  }
237  // no chance to have a peak at all, continue looping
238  else
239  {
240  //++i;
241  previous_rt = current_rt;
242  current_rt++;
243  current_rt_on_maxmin++;
244  }
245  } // end loop for peaks
246  qDebug();
247 }
248 } // namespace pappso
transform the trace with the maximum of the minimum equivalent of the erode filter for pictures
Definition: filtermorpho.h:138
Trace & filter(Trace &data_points) const override
std::size_t getMaxMinHalfEdgeWindows() const
mean filter apply mean of y values inside the window : this results in a kind of smoothing
Definition: filtermorpho.h:210
std::size_t getMeanHalfEdgeWindows() const
transform the trace with the minimum of the maximum equivalent of the dilate filter for pictures
Definition: filtermorpho.h:118
std::size_t getMinMaxHalfEdgeWindows() const
Trace & filter(Trace &data_points) const override
virtual Trace & filter(Trace &data_points) const override
virtual const QString & qwhat() const
virtual void setTracePeak(TracePeak &xic_peak)=0
void setFilterMorphoMinMax(const FilterMorphoMinMax &m_minMax)
void setDetectionThresholdOnMaxmin(double detectionThresholdOnMaxMin)
unsigned int getMinMaxHalfEdgeWindows() const
void detect(const Trace &xic, TraceDetectionSinkInterface &sink, bool remove_peak_base) const override
detect peaks on a trace
unsigned int getSmoothingHalfEdgeWindows() const
pappso_double getDetectionThresholdOnMinmax() const
pappso_double getDetectionThresholdOnMaxmin() const
void setFilterMorphoMaxMin(const FilterMorphoMaxMin &maxMin)
TraceDetectionZivy(unsigned int smoothing_half_window_length, unsigned int minmax_half_window_length, unsigned int maxmin_half_window_length, pappso_double detection_threshold_on_minmax, pappso_double detection_threshold_on_maxmin)
pappso_double m_detectionThresholdOnMaxMin
void setFilterMorphoMean(const FilterMorphoMean &smooth)
pappso_double m_detectionThresholdOnMinMax
unsigned int getMaxMinHalfEdgeWindows() const
void setDetectionThresholdOnMinmax(double detectionThresholdOnMinMax)
A simple container of DataPoint instances.
Definition: trace.h:148
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition: aa.cpp:39
std::vector< DataPoint >::iterator findFirstEqualOrGreaterX(std::vector< DataPoint >::iterator begin, std::vector< DataPoint >::iterator end, const double &value)
find the first element in which X is equal or greater than the value searched important : it implies ...
Definition: trace.cpp:69
std::vector< DataPoint >::const_iterator moveLowerYLeftDataPoint(const Trace &trace, std::vector< DataPoint >::const_iterator begin)
Move left to the lower value.
Definition: trace.cpp:223
double pappso_double
A type definition for doubles.
Definition: types.h:49
std::vector< DataPoint >::const_iterator moveLowerYRigthDataPoint(const Trace &trace, std::vector< DataPoint >::const_iterator begin)
Move right to the lower value.
Definition: trace.cpp:205