Mantid
Loading...
Searching...
No Matches
RadiusSum.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 +
18#include "MantidKernel/Unit.h"
21
22#include "MantidHistogramData/LinearGenerator.h"
23
24#include <array>
25#include <cmath>
26#include <limits>
27#include <numeric>
28#include <sstream>
29
30using namespace Mantid::Kernel;
31
32namespace Mantid::Algorithms {
33
34// Register the algorithm into the AlgorithmFactory
35DECLARE_ALGORITHM(RadiusSum)
36
37
38const std::string RadiusSum::name() const { return "RadiusSum"; }
39
41int RadiusSum::version() const { return 1; }
42
44const std::string RadiusSum::category() const { return "Transforms\\Grouping"; }
45
50 std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
51 "An input workspace.");
52 declareProperty(std::make_unique<API::WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
53 "An output workspace.");
54
55 auto twoOrThreeElements = std::make_shared<ArrayLengthValidator<double>>(2, 3);
56 std::vector<double> myInput(3, 0);
57 declareProperty(std::make_unique<ArrayProperty<double>>("Centre", std::move(myInput), std::move(twoOrThreeElements)),
58 "Coordinate of the centre of the ring");
59
60 auto nonNegative = std::make_shared<BoundedValidator<double>>();
61 nonNegative->setLower(0);
62 declareProperty("MinRadius", 0.0, nonNegative->clone(), "Length of the inner ring. Default=0");
63 declareProperty(std::make_unique<PropertyWithValue<double>>("MaxRadius", std::numeric_limits<double>::max(),
64 std::move(nonNegative)),
65 "Length of the outer ring. Default=ImageSize.");
66
67 auto nonNegativeInt = std::make_shared<BoundedValidator<int>>();
68 nonNegativeInt->setLower(1);
69 declareProperty("NumBins", 100, std::move(nonNegativeInt), "Number of slice bins for the output. Default=100");
70
71 const char *normBy = "NormalizeByRadius";
72 const char *normOrder = "NormalizationOrder";
73
74 declareProperty(normBy, false,
75 "Divide the sum of each ring by the radius "
76 "powered by Normalization Order");
77 declareProperty(normOrder, 1.0,
78 "If 2, the normalization will be divided by "
79 "the quadratic value of the ring for each "
80 "radius.");
81 setPropertySettings(normOrder, std::make_unique<VisibleWhenProperty>(normBy, IS_EQUAL_TO, "1"));
82
83 const char *groupNorm = "Normalization";
84 setPropertyGroup(normBy, groupNorm);
85 setPropertyGroup(normOrder, groupNorm);
86}
87
93
94 // this is the main method, that executes the algorithm of radius sum
95 std::vector<double> output;
98 else
100
101 if (getProperty("NormalizeByRadius"))
102 normalizeOutputByRadius(output, getProperty("NormalizationOrder"));
103
104 setUpOutputWorkspace(output);
105}
106
108 g_log.debug() << "Process Instrument related image\n";
109
110 std::vector<double> accumulator(num_bins, 0);
111
112 g_log.debug() << "For every detector in the image get its position "
113 << " and sum up all the counts inside the related spectrum\n";
114
115 const auto &spectrumInfo = inputWS->spectrumInfo();
116 for (size_t i = 0; i < inputWS->getNumberHistograms(); i++) {
117 if (!spectrumInfo.hasDetectors(i)) {
118 // it may occur because there is no detector assigned, but it does not
119 // cause problem. Hence, continue the loop.
120 g_log.information() << "Spectrum " << i << " has no detector assigned.\n";
121 continue;
122 }
123 if (spectrumInfo.isMonitor(i))
124 continue;
125
126 int bin_n = getBinForPixelPos(spectrumInfo.position(i));
127
128 if (bin_n < 0)
129 continue; // not in the limits of min_radius and max_radius
130
131 auto &refY = inputWS->getSpectrum(i).y();
132 accumulator[bin_n] += std::accumulate(refY.cbegin(), refY.cend(), 0.0);
133 }
134 return accumulator;
135}
136
138 g_log.debug() << "Process Numeric Image\n";
139
140 std::vector<double> accumulator(num_bins, 0);
141
142 // the position of the pixel in the vertical axis comes from the axis(1)
143 // (NumericAxis)
144 API::NumericAxis *verticalAxis = dynamic_cast<API::NumericAxis *>(inputWS->getAxis(1));
145
146 // assumption: in a numeric image, all the bins position for every row in the
147 // image
148 // will be in the same place.
149 // The positions in the horizontal axis is the position of X if it is not a
150 // histogram,
151 // or the center of the bin, if it is a histogram.
152
153 g_log.debug() << "Define the X positions of the pixels\n";
154 auto x_pos = inputWS->points(0);
155
156 g_log.debug() << "For every pixel define its bin position and sum them up\n";
157 // for each row in the image
158 for (size_t i = 0; i < inputWS->getNumberHistograms(); i++) {
159
160 // pixel values
161 auto &refY = inputWS->getSpectrum(i).y();
162
163 // for every pixel
164 for (size_t j = 0; j < refY.size(); j++) {
165
166 // the position of the pixel is the horizontal position of the pixel
167 // and the vertical row position.
168 V3D pixelPos(x_pos[j], (*verticalAxis)(i), 0);
169
170 int bin_pos = getBinForPixelPos(pixelPos);
171
172 if (bin_pos < 0)
173 continue; // not in the region [min_radius, max_radius]
174
175 accumulator[bin_pos] += refY[j];
176 }
177 }
178
179 return accumulator;
180}
181
183 double radio, theta, phi;
184 V3D diff_vector = pos - centre;
185 diff_vector.getSpherical(radio, theta, phi);
186
187 // the distance is the projection of the radio in the plane of the image
188 // which is given by radio * sin (theta)
189 double effect_distance = radio * sin(theta * M_PI / 180);
190
191 if (effect_distance < min_radius || effect_distance > max_radius)
192 return -1; // outside the limits [min_radius, max_radius]
193
194 return fromDistanceToBin(effect_distance);
195}
196
198 g_log.debug() << "Copy the input property values\n";
199 inputWS = getProperty("InputWorkspace");
200
201 g_log.debug() << "Extract the center and make it a V3D object\n";
202 std::vector<double> centre_aux = getProperty("Centre");
203 if (centre_aux.size() == 2)
204 centre = V3D(centre_aux[0], centre_aux[1], 0);
205 else
206 centre = V3D(centre_aux[0], centre_aux[1], centre_aux[2]);
207
208 g_log.debug() << "Copy the remaning properties: MinRadius, MaxRadius and NumBins\n";
209 min_radius = getProperty("MinRadius");
210 max_radius = getProperty("MaxRadius");
211 num_bins = getProperty("NumBins");
212}
213
215 g_log.debug() << "Sanity check\n";
216
217 g_log.debug() << "Check MinRadius < MaxRadius\n";
218 if (min_radius >= max_radius) {
219 std::stringstream s;
220 s << "Wrong definition of the radius min and max. The minimum radius can "
221 "not be bigger than maximum. "
222 << "\nInputs (" << min_radius << ", " << max_radius << ").\n";
223 throw std::invalid_argument(s.str());
224 }
225
226 const std::vector<double> boundary_limits = getBoundariesOfInputWorkspace();
227 std::stringstream s;
228 std::copy(boundary_limits.begin(), std::prev(boundary_limits.end()), std::ostream_iterator<double>(s, ", "));
229 s << boundary_limits.back();
230 g_log.information() << "Boundary limits are: " << s.str() << '\n';
231
232 g_log.debug() << "Check: centre is defined inside the region defined by the "
233 "image or instrument\n";
234 centerIsInsideLimits(getProperty("centre"), boundary_limits);
235
236 g_log.debug() << "Recalculate MaxRadius if default value is given\n";
237 if (max_radius > 0.9 * std::numeric_limits<double>::max()) {
238 max_radius = getMaxDistance(centre, boundary_limits);
239 g_log.notice() << "RadiusMax automatically calculated and set to " << max_radius << '\n';
240 }
241
242 g_log.debug() << "Check number of bins to alert user if many bins will end up empty\n";
244}
245
260 return inWS->getAxis(1)->isSpectra();
261}
262
266 } else {
268 }
269}
270
289
290 // horizontal axis
291
292 // get the pixel position in the horizontal axis from the middle of the image.
293 const auto &refX = inWS->x(inWS->getNumberHistograms() / 2);
294
295 double min_x, max_x;
296
297 const double &first_x = refX.front();
298 const double &last_x = refX.back();
299 if (first_x < last_x) {
300 min_x = first_x;
301 max_x = last_x;
302 } else {
303 min_x = last_x;
304 max_x = first_x;
305 }
306
307 // vertical axis
308 API::NumericAxis *verticalAxis = dynamic_cast<API::NumericAxis *>(inWS->getAxis(1));
309 if (!verticalAxis)
310 throw std::invalid_argument("Vertical axis is not a numeric axis. Can not "
311 "find the limits of the image.");
312
313 double min_y, max_y;
314 min_y = verticalAxis->getMin();
315 max_y = verticalAxis->getMax();
316
317 // check the assumption made that verticalAxis will provide the correct
318 // answer.
319 if (min_y > max_y) {
320 throw std::logic_error("Failure to get the boundaries of this image. "
321 "Internal logic error. Please, inform MantidHelp");
322 }
323
324 std::vector<double> output(4); // output = {min_x, max_x, min_y, max_y}; not
325 // supported in all compilers
326 output[0] = min_x;
327 output[1] = max_x;
328 output[2] = min_y;
329 output[3] = max_y;
330 return output;
331}
332
344
345 // This function is implemented based in the following assumption:
346 // - The workspace is composed by spectrum with associated spectrum No which
347 // is associated to one detector or monitor
348 // - The first spectrum No (non monitor) is associated with one detector
349 // while the last spectrum No (non monitor)
350 // is associated with one detector.
351 // - They are in complete oposite direction.
352 //
353 // Consider the following 'image' (where the ID is the number and the
354 // position is where it is displayed)
355 //
356 // 1 2 3
357 // 4 5 6
358 // 7 8 9
359 // 10 11 12
360 //
361 // In this image, the assumption is true, because, we can derive the
362 // boundaries of the image looking just to the
363 // ids 1 and 12.
364 //
365 // But the following image:
366 //
367 // 1 2 3 6 5 4
368 // 6 5 4 1 2 3
369 // 7 8 9 12 11 10
370 // 12 11 12 7 8 9
371 //
372 // Although valid 'IDF' instrument, fail the assumption, and will return
373 // wrong values.
374 // Bear in mind these words if you face problems with the return of the
375 // boundaries of one instrument
376 //
377
378 double first_x, first_y, first_z;
379 size_t i = 0;
380 const auto &spectrumInfo = inWS->spectrumInfo();
381 while (true) {
382 i++;
383 if (i >= inWS->getNumberHistograms())
384 throw std::invalid_argument("Did not find any non monitor detector. "
385 "Failed to identify the boundaries of this "
386 "instrument.");
387
388 if (spectrumInfo.isMonitor(i))
389 continue;
390
391 const auto pos = spectrumInfo.position(i);
392 // get the position of the first valid (non-monitor) detector.
393 first_x = pos.X();
394 first_y = pos.Y();
395 first_z = pos.Z();
396 break;
397 }
398
399 double last_x, last_y, last_z;
400 i = inWS->getNumberHistograms() - 1;
401 while (true) {
402 i--;
403 if (i == 0)
404 throw std::invalid_argument("There is no region defined for the "
405 "instrument of this workspace. Failed to "
406 "identify the boundaries of this instrument");
407
408 if (spectrumInfo.isMonitor(i))
409 continue;
410
411 const auto pos = spectrumInfo.position(i);
412 // get the last valid detector position
413 last_x = pos.X();
414 last_y = pos.Y();
415 last_z = pos.Z();
416 break;
417 }
418
419 // order the values
420 double xMax, yMax, zMax;
421 double xMin, yMin, zMin;
422 xMax = std::max(first_x, last_x);
423 yMax = std::max(first_y, last_y);
424 zMax = std::max(first_z, last_z);
425 xMin = std::min(first_x, last_x);
426 yMin = std::min(first_y, last_y);
427 zMin = std::min(first_z, last_z);
428
429 std::vector<double> output(6); // output = {xMin, xMax, yMin, yMax, zMin,
430 // zMax }; not supported in all compilers
431 output[0] = xMin;
432 output[1] = xMax;
433 output[2] = yMin;
434 output[3] = yMax;
435 output[4] = zMin;
436 output[5] = zMax;
437
438 return output;
439}
440
462void RadiusSum::centerIsInsideLimits(const std::vector<double> &centre, const std::vector<double> &boundaries) {
463 // sanity check
464 if ((2 * centre.size()) != boundaries.size())
465 throw std::invalid_argument("CenterIsInsideLimits: The centre and boundaries were not given in the "
466 "correct as {x1,x2,...} and {x1_min, x1_max, x2_min, x2_max, ... }");
467
468 bool check_passed = true;
469 std::stringstream s;
470
471 for (size_t i = 0; i < 2; i++) { // TODO: fix this
472 if (centre[i] < boundaries[2 * i] || centre[i] > boundaries[2 * i + 1]) {
473 check_passed = false;
474 s << "The position for axis " << i + 1 << " (" << centre[i] << ") is outside the limits [" << boundaries[2 * i]
475 << ", " << boundaries[2 * i + 1] << "]. \n";
476 }
477 }
478
479 if (!check_passed)
480 throw std::invalid_argument(s.str());
481}
482
484 double size_bins = (max_radius - min_radius) / num_bins;
485
486 double min_bin_size;
488 min_bin_size = getMinBinSizeForInstrument(inputWS);
489 else
491
492 if (size_bins < min_bin_size)
493 g_log.warning() << "The number of bins provided is too big. "
494 << "It corresponds to a separation smaller than the image "
495 "resolution (detector size). "
496 << "A resonable number is smaller than "
497 << static_cast<int>((max_radius - min_radius) / min_bin_size) << '\n';
498}
499
501 // Assumption made: the detectors are placed one after the other, so the
502 // minimum
503 // reasonalbe size for the bin is the width of one detector.
504
505 const auto &spectrumInfo = inWS->spectrumInfo();
506 for (size_t i = 0; i < inWS->getNumberHistograms(); ++i) {
507 if (spectrumInfo.isMonitor(i))
508 continue;
510 spectrumInfo.detector(i).getBoundingBox(bbox);
511 return bbox.width().norm();
512 }
513 // this should never happen because it was done in getBoundariesOfInstrument,
514 throw std::invalid_argument("Did not find any non monitor detector position");
515}
516
518 // The pixel dimensions:
519 // - width: image width/ number of pixels in one row
520 // - height: image height/ number of pixels in one column
521 // The minimum bin size is the smallest value between this two values.
522
523 std::vector<double> boundaries = getBoundariesOfNumericImage(inWS);
524 const auto &refX = inWS->x(inputWS->getNumberHistograms() / 2);
525 auto nX = static_cast<int>(refX.size());
526 int nY = static_cast<int>(inWS->getAxis(1)->length());
527
528 // remembering boundaries is defined as { xMin, xMax, yMin, yMax}
529 return std::min(((boundaries[1] - boundaries[0]) / nX), ((boundaries[3] - boundaries[2]) / nY));
530}
531
532void RadiusSum::normalizeOutputByRadius(std::vector<double> &values, double exp_power) {
533 g_log.debug() << "Normalization of the output in relation to the 'radius' (distance)\n";
534
535 // the radius can be defined as:
536 // radius_min + bin_size/2 + n * bin_size ; for 0<= n <= num_bins
537 double bin_size = (max_radius - min_radius) / num_bins;
538 double first_radius = min_radius + bin_size / 2;
539
540 g_log.debug() << "Calculate Output[i] = Counts[i] / (Radius[i] ^ " << exp_power << ") << \n";
541 if (exp_power > 1.00001 || exp_power < 0.99999) {
542 for (int i = 0; i < static_cast<int>(values.size()); i++) {
543 values[i] = values[i] / pow(first_radius + i * bin_size, exp_power);
544 }
545 } else { // avoid calling pow because exp_power == 1 (for performance)
546 for (int i = 0; i < static_cast<int>(values.size()); i++) {
547 values[i] = values[i] / (first_radius + i * bin_size);
548 }
549 }
550}
551
552double RadiusSum::getMaxDistance(const V3D &centre, const std::vector<double> &boundary_limits) {
553
554 std::array<double, 2> Xs = {{boundary_limits[0], boundary_limits[1]}};
555 std::array<double, 2> Ys = {{boundary_limits[2], boundary_limits[3]}};
556 std::array<double, 2> Zs = {{0., 0.}};
557
558 if (boundary_limits.size() == 6) {
559 Zs[0] = boundary_limits[4];
560 Zs[1] = boundary_limits[5];
561 }
562
563 double max_distance = 0;
564 for (auto &x : Xs)
565 for (auto &y : Ys)
566 for (auto &z : Zs) {
567 // define all the possible combinations for the limits
568
569 double curr_distance = centre.distance(V3D(x, y, z));
570
571 if (curr_distance > max_distance)
572 max_distance = curr_distance; // keep the maximum distance.
573 }
574 return max_distance;
575}
576
577void RadiusSum::setUpOutputWorkspace(const std::vector<double> &values) {
578
579 g_log.debug() << "Output calculated, setting up the output workspace\n";
580
582 API::WorkspaceFactory::Instance().create(inputWS, 1, values.size() + 1, values.size());
583
584 g_log.debug() << "Set the data\n";
585 outputWS->mutableY(0) = values;
586
587 g_log.debug() << "Set the bins limits\n";
588
589 auto xSize = outputWS->mutableX(0).size();
590 double bin_size = (max_radius - min_radius) / num_bins;
591 outputWS->setBinEdges(0, xSize, HistogramData::LinearGenerator(min_radius, bin_size));
592
593 // configure the axis:
594 // for numeric images, the axis are the same as the input workspace, and are
595 // copied in the creation.
596
597 // for instrument related, the axis Y (1) continues to be the same.
598 // it is necessary to change only the axis X. We have to change it to radius.
600 auto horizontal = std::make_unique<API::NumericAxis>(xSize);
601 auto labelX = UnitFactory::Instance().create("Label");
602 std::dynamic_pointer_cast<Units::Label>(labelX)->setLabel("Radius");
603 horizontal->unit() = labelX;
604 outputWS->replaceAxis(0, std::move(horizontal));
605 }
606
607 setProperty("OutputWorkspace", outputWS);
608}
609
610} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
Kernel::Logger & g_log
Definition: Algorithm.h:451
Class to represent a numeric axis of a workspace.
Definition: NumericAxis.h:29
double getMax() const override
returns max value defined on axis
Definition: NumericAxis.h:54
double getMin() const override
returns min value defined on axis
Definition: NumericAxis.h:52
A property class for workspaces.
int fromDistanceToBin(double distance)
Return the bin position for a given distance From the input, it is defined the limits of distances as...
Definition: RadiusSum.h:94
double getMaxDistance(const Kernel::V3D &centre, const std::vector< double > &boundary_limits)
Definition: RadiusSum.cpp:552
static void centerIsInsideLimits(const std::vector< double > &centre, const std::vector< double > &boundaries)
Check if a given position is inside the limits defined by the boundaries.
Definition: RadiusSum.cpp:462
std::vector< double > processNumericImageRadiusSum()
Definition: RadiusSum.cpp:137
void normalizeOutputByRadius(std::vector< double > &values, double exp_power)
Definition: RadiusSum.cpp:532
static bool inputWorkspaceHasInstrumentAssociated(const API::MatrixWorkspace_sptr &)
Differentiate between Instrument related image (where the position of the pixels depend on the instru...
Definition: RadiusSum.cpp:259
void init() override
Initialize the algorithm's properties.
Definition: RadiusSum.cpp:48
double getMinBinSizeForInstrument(const API::MatrixWorkspace_sptr &)
Definition: RadiusSum.cpp:500
API::MatrixWorkspace_sptr inputWS
Definition: RadiusSum.h:65
double getMinBinSizeForNumericImage(const API::MatrixWorkspace_sptr &)
Definition: RadiusSum.cpp:517
const std::string category() const override
Algorithm's category for identification.
Definition: RadiusSum.cpp:44
int version() const override
Algorithm's version for identification.
Definition: RadiusSum.cpp:41
static std::vector< double > getBoundariesOfNumericImage(const API::MatrixWorkspace_sptr &)
Assuming that the input workspace is a Numeric Image where the pixel positions depend on their relati...
Definition: RadiusSum.cpp:288
void setUpOutputWorkspace(const std::vector< double > &values)
Definition: RadiusSum.cpp:577
std::vector< double > processInstrumentRadiusSum()
Definition: RadiusSum.cpp:107
void exec() override
Execute the algorithm.
Definition: RadiusSum.cpp:90
int getBinForPixelPos(const Kernel::V3D &pos)
Definition: RadiusSum.cpp:182
static std::vector< double > getBoundariesOfInstrument(const API::MatrixWorkspace_sptr &)
Assuming that the workspace has an instrument associated with it from which the pixel positions has t...
Definition: RadiusSum.cpp:343
std::vector< double > getBoundariesOfInputWorkspace()
Definition: RadiusSum.cpp:263
A simple structure that defines an axis-aligned cuboid shaped bounding box for a geometrical object.
Definition: BoundingBox.h:34
Kernel::V3D width() const
Returns the width of the box.
Definition: BoundingBox.h:98
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 setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
void setPropertyGroup(const std::string &name, const std::string &group)
Set the group for a given property.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
The concrete, templated class for properties.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Class for 3D vectors.
Definition: V3D.h:34
double distance(const V3D &v) const noexcept
Calculates the distance between two vectors.
Definition: V3D.h:287
double norm() const noexcept
Definition: V3D.h:263
void getSpherical(double &R, double &theta, double &phi) const noexcept
Return the vector's position in spherical coordinates.
Definition: V3D.cpp:117
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54