Home  · Classes  · Annotated Classes  · Modules  · Members  · Namespaces  · Related Pages
SavitzkyGolayFilter.h
Go to the documentation of this file.
1 // --------------------------------------------------------------------------
2 // OpenMS -- Open-Source Mass Spectrometry
3 // --------------------------------------------------------------------------
4 // Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
5 // ETH Zurich, and Freie Universitaet Berlin 2002-2013.
6 //
7 // This software is released under a three-clause BSD license:
8 // * Redistributions of source code must retain the above copyright
9 // notice, this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright
11 // notice, this list of conditions and the following disclaimer in the
12 // documentation and/or other materials provided with the distribution.
13 // * Neither the name of any author or any participating institution
14 // may be used to endorse or promote products derived from this software
15 // without specific prior written permission.
16 // For a full list of authors, refer to the file AUTHORS.
17 // --------------------------------------------------------------------------
18 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 // ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
22 // INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
23 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
24 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
25 // OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
26 // WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
27 // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
28 // ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29 //
30 // --------------------------------------------------------------------------
31 // $Maintainer: Alexandra Zerck $
32 // $Authors: Eva Lange $
33 // --------------------------------------------------------------------------
34 
35 #ifndef OPENMS_FILTERING_SMOOTHING_SAVITZKYGOLAYFILTER_H
36 #define OPENMS_FILTERING_SMOOTHING_SAVITZKYGOLAYFILTER_H
37 
41 
42 #include <gsl/gsl_vector.h>
43 #include <gsl/gsl_matrix.h>
44 #include <gsl/gsl_linalg.h>
45 #include <gsl/gsl_permutation.h>
46 #include <gsl/gsl_pow_int.h>
47 
48 namespace OpenMS
49 {
107  class OPENMS_DLLAPI SavitzkyGolayFilter :
108  public ProgressLogger,
109  public DefaultParamHandler
110  {
111 public:
114 
116  virtual ~SavitzkyGolayFilter();
117 
121  template <typename PeakType>
122  void filter(MSSpectrum<PeakType> & spectrum)
123  {
124  UInt n = (UInt)spectrum.size();
125 
126  typename MSSpectrum<PeakType>::iterator first = spectrum.begin();
127  typename MSSpectrum<PeakType>::iterator last = spectrum.end();
128 
129  //copy the data AND META DATA to the output container
130  MSSpectrum<PeakType> output = spectrum;
131 
132  if (frame_size_ > n)
133  {
134  return;
135  }
136 
137  int i;
138  UInt j;
139  int mid = (frame_size_ / 2);
140  double help;
141 
142  typename MSSpectrum<PeakType>::iterator it_forward;
143  typename MSSpectrum<PeakType>::iterator it_help;
144  typename MSSpectrum<PeakType>::iterator out_it = output.begin();
145 
146  // compute the transient on
147  for (i = 0; i <= mid; ++i)
148  {
149  it_forward = (first - i);
150  help = 0;
151 
152  for (j = 0; j < frame_size_; ++j)
153  {
154  help += it_forward->getIntensity() * coeffs_[(i + 1) * frame_size_ - 1 - j];
155  ++it_forward;
156  }
157 
158 
159  out_it->setPosition(first->getPosition());
160  out_it->setIntensity(std::max(0.0, help));
161  ++out_it;
162  ++first;
163  }
164 
165  // compute the steady state output
166  it_help = (last - mid);
167  while (first != it_help)
168  {
169  it_forward = (first - mid);
170  help = 0;
171 
172  for (j = 0; j < frame_size_; ++j)
173  {
174  help += it_forward->getIntensity() * coeffs_[mid * frame_size_ + j];
175  ++it_forward;
176  }
177 
178 
179  out_it->setPosition(first->getPosition());
180  out_it->setIntensity(std::max(0.0, help));
181  ++out_it;
182  ++first;
183  }
184 
185  // compute the transient off
186  for (i = (mid - 1); i >= 0; --i)
187  {
188  it_forward = (first - (frame_size_ - i - 1));
189  help = 0;
190 
191  for (j = 0; j < frame_size_; ++j)
192  {
193  help += it_forward->getIntensity() * coeffs_[i * frame_size_ + j];
194  ++it_forward;
195  }
196 
197  out_it->setPosition(first->getPosition());
198  out_it->setIntensity(std::max(0.0, help));
199  ++out_it;
200  ++first;
201  }
202 
203  spectrum = output;
204  }
205 
206  template <typename PeakType>
207  void filter(MSChromatogram<PeakType> & chromatogram)
208  {
209 
210  MSSpectrum<PeakType> filter_spectra;
211  for (typename MSChromatogram<PeakType>::const_iterator it = chromatogram.begin(); it != chromatogram.end(); ++it)
212  {
213  filter_spectra.push_back(*it);
214  }
215  filter(filter_spectra);
216  chromatogram.clear(false);
217  for (typename MSSpectrum<PeakType>::const_iterator it = filter_spectra.begin(); it != filter_spectra.end(); ++it)
218  {
219  chromatogram.push_back(*it);
220  }
221 
222  }
223 
227  template <typename PeakType>
229  {
230  Size progress = 0;
231  startProgress(0, map.size() + map.getChromatograms().size(), "smoothing data");
232  for (Size i = 0; i < map.size(); ++i)
233  {
234  filter(map[i]);
235  setProgress(++progress);
236  }
237  for (Size i = 0; i < map.getChromatograms().size(); ++i)
238  {
239  filter(map.getChromatogram(i));
240  setProgress(++progress);
241  }
242  endProgress();
243  }
244 
245 protected:
247  std::vector<DoubleReal> coeffs_;
252  // Docu in base class
253  virtual void updateMembers_();
254 
255  };
256 
257 } // namespace OpenMS
258 #endif // OPENMS_FILTERING_SMOOTHING_SAVITZKYGOLAYFILTER_H
Size size() const
Definition: MSExperiment.h:117
The representation of a chromatogram.
Definition: MSChromatogram.h:53
unsigned int UInt
Unsigned integer type.
Definition: Types.h:92
UInt frame_size_
UInt of the filter kernel (number of pre-tabulated coefficients)
Definition: SavitzkyGolayFilter.h:249
void filter(MSSpectrum< PeakType > &spectrum)
Removed the noise from an MSSpectrum containing profile data.
Definition: SavitzkyGolayFilter.h:122
UInt order_
The order of the smoothing polynomial.
Definition: SavitzkyGolayFilter.h:251
Computes the Savitzky-Golay filter coefficients using QR decomposition.
Definition: SavitzkyGolayFilter.h:107
MSChromatogram< ChromatogramPeakType > & getChromatogram(Size id)
returns a single chromatogram
Definition: MSExperiment.h:774
void filter(MSChromatogram< PeakType > &chromatogram)
Definition: SavitzkyGolayFilter.h:207
Representation of a mass spectrometry experiment.
Definition: MSExperiment.h:68
void filterExperiment(MSExperiment< PeakType > &map)
Removed the noise from an MSExperiment containing profile data.
Definition: SavitzkyGolayFilter.h:228
std::vector< DoubleReal > coeffs_
Coefficients.
Definition: SavitzkyGolayFilter.h:247
size_t Size
Size type e.g. used as variable which can hold result of size()
Definition: Types.h:144
Base class for all classes that want to report their progess.
Definition: ProgressLogger.h:56
A base class for all classes handling default parameters.
Definition: DefaultParamHandler.h:90
void clear(bool clear_meta_data)
Clears all data and meta data.
Definition: MSChromatogram.h:564
const std::vector< MSChromatogram< ChromatogramPeakType > > & getChromatograms() const
returns the chromatogram list
Definition: MSExperiment.h:768

OpenMS / TOPP release 1.11.1 Documentation generated on Thu Nov 14 2013 11:19:20 using doxygen 1.8.5