Mantid
Loading...
Searching...
No Matches
WienerSmooth.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
8
11#include "MantidAPI/TextAxis.h"
14
15#include <numeric>
16
17namespace Mantid::Algorithms {
18
21
22// Register the algorithm into the AlgorithmFactory
23DECLARE_ALGORITHM(WienerSmooth)
24
25namespace {
26// Square values
27struct PowerSpectrum {
28 double operator()(double x) const { return x * x; }
29};
30
31// To be used when actual noise level cannot be estimated
32const double guessSignalToNoiseRatio = 1e15;
33} // namespace
34
35//----------------------------------------------------------------------------------------------
36
38int WienerSmooth::version() const { return 1; }
39
41const std::string WienerSmooth::category() const { return "Arithmetic\\FFT;Transforms\\Smoothing"; }
42
44const std::string WienerSmooth::summary() const { return "Smooth spectra using Wiener filter."; }
45
46//----------------------------------------------------------------------------------------------
50 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input), "An input workspace.");
51 declareProperty(std::make_unique<Kernel::ArrayProperty<int>>("WorkspaceIndexList"),
52 "Workspace indices for spectra to process. "
53 "If empty smooth all spectra.");
54 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
55 "An output workspace.");
56}
57
58//----------------------------------------------------------------------------------------------
62 // Get the data to smooth
63 API::MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
64 std::vector<int> wsIndexList = getProperty("WorkspaceIndexList");
65
66 // number of spectra in the input workspace
67 const size_t nInputSpectra = inputWS->getNumberHistograms();
68
69 // Validate the input
70 if (wsIndexList.size() > nInputSpectra) {
71 throw std::invalid_argument("Workspace index list has more indices than "
72 "there are spectra in the input workspace.");
73 }
74
75 // if empty do whole workspace
76 if (wsIndexList.empty()) {
77 // fill wsIndexList with consecutive integers from 0 to nSpectra - 1
78 wsIndexList.resize(nInputSpectra);
79 wsIndexList.front() = 0;
80 for (auto index = wsIndexList.begin() + 1; index != wsIndexList.end(); ++index) {
81 *index = *(index - 1) + 1;
82 }
83 }
84
85 // number of spectra in the output workspace
86 const size_t nOutputSpectra = wsIndexList.size();
87
88 // smooth the first spectrum to find out the output blocksize
89 auto wsIndex = static_cast<size_t>(wsIndexList.front());
90 auto first = smoothSingleSpectrum(inputWS, wsIndex);
91
92 // create full output workspace by copying all settings from tinputWS
93 // blocksize is taken form first
95 API::WorkspaceFactory::Instance().create(inputWS, nOutputSpectra, first->x(0).size(), first->y(0).size());
96
97 // TODO: ideally axis cloning should be done via API::Axis interface but it's
98 // not possible
99 // at he moment and as it turned out not straight-forward to implement
100 auto inAxis = inputWS->getAxis(1);
101 auto outAxis = std::unique_ptr<API::Axis>(inAxis->clone(nOutputSpectra, outputWS.get()));
102
103 bool isSpectra = outAxis->isSpectra();
104 bool isNumeric = outAxis->isNumeric();
105 auto inTextAxis = dynamic_cast<API::TextAxis *>(inAxis);
106 auto outTextAxis = dynamic_cast<API::TextAxis *>(outAxis.get());
107
108 // Initialise the progress reporting object
109 API::Progress progress(this, 0.0, 1.0, nOutputSpectra);
110
111 // smooth the rest of the input
112 for (size_t outIndex = 0; outIndex < nOutputSpectra; ++outIndex) {
113 auto inIndex = wsIndexList[outIndex];
114 auto next = outIndex == 0 ? first : smoothSingleSpectrum(inputWS, inIndex);
115
116 // copy the values
117 outputWS->setHistogram(outIndex, next->histogram(0));
118
119 // set the axis value
120 if (isSpectra) {
121 auto &inSpectrum = inputWS->getSpectrum(inIndex);
122 auto &outSpectrum = outputWS->getSpectrum(outIndex);
123 outSpectrum.setSpectrumNo(inSpectrum.getSpectrumNo());
124 outSpectrum.setDetectorIDs(inSpectrum.getDetectorIDs());
125 } else if (isNumeric) {
126 outAxis->setValue(outIndex, inAxis->getValue(inIndex));
127 } else if (inTextAxis && outTextAxis) {
128 outTextAxis->setLabel(outIndex, inTextAxis->label(inIndex));
129 }
130 progress.report();
131 }
132 outputWS->replaceAxis(1, std::move(outAxis));
133 // set the output
134 setProperty("OutputWorkspace", outputWS);
135}
136
137//----------------------------------------------------------------------------------------------
145 size_t dataSize = inputWS->blocksize();
146
147 // it won't work for very small workspaces
148 if (dataSize < 4) {
149 g_log.debug() << "No smoothing, spectrum copied.\n";
150 return copyInput(inputWS, wsIndex);
151 }
152
153 // Due to the way RealFFT works the input should be even-sized
154 const bool isOddSize = dataSize % 2 != 0;
155 if (isOddSize) {
156 // add a fake value to the end to make size even
157 inputWS = copyInput(inputWS, wsIndex);
158 wsIndex = 0;
159
160 auto &X = inputWS->x(wsIndex);
161 auto &Y = inputWS->y(wsIndex);
162 auto &E = inputWS->e(wsIndex);
163 double dx = X[dataSize - 1] - X[dataSize - 2];
164 auto histogram = inputWS->histogram(wsIndex);
165 histogram.resize(histogram.size() + 1);
166 histogram.mutableX().back() = X.back() + dx;
167 histogram.mutableY().back() = Y.back();
168 histogram.mutableE().back() = E.back();
169 inputWS->setHistogram(wsIndex, histogram);
170 }
171
172 // the input vectors
173 auto &X = inputWS->x(wsIndex);
174 auto &Y = inputWS->y(wsIndex);
175
176 // Digital fourier transform works best for data oscillating around 0.
177 // Fit a spline with a small number of break points to the data.
178 // Make sure that the spline passes through the first and the last points
179 // of the data.
180 // The fitted spline will be subtracted from the data and the difference
181 // will be smoothed with the Wiener filter. After that the spline will be
182 // added to the smoothed data to produce the output.
183
184 // number of spline break points, must be smaller than the data size but
185 // between 2 and 10
186 size_t nbreak = 10;
187 if (nbreak * 3 > dataSize)
188 nbreak = dataSize / 3;
189
190 // NB. The spline mustn't fit too well to the data. If it does smoothing
191 // doesn't happen.
192 // TODO: it's possible that the spline is unnecessary and a simple linear
193 // function will
194 // do a better job.
195
196 g_log.debug() << "Spline break points " << nbreak << '\n';
197
198 // define the spline
199 API::IFunction_sptr spline = API::FunctionFactory::Instance().createFunction("BSpline");
200 auto xInterval = getStartEnd(X, inputWS->isHistogramData());
201 setSplineDataBounds(spline, xInterval.first, xInterval.second);
202 spline->setAttributeValue("NBreak", static_cast<int>(nbreak));
203 // fix the first and last parameters to the first and last data values
204 spline->setParameter(0, Y.front());
205 spline->fix(0);
206 size_t lastParamIndex = spline->nParams() - 1;
207 spline->setParameter(lastParamIndex, Y.back());
208 spline->fix(lastParamIndex);
209
210 // fit the spline to the data
211 auto fit = createChildAlgorithm("Fit");
212 fit->initialize();
213 fit->setProperty("Function", spline);
214 fit->setProperty("InputWorkspace", inputWS);
215 fit->setProperty("WorkspaceIndex", static_cast<int>(wsIndex));
216 fit->setProperty("CreateOutput", true);
217 fit->execute();
218
219 // get the fit output workspace; spectrum 2 contains the difference that is to
220 // be smoothed
221 API::MatrixWorkspace_sptr fitOut = fit->getProperty("OutputWorkspace");
222
223 // Fourier transform the difference spectrum
224 auto fourier = createChildAlgorithm("RealFFT");
225 fourier->initialize();
226 fourier->setProperty("InputWorkspace", fitOut);
227 fourier->setProperty("WorkspaceIndex", 2);
228 // we don't require bin linearity as we don't need the exact transform
229 fourier->setProperty("IgnoreXBins", true);
230 fourier->execute();
231
232 API::MatrixWorkspace_sptr fourierOut = fourier->getProperty("OutputWorkspace");
233
234 // spectrum 2 of the transformed workspace has the transform modulus which is
235 // a square
236 // root of the power spectrum
237 auto &powerSpec = fourierOut->mutableY(2);
238 // convert the modulus to power spectrum wich is the base of the Wiener filter
239 std::transform(powerSpec.begin(), powerSpec.end(), powerSpec.begin(), PowerSpectrum());
240
241 // estimate power spectrum's noise as the average of its high frequency half
242 size_t n2 = powerSpec.size();
243 double noise = std::accumulate(powerSpec.begin() + n2 / 2, powerSpec.end(), 0.0);
244 noise /= static_cast<double>(n2);
245
246 // index of the maximum element in powerSpec
247 const size_t imax =
248 static_cast<size_t>(std::distance(powerSpec.begin(), std::max_element(powerSpec.begin(), powerSpec.end())));
249
250 if (noise == 0.0) {
251 noise = powerSpec[imax] / guessSignalToNoiseRatio;
252 }
253
254 g_log.debug() << "Maximum signal " << powerSpec[imax] << '\n';
255 g_log.debug() << "Noise " << noise << '\n';
256
257 // storage for the Wiener filter, initialized with 0.0's
258 std::vector<double> wf(n2);
259
260 // The filter consists of two parts:
261 // 1) low frequency region, from 0 until the power spectrum falls to the
262 // noise level, filter is calculated
263 // from the power spectrum
264 // 2) high frequency noisy region, filter is a smooth function of frequency
265 // decreasing to 0
266
267 // the following code is an adaptation of a fortran routine
268 // noise starting index
269 size_t i0 = 0;
270 // intermediate variables
271 double xx = 0.0;
272 double xy = 0.0;
273 double ym = 0.0;
274 // low frequency filter values: the higher the power spectrum the closer the
275 // filter to 1.0
276 for (size_t i = 0; i < n2; ++i) {
277 double cd1 = powerSpec[i] / noise;
278 if (cd1 < 1.0 && i > imax) {
279 i0 = i;
280 break;
281 }
282 double cd2 = log(cd1);
283 wf[i] = cd1 / (1.0 + cd1);
284 auto j = static_cast<double>(i + 1);
285 xx += j * j;
286 xy += j * cd2;
287 ym += cd2;
288 }
289
290 // i0 should always be > 0 but in case something goes wrong make a check
291 if (i0 > 0) {
292 g_log.debug() << "Noise start index " << i0 << '\n';
293
294 // high frequency filter values: smooth decreasing function
295 auto ri0f = static_cast<double>(i0 + 1);
296 double xm = (1.0 + ri0f) / 2;
297 ym /= ri0f;
298 double a1 = (xy - ri0f * xm * ym) / (xx - ri0f * xm * xm);
299 double b1 = ym - a1 * xm;
300
301 g_log.debug() << "(a1,b1) = (" << a1 << ',' << b1 << ')' << '\n';
302
303 const double dblev = -20.0;
304 // cut-off index
305 double ri1 = floor((dblev / 4 - b1) / a1);
306 if (ri1 < static_cast<double>(i0)) {
307 g_log.warning() << "Failed to build Wiener filter: no smoothing.\n";
308 ri1 = static_cast<double>(i0);
309 }
310 auto i1 = static_cast<size_t>(ri1);
311 if (i1 > n2)
312 i1 = n2;
313 for (size_t i = i0; i < i1; ++i) {
314 double s = exp(a1 * static_cast<double>(i + 1) + b1);
315 wf[i] = s / (1.0 + s);
316 }
317 // wf[i] for i1 <= i < n2 are 0.0
318
319 g_log.debug() << "Cut-off index " << i1 << '\n';
320 } else {
321 g_log.warning() << "Power spectrum has an unexpected shape: no smoothing\n";
322 return copyInput(inputWS, wsIndex);
323 }
324
325 // multiply the fourier transform by the filter
326 auto &re = fourierOut->mutableY(0);
327 auto &im = fourierOut->mutableY(1);
328
329 std::transform(re.begin(), re.end(), wf.begin(), re.begin(), std::multiplies<double>());
330 std::transform(im.begin(), im.end(), wf.begin(), im.begin(), std::multiplies<double>());
331
332 // inverse fourier transform
333 fourier = createChildAlgorithm("RealFFT");
334 fourier->initialize();
335 fourier->setProperty("InputWorkspace", fourierOut);
336 fourier->setProperty("IgnoreXBins", true);
337 fourier->setPropertyValue("Transform", "Backward");
338 fourier->execute();
339
340 API::MatrixWorkspace_sptr out = fourier->getProperty("OutputWorkspace");
341 auto &background = fitOut->y(1);
342 auto &y = out->mutableY(0);
343
344 if (y.size() != background.size()) {
345 throw std::logic_error("Logic error: inconsistent arrays");
346 }
347
348 // add the spline "background" to the smoothed data
349 std::transform(y.begin(), y.end(), background.begin(), y.begin(), std::plus<double>());
350
351 // copy the x-values and errors from the original spectrum
352 // remove the last values for odd-sized inputs
353
354 if (isOddSize) {
355 auto histogram = out->histogram(0);
356 histogram.setSharedX(inputWS->sharedX(wsIndex));
357 histogram.setSharedE(inputWS->sharedE(wsIndex));
358
359 auto newSize = histogram.y().size() - 1;
360 histogram.resize(newSize);
361
362 out->setHistogram(0, histogram);
363 } else {
364 out->setSharedX(0, inputWS->sharedX(wsIndex));
365 out->setSharedE(0, inputWS->sharedE(wsIndex));
366 }
367
368 return out;
369}
370
378std::pair<double, double> WienerSmooth::getStartEnd(const Mantid::HistogramData::HistogramX &X,
379 bool isHistogram) const {
380 const size_t n = X.size();
381 if (n < 3) {
382 // 3 is the smallest number for this method to work without breaking
383 throw std::runtime_error("Number of bins/data points cannot be smaller than 3.");
384 }
385 if (isHistogram) {
386 return std::make_pair((X[0] + X[1]) / 2, (X[n - 1] + X[n - 2]) / 2);
387 }
388 // else
389 return std::make_pair(X.front(), X.back());
390}
391
399 auto alg = createChildAlgorithm("ExtractSingleSpectrum");
400 alg->initialize();
401 alg->setProperty("InputWorkspace", inputWS);
402 alg->setProperty("WorkspaceIndex", static_cast<int>(wsIndex));
403 alg->execute();
404 API::MatrixWorkspace_sptr ws = alg->getProperty("OutputWorkspace");
405 return ws;
406}
407
408void WienerSmooth::setSplineDataBounds(API::IFunction_sptr &spline, const double startX, const double endX) {
409 const double currentEndX = spline->getAttribute("EndX").asDouble();
410 if (startX >= currentEndX) {
411 spline->setAttributeValue("EndX", endX);
412 spline->setAttributeValue("StartX", startX);
413 } else {
414 spline->setAttributeValue("StartX", startX);
415 spline->setAttributeValue("EndX", endX);
416 }
417}
418
419} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
Class to represent a text axis of a workspace.
Definition: TextAxis.h:36
A property class for workspaces.
void init() override
Initialize the algorithm's properties.
void setSplineDataBounds(API::IFunction_sptr &spline, const double startX, const double endX)
std::pair< double, double > getStartEnd(const Mantid::HistogramData::HistogramX &X, bool isHistogram) const
Get the start and end of the x-interval.
void exec() override
Execute the algorithm.
const std::string category() const override
Algorithm's category for identification.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
API::MatrixWorkspace_sptr copyInput(const API::MatrixWorkspace_sptr &inputWS, size_t wsIndex)
Exctract the input spectrum into a separate workspace.
API::MatrixWorkspace_sptr smoothSingleSpectrum(API::MatrixWorkspace_sptr inputWS, size_t wsIndex)
Execute smoothing of a single spectrum.
int version() const override
Algorithm's version for identification.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
Definition: IFunction.h:732
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
Describes the direction (within an algorithm) of a Property.
Definition: Property.h:50
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54