Mantid
Loading...
Searching...
No Matches
LineProfile.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
15#include "MantidHistogramData/HistogramBuilder.h"
16#include "MantidHistogramData/HistogramIterator.h"
21#include "MantidKernel/Unit.h"
22
23#include "boost/make_shared.hpp"
24#include <algorithm>
25
26namespace Mantid::Algorithms {
27
39using Mantid::HistogramData::HistogramBuilder;
45
46namespace {
48enum class LineDirection { horizontal, vertical };
49
51namespace DirectionChoices {
52const static std::string HORIZONTAL{"Horizontal"};
53const static std::string VERTICAL{"Vertical"};
54} // namespace DirectionChoices
55
57namespace ModeChoices {
58const static std::string AVERAGE{"Average"};
59const static std::string SUM{"Sum"};
60} // namespace ModeChoices
61
63namespace PropertyNames {
64const static std::string CENTRE{"Centre"};
65const static std::string DIRECTION{"Direction"};
66const static std::string END{"End"};
67const static std::string INPUT_WORKSPACE{"InputWorkspace"};
68const static std::string HALF_WIDTH{"HalfWidth"};
69const static std::string IGNORE_INFS{"IgnoreInfs"};
70const static std::string IGNORE_NANS{"IgnoreNans"};
71const static std::string MODE{"Mode"};
72const static std::string OUTPUT_WORKSPACE{"OutputWorkspace"};
73const static std::string START{"Start"};
74} // namespace PropertyNames
75
77struct Box {
78 double top;
79 double bottom;
80 double left;
81 double right;
82};
83
85struct IndexLimits {
86 size_t lineStart;
87 size_t lineEnd;
88 size_t widthStart;
89 size_t widthEnd;
90};
91
100Workspace2D_sptr makeOutput(const MatrixWorkspace &parent, const LineDirection direction, std::vector<double> &&Xs,
101 std::vector<double> &&Ys, std::vector<double> &&Es) {
102 HistogramBuilder builder;
103 builder.setX(std::move(Xs));
104 builder.setY(std::move(Ys));
105 builder.setE(std::move(Es));
106 builder.setDistribution(direction == LineDirection::horizontal && parent.isDistribution());
107 return create<Workspace2D>(parent, 1, builder.build());
108}
109
117void setAxesAndUnits(Workspace2D &outWS, const MatrixWorkspace &ws, const Box &box, const LineDirection dir) {
118 // Y units.
119 outWS.setYUnit(ws.YUnit());
120 outWS.setYUnitLabel(ws.YUnitLabel());
121 // Horizontal axis.
122 auto axisIndex = dir == LineDirection::horizontal ? 0 : 1;
123 if (ws.getAxis(axisIndex)->isSpectra()) {
124 outWS.getAxis(axisIndex)->setUnit("Empty");
125 } else {
126 outWS.getAxis(0)->setUnit(ws.getAxis(axisIndex)->unit()->unitID());
127 }
128 // Vertical axis. We'll use bin edges set to Centre +/- HalfWidth.
129 std::vector<double> vertBins(2);
130 vertBins.front() = dir == LineDirection::horizontal ? box.top : box.left;
131 vertBins.back() = dir == LineDirection::horizontal ? box.bottom : box.right;
132 auto outVertAxis = std::make_unique<BinEdgeAxis>(vertBins);
133 axisIndex = dir == LineDirection::horizontal ? 1 : 0;
134 if (ws.getAxis(axisIndex)->isSpectra()) {
135 outVertAxis->setUnit("Empty");
136 } else {
137 outVertAxis->setUnit(ws.getAxis(axisIndex)->unit()->unitID());
138 }
139 outWS.replaceAxis(1, std::move(outVertAxis));
140}
141
151template <typename Container>
152std::pair<size_t, size_t> startAndEnd(const Container &bins, const bool isBinEdges, const double lowerLimit,
153 const double upperLimit) {
154 auto lowerIt = std::upper_bound(bins.cbegin(), bins.cend(), lowerLimit);
155 if (lowerIt == bins.cend()) {
156 throw std::runtime_error("Profile completely outside input workspace.");
157 }
158 if (lowerIt != bins.cbegin()) {
159 --lowerIt;
160 }
161 auto upperIt = std::upper_bound(lowerIt, bins.cend(), upperLimit);
162 if (upperIt == bins.cbegin()) {
163 throw std::runtime_error("Profile completely outside input workspace.");
164 }
165 if (isBinEdges && upperIt == bins.cend()) {
166 --upperIt;
167 }
168 const auto start = std::distance(bins.cbegin(), lowerIt);
169 const auto end = std::distance(bins.cbegin(), upperIt);
170 return std::pair<size_t, size_t>{start, end};
171}
172
180std::vector<double> extractVerticalBins(const Axis &axis, const size_t numberHistograms) {
181 if (axis.isSpectra()) {
182 std::vector<double> spectrumNumbers(numberHistograms);
183 std::iota(spectrumNumbers.begin(), spectrumNumbers.end(), 1.0);
184 return spectrumNumbers;
185 }
186 std::vector<double> bins(axis.length());
187 for (size_t i = 0; i < bins.size(); ++i) {
188 bins[i] = axis.getValue(i);
189 }
190 return bins;
191}
192
207template <typename Container, typename Function>
208void profile(std::vector<double> &Xs, std::vector<double> &Ys, std::vector<double> &Es, const MatrixWorkspace &ws,
209 const LineDirection dir, const IndexLimits &limits, const Container &lineBins, const bool isBinEdges,
210 Function modeFunction, const bool ignoreNans, const bool ignoreInfs) {
211 const auto lineSize = limits.lineEnd - limits.lineStart;
212 Xs.resize(lineSize + (isBinEdges ? 1 : 0));
213 Ys.resize(lineSize);
214 Es.resize(lineSize);
215 for (size_t i = limits.lineStart; i < limits.lineEnd; ++i) {
216 Xs[i - limits.lineStart] = lineBins[i];
217 double ySum = 0;
218 double eSqSum = 0;
219 int n = 0;
220 for (size_t j = limits.widthStart; j < limits.widthEnd; ++j) {
221 const size_t iHor = dir == LineDirection::horizontal ? i : j;
222 const size_t iVert = dir == LineDirection::horizontal ? j : i;
223 auto histogram = ws.histogram(iVert);
224 auto iter = histogram.begin();
225 std::advance(iter, iHor);
226 const double y = iter->counts();
227 if ((ignoreNans && std::isnan(y)) || (ignoreInfs && std::isinf(y))) {
228 continue;
229 }
230 const double e = iter->countStandardDeviation();
231 ySum += y;
232 eSqSum += e * e;
233 ++n;
234 }
235 const size_t nTotal = limits.widthEnd - limits.widthStart;
236 Ys[i - limits.lineStart] = n == 0 ? std::nan("") : modeFunction(ySum, n, nTotal);
237 const double e = modeFunction(std::sqrt(eSqSum), n, nTotal);
238 Es[i - limits.lineStart] = std::isnan(e) ? 0 : e;
239 }
240 if (isBinEdges) {
241 Xs.back() = lineBins[limits.lineEnd];
242 }
243}
244
251double averageMode(const double sum, const size_t n, const size_t /*unused*/) noexcept {
252 return sum / static_cast<double>(n);
253}
254
263double sumMode(const double sum, const size_t n, const size_t nTot) noexcept {
264 return static_cast<double>(nTot) / static_cast<double>(n) * sum;
265}
266
271auto createMode(const std::string &modeName) noexcept {
272 if (modeName == ModeChoices::AVERAGE) {
273 return averageMode;
274 }
275 return sumMode;
276}
277
278void divideByBinHeight(MatrixWorkspace &ws) {
279 const BinEdgeAxis &axis = *static_cast<BinEdgeAxis *>(ws.getAxis(1));
280 const auto height = axis.getMax() - axis.getMin();
281 ws.mutableY(0) /= height;
282 ws.mutableE(0) /= height;
283}
284} // namespace
285
286// Register the algorithm into the AlgorithmFactory
287DECLARE_ALGORITHM(LineProfile)
288
289//----------------------------------------------------------------------------------------------
290
291
292const std::string LineProfile::name() const { return "LineProfile"; }
293
295int LineProfile::version() const { return 1; }
296
298const std::string LineProfile::category() const { return "Utility"; }
299
301const std::string LineProfile::summary() const { return "Calculates a line profile over a MatrixWorkspace."; }
302
303//----------------------------------------------------------------------------------------------
307 const auto mandatoryDouble = std::make_shared<MandatoryValidator<double>>();
308 const auto positiveDouble = std::make_shared<BoundedValidator<double>>();
309 positiveDouble->setLower(0.0);
310 positiveDouble->setLowerExclusive(true);
311 const auto mandatoryPositiveDouble = std::make_shared<CompositeValidator>();
312 mandatoryPositiveDouble->add(mandatoryDouble);
313 mandatoryPositiveDouble->add(positiveDouble);
314 const auto inputWorkspaceValidator = std::make_shared<CompositeValidator>();
315 inputWorkspaceValidator->add(std::make_shared<CommonBinsValidator>());
316 inputWorkspaceValidator->add(std::make_shared<IncreasingAxisValidator>());
318 Direction::Input, inputWorkspaceValidator),
319 "An input workspace.");
322 "A single histogram workspace containing the profile.");
323 declareProperty(PropertyNames::CENTRE, EMPTY_DBL(), mandatoryDouble, "Centre of the line.");
324 declareProperty(PropertyNames::HALF_WIDTH, EMPTY_DBL(), mandatoryPositiveDouble,
325 "Half of the width over which to calcualte the profile.");
326 const std::array<std::string, 2> directions = {{DirectionChoices::HORIZONTAL, DirectionChoices::VERTICAL}};
327 declareProperty(PropertyNames::DIRECTION, DirectionChoices::HORIZONTAL,
328 std::make_shared<ListValidator<std::string>>(directions), "Orientation of the profile line.");
329 declareProperty(PropertyNames::START, EMPTY_DBL(), "Starting point of the line.");
330 declareProperty(PropertyNames::END, EMPTY_DBL(), "End point of the line.");
331 const std::array<std::string, 2> modes = {{ModeChoices::AVERAGE, ModeChoices::SUM}};
332 declareProperty(PropertyNames::MODE, ModeChoices::AVERAGE, std::make_shared<ListValidator<std::string>>(modes),
333 "How the profile is calculated over the line width.");
334 declareProperty(PropertyNames::IGNORE_INFS, false, "If true, ignore infinities when calculating the profile.");
335 declareProperty(PropertyNames::IGNORE_NANS, true, "If true, ignore not-a-numbers when calculating the profile.");
336}
337
338//----------------------------------------------------------------------------------------------
342 // Extract properties.
344 const bool ignoreNans = getProperty(PropertyNames::IGNORE_NANS);
345 const bool ignoreInfs = getProperty(PropertyNames::IGNORE_INFS);
346 const auto &horizontalBins = ws->x(0);
347 const auto horizontalIsBinEdges = ws->isHistogramData();
348 const auto vertAxis = ws->getAxis(1);
349 // It is easier if the vertical axis values are in a vector.
350 const auto verticalBins = extractVerticalBins(*vertAxis, ws->getNumberHistograms());
351 const auto verticalIsBinEdges = verticalBins.size() > ws->getNumberHistograms();
352 const std::string directionString = getProperty(PropertyNames::DIRECTION);
353 LineDirection dir{LineDirection::horizontal};
354 if (directionString == DirectionChoices::VERTICAL) {
355 dir = LineDirection::vertical;
356 }
357 const double centre = getProperty(PropertyNames::CENTRE);
358 const double halfWidth = getProperty(PropertyNames::HALF_WIDTH);
359 double start = getProperty(PropertyNames::START);
360 if (start == EMPTY_DBL()) {
361 start = std::numeric_limits<double>::lowest();
362 }
363 double end = getProperty(PropertyNames::END);
364 if (end == EMPTY_DBL()) {
365 end = std::numeric_limits<double>::max();
366 }
367 // Define a box in workspace's units to have a standard representation
368 // of the profile's dimensions.
369 Box bounds;
370 if (dir == LineDirection::horizontal) {
371 bounds.top = centre - halfWidth;
372 bounds.bottom = centre + halfWidth;
373 bounds.left = start;
374 bounds.right = end;
375 } else {
376 bounds.top = start;
377 bounds.bottom = end;
378 bounds.left = centre - halfWidth;
379 bounds.right = centre + halfWidth;
380 }
381 // Convert the bounds from workspace units to indices.
382 const auto vertInterval = startAndEnd(verticalBins, verticalIsBinEdges, bounds.top, bounds.bottom);
383 const auto horInterval = startAndEnd(horizontalBins, horizontalIsBinEdges, bounds.left, bounds.right);
384 // Choose mode.
385 auto mode = createMode(getProperty(PropertyNames::MODE));
386 // Build the actual profile.
387 std::vector<double> profileYs;
388 std::vector<double> profileEs;
389 std::vector<double> Xs;
390 if (dir == LineDirection::horizontal) {
391 IndexLimits limits;
392 limits.lineStart = horInterval.first;
393 limits.lineEnd = horInterval.second;
394 limits.widthStart = vertInterval.first;
395 limits.widthEnd = vertInterval.second;
396 profile(Xs, profileYs, profileEs, *ws, dir, limits, horizontalBins, horizontalIsBinEdges, mode, ignoreNans,
397 ignoreInfs);
398 } else {
399 IndexLimits limits;
400 limits.lineStart = vertInterval.first;
401 limits.lineEnd = vertInterval.second;
402 limits.widthStart = horInterval.first;
403 limits.widthEnd = horInterval.second;
404 profile(Xs, profileYs, profileEs, *ws, dir, limits, verticalBins, verticalIsBinEdges, mode, ignoreNans, ignoreInfs);
405 }
406 // Prepare and set output.
407 auto outWS = makeOutput(*ws, dir, std::move(Xs), std::move(profileYs), std::move(profileEs));
408 // The actual profile might be of different size than what user
409 // specified.
410 Box actualBounds;
411 actualBounds.top = verticalBins[vertInterval.first];
412 actualBounds.bottom =
413 vertInterval.second < verticalBins.size() ? verticalBins[vertInterval.second] : verticalBins.back();
414 actualBounds.left = horizontalBins[horInterval.first];
415 actualBounds.right =
416 horInterval.second < horizontalBins.size() ? horizontalBins[horInterval.second] : horizontalBins.back();
417 setAxesAndUnits(*outWS, *ws, actualBounds, dir);
418 if (dir == LineDirection::vertical && ws->isDistribution()) {
419 divideByBinHeight(*outWS);
420 }
422}
423
426std::map<std::string, std::string> LineProfile::validateInputs() {
427 std::map<std::string, std::string> issues;
428 const double start = getProperty(PropertyNames::START);
429 const double end = getProperty(PropertyNames::END);
430 if (start > end) {
431 issues[PropertyNames::START] = PropertyNames::START + " greater than " + PropertyNames::END + ".";
432 }
434 if (ws->getAxis(1)->isText()) {
435 issues[PropertyNames::INPUT_WORKSPACE] = "The vertical axis in " + PropertyNames::INPUT_WORKSPACE + " is text.";
436 }
437 return issues;
438}
439
440} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double height
Definition: GetAllEi.cpp:155
size_t lineEnd
Definition: LineProfile.cpp:87
double bottom
Definition: LineProfile.cpp:79
double top
Definition: LineProfile.cpp:78
double left
Definition: LineProfile.cpp:80
size_t widthStart
Definition: LineProfile.cpp:88
size_t widthEnd
Definition: LineProfile.cpp:89
double right
Definition: LineProfile.cpp:81
size_t lineStart
Definition: LineProfile.cpp:86
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
Class to represent the axis of a workspace.
Definition: Axis.h:30
Stores numeric values that are assumed to be bin edge values.
Definition: BinEdgeAxis.h:20
A validator which provides a TENTATIVE check that a workspace contains common bins in each spectrum.
A validator which checks that the X axis of a workspace is increasing from left to right.
Base MatrixWorkspace Abstract Class.
A property class for workspaces.
LineProfile : Calculates a horizontal or vertical line profile over a MatrixWorkspace.
Definition: LineProfile.h:18
void init() override
Initialize the algorithm's properties.
int version() const override
Algorithm's version for identification.
const std::string category() const override
Algorithm's category for identification.
void exec() override
Execute the algorithm.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
std::map< std::string, std::string > validateInputs() override
Validate the algorithm's inputs.
Concrete workspace implementation.
Definition: Workspace2D.h:29
BoundedValidator is a validator that requires the values to be between upper or lower bounds,...
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
ListValidator is a validator that requires the value of a property to be one of a defined list of pos...
Definition: ListValidator.h:29
Validator to check that a property is not left empty.
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
static const std::string OUTPUT_WORKSPACE
static const std::string INPUT_WORKSPACE
std::string CENTRE("Centre")
std::shared_ptr< Workspace2D > Workspace2D_sptr
shared pointer to Mantid::DataObjects::Workspace2D
std::unique_ptr< T > create(const P &parent, const IndexArg &indexArg, const HistArg &histArg)
This is the create() method that all the other create() methods call.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
STL namespace.
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