Mantid
Loading...
Searching...
No Matches
RingProfile.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 +
11#include "MantidAPI/TextAxis.h"
19#include <climits>
20#include <cmath>
21
22namespace Mantid::Algorithms {
23
24// Register the algorithm into the AlgorithmFactory
25DECLARE_ALGORITHM(RingProfile)
26
27
30 : min_radius(-1.), max_radius(-1), start_angle(-1), clockwise(false), num_bins(-1), centre_x(-1.), centre_y(-1.),
31 centre_z(-1), bin_size(-1) {}
32
47 std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>("InputWorkspace", "", Kernel::Direction::Input),
48 "An input workspace.");
49 declareProperty(std::make_unique<API::WorkspaceProperty<>>("OutputWorkspace", "", Kernel::Direction::Output),
50 "An output workspace.");
51
52 auto twoOrThree = std::make_shared<Kernel::ArrayLengthValidator<double>>(2, 3);
53 std::vector<double> myInput(3, 0);
54 declareProperty(std::make_unique<Kernel::ArrayProperty<double>>("Centre", std::move(myInput), std::move(twoOrThree)),
55 "Coordinate of the centre of the ring");
56 auto nonNegative = std::make_shared<Kernel::BoundedValidator<double>>();
57 nonNegative->setLower(0);
58
59 declareProperty<double>("MinRadius", 0, nonNegative->clone(), "Radius of the inner ring(m)");
60 declareProperty(std::make_unique<Kernel::PropertyWithValue<double>>("MaxRadius", std::numeric_limits<double>::max(),
61 std::move(nonNegative)),
62 "Radius of the outer ring(m)");
63 auto nonNegativeInt = std::make_shared<Kernel::BoundedValidator<int>>();
64 nonNegativeInt->setLower(1);
65 declareProperty<int>("NumBins", 100, std::move(nonNegativeInt), "Number of slice bins for the output");
66 auto degreesLimits = std::make_shared<Kernel::BoundedValidator<double>>();
67 degreesLimits->setLower(-360);
68 degreesLimits->setUpper(360);
69 declareProperty<double>("StartAngle", 0, degreesLimits, "The angle to start from.");
70
71 std::vector<std::string> op(2);
72 op[0] = "ClockWise";
73 op[1] = "Anti-ClockWise";
74 declareProperty("Sense", "Anti-ClockWise", std::make_shared<Kernel::StringListValidator>(op),
75 "The direction of the integration around the ring");
76}
77
78//----------------------------------------------------------------------------------------------
102 // get input workspace
103 API::MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
104
105 // the RingProfile does not support eventworkspace
106 auto checkEvent = std::dynamic_pointer_cast<API::IEventWorkspace>(inputWS);
107 if (checkEvent) {
108 throw std::invalid_argument("RingProfile is not defined for EventWorkspaces.");
109 }
110 g_log.debug() << "Get the input parameters \n";
111 // get the algorithm parameters
112 std::vector<double> centre = getProperty("Centre");
113 centre_x = centre[0];
114 centre_y = centre[1];
115 if (centre.size() == 3)
116 centre_z = centre[2];
117 min_radius = getProperty("MinRadius");
118 max_radius = getProperty("MaxRadius");
119 start_angle = getProperty("StartAngle");
120 num_bins = getProperty("NumBins");
121 bin_size = 360.0 / num_bins;
122 clockwise = (getPropertyValue("Sense") == "ClockWise");
123
124 g_log.debug() << "Check the inputs of the algorithm\n";
125 // VALIDATE THE INPUTS
126 if (inputWS->getAxis(1)->isSpectra()) {
128 } else {
130 }
131
132 m_progress = std::shared_ptr<API::Progress>(new API::Progress(this, 0.0, 1.0, inputWS->getNumberHistograms() + 1));
133
134 // prepare the vector to hold the output
135 std::vector<double> output_bins(num_bins, 0);
136
137 g_log.debug() << "Execute the ring profile calculation\n";
138 // perform the ring profile calculation
139 if (inputWS->getAxis(1)->isSpectra()) {
140 processInstrumentRingProfile(inputWS, output_bins);
141
142 } else {
143 processNumericImageRingProfile(inputWS, output_bins);
144 }
145
146 g_log.debug() << "Prepare the output\n";
147 // create the output
149 API::WorkspaceFactory::Instance().create(inputWS, 1, output_bins.size() + 1, output_bins.size());
150 m_progress->report("Preparing the output");
151 // populate Y data getting the values from the output_bins
152 MantidVec &refY = outputWS->dataY(0);
153 if (clockwise) {
154 for (size_t j = 0; j < output_bins.size(); j++)
155 refY[j] = output_bins[output_bins.size() - j - 1];
156 } else {
157 for (size_t j = 0; j < output_bins.size(); j++)
158 refY[j] = output_bins[j];
159 }
160
161 // prepare the output for the angles:
162 // populate X data
163 MantidVec &refX = outputWS->dataX(0);
164
165 std::vector<double> angles(output_bins.size() + 1);
166 // we always keep the angle in relative from where it starts and growing in
167 // its sense.
168 for (int j = 0; j < static_cast<int>(angles.size()); j++)
169 refX[j] = bin_size * j;
170
171 // configure the axis
172
173 // the horizontal axis is configured as degrees and copy the values of X
174 auto horizontal = std::make_unique<API::NumericAxis>(refX.size());
175 horizontal->unit() = std::make_shared<Kernel::Units::Phi>();
176 horizontal->title() = "Ring Angle";
177 for (size_t j = 0; j < refX.size(); j++)
178 horizontal->setValue(j, refX[j]);
179 outputWS->replaceAxis(0, std::move(horizontal));
180
181 // the vertical axis get the same unit and information from the input
182 // workspace
183 auto verticalAxis = std::make_unique<API::TextAxis>(1);
184 verticalAxis->unit() = inputWS->getAxis(1)->unit();
185 verticalAxis->title() = inputWS->getAxis(1)->title();
186 outputWS->replaceAxis(1, std::move(verticalAxis));
187
188 // set up the output
189 setProperty("OutputWorkspace", outputWS);
190}
191
206 try {
207 // finding the limits of the instrument
208 const auto &spectrumInfo = inputWS->spectrumInfo();
209 double first_x, first_y, first_z;
210 size_t i = 0;
211 while (true) {
212 i++;
213 if (i >= inputWS->getNumberHistograms())
214 throw std::invalid_argument("Did not find any non monitor detector position");
215
216 if (spectrumInfo.isMonitor(i))
217 continue;
218 const auto pos = spectrumInfo.position(i);
219 first_x = pos.X();
220 first_y = pos.Y();
221 first_z = pos.Z();
222 break;
223 }
224
225 double last_x, last_y, last_z;
226 i = inputWS->getNumberHistograms() - 1;
227 while (true) {
228 i--;
229 if (i == 0)
230 throw std::invalid_argument("There is no region defined for the instrument of this workspace");
231
232 if (spectrumInfo.isMonitor(i))
233 continue;
234 const auto pos = spectrumInfo.position(i);
235 last_x = pos.X();
236 last_y = pos.Y();
237 last_z = pos.Z();
238 break;
239 }
240
241 double xMax, yMax, zMax;
242 double xMin, yMin, zMin;
243 xMax = std::max(first_x, last_x);
244 yMax = std::max(first_y, last_y);
245 zMax = std::max(first_z, last_z);
246 xMin = std::min(first_x, last_x);
247 yMin = std::min(first_y, last_y);
248 zMin = std::min(first_z, last_z);
249
250 std::stringstream limits_s;
251 limits_s << "([" << xMin << ", " << xMax << "], [" << yMin << ", " << yMax << "], [" << zMin << ", " << zMax
252 << "])";
253 g_log.debug() << "The limits for the instrument is : " << limits_s.str() << '\n';
254 int xOutside = 0, yOutside = 0, zOutside = 0;
255 if (centre_x < xMin || centre_x > xMax)
256 xOutside = 1;
257 if (centre_y < yMin || centre_y > yMax)
258 yOutside = 1;
259 if (centre_z < zMin || centre_z > zMax)
260 zOutside = 1;
261 int summed = xOutside + yOutside + zOutside;
262 // if at least 2 are outside, the centre is considered outside the box.
263 if (summed >= 2) {
264 std::stringstream s;
265 s << "The defined centre (" << centre_x << ", " << centre_y << ", " << centre_z
266 << ") is outside the limits of the detectors inside this instrument: " << limits_s.str();
267 throw std::invalid_argument(s.str());
268 }
269
270 xOutside = yOutside = zOutside = 0;
271 if (centre_x - min_radius > xMax || centre_x + min_radius < xMin)
272 xOutside = 1;
273 if (centre_y - min_radius > yMax || centre_y + min_radius < yMin)
274 yOutside = 1;
275 if (centre_z - min_radius > zMax || centre_z + min_radius < zMin)
276 zOutside = 1;
277
278 summed = xOutside + yOutside + zOutside;
279
280 if (summed >= 2) {
281 std::stringstream s;
282 s << "The defined minRadius make the inner ring outside the limits of "
283 "the detectors inside this instrument: "
284 << limits_s.str();
285 throw std::invalid_argument(s.str());
286 }
287
289 throw std::invalid_argument("Invalid input workspace. This workspace does "
290 "not has detectors to get the positions "
291 "from. ");
292 }
293}
294
310 g_log.notice() << "CheckingInputs For Numeric Workspace\n";
311
312 // The Axis0 is defined by the values of readX inside the spectra of the
313 // workspace.
314 // The limits of this axis will be get by inspection of the readX vector
315 // taking the first
316 // and the last value.
317
318 // check that centre is inside the range available for the instrument
319 const MantidVec &refX = inputWS->readX(inputWS->getNumberHistograms() / 2);
320 // get the limits of the axis 0 (X)
321 double min_v_x, max_v_x;
322 min_v_x = std::min(refX[0], refX.back());
323 max_v_x = std::max(refX[0], refX.back());
324 g_log.notice() << "Limits X = " << min_v_x << " " << max_v_x << '\n';
325 // check centre is inside the X domain
326 if (centre_x < min_v_x || centre_x > max_v_x) {
327 std::stringstream s;
328 s << "The input value for centre (X=" << centre_x << ") is outside the limits of the instrument [" << min_v_x
329 << ", " << max_v_x << "]";
330 throw std::invalid_argument(s.str());
331 }
332
333 // The Axis1 is defined by the spectra inside the workspace. Its limits and
334 // values are given by the
335 // ws->getAxis(1)
336
337 // get the limits of the axis1 (Y)
338 API::NumericAxis *oldAxis2 = dynamic_cast<API::NumericAxis *>(inputWS->getAxis(1));
339 // we cannot have the positions in Y direction without a NumericAxis
340 if (!oldAxis2)
341 throw std::invalid_argument("Vertical axis is not a numeric axis. If it is "
342 "a spectra axis try running "
343 "ConvertSpectrumAxis first.");
344 double min_v_y = std::min(oldAxis2->getMin(), oldAxis2->getMax());
345 double max_v_y = std::max(oldAxis2->getMin(), oldAxis2->getMax());
346 g_log.notice() << "Limits Y = " << min_v_y << " " << max_v_y << '\n';
347 // check centre is inside the Y domain
348 if (centre_y < min_v_y || centre_y > max_v_y) {
349 std::stringstream s;
350 s << "The input value for centre (Y=" << centre_y << ") is outside the limits of the instrument [" << min_v_y
351 << ", " << max_v_y << "]";
352 throw std::invalid_argument(s.str());
353 }
354 g_log.notice() << "Centre: " << centre_x << " " << centre_y << '\n';
355 // check minradius is inside the limits of the region of the instrument
356 if (centre_x - min_radius > max_v_x || centre_x + min_radius < min_v_x || centre_y - min_radius > max_v_y ||
357 centre_y + min_radius < min_v_y)
358 throw std::invalid_argument("The minimun radius is outside the region of the instrument");
359}
360
376 std::vector<double> &output_bins) {
377
378 const auto &spectrumInfo = inputWS->spectrumInfo();
379 for (int i = 0; i < static_cast<int>(inputWS->getNumberHistograms()); i++) {
380 m_progress->report("Computing ring bins positions for detectors");
381 // for the detector based, the positions will be taken from the detector
382 // itself.
383 if (!spectrumInfo.hasDetectors(i)) {
384 g_log.information() << "Spectrum " << i << " has no detector assigned.\n";
385 continue;
386 }
387
388 if (spectrumInfo.isMonitor(i))
389 continue;
390
391 // this part will be executed if the instrument is attached to the workspace
392
393 // get the bin position
394 int bin_n = getBinForPixel(spectrumInfo.position(i));
395
396 if (bin_n < 0) // -1 is the agreement for an invalid bin, or outside the
397 // ring being integrated
398 continue;
399
400 g_log.debug() << "Bin for the index " << i << " = " << bin_n << " Pos = " << spectrumInfo.position(i) << '\n';
401
402 const MantidVec &refY = inputWS->getSpectrum(i).dataY();
403 // accumulate the values of this spectrum inside this bin
404 for (double sp_ind : refY)
405 output_bins[bin_n] += sp_ind;
406 }
407}
408
427
429 V3D origin(centre_x, centre_y, centre_z);
430 V3D diff_vector = position - origin;
431 double radio, theta, phi;
432 // get the spherical values of the vector from center to detector position
433 diff_vector.getSpherical(radio, theta, phi);
434
435 // the distance from the centre to the ring will be calculated as radio *
436 // sin(theta).
437 double effect_distance = radio * sin(theta * M_PI / 180);
438
439 // g_log.debug() << "effect Distance = " << effect_distance << '\n';
440
441 // check if it is inside the ring defined by min_radius, max_radius
442 if (effect_distance < min_radius || effect_distance > max_radius || effect_distance == 0)
443 return -1;
444
445 // get the angle
446 // g_log.debug() << "The real angle is " << phi << '\n';
447 return fromAngleToBin(phi);
448}
449
463 std::vector<double> &output_bins) {
464 // allocate the bin positions vector
465 std::vector<int> bin_n(inputWS->dataY(0).size(), -1);
466
467 // consider that each spectrum is a row in the image
468 for (int i = 0; i < static_cast<int>(inputWS->getNumberHistograms()); i++) {
469 m_progress->report("Computing ring bins positions for pixels");
470 // get bin for the pixels inside this spectrum
471 // for each column of the image
472 getBinForPixel(inputWS, i, bin_n);
473
474 // accumulate the values from the spectrum to the target bin
475 // each column has it correspondend bin_position inside bin_n
476 const MantidVec &refY = inputWS->dataY(i);
477 for (size_t j = 0; j < bin_n.size(); j++) {
478
479 // is valid bin? No -> skip
480 if (bin_n[j] < 0)
481 continue;
482
483 // accumulate the values of this spectrum inside this bin
484 output_bins[bin_n[j]] += refY[j];
485 }
486 }
487}
488
516void RingProfile::getBinForPixel(const API::MatrixWorkspace_sptr &ws, int spectrum_index, std::vector<int> &bins_pos) {
517
518 if (bins_pos.size() != ws->dataY(spectrum_index).size())
519 throw std::runtime_error("Invalid bin positions vector");
520
521 API::NumericAxis *oldAxis2 = dynamic_cast<API::NumericAxis *>(ws->getAxis(1));
522 // assumption y position is the ws->getAxis(1)(spectrum_index)
523 if (!oldAxis2) {
524 throw std::logic_error("Failed to cast workspace axis to NumericAxis");
525 }
526
527 // calculate ypos, the difference of y - centre and the square of this
528 // difference
529 double ypos = (*oldAxis2)(spectrum_index);
530 double diffy = ypos - centre_y;
531 double diffy_quad = pow(diffy, 2.0);
532
533 // the reference to X bins (the limits for each pixel in the horizontal
534 // direction)
535 auto xvec = ws->dataX(spectrum_index);
536
537 // for each pixel inside this row
538 for (size_t i = 0; i < xvec.size() - 1; i++) {
539
540 double xpos = (xvec[i] + xvec[i + 1]) / 2.0; // the x position is the centre of the bins boundaries
541 double diffx = xpos - centre_x;
542 // calculate the distance => norm of pixel position - centre
543 double distance = sqrt(pow(diffx, 2.0) + diffy_quad);
544
545 // check if the distance is inside the ring
546 if (distance < min_radius || distance > max_radius || distance == 0) {
547 bins_pos[i] = -1;
548 continue;
549 }
550
551 double angle = atan2(diffy, diffx);
552
553 // call fromAngleToBin (radians)
554 bins_pos[i] = fromAngleToBin(angle, false);
555 }
556}
557
558/* Return the bin position for a given angle.
559 *
560 * The whole ring has 360 degree which is divided in NumBins bins.
561 * Hence, defining the bin_size as 360/NumBins, you have the dimension
562 * of each bin. This means, that the bins will follow the following rule:
563 *
564 * Bins(n) = [startAngle + n * bin_size, startAngle + (n+1) * bin_size]
565 *
566 * Having a given angle x, we need to find n so that:
567 *
568 * startAngle + n bin_size < x < startAngle + (n+1) bin_size
569 * n bin_size < x-startAngle < (n+1) bin_size
570 * n < (x-startAngle)/bin_size < n+1
571 *
572 * So, n = truncate (x-startAngle)/bin_size
573 *
574 * @param angle: the input angle
575 * @param degree: flag that indicates that the input angle is in degree (default
576 *= true)
577 * @return n
578 */
579int RingProfile::fromAngleToBin(double angle, bool degree) {
580 if (!degree)
581 angle *= (180 / M_PI);
582
583 angle -= start_angle;
584
585 if (angle < 0) {
586 while (angle < 0)
587 angle += 360;
588 }
589
590 angle /= bin_size;
591 return static_cast<int>(angle);
592}
593
594} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double position
Definition: GetAllEi.cpp:154
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
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
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
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
Calculates the sum of the counts against a circular ring.
Definition: RingProfile.h:23
double centre_x
copy of centre for axis(0)
Definition: RingProfile.h:66
bool clockwise
flag that indicate if Sense input was configure to clockwise
Definition: RingProfile.h:62
void exec() override
Execute the algorithm.
double centre_y
copy of centre for axis(1)
Definition: RingProfile.h:68
void init() override
Initialize the algorithm's properties.
Definition: RingProfile.cpp:45
void checkInputsForNumericWorkspace(const API::MatrixWorkspace_sptr &)
validate the inputs of the algorithm for 2d matrix based instrument
void checkInputsForSpectraWorkspace(const API::MatrixWorkspace_sptr &)
validate the inputs of the algorithm for instrument based workspace
void processInstrumentRingProfile(const API::MatrixWorkspace_sptr &inputWS, std::vector< double > &output_bins)
process ring profile for instrument based workspace
double min_radius
copy of the minRadius input
Definition: RingProfile.h:56
double bin_size
The size of the bins in angle that is equal to 360 / NumBins.
Definition: RingProfile.h:72
double max_radius
copy of the maxRadius input
Definition: RingProfile.h:58
int fromAngleToBin(double angle, bool degree=true)
get the bin position for the given angle
void getBinForPixel(const API::MatrixWorkspace_sptr &, int, std::vector< int > &)
identify the bin position for the given pixel in the image based workspace
double start_angle
copy of the StartAngle input
Definition: RingProfile.h:60
std::shared_ptr< API::Progress > m_progress
Definition: RingProfile.h:36
double centre_z
copy of centre for axis(2)
Definition: RingProfile.h:70
void processNumericImageRingProfile(const API::MatrixWorkspace_sptr &inputWS, std::vector< double > &output_bins)
process ring profile for image based workspace
int num_bins
copy of NumBins input
Definition: RingProfile.h:64
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
Exception for when an item is not found in a collection.
Definition: Exception.h:145
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 notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
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
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
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54