Mantid
Loading...
Searching...
No Matches
SlicingAlgorithm.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#include "MantidAPI/Run.h"
18#include "MantidKernel/System.h"
20
21#include <boost/regex.hpp>
22#include <utility>
23
24using namespace Mantid::Kernel;
25using namespace Mantid::API;
26using namespace Mantid::Geometry;
28
29namespace Mantid::MDAlgorithms {
30
34 : m_transform(), m_transformFromOriginal(), m_transformToOriginal(), m_transformFromIntermediate(),
35 m_transformToIntermediate(), m_axisAligned(true), m_outD(0), // unititialized and should be invalid
36 m_NormalizeBasisVectors(false) {}
37
41 std::string dimChars = getDimensionChars();
42
43 // --------------- Axis-aligned properties
44 // ---------------------------------------
45 declareProperty("AxisAligned", true, "Perform binning aligned with the axes of the input MDEventWorkspace?");
46 setPropertyGroup("AxisAligned", "Axis-Aligned Binning");
47 for (size_t i = 0; i < dimChars.size(); i++) {
48 std::string dim(" ");
49 dim[0] = dimChars[i];
50 std::string propName = "AlignedDim" + dim;
52 "Binning parameters for the " + Strings::toString(i) +
53 "th dimension.\n"
54 "Enter it as a comma-separated list of values with the format: "
55 "'name,minimum,maximum,number_of_bins'. Leave blank for NONE.");
56 setPropertySettings(propName, std::make_unique<VisibleWhenProperty>("AxisAligned", IS_EQUAL_TO, "1"));
57 setPropertyGroup(propName, "Axis-Aligned Binning");
58 }
59
60 // --------------- NON-Axis-aligned properties
61 // ---------------------------------------
62 std::string grpName = "Non-Aligned Binning";
63
64 auto ps = [] {
65 std::unique_ptr<IPropertySettings> settings =
66 std::make_unique<VisibleWhenProperty>("AxisAligned", IS_EQUAL_TO, "0");
67 return settings;
68 };
69 for (size_t i = 0; i < dimChars.size(); i++) {
70 std::string dim(" ");
71 dim[0] = dimChars[i];
72 std::string propName = "BasisVector" + dim;
74 "Description of the basis vector of the " + Strings::toString(i) +
75 "th output dimension."
76 "Format: 'name, units, x,y,z,..'.\n"
77 " name : string for the name of the output dimension.\n"
78 " units : string for the units of the output dimension.\n"
79 " x,y,z,...: vector definining the basis in the input dimensions "
80 "space.\n"
81 "Leave blank for NONE.");
82 setPropertySettings(propName, ps());
83 setPropertyGroup(propName, grpName);
84 }
85 declareProperty(std::make_unique<ArrayProperty<double>>("Translation", Direction::Input),
86 "Coordinates in the INPUT workspace that corresponds to "
87 "(0,0,0) in the OUTPUT workspace.\n"
88 "Enter as a comma-separated string.\n"
89 "Default: 0 in all dimensions (no translation).");
90
91 declareProperty(std::make_unique<ArrayProperty<double>>("OutputExtents", Direction::Input),
92 "The minimum, maximum edges of space of each dimension of "
93 "the OUTPUT workspace, as a comma-separated list");
94
95 declareProperty(std::make_unique<ArrayProperty<int>>("OutputBins", Direction::Input),
96 "The number of bins for each dimension of the OUTPUT workspace.");
97
98 declareProperty(std::make_unique<PropertyWithValue<bool>>("NormalizeBasisVectors", true, Direction::Input),
99 "Normalize the given basis vectors to unity. \n"
100 "If true, then a distance of 1 in the INPUT dimensions = 1 "
101 "in the OUTPUT dimensions.\n"
102 "If false, then a distance of norm(basis_vector) in the "
103 "INPUT dimension = 1 in the OUTPUT dimensions.");
104
105 declareProperty(std::make_unique<PropertyWithValue<bool>>("ForceOrthogonal", false, Direction::Input),
106 "Force the input basis vectors to form an orthogonal coordinate system. "
107 "Only works in 3 dimension!");
108
109 // For GUI niceness
110 setPropertyGroup("Translation", grpName);
111 setPropertyGroup("OutputExtents", grpName);
112 setPropertyGroup("OutputBins", grpName);
113 setPropertyGroup("NormalizeBasisVectors", grpName);
114 setPropertyGroup("ForceOrthogonal", grpName);
115 setPropertySettings("Translation", ps());
116 setPropertySettings("OutputExtents", ps());
117 setPropertySettings("OutputBins", ps());
118 setPropertySettings("NormalizeBasisVectors", ps());
119 setPropertySettings("ForceOrthogonal", ps());
120}
121
122//----------------------------------------------------------------------------------------------
136 std::string input = Strings::strip(str);
137 if (input.empty())
138 return;
139 if (input.size() < 3)
140 throw std::invalid_argument("Dimension string is too short to be valid: " + str);
141
142 // The current dimension index.
143 size_t dim = m_binDimensions.size();
144
145 size_t n_first_comma;
146
147 // Special case: accept dimension names [x,y,z]
148 if (input[0] == '[') {
149 // Find the name at the closing []
150 size_t n = input.find_first_of(']', 1);
151 if (n == std::string::npos)
152 throw std::invalid_argument("No closing ] character in the dimension name of : " + str);
153 // Find the comma after the name
154 n_first_comma = input.find_first_of(',', n);
155 if (n_first_comma == std::string::npos)
156 throw std::invalid_argument("No comma after the closing ] character in the dimension string: " + str);
157 } else
158 // Find the comma after the name
159 n_first_comma = input.find_first_of(',');
160
161 if (n_first_comma == std::string::npos)
162 throw std::invalid_argument("No comma in the dimension string: " + str);
163 if (n_first_comma == input.size() - 1)
164 throw std::invalid_argument("Dimension string ends in a comma: " + str);
165
166 // Get the entire name
167 std::string name = Strings::strip(input.substr(0, n_first_comma));
168 if (name.empty())
169 throw std::invalid_argument("name should not be blank.");
170
171 // Now remove the name and comma
172 // And split the rest of it
173 input = input.substr(n_first_comma + 1);
174 std::vector<std::string> strs;
175 boost::split(strs, input, boost::is_any_of(","));
176 if (strs.size() != this->m_inWS->getNumDims() + 1)
177 throw std::invalid_argument("Wrong number of values (expected 2 + # of "
178 "input dimensions) in the dimensions string: " +
179 str);
180
181 // Get the number of bins from
182 int numBins = m_numBins[dim];
183 if (numBins < 1)
184 throw std::invalid_argument("Number of bins for output dimension " + Strings::toString(dim) + " should be >= 1.");
185
186 // Get the min/max extents in this OUTPUT dimension
187 double min = m_minExtents[dim];
188 double max = m_maxExtents[dim];
189 double lengthInOutput = max - min;
190 if (lengthInOutput <= 0)
191 throw std::invalid_argument("The maximum extents for dimension " + Strings::toString(dim) + " should be > 0.");
192
193 // Create the basis vector with the right # of dimensions
194 VMD basis(this->m_inWS->getNumDims());
195 for (size_t d = 0; d < this->m_inWS->getNumDims(); d++)
196 if (!Strings::convert(strs[d + 1], basis[d]))
197 throw std::invalid_argument("Error converting argument '" + strs[d + 1] + "' in the dimensions string '" + str +
198 "' to a number.");
199
200 // If B was binned from A (m_originalWS), and we are binning C from B,
201 // convert the basis vector from B space -> A space
202 if (m_originalWS) {
203 // Turn basis vector into two points
204 VMD basis0(this->m_inWS->getNumDims());
205 VMD basis1 = basis;
206 // Convert the points to the original coordinates (from inWS to originalWS)
207 CoordTransform const *toOrig = m_inWS->getTransformToOriginal();
208 VMD origBasis0 = toOrig->applyVMD(basis0);
209 VMD origBasis1 = toOrig->applyVMD(basis1);
210 // New basis vector, now in the original workspace
211 basis = origBasis1 - origBasis0;
212 }
213
214 // Check on the length of the basis vector
215 double basisLength = basis.norm();
216 if (basisLength <= 0)
217 throw std::invalid_argument("direction should not be 0-length.");
218
219 // Normalize it to unity, if desized
220 double transformScaling = 1.0;
222 // A distance of 1 in the INPUT space = a distance of 1.0 in the OUTPUT
223 // space
224 basis.normalize();
225 transformScaling = 1.0;
226 } else
227 // A distance of |basisVector| in the INPUT space = a distance of 1.0 in the
228 // OUTPUT space
229 transformScaling = (1.0 / basisLength);
230
231 // This is the length of this dimension as measured in the INPUT space
232 double lengthInInput = lengthInOutput / transformScaling;
233
234 // Scaling factor, to convert from units in the INPUT dimensions to the output
235 // BIN number
236 double binningScaling = double(numBins) / (lengthInInput);
237
238 // Extract the arguments
239 std::string units = Strings::strip(strs[0]);
240
241 // Create the appropriate frame
242 auto frame = createMDFrameForNonAxisAligned(units, basis);
243
244 // Create the output dimension
245 auto out = std::make_shared<MDHistoDimension>(name, name, *frame, static_cast<coord_t>(min),
246 static_cast<coord_t>(max), numBins);
247
248 // Put both in the algo for future use
249 m_bases.emplace_back(basis);
250 m_binDimensions.emplace_back(std::move(out));
251 m_binningScaling.emplace_back(binningScaling);
252 m_transformScaling.emplace_back(transformScaling);
253}
254
255//----------------------------------------------------------------------------------------------
262 // Count the number of output dimensions
263 m_outD = 0;
264 std::string dimChars = this->getDimensionChars();
265 for (char dimChar : dimChars) {
266 std::string propName = "BasisVector0";
267 propName[11] = dimChar;
268 if (!Strings::strip(this->getPropertyValue(propName)).empty())
269 m_outD++;
270 }
271
272 std::vector<double> extents = this->getProperty("OutputExtents");
273 if (extents.size() != m_outD * 2)
274 throw std::invalid_argument("The OutputExtents parameter must have " + Strings::toString(m_outD * 2) +
275 " entries "
276 "(2 for each dimension in the OUTPUT workspace).");
277
278 m_minExtents.clear();
279 m_maxExtents.clear();
280 for (size_t d = 0; d < m_outD; d++) {
281 m_minExtents.emplace_back(extents[d * 2]);
282 m_maxExtents.emplace_back(extents[d * 2 + 1]);
283 }
284
285 m_numBins = this->getProperty("OutputBins");
286 if (m_numBins.size() != m_outD)
287 throw std::invalid_argument("The OutputBins parameter must have 1 entry "
288 "for each dimension in the OUTPUT workspace.");
289
290 m_NormalizeBasisVectors = this->getProperty("NormalizeBasisVectors");
291 m_transformScaling.clear();
292
293 // Create the dimensions based on the strings from the user
294 for (char dimChar : dimChars) {
295 std::string propName = "BasisVector0";
296 propName[11] = dimChar;
297 try {
299 } catch (std::exception &e) {
300 throw std::invalid_argument("Error parsing the " + propName + " parameter: " + std::string(e.what()));
301 }
302 }
303
304 // Number of output binning dimensions found
305 m_outD = m_binDimensions.size();
306 if (m_outD == 0)
307 throw std::runtime_error("No output dimensions were found in the MDEventWorkspace. Cannot bin!");
308
309 // Get the Translation parameter
310 std::vector<double> translVector;
311 try {
312 translVector = getProperty("Translation");
313 } catch (std::exception &e) {
314 throw std::invalid_argument("Error parsing the Translation parameter: " + std::string(e.what()));
315 }
316
317 // Default to 0,0,0 when not specified
318 if (translVector.empty())
319 translVector.resize(m_inWS->getNumDims(), 0);
320 m_translation = VMD(translVector);
321
322 if (m_translation.getNumDims() != m_inWS->getNumDims())
323 throw std::invalid_argument("The number of dimensions in the Translation parameter is "
324 "not consistent with the number of dimensions in the input workspace.");
325
326 // Validate
327 if (m_outD > m_inWS->getNumDims())
328 throw std::runtime_error("More output dimensions were specified than input dimensions "
329 "exist in the MDEventWorkspace. Cannot bin!");
330 if (m_binningScaling.size() != m_outD)
331 throw std::runtime_error("Inconsistent number of entries in scaling vector.");
332}
333
334//----------------------------------------------------------------------------------------------
339 // Process all the input properties
341
342 // Number of input dimensions
343 size_t inD = m_inWS->getNumDims();
344
345 // ----- Make the basis vectors orthogonal -------------------------
346 bool ForceOrthogonal = getProperty("ForceOrthogonal");
347 if (ForceOrthogonal && m_bases[0].getNumDims() == 3 && m_bases.size() >= 2) {
348 std::vector<VMD> firstTwo = m_bases;
349 firstTwo.resize(2, VMD(3));
350 std::vector<VMD> ortho = VMD::makeVectorsOrthogonal(firstTwo);
351 // Set the bases back
352 ortho.resize(m_bases.size(), VMD(3));
353 m_bases = ortho;
354 g_log.information() << "Basis vectors forced to be orthogonal: ";
355 for (auto &base : m_bases)
356 g_log.information() << base.toString(",") << "; ";
357 g_log.information() << '\n';
358 }
359
360 // Now, convert the original vector to the coordinates of the ORIGNAL ws, if
361 // any
362 if (m_originalWS) {
363 CoordTransform const *toOrig = m_inWS->getTransformToOriginal();
365 }
366
367 // OK now find the min/max coordinates of the edges in the INPUT workspace
369 for (size_t d = 0; d < m_outD; d++) {
370 // Translate from the outCoords=(0,0,0) to outCoords=(min,min,min)
371 m_inputMinPoint += (m_bases[d] * m_binDimensions[d]->getMinimum());
372 }
373
374 // Create the CoordTransformAffine for BINNING with these basis vectors
375 auto ct = std::make_unique<DataObjects::CoordTransformAffine>(inD, m_outD);
376 // Note: the scaling makes the coordinate correspond to a bin index
377 ct->buildNonOrthogonal(m_inputMinPoint, this->m_bases, VMD(this->m_binningScaling) / VMD(m_transformScaling));
378 this->m_transform = std::move(ct);
379
380 // Transformation original->binned
381 auto ctFrom = std::make_unique<DataObjects::CoordTransformAffine>(inD, m_outD);
382 ctFrom->buildNonOrthogonal(m_translation, this->m_bases, VMD(m_transformScaling) / VMD(m_transformScaling));
383 m_transformFromOriginal = std::move(ctFrom);
384
385 // Validate
386 if (m_transform->getInD() != inD)
387 throw std::invalid_argument("The number of input dimensions in the CoordinateTransform "
388 "object is not consistent with the number of dimensions in the input "
389 "workspace.");
390 if (m_transform->getOutD() != m_outD)
391 throw std::invalid_argument("The number of output dimensions in the CoordinateTransform "
392 "object is not consistent with the number of dimensions specified in "
393 "the OutDimX, etc. properties.");
394
395 // Now the reverse transformation
396 m_transformToOriginal = nullptr;
397 if (m_outD == inD) {
398 // Can't reverse transform if you lost dimensions.
399 auto ctTo = std::make_unique<DataObjects::CoordTransformAffine>(inD, m_outD);
400 auto toMatrix = static_cast<DataObjects::CoordTransformAffine *>(m_transformFromOriginal.get())->getMatrix();
401 // Invert the affine matrix to get the reverse transformation
402 toMatrix.Invert();
403 ctTo->setMatrix(toMatrix);
404 m_transformToOriginal = std::move(ctTo);
405 }
406}
407
408//----------------------------------------------------------------------------------------------
416 if (str.empty()) {
417 throw std::runtime_error("Empty string passed to one of the AlignedDim0 parameters.");
418 } else {
419 // Strip spaces
420 std::string input = Strings::strip(str);
421 if (input.size() < 4)
422 throw std::invalid_argument("Dimensions string is too short to be valid: " + str);
423
424 // Find the 3rd comma from the end
425 size_t n = std::string::npos;
426 for (size_t i = 0; i < 3; i++) {
427 n = input.find_last_of(',', n);
428 if (n == std::string::npos)
429 throw std::invalid_argument("Wrong number of values (4 are expected) "
430 "in the dimensions string: " +
431 str);
432 if (n == 0)
433 throw std::invalid_argument("Dimension string starts with a comma: " + str);
434 n--;
435 }
436 // The name is everything before the 3rd comma from the end
437 std::string name = input.substr(0, n + 1);
439
440 // And split the rest of it
441 input = input.substr(n + 2);
442 std::vector<std::string> strs;
443 boost::split(strs, input, boost::is_any_of(","));
444 if (strs.size() != 3)
445 throw std::invalid_argument("Wrong number of values (3 are expected) after the name "
446 "in the dimensions string: " +
447 str);
448
449 // Extract the arguments
450 coord_t min, max;
451 int numBins = 0;
452 Strings::convert(strs[0], min);
453 Strings::convert(strs[1], max);
454 Strings::convert(strs[2], numBins);
455 if (name.empty())
456 throw std::invalid_argument("Name should not be blank.");
457 if (min >= max)
458 throw std::invalid_argument("Min should be > max.");
459 if (numBins < 1)
460 throw std::invalid_argument("Number of bins should be >= 1.");
461
462 // Find the named axis in the input workspace
463 size_t dim_index = 0;
464 try {
465 dim_index = m_inWS->getDimensionIndexByName(name);
466 } catch (std::runtime_error &) {
467 // The dimension was not found by name. Try using the ID
468 try {
469 dim_index = m_inWS->getDimensionIndexById(name);
470 } catch (std::runtime_error &) {
471 throw std::runtime_error("Dimension " + name +
472 " was not found in the "
473 "MDEventWorkspace! Cannot continue.");
474 }
475 }
476
477 // Copy the dimension name, ID and units
478 IMDDimension_const_sptr inputDim = m_inWS->getDimension(dim_index);
479 const auto &frame = inputDim->getMDFrame();
481 new MDHistoDimension(inputDim->getName(), inputDim->getDimensionId(), frame, min, max, numBins)));
482
483 // Add the index from which we're binning to the vector
484 this->m_dimensionToBinFrom.emplace_back(dim_index);
485 }
486}
487//----------------------------------------------------------------------------------------------
492 std::string dimChars = this->getDimensionChars();
493
494 // Validate inputs
495 bool previousWasEmpty = false;
496 size_t numDims = 0;
497 for (char dimChar : dimChars) {
498 std::string propName = "AlignedDim0";
499 propName[10] = dimChar;
500 std::string prop = Strings::strip(getPropertyValue(propName));
501 if (!prop.empty())
502 numDims++;
503 if (!prop.empty() && previousWasEmpty)
504 throw std::invalid_argument("Please enter the AlignedDim parameters in the order 0,1,2, etc.,"
505 "without skipping any entries.");
506 previousWasEmpty = prop.empty();
507 }
508
509 // Number of input dimension
510 size_t inD = m_inWS->getNumDims();
511 // Validate
512 if (numDims == 0)
513 throw std::runtime_error("No output dimensions specified.");
514 if (numDims > inD)
515 throw std::runtime_error("More output dimensions were specified than input dimensions "
516 "exist in the MDEventWorkspace.");
517
518 // Create the dimensions based on the strings from the user
519 for (size_t i = 0; i < numDims; i++) {
520 std::string propName = "AlignedDim0";
521 propName[10] = dimChars[i];
523 }
524
525 // Number of output binning dimensions found
526 m_outD = m_binDimensions.size();
527
528 // Now we build the coordinate transformation object
529 m_translation = VMD(inD);
530 m_bases.clear();
531 std::vector<coord_t> origin(m_outD), scaling(m_outD);
532 for (size_t d = 0; d < m_outD; d++) {
533 origin[d] = m_binDimensions[d]->getMinimum();
534 scaling[d] = 1.0f / m_binDimensions[d]->getBinWidth();
535 // Origin in the input
537 // Create a unit basis vector that corresponds to this
538 VMD basis(inD);
539 basis[m_dimensionToBinFrom[d]] = 1.0;
540 m_bases.emplace_back(basis);
541 }
542
543 // Transform for binning
544 m_transform = std::make_unique<DataObjects::CoordTransformAligned>(m_inWS->getNumDims(), m_outD, m_dimensionToBinFrom,
545 origin, scaling);
546
547 // Transformation original->binned. There is no offset or scaling!
548 std::vector<coord_t> unitScaling(m_outD, 1.0);
549 std::vector<coord_t> zeroOrigin(m_outD, 0.0);
551 std::make_unique<DataObjects::CoordTransformAligned>(inD, m_outD, m_dimensionToBinFrom, zeroOrigin, unitScaling);
552
553 // Now the reverse transformation.
554 if (m_outD == inD) {
555 // Make the reverse map = if you're in the output dimension "od", what INPUT
556 // dimension index is that?
557 Matrix<coord_t> mat = m_transformFromOriginal->makeAffineMatrix();
558 mat.Invert();
559 auto tmp = std::make_unique<DataObjects::CoordTransformAffine>(inD, m_outD);
560 tmp->setMatrix(mat);
561 m_transformToOriginal = std::move(tmp);
562 } else {
563 // Changed # of dimensions - can't reverse the transform
564 m_transformToOriginal = nullptr;
565 g_log.warning("SlicingAlgorithm: Your slice will cause the output "
566 "workspace to have fewer dimensions than the input. This will "
567 "affect your ability to create subsequent slices.");
568 }
569}
570
571//-----------------------------------------------------------------------------------------------
584 if (!m_inWS)
585 throw std::runtime_error("SlicingAlgorithm::createTransform(): input "
586 "MDWorkspace must be set first!");
587 if (std::dynamic_pointer_cast<MatrixWorkspace>(m_inWS))
588 throw std::runtime_error(this->name() + " cannot be run on a MatrixWorkspace!");
589
590 // Is the transformation aligned with axes?
591 m_axisAligned = getProperty("AxisAligned");
592
593 // Refer to the original workspace. Make sure that is possible
594 if (m_inWS->numOriginalWorkspaces() > 0)
595 m_originalWS = std::dynamic_pointer_cast<IMDWorkspace>(m_inWS->getOriginalWorkspace());
596 if (m_originalWS) {
597 if (m_axisAligned)
598 throw std::runtime_error("Cannot perform axis-aligned binning on a MDHistoWorkspace. "
599 "Please use non-axis aligned binning.");
600
601 if (m_originalWS->getNumDims() != m_inWS->getNumDims())
602 throw std::runtime_error("SlicingAlgorithm::createTransform(): Cannot propagate "
603 "a transformation if the number of dimensions has changed.");
604
605 if (!m_inWS->getTransformToOriginal())
606 throw std::runtime_error("SlicingAlgorithm::createTransform(): Cannot propagate "
607 "a transformation. There is no transformation saved from " +
608 m_inWS->getName() + " back to " + m_originalWS->getName() + ".");
609
610 // Fail if the MDHistoWorkspace was modified by binary operation
611 DataObjects::MDHistoWorkspace_sptr inHisto = std::dynamic_pointer_cast<DataObjects::MDHistoWorkspace>(m_inWS);
612 if (inHisto) {
613 if (inHisto->getNumExperimentInfo() > 0) {
614 const Run &run = inHisto->getExperimentInfo(0)->run();
615 if (run.hasProperty("mdhisto_was_modified")) {
616 Property *prop = run.getProperty("mdhisto_was_modified");
617 if (prop) {
618 if (prop->value() == "1") {
619 throw std::runtime_error("This MDHistoWorkspace was modified by a binary operation "
620 "(e.g. Plus, Minus). "
621 "It is not currently possible to rebin a modified "
622 "MDHistoWorkspace because that requires returning to the "
623 "original "
624 "(unmodified) MDEventWorkspace, and so would give incorrect "
625 "results. "
626 "Instead, you can use SliceMD and perform operations on the "
627 "resulting "
628 "MDEventWorkspaces, which preserve all events. "
629 "You can override this check by removing the "
630 "'mdhisto_was_modified' sample log.");
631 }
632 }
633 }
634 }
635 }
636
637 g_log.notice() << "Performing " << this->name() << " on the original workspace, '" << m_originalWS->getName()
638 << "'\n";
639 }
640
641 // Create the coordinate transformation
642 m_transform = nullptr;
643 if (m_axisAligned)
645 else
647
648 // Finalize, for binnign MDHistoWorkspace
649 if (m_originalWS) {
650 // The intermediate workspace is the MDHistoWorkspace being BINNED
652 CoordTransform const *originalToIntermediate = m_intermediateWS->getTransformFromOriginal();
653 if (originalToIntermediate && (m_originalWS->getNumDims() == m_intermediateWS->getNumDims())) {
654 try {
655 // The transform from the INPUT to the INTERMEDIATE ws
656 // intermediate_coords = [OriginalToIntermediate] * [thisToOriginal] *
657 // these_coords
658 Matrix<coord_t> matToOriginal = m_transformToOriginal->makeAffineMatrix();
659 Matrix<coord_t> matOriginalToIntermediate = originalToIntermediate->makeAffineMatrix();
660 Matrix<coord_t> matToIntermediate = matOriginalToIntermediate * matToOriginal;
661
662 m_transformToIntermediate = std::make_unique<DataObjects::CoordTransformAffine>(m_originalWS->getNumDims(),
663 m_intermediateWS->getNumDims());
664 m_transformToIntermediate->setMatrix(matToIntermediate);
665 // And now the reverse
666 matToIntermediate.Invert();
667 m_transformFromIntermediate = std::make_unique<DataObjects::CoordTransformAffine>(
668 m_intermediateWS->getNumDims(), m_originalWS->getNumDims());
669 m_transformFromIntermediate->setMatrix(matToIntermediate);
670 } catch (std::runtime_error &) {
671 // Ignore error. Have no transform.
672 }
673 }
674
675 // Replace the input workspace with the original MDEventWorkspace
676 // for future binning
678 }
679}
680
681//----------------------------------------------------------------------------------------------
698std::unique_ptr<MDImplicitFunction> SlicingAlgorithm::getGeneralImplicitFunction(const size_t *const chunkMin,
699 const size_t *const chunkMax) {
700 size_t nd = m_inWS->getNumDims();
701
702 // General implicit function
703 auto func = std::make_unique<MDImplicitFunction>();
704
705 // First origin = min of each basis vector
706 VMD o1 = m_translation;
707 // Second origin = max of each basis vector
708 VMD o2 = m_translation;
709 // And this the list of basis vectors. Each vertex is given by o1+bases[i].
710 std::vector<VMD> bases;
711 VMD x;
712
713 for (size_t d = 0; d < m_bases.size(); d++) {
714 double xMin = m_binDimensions[d]->getMinimum();
715 double xMax = m_binDimensions[d]->getMaximum();
716 // Move the position if you're using a chunk
717 if (chunkMin != nullptr)
718 xMin = m_binDimensions[d]->getX(chunkMin[d]);
719 if (chunkMax != nullptr)
720 xMax = m_binDimensions[d]->getX(chunkMax[d]);
721 // Offset the origin by the position along the basis vector
722 o1 += (m_bases[d] * xMin);
723 o2 += (m_bases[d] * xMax);
724
725 VMD thisBase = m_bases[d] * (xMax - xMin);
726 bases.emplace_back(thisBase);
727 if (d == 0)
728 x = thisBase;
729 }
730
731 // Dimensionality of the box
732 size_t boxDim = bases.size();
733
734 // Point that is sure to be inside the volume of interest
735 VMD insidePoint = (o1 + o2) / 2.0;
736
737 if (boxDim == 1) {
738 // 2 planes defined by 1 basis vector
739 // Your normal = the x vector
740 func->addPlane(MDPlane(x, o1));
741 func->addPlane(MDPlane(x * -1.0, o2));
742 } else if (boxDim == nd || boxDim == nd - 1) {
743 // Create a pair of planes for each base supplied. This is general to non-
744 // orthogonal bases. If we have bases (x y z t) then we create the planes
745 //
746 // y z t
747 // x z t
748 // x y t
749 // x y z
750 //
751 // Note: the last plane may or may not be created depending on the number
752 // of bases supplied to the slicing algorithm relative to the number
753 // of dimensions. i.e. if 3 bases were supplied and we have 4 dimensions
754 // then 6 planes are created instead of 8.
755 std::vector<VMD> vectors;
756
757 for (size_t ignoreIndex = 0; ignoreIndex < boxDim; ++ignoreIndex) {
758 vectors.clear();
759 // Create a list of vectors that excludes the "current" basis
760 for (size_t baseIndex = 0; baseIndex < boxDim; ++baseIndex) {
761 if (baseIndex != ignoreIndex)
762 vectors.emplace_back(bases[baseIndex]);
763 }
764
765 // if we have fewer basis vectors than dimensions
766 // create a normal for the final dimension
767 if (boxDim == nd - 1)
768 vectors.emplace_back(VMD::getNormalVector(bases));
769
770 // Add two planes for each set of vectors
771 func->addPlane(MDPlane(vectors, o1, insidePoint));
772 func->addPlane(MDPlane(vectors, o2, insidePoint));
773 }
774 } else {
775 // Last-resort, totally general case
776 // 2*N planes defined by N basis vectors, in any dimensionality workspace.
777 // Assumes orthogonality!
778 g_log.warning("SlicingAlgorithm given " + std::to_string(boxDim) + " bases and " + std::to_string(nd) +
779 " dimensions and " + "therefore will assume orthogonality");
780 for (auto &base : bases) {
781 // For each basis vector, make two planes, perpendicular to it and facing
782 // inwards
783 func->addPlane(MDPlane(base, o1));
784 func->addPlane(MDPlane(base * -1.0, o2));
785 }
786 }
787
788 return func;
789}
790
791//----------------------------------------------------------------------------------------------
804std::unique_ptr<MDImplicitFunction> SlicingAlgorithm::getImplicitFunctionForChunk(const size_t *const chunkMin,
805 const size_t *const chunkMax) {
806 size_t nd = m_inWS->getNumDims();
807 if (m_axisAligned) {
808 std::vector<coord_t> function_min(nd, -1e30f); // default to all space if the dimension is not specified
809 std::vector<coord_t> function_max(nd, +1e30f); // default to all space if the dimension is not specified
810 for (size_t bd = 0; bd < m_outD; bd++) {
811 // Dimension in the MDEventWorkspace
812 size_t d = m_dimensionToBinFrom[bd];
813 if (chunkMin)
814 function_min[d] = m_binDimensions[bd]->getX(chunkMin[bd]);
815 else
816 function_min[d] = m_binDimensions[bd]->getX(0);
817 if (chunkMax)
818 function_max[d] = m_binDimensions[bd]->getX(chunkMax[bd]);
819 else
820 function_max[d] = m_binDimensions[bd]->getX(m_binDimensions[bd]->getNBins());
821 }
822 return std::make_unique<MDBoxImplicitFunction>(function_min, function_max);
823 } else {
824 // General implicit function
825 return getGeneralImplicitFunction(chunkMin, chunkMax);
826 }
827}
828
839 const Mantid::Kernel::VMD &basisVector) const {
840 // Get set of basis vectors
841 auto oldBasis = getOldBasis(m_inWS->getNumDims());
842
843 // Get indices onto which the vector projects
844 auto indicesWithProjection = getIndicesWithProjection(basisVector, oldBasis);
845
846 // Extract MDFrame
847 return extractMDFrameForNonAxisAligned(indicesWithProjection, units);
848}
849
850std::vector<Mantid::Kernel::VMD> SlicingAlgorithm::getOldBasis(size_t dimension) const {
851 std::vector<Mantid::Kernel::VMD> oldBasis;
852 for (size_t i = 0; i < dimension; ++i) {
853 Mantid::Kernel::VMD basisVector(dimension);
854 basisVector[i] = 1.0;
855 oldBasis.emplace_back(basisVector);
856 }
857 return oldBasis;
858}
859
867 const Mantid::Kernel::VMD &basisVector) const {
868 return std::fabs(oldVector.scalar_prod(basisVector)) > 0.0;
869}
870
878 const std::vector<Mantid::Kernel::VMD> &oldBasis) const {
879 std::vector<size_t> indexWithProjection;
880 for (size_t index = 0; index < oldBasis.size(); ++index) {
881 if (isProjectingOnFrame(oldBasis[index], basisVector)) {
882 indexWithProjection.emplace_back(index);
883 }
884 }
885 return indexWithProjection;
886}
887
896SlicingAlgorithm::extractMDFrameForNonAxisAligned(std::vector<size_t> indicesWithProjection,
897 const std::string &units) const {
898 if (indicesWithProjection.empty()) {
899 g_log.warning() << "Slicing Algorithm: Chosen vector does not "
900 "project on any vector of the old basis.";
901 }
902 // Get a reference frame to perform pairwise comparison
903 const auto &referenceMDFrame = m_inWS->getDimension(indicesWithProjection[0])->getMDFrame();
904
905 for (auto &index : indicesWithProjection) {
906 const auto &toCheckMDFrame = m_inWS->getDimension(index)->getMDFrame();
907 if (!referenceMDFrame.isSameType(toCheckMDFrame)) {
908 g_log.warning() << "Slicing Algorithm: New basis vector tries to "
909 "mix un-mixable MDFrame types.";
910 }
911 }
912
913 Mantid::Geometry::MDFrame_uptr mdFrame(referenceMDFrame.clone());
914 setTargetUnits(mdFrame, units);
915
916 return mdFrame;
917}
918
919/*
920 * Set units of the output workspace
921 * @param mdFrame: MDFrame to be added to the output workspace
922 * @param unit: the unit to use in mdFrame
923 */
924void SlicingAlgorithm::setTargetUnits(Mantid::Geometry::MDFrame_uptr &mdFrame, const std::string &unit) const {
925 boost::regex pattern("in.*A.*\\^-1");
926
927 if (boost::regex_match(unit, pattern)) {
928 // RLU with special label
929 auto md_unit = ReciprocalLatticeUnit(unit);
930 mdFrame->setMDUnit(md_unit);
931 } else if (unit == "r") {
932 // RLU
933 auto md_unit = ReciprocalLatticeUnit();
934 mdFrame->setMDUnit(md_unit);
935 } else if (unit == "a") {
936 // Inverse angstroms
937 auto md_unit = InverseAngstromsUnit();
938 mdFrame->setMDUnit(md_unit);
939 }
940 // else leave the unit the same as the input workspace
941}
942
943} // namespace Mantid::MDAlgorithms
gsl_vector * tmp
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
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
const std::string name() const override=0
function to return a name of the algorithm, must be overridden in all algorithms
Unique SingleValueParameter Declaration for InputNDimensions.
Mantid::Kernel::VMD applyVMD(const Mantid::Kernel::VMD &inputVector) const
Wrapper for VMD.
virtual Mantid::Kernel::Matrix< coord_t > makeAffineMatrix() const
bool hasProperty(const std::string &name) const
Does the property exist on the object.
Definition: LogManager.cpp:265
Kernel::Property * getProperty(const std::string &name) const
Returns the named property as a pointer.
Definition: LogManager.cpp:404
This class stores information regarding an experimental run as a series of log entries.
Definition: Run.h:38
Generic class to transform from M input dimensions to N output dimensions.
void setMatrix(const Mantid::Kernel::Matrix< coord_t > &newMatrix)
Directly set the affine matrix to use.
A generalized description of a N-dimensional hyperplane.
Definition: MDPlane.h:41
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
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.
Inverse Angstroms unit.
Definition: MDUnit.h:51
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
Numerical Matrix class.
Definition: Matrix.h:42
T Invert()
LU inversion routine.
Definition: Matrix.cpp:924
The concrete, templated class for properties.
Base class for properties.
Definition: Property.h:94
virtual std::string value() const =0
Returns the value of the property as a string.
static std::vector< VMDBase > makeVectorsOrthogonal(std::vector< VMDBase > &vectors)
Make an orthogonal system with 2 input 3D vectors.
Definition: VMD.cpp:450
TYPE normalize()
Normalize this vector to unity length.
Definition: VMD.cpp:427
TYPE scalar_prod(const VMDBase &v) const
Scalar product of two vectors.
Definition: VMD.cpp:391
size_t size() const
Definition: VMD.cpp:239
size_t getNumDims() const
Definition: VMD.cpp:236
static VMDBase getNormalVector(const std::vector< VMDBase > &vectors)
Given N-1 vectors defining a N-1 dimensional hyperplane in N dimensions, returns a vector that is nor...
Definition: VMD.cpp:494
TYPE norm() const
Definition: VMD.cpp:420
std::unique_ptr< API::CoordTransform > m_transformFromOriginal
Coordinate transformation to save in the output workspace (original->binned)
Mantid::API::IMDWorkspace_sptr m_intermediateWS
Intermediate original workspace.
std::vector< size_t > m_dimensionToBinFrom
Index of the dimension in the MDEW for the dimension in the output.
std::unique_ptr< Mantid::Geometry::MDImplicitFunction > getImplicitFunctionForChunk(const size_t *const chunkMin, const size_t *const chunkMax)
Create an implicit function for picking boxes, based on the indexes in the output MDHistoWorkspace.
bool isProjectingOnFrame(const Mantid::Kernel::VMD &oldVector, const Mantid::Kernel::VMD &basisVector) const
Check if the two vectors are orthogonal or not.
std::unique_ptr< API::CoordTransform > m_transform
Coordinate transformation to apply.
std::vector< size_t > getIndicesWithProjection(const Mantid::Kernel::VMD &basisVector, const std::vector< Mantid::Kernel::VMD > &oldBasis) const
Get indices which have a projection contribution.
Mantid::API::IMDWorkspace_sptr m_originalWS
Original (MDEventWorkspace) that inWS was based on.
void createGeneralTransform()
Loads the dimensions and create the coordinate transform, using the inputs.
void createTransform()
Read the algorithm properties and creates the appropriate transforms for slicing the MDEventWorkspace...
void createAlignedTransform()
Using the parameters, create a coordinate transformation for aligned cuts.
void setTargetUnits(Mantid::Geometry::MDFrame_uptr &frame, const std::string &units) const
std::vector< Mantid::Kernel::VMD > getOldBasis(size_t dimension) const
Mantid::Kernel::VMD m_inputMinPoint
Coordinates in the INPUT workspace corresponding to the minimum edge in all dimensions.
std::vector< double > m_maxExtents
For non-aligned, the maximum coordinate extents in each OUTPUT dimension.
Mantid::API::IMDWorkspace_sptr m_inWS
Input workspace.
std::unique_ptr< API::CoordTransform > m_transformToOriginal
Coordinate transformation to save in the output workspace (binned->original)
std::unique_ptr< DataObjects::CoordTransformAffine > m_transformFromIntermediate
Coordinate transformation to save in the output WS, from the intermediate WS.
std::vector< double > m_binningScaling
Scaling factor to apply for each basis vector (to map to the bins).
size_t m_outD
Number of dimensions in the output (binned) workspace.
void makeBasisVectorFromString(const std::string &str)
Generate the MDHistoDimension and basis vector for a given string from BasisVector0 etc.
Mantid::Geometry::MDFrame_uptr createMDFrameForNonAxisAligned(const std::string &units, const Mantid::Kernel::VMD &basisVector) const
Create an MDFrame for the Non-Axis-Aligned case.
bool m_NormalizeBasisVectors
The NormalizeBasisVectors option.
void processGeneralTransformProperties()
Reads the various Properties for the general (non-aligned) case and fills in members on the Algorithm...
std::vector< Mantid::Kernel::VMD > m_bases
Basis vectors of the output dimensions, normalized to unity length.
std::unique_ptr< Mantid::Geometry::MDImplicitFunction > getGeneralImplicitFunction(const size_t *const chunkMin, const size_t *const chunkMax)
Create an implicit function for picking boxes, based on the indexes in the output MDHistoWorkspace.
Mantid::Geometry::MDFrame_uptr extractMDFrameForNonAxisAligned(std::vector< size_t > indicesWithProjection, const std::string &units) const
Extract the MDFrame.
bool m_axisAligned
Set to true if the cut is aligned with the axes.
void initSlicingProps()
Initialise the properties.
std::vector< int > m_numBins
For non-aligned, the number of bins in each OUTPUT dimension.
Mantid::Kernel::VMD m_translation
Translation from the OUTPUT to the INPUT workspace i.e.
std::vector< double > m_transformScaling
Scaling factor to apply for each basis vector to transfor to the output dimensions.
std::vector< double > m_minExtents
For non-aligned, the minimum coordinate extents in each OUTPUT dimension.
std::unique_ptr< DataObjects::CoordTransformAffine > m_transformToIntermediate
Coordinate transformation to save in the intermediate WS.
void makeAlignedDimensionFromString(const std::string &str)
Generate a MDHistoDimension_sptr from a comma-sep string (for AlignedDim0, etc.) Must be called in or...
std::vector< Mantid::Geometry::MDHistoDimension_sptr > m_binDimensions
Bin dimensions to actually use.
std::shared_ptr< MDHistoWorkspace > MDHistoWorkspace_sptr
A shared pointer to a MDHistoWorkspace.
std::unique_ptr< MDFrame > MDFrame_uptr
Definition: MDFrame.h:36
std::shared_ptr< MDHistoDimension > MDHistoDimension_sptr
Shared pointer to a MDHistoDimension.
std::shared_ptr< const IMDDimension > IMDDimension_const_sptr
Shared Pointer to const IMDDimension.
Definition: IMDDimension.h:101
MANTID_KERNEL_DLL std::string strip(const std::string &A)
strip pre/post spaces
Definition: Strings.cpp:397
int convert(const std::string &A, T &out)
Convert a string into a number.
Definition: Strings.cpp:665
std::string toString(const T &value)
Convert a number to a string.
Definition: Strings.cpp:703
VMDBase< VMD_t > VMD
Define the VMD as using the double or float data type.
Definition: VMD.h:86
float coord_t
Typedef for the data type to use for coordinate axes in MD objects such as MDBox, MDEventWorkspace,...
Definition: MDTypes.h:27
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
Definition: Property.h:53