-
Notifications
You must be signed in to change notification settings - Fork 122
/
FindPeaks.h
259 lines (210 loc) · 9.87 KB
/
FindPeaks.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
#ifndef MANTID_ALGORITHMS_FINDPEAKS_H_
#define MANTID_ALGORITHMS_FINDPEAKS_H_
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidKernel/System.h"
#include "MantidAPI/Algorithm.h"
// #include "MantidAPI/IFunction.h"
#include "MantidAPI/IPeakFunction.h"
#include "MantidAPI/IBackgroundFunction.h"
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidKernel/cow_ptr.h"
namespace Mantid {
namespace HistogramData {
class HistogramX;
class HistogramY;
}
namespace Algorithms {
/** This algorithm searches for peaks in a dataset.
The method used is detailed in: M.A.Mariscotti, NIM 50 (1967) 309.
Required Properties:
<UL>
<LI> InputWorkspace - The name of the Workspace to search for peaks. </LI>
<LI> PeaksList - The name of the TableWorkspace in which to store the
list of peaks found. </LI>
</UL>
Optional Properties:
<UL>
<LI> fwhm - The number of points covered on average by the fwhm of a peak
(default 7) </LI>
<LI> Tolerance - Sets the strictness desired in meeting the conditions on
peak candidates (default 4, Mariscotti recommended 2) </LI>
<LI> WorkspaceIndex - The spectrum to search for peaks. Will search all
spectra if absent. </LI>
</UL>
@author Russell Taylor, Tessella Support Services plc
@date 25/11/2008
Copyright © 2008-9 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
National Laboratory & European Spallation Source
This file is part of Mantid.
Mantid 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.
Mantid 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 this program. If not, see <http://www.gnu.org/licenses/>.
File change history is stored at: <https://github.com/mantidproject/mantid>
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
class DLLExport FindPeaks : public API::Algorithm {
public:
/// Constructor
FindPeaks();
/// Virtual destructor
~FindPeaks() override {
if (m_progress)
delete m_progress;
m_progress = nullptr;
}
/// Algorithm's name
const std::string name() const override { return "FindPeaks"; }
/// Summary of algorithms purpose
const std::string summary() const override {
return "Searches for peaks in a dataset.";
}
/// Algorithm's version
int version() const override { return (1); }
/// Algorithm's category for identification
const std::string category() const override {
return "Optimization\\PeakFinding";
}
/// needed by FindPeaksBackground
int getIndex(const HistogramData::HistogramX &vecX, double x);
private:
void init() override;
void exec() override;
/// Process algorithm's properties
void processAlgorithmProperties();
/// Find peaks by searching peak position using Mariscotti algorithm
void findPeaksUsingMariscotti();
/// Find peaks according to given peak positions
void findPeaksGivenStartingPoints(const std::vector<double> &peakcentres,
const std::vector<double> &fitwindows);
/// Methods searving for findPeaksUsingMariscotti()
API::MatrixWorkspace_sptr
calculateSecondDifference(const API::MatrixWorkspace_const_sptr &input);
void smoothData(API::MatrixWorkspace_sptr &WS, const int &w);
void calculateStandardDeviation(const API::MatrixWorkspace_const_sptr &input,
const API::MatrixWorkspace_sptr &smoothed,
const int &w);
long long computePhi(const int &w) const;
/// Fit peak confined in a given window (x-min, x-max)
void fitPeakInWindow(const API::MatrixWorkspace_sptr &input,
const int spectrum, const double centre_guess,
const double xmin, const double xmax);
/// Fit peak by given/guessed FWHM
void fitPeakGivenFWHM(const API::MatrixWorkspace_sptr &input,
const int spectrum, const double center_guess,
const int fitWidth, const bool hasleftpeak,
const double leftpeakcentre, const bool hasrightpeak,
const double rightpeakcentre);
/// Fit peak: this is a basic peak fit function as a root function for all
/// different type of user input
void fitSinglePeak(const API::MatrixWorkspace_sptr &input, const int spectrum,
const int i_min, const int i_max, const int i_centre);
void fitPeakHighBackground(const API::MatrixWorkspace_sptr &input,
const size_t spectrum, int i_centre, int i_min,
int i_max, double &in_bg0, double &in_bg1,
double &in_bg2, int i_peakmin, int i_peakmax);
void fitPeakOneStep(const API::MatrixWorkspace_sptr &input,
const int spectrum, const int &i0, const int &i2,
const int &i4, const double &in_bg0, const double &in_bg1,
const double &in_bg2);
/// Add a new row in output TableWorkspace containing information of the
/// fitted peak+background
void addInfoRow(const size_t spectrum,
const API::IPeakFunction_const_sptr &peakfunction,
const API::IBackgroundFunction_sptr &bkgdfunction,
const bool isoutputraw, const double mincost);
/// Add the fit record (failure) to output workspace
void addNonFitRecord(const size_t spectrum, const double centre);
/// Create peak and background functions
void createFunctions();
/// Find peak background
int findPeakBackground(const API::MatrixWorkspace_sptr &input, int spectrum,
size_t i_min, size_t i_max,
std::vector<double> &vecBkgdParamValues,
std::vector<double> &vecpeakrange);
/// Estimate background of a given range
void estimateBackground(const HistogramData::HistogramX &X,
const HistogramData::HistogramY &Y,
const size_t i_min, const size_t i_max,
std::vector<double> &vecbkgdparvalues);
/// Estimate peak range based on background peak parameter
void estimatePeakRange(const HistogramData::HistogramX &vecX, size_t i_centre,
size_t i_min, size_t i_max, const double &leftfwhm,
const double &rightfwhm,
std::vector<double> &vecpeakrange);
/// Estimate peak parameters
std::string estimatePeakParameters(
const HistogramData::HistogramX &vecX,
const HistogramData::HistogramY &vecY, size_t i_min, size_t i_max,
const std::vector<double> &vecbkgdparvalues, size_t &iobscentre,
double &height, double &fwhm, double &leftfwhm, double &rightfwhm);
/// Generate a table workspace for output peak parameters
void generateOutputPeakParameterTable();
std::vector<double> getStartingPeakValues();
std::vector<double> getStartingBkgdValues();
/// Fit peak by calling 'FitPeak'
double callFitPeak(const API::MatrixWorkspace_sptr &dataws, int wsindex,
const API::IPeakFunction_sptr peakfunction,
const API::IBackgroundFunction_sptr backgroundfunction,
const std::vector<double> &vec_fitwindow,
const std::vector<double> &vec_peakrange, int minGuessFWHM,
int maxGuessFWHM, int guessedFWHMStep,
int estBackResult = 0);
std::vector<std::string> m_peakParameterNames;
std::vector<std::string> m_bkgdParameterNames;
size_t m_bkgdOrder;
/// The number of smoothing iterations. Set to 5, the optimum value according
/// to Mariscotti.
static const int g_z = 5;
/// Storage of the peak data
API::ITableWorkspace_sptr m_outPeakTableWS;
/// Progress reporting
API::Progress *m_progress;
// Properties saved in the algo.
API::MatrixWorkspace_sptr m_dataWS; ///<workspace to check for peaks
int m_inputPeakFWHM; ///<holder for the requested peak FWHM
int m_wsIndex; ///<list of workspace indicies to check
bool singleSpectrum; ///<flag for if only a single spectrum is present
bool m_highBackground; ///<flag for find relatively weak peak in high
/// background
bool m_rawPeaksTable; ///<flag for whether the output is the raw peak
/// parameters or effective (centre, width, height)
std::size_t
m_numTableParams; //<Number of parameters in the output table workspace
std::size_t m_centreIndex; //< Column in output table of peak centre
std::string m_peakFuncType; //< The name of the peak function to fit
std::string m_backgroundType; //< The type of background to fit
// Peaks positions
std::vector<double> m_vecPeakCentre;
std::vector<double> m_vecFitWindows;
// Functions for reused
API::IBackgroundFunction_sptr m_backgroundFunction;
API::IPeakFunction_sptr m_peakFunction;
int m_minGuessedPeakWidth;
int m_maxGuessedPeakWidth;
int m_stepGuessedPeakWidth;
bool m_usePeakPositionTolerance;
double m_peakPositionTolerance;
std::vector<API::IFunction_sptr> m_fitFunctions;
std::vector<size_t> m_peakLeftIndexes;
std::vector<size_t> m_peakRightIndexes;
std::string m_minimizer;
std::string m_costFunction;
/// Minimum peak height
double m_minHeight;
/// Minimum value of peak's observed maximum Y value
double m_leastMaxObsY;
/// Start values
bool m_useObsCentre;
};
} // namespace Algorithms
} // namespace Mantid
#endif /*MANTID_ALGORITHMS_FINDPEAKS_H_*/