libpappsomspp
Library for mass spectrometry
filterlowintensitysignalremoval.cpp
Go to the documentation of this file.
1 /* BEGIN software license
2  *
3  * msXpertSuite - mass spectrometry software suite
4  * -----------------------------------------------
5  * Copyright(C) 2009,...,2021 Filippo Rusconi
6  *
7  * http://www.msxpertsuite.org
8  *
9  * This file is part of the msXpertSuite project.
10  *
11  * The msXpertSuite project is the successor of the massXpert project. This
12  * project now includes various independent modules:
13  *
14  * - massXpert, model polymer chemistries and simulate mass spectrometric data;
15  * - mineXpert, a powerful TIC chromatogram/mass spectrum viewer/miner;
16  *
17  * This program is free software: you can redistribute it and/or modify
18  * it under the terms of the GNU General Public License as published by
19  * the Free Software Foundation, either version 3 of the License, or
20  * (at your option) any later version.
21  *
22  * This program is distributed in the hope that it will be useful,
23  * but WITHOUT ANY WARRANTY; without even the implied warranty of
24  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25  * GNU General Public License for more details.
26  *
27  * You should have received a copy of the GNU General Public License
28  * along with this program. If not, see <http://www.gnu.org/licenses/>.
29  *
30  * END software license
31  */
32 
33 
34 #include <qmath.h>
35 
36 #include <limits>
37 #include <iterator>
38 
39 #include <QVector>
40 #include <QDebug>
41 
42 #include "../../exception/exceptionnotrecognized.h"
43 
45 
46 
47 namespace pappso
48 {
49 
50 
52  double mean, double std_dev, double threshold)
53 {
54  m_noiseMean = mean;
55  m_noiseStdDev = std_dev;
56  m_threshold = threshold;
57 }
58 
59 
61  const QString &parameters)
62 {
63  buildFilterFromString(parameters);
64 }
65 
66 
69 {
70  m_noiseMean = other.m_noiseMean;
72  m_threshold = other.m_threshold;
73 }
74 
75 
77 {
78 }
79 
80 
84 {
85  if(&other == this)
86  return *this;
87 
88  m_noiseMean = other.m_noiseMean;
90  m_threshold = other.m_threshold;
91 
92  return *this;
93 }
94 
95 
96 void
98  const QString &parameters)
99 {
100  // Typical string: "FloorAmplitudePercentage|15"
101  if(parameters.startsWith(QString("%1|").arg(name())))
102  {
103  QStringList params = parameters.split("|").back().split(";");
104 
105  m_noiseMean = params.at(0).toDouble();
106  m_noiseStdDev = params.at(1).toDouble();
107  m_threshold = params.at(2).toDouble();
108  }
109  else
110  {
112  QString(
113  "Building of FilterLowIntensitySignalRemoval from string %1 failed")
114  .arg(parameters));
115  }
116 }
117 
118 
119 std::size_t
121 {
122  // We want to iterate in the trace and check if we can get the clusters
123  // isolated one by one.
124 
125  m_isotopicClusters.clear();
126 
127  std::size_t trace_size = trace.size();
128  if(trace_size <= 2)
129  {
130  qDebug() << "The original trace has less than 3 points. Returning it "
131  "without modification.";
132  return m_isotopicClusters.size();
133  }
134  else
135  qDebug() << "The original trace has" << trace_size << "data points";
136 
137  std::size_t index = 0;
138 
139  // Seed the system with the first point of the trace.
140  m_prev = IndexIntensity(index, trace[index].y);
141 
142  qDebug() << "First trace point: "
143  << "(" << trace[m_prev.first].x << "," << m_prev.second << ");";
144 
145  // We still have not found any apex!
146  m_prevApex = std::pair(0, nan);
147 
148  // We will need to know if we were ascending to an apex.
149  bool was_ascending_to_apex = false;
150 
151  ClusterSPtr cluster_sp = std::make_shared<Cluster>();
152 
153  // Now that we have seeded the system with the first point of the original
154  // trace, go to the next one.
155  ++index;
156 
157  // And now iterate in the original trace and make sure we detect all the
158  // clusters.
159 
160  while(index < trace_size)
161  {
162  // Get the current trace point as an IndexIntensity.
163  m_cur = IndexIntensity(index, trace[index].y);
164 
165  qDebug() << "Current trace point: "
166  << "(" << trace[m_cur.first].x << "," << m_cur.second << ");";
167 
168  // We monitor if we are going up or down a peak.
169 
170  if(m_cur.second > m_prev.second)
171  // We are ascending a peak. We do not know if we are at the apex.
172  {
173  qDebug().noquote() << "We are ascending to an apex.\n";
174  was_ascending_to_apex = true;
175  }
176  else
177  // We are descending a peak.
178  {
179  qDebug().noquote() << "Descending a peak. ";
180 
181  // There are two situations:
182  //
183  // 1. Either we were ascending to an apex, and m_prev is that apex,
184  //
185  // 2. Or we were not ascending to an apex and in fact all we are doing
186  // is going down an apex that occurred more than one trace point ago.
187 
188  if(!was_ascending_to_apex)
189  // We were not ascending to an apex.
190  {
191  // Do nothing here.
192 
193  qDebug().noquote()
194  << "But, we were not ascending to an apex, so do nothing.\n";
195  }
196  else
197  // We are effectively descending a peak right after the apex was
198  // reached.
199  {
200  qDebug().noquote()
201  << "And, we were ascending to an apex, so "
202  "m_curApeex becomes m_prev: ("
203  << trace[m_curApex.first].x << "," << m_curApex.second << ")\n";
204 
205  m_curApex = m_prev;
206 
207  // We might have two situations:
208 
209  // 1. We had already encountered an apex.
210  // 2. We had not yet encountered an apex.
211 
212  if(!std::isnan(m_prevApex.second))
213  // We had already encountered an apex.
214  {
215  // Was that apex far on the left of the current apex ?
216 
217  if(trace[m_curApex.first].x - trace[m_prevApex.first].x >
219  // The distance is not compatible with both apices to belong
220  // to the same cluster.
221  {
222  // We are creating a new isotopic cluster.
223 
224  // But, since another apex had been encountered already,
225  // that means that an isotopic cluster was cooking
226  // already. We must store it.
227 
228  if(!cluster_sp->size())
229  qFatal("Cannot be that the cluster has no apex.");
230 
231  m_isotopicClusters.push_back(cluster_sp);
232 
233  qDebug().noquote()
234  << "There was a previous apex already, BUT "
235  "outside of cluster range. "
236  "Pushing the cooking cluster that has size:"
237  << cluster_sp->size();
238 
239  // Now create a brand new cluster.
240  cluster_sp = std::make_shared<Cluster>();
241 
242  qDebug() << "Created a brand new cluster.";
243 
244  // We only start the new cluster with the current apex if
245  // that apex is above the threshold.
246 
247  if(m_curApex.second > m_threshold)
248  {
249  // And start it with the current apex.
250  cluster_sp->push_back(
251  std::make_shared<IndexIntensity>(m_curApex));
252 
253  qDebug()
254  << "Since the current apex is above the threshold, "
255  "we PUSH it to the newly created cluster: ("
256  << trace[m_curApex.first].x << ","
257  << m_curApex.second << ")";
258 
260 
261  qDebug() << "Set prev apex to be cur apex.";
262  }
263  else
264  {
265  qDebug()
266  << "Since the current apex is below the threshold, "
267  "we DO NOT push it to the newly created "
268  "cluster: ("
269  << trace[m_curApex.first].x << ","
270  << m_curApex.second << ")";
271 
272  // Since previous apex went to the closed cluster, we
273  // need to reset it.
274 
275  m_prevApex = std::pair(0, nan);
276 
277  qDebug() << "Since the previous apex went to the "
278  "closed cluster, and cur apex has too "
279  "small an intensity, we reset prev apex "
280  "to std::pair(0, nan).";
281  }
282  }
283  else
284  // The distance is compatible with both apices to belong to
285  // the same isotopic cluster.
286  {
287 
288  // But we only push back the current apex to the cluster
289  // if its intensity is above the threshold.
290 
291  if(m_curApex.second > m_threshold)
292  // The current apex was above the threshold
293  {
294  cluster_sp->push_back(
295  std::make_shared<IndexIntensity>(m_curApex));
296 
297  qDebug().noquote()
298  << "There was an apex already inside of cluster "
299  "range. "
300  "AND, since the current apex was above the "
301  "threshold, we indeed push it to cluster.\n";
302 
303  qDebug().noquote() << "Current apex PUSHED: "
304  << trace[m_curApex.first].x << ", "
305  << m_curApex.second;
306 
308 
309  qDebug() << "We set prev apex to be cur apex.";
310  }
311  else
312  {
313  qDebug().noquote()
314  << "There was an apex already inside of cluster "
315  "range. "
316  "BUT, since the current apex was below the "
317  "threshold, we do not push it to cluster.\n";
318 
319  qDebug().noquote() << "Current apex NOT pushed: "
320  << trace[m_curApex.first].x << ", "
321  << m_curApex.second;
322  }
323  }
324  }
325  else
326  // No apex was previously found. We are fillin-up a new isotopic
327  // cluster.
328  {
329  if(m_curApex.second > m_threshold)
330  // We can actually add that apex to a new isotopic cluster.
331  {
332  if(cluster_sp->size() != 0)
333  qCritical(
334  "At this point, the cluster should be new and "
335  "empty.");
336 
337  cluster_sp->push_back(
338  std::make_shared<IndexIntensity>(m_curApex));
339 
340  qDebug().noquote()
341  << "No previous apex was found. Since current apex' "
342  "intensity is above threshold, we push it back to "
343  "the "
344  "cluster.\n";
345 
346  // Store current apex as previous apex for next rounds.
348 
349  qDebug().noquote() << "And thus we store the current "
350  "apex as previous apex.\n";
351  }
352  else
353  {
354  qDebug().noquote()
355  << "No previous apex was found. Since current apex' "
356  "intensity is below threshold, we do nothing.\n";
357  }
358  }
359  }
360  // End of
361  // ! if(!was_ascending_to_apex)
362  // That is, we were ascending to an apex.
363 
364  // Tell what it is: we were not ascending to an apex.
365  was_ascending_to_apex = false;
366  }
367  // End of
368  // ! if(m_cur.second > m_prev.second)
369  // That is, we are descending a peak.
370 
371  // At this point, prepare next round.
372 
373  qDebug().noquote()
374  << "Preparing next round, with m_prev = m_cur and ++index.\n";
375 
376  m_prev = m_cur;
377  ++index;
378  }
379  // End of
380  // while(index < trace_size)
381 
382  // At this point, if a cluster had been cooking, add it.
383 
384  if(cluster_sp->size())
385  m_isotopicClusters.push_back(cluster_sp);
386 
387  return m_isotopicClusters.size();
388 }
389 
390 
391 Trace &
393 {
394  qDebug();
395 
396  // Horrible hack to have a non const filtering process.
397  return const_cast<FilterLowIntensitySignalRemoval *>(this)->nonConstFilter(
398  trace);
399 }
400 
401 
402 Trace &
404 {
405  qDebug();
406 
407  // Make a copy of the trace.
408  m_helperTraceCopy = Trace(trace);
409 
410  if(trace.size() <= 2)
411  {
412  qDebug() << "The original trace has less than 3 points. Returning it "
413  "without modification.";
414  return trace;
415  }
416  else
417  qDebug() << "The original trace had" << trace.size() << "data points";
418 
419  std::size_t isotopic_cluster_count = detectIsotopicClusters(trace);
420 
421  qDebug() << "Number of isotopic clusters: " << isotopic_cluster_count;
422 
423  Trace final_trace;
424 
425  for(auto &&apices_p : m_isotopicClusters)
426  {
427  qDebug() << "iterating in new cluster.";
428 
429  for(auto &&apex_p : *apices_p)
430  {
431  DataPoint data_point =
432  DataPoint(trace[apex_p->first].x, apex_p->second);
433 
434  qDebug() << "Current cluster data point:" << data_point.toString();
435 
436  final_trace.push_back(data_point);
437  }
438  }
439 
440  trace = std::move(final_trace);
441 
442  return trace;
443 }
444 
445 
446 double
448 {
449  return m_threshold;
450 }
451 
452 
453 //! Return a string with the textual representation of the configuration data.
454 QString
456 {
457  return QString("%1|%2").arg(name()).arg(QString::number(m_threshold, 'f', 2));
458 }
459 
460 
461 QString
463 {
464  return "FilterLowIntensitySignalRemoval";
465 }
466 
467 } // namespace pappso
excetion to use when an item type is not recognized
Redefines the floor intensity of the Trace.
std::shared_ptr< std::vector< ApexSPtr > > ClusterSPtr
FilterLowIntensitySignalRemoval(double mean, double std_dev, double threshold)
Trace & filter(Trace &data_points) const override
FilterLowIntensitySignalRemoval & operator=(const FilterLowIntensitySignalRemoval &other)
void buildFilterFromString(const QString &strBuildParams) override
build this filter using a string
QString toString() const override
Return a string with the textual representation of the configuration data.
A simple container of DataPoint instances.
Definition: trace.h:132
tries to keep as much as possible monoisotopes, removing any possible C13 peaks and changes multichar...
Definition: aa.cpp:39
QString toString() const
Definition: datapoint.cpp:139