#include "MantidGeometry/MDGeometry/MDBoxImplicitFunction.h"
#include "MantidGeometry/MDGeometry/MDHistoDimension.h"
#include "MantidGeometry/MDGeometry/MDPlane.h"
#include "MantidKernel/EnabledWhenProperty.h"
#include "MantidKernel/Strings.h"
#include "MantidKernel/System.h"
#include "MantidDataObjects/CoordTransformAffine.h"
#include "MantidDataObjects/CoordTransformAligned.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/VisibleWhenProperty.h"
#include "MantidKernel/ArrayLengthValidator.h"
#include "MantidMDAlgorithms/SlicingAlgorithm.h"
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;
using Mantid::Kernel::Strings::strip;
namespace Mantid {
namespace MDAlgorithms {
//----------------------------------------------------------------------------------------------
/** Constructor
*/
SlicingAlgorithm::SlicingAlgorithm()
: m_transform(), m_transformFromOriginal(), m_transformToOriginal(),
m_transformFromIntermediate(), m_transformToIntermediate(),
m_axisAligned(true), m_outD(0), // unititialized and should be invalid
m_NormalizeBasisVectors(false) {}
//----------------------------------------------------------------------------------------------
/** Destructor
*/
SlicingAlgorithm::~SlicingAlgorithm() {}
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void SlicingAlgorithm::initSlicingProps() {
std::string dimChars = getDimensionChars();
// --------------- Axis-aligned properties
// ---------------------------------------
declareProperty(
"AxisAligned", true,
"Perform binning aligned with the axes of the input MDEventWorkspace?");
setPropertyGroup("AxisAligned", "Axis-Aligned Binning");
for (size_t i = 0; i < dimChars.size(); i++) {
std::string dim(" ");
dim[0] = dimChars[i];
std::string propName = "AlignedDim" + dim;
declareProperty(
new PropertyWithValue<std::string>(propName, "", Direction::Input),
"Binning parameters for the " + Strings::toString(i) +
"th dimension.\n"
"Enter it as a comma-separated list of values with the format: "
"'name,minimum,maximum,number_of_bins'. Leave blank for NONE.");
setPropertySettings(
propName, new VisibleWhenProperty("AxisAligned", IS_EQUAL_TO, "1"));
setPropertyGroup(propName, "Axis-Aligned Binning");
}
// --------------- NON-Axis-aligned properties
// ---------------------------------------
std::string grpName = "Non-Aligned Binning";
IPropertySettings *ps =
new VisibleWhenProperty("AxisAligned", IS_EQUAL_TO, "0");
for (size_t i = 0; i < dimChars.size(); i++) {
std::string dim(" ");
dim[0] = dimChars[i];
std::string propName = "BasisVector" + dim;
declareProperty(
new PropertyWithValue<std::string>(propName, "", Direction::Input),
"Description of the basis vector of the " + Strings::toString(i) +
"th output dimension."
"Format: 'name, units, x,y,z,..'.\n"
" name : string for the name of the output dimension.\n"
" units : string for the units of the output dimension.\n"
" x,y,z,...: vector definining the basis in the input dimensions "
"space.\n"
"Leave blank for NONE.");
setPropertySettings(propName, ps->clone());
setPropertyGroup(propName, grpName);
}
declareProperty(new ArrayProperty<double>("Translation", Direction::Input),
"Coordinates in the INPUT workspace that corresponds to "
"(0,0,0) in the OUTPUT workspace.\n"
"Enter as a comma-separated string.\n"
"Default: 0 in all dimensions (no translation).");
declareProperty(new ArrayProperty<double>("OutputExtents", Direction::Input),
"The minimum, maximum edges of space of each dimension of "
"the OUTPUT workspace, as a comma-separated list");
declareProperty(
new ArrayProperty<int>("OutputBins", Direction::Input),
"The number of bins for each dimension of the OUTPUT workspace.");
declareProperty(new PropertyWithValue<bool>("NormalizeBasisVectors", true,
Direction::Input),
"Normalize the given basis vectors to unity. \n"
"If true, then a distance of 1 in the INPUT dimensions = 1 "
"in the OUTPUT dimensions.\n"
"If false, then a distance of norm(basis_vector) in the "
"INPUT dimension = 1 in the OUTPUT dimensions.");
declareProperty(
new PropertyWithValue<bool>("ForceOrthogonal", false, Direction::Input),
"Force the input basis vectors to form an orthogonal coordinate system. "
"Only works in 3 dimension!");
// For GUI niceness
setPropertyGroup("Translation", grpName);
setPropertyGroup("OutputExtents", grpName);
setPropertyGroup("OutputBins", grpName);
setPropertyGroup("NormalizeBasisVectors", grpName);
setPropertyGroup("ForceOrthogonal", grpName);
setPropertySettings("Translation", ps->clone());
setPropertySettings("OutputExtents", ps->clone());
setPropertySettings("OutputBins", ps->clone());
setPropertySettings("NormalizeBasisVectors", ps->clone());
setPropertySettings("ForceOrthogonal", ps);
}
//----------------------------------------------------------------------------------------------
/** Generate the MDHistoDimension and basis vector for a given string from
*BasisVector0 etc.
*
* If the workspace being binned has an original workspace, then the vector
* is transformed to THOSE coordinates.
*
* "Format: 'name, units, x,y,z,..'.\n"
* Adds values to m_bases, m_binDimensions,
* m_binningScaling and m_transformScaling
*
* @param str :: name,number_of_bins
*/
void SlicingAlgorithm::makeBasisVectorFromString(const std::string &str) {
std::string input = Strings::strip(str);
if (input.empty())
return;
if (input.size() < 3)
throw std::invalid_argument("Dimension string is too short to be valid: " +
str);
// The current dimension index.
size_t dim = m_binDimensions.size();
size_t n_first_comma = std::string::npos;
// Special case: accept dimension names [x,y,z]
if (input[0] == '[') {
// Find the name at the closing []
size_t n = input.find_first_of("]", 1);
if (n == std::string::npos)
throw std::invalid_argument(
"No closing ] character in the dimension name of : " + str);
// Find the comma after the name
n_first_comma = input.find_first_of(",", n);
if (n_first_comma == std::string::npos)
throw std::invalid_argument(
"No comma after the closing ] character in the dimension string: " +
str);
} else
// Find the comma after the name
n_first_comma = input.find_first_of(",");
if (n_first_comma == std::string::npos)
throw std::invalid_argument("No comma in the dimension string: " + str);
if (n_first_comma == input.size() - 1)
throw std::invalid_argument("Dimension string ends in a comma: " + str);
// Get the entire name
std::string name = Strings::strip(input.substr(0, n_first_comma));
if (name.size() == 0)
throw std::invalid_argument("name should not be blank.");
// Now remove the name and comma
// And split the rest of it
input = input.substr(n_first_comma + 1);
std::vector<std::string> strs;
boost::split(strs, input, boost::is_any_of(","));
if (strs.size() != this->m_inWS->getNumDims() + 1)
throw std::invalid_argument("Wrong number of values (expected 2 + # of "
"input dimensions) in the dimensions string: " +
str);
// Get the number of bins from
int numBins = m_numBins[dim];
if (numBins < 1)
throw std::invalid_argument("Number of bins for output dimension " +
Strings::toString(dim) + " should be >= 1.");
// Get the min/max extents in this OUTPUT dimension
double min = m_minExtents[dim];
double max = m_maxExtents[dim];
double lengthInOutput = max - min;
if (lengthInOutput <= 0)
throw std::invalid_argument("The maximum extents for dimension " +
Strings::toString(dim) + " should be > 0.");
// Create the basis vector with the right # of dimensions
VMD basis(this->m_inWS->getNumDims());
for (size_t d = 0; d < this->m_inWS->getNumDims(); d++)
if (!Strings::convert(strs[d + 1], basis[d]))
throw std::invalid_argument("Error converting argument '" + strs[d + 1] +
"' in the dimensions string '" + str +
"' to a number.");
// If B was binned from A (m_originalWS), and we are binning C from B,
// convert the basis vector from B space -> A space
if (m_originalWS) {
// Turn basis vector into two points
VMD basis0(this->m_inWS->getNumDims());
VMD basis1 = basis;
// Convert the points to the original coordinates (from inWS to originalWS)
CoordTransform const *toOrig = m_inWS->getTransformToOriginal();
VMD origBasis0 = toOrig->applyVMD(basis0);
VMD origBasis1 = toOrig->applyVMD(basis1);
// New basis vector, now in the original workspace
basis = origBasis1 - origBasis0;
}
// Check on the length of the basis vector
double basisLength = basis.norm();
if (basisLength <= 0)
throw std::invalid_argument("direction should not be 0-length.");
// Normalize it to unity, if desized
double transformScaling = 1.0;
if (m_NormalizeBasisVectors) {
// A distance of 1 in the INPUT space = a distance of 1.0 in the OUTPUT
// space
basis.normalize();
transformScaling = 1.0;
} else
// A distance of |basisVector| in the INPUT space = a distance of 1.0 in the
// OUTPUT space
transformScaling = (1.0 / basisLength);
// This is the length of this dimension as measured in the INPUT space
double lengthInInput = lengthInOutput / transformScaling;
// Scaling factor, to convert from units in the INPUT dimensions to the output
// BIN number
double binningScaling = double(numBins) / (lengthInInput);
// Extract the arguments
std::string id = name;
std::string units = Strings::strip(strs[0]);
// Create the appropriate frame
//auto frame = createMDFrame(units);
// Create the output dimension
MDHistoDimension_sptr out(
new MDHistoDimension(name, id, units, static_cast<coord_t>(min),
static_cast<coord_t>(max), numBins));
// Put both in the algo for future use
m_bases.push_back(basis);
m_binDimensions.push_back(out);
m_binningScaling.push_back(binningScaling);
m_transformScaling.push_back(transformScaling);
}
//----------------------------------------------------------------------------------------------
/** Reads the various Properties for the general (non-aligned) case
* and fills in members on the Algorithm for later use
*
* @throw if some of the inputs are invalid
*/
void SlicingAlgorithm::processGeneralTransformProperties() {
// Count the number of output dimensions
m_outD = 0;
std::string dimChars = this->getDimensionChars();
for (size_t i = 0; i < dimChars.size(); i++) {
std::string propName = "BasisVector0";
propName[11] = dimChars[i];
if (!Strings::strip(this->getPropertyValue(propName)).empty())
m_outD++;
}
std::vector<double> extents = this->getProperty("OutputExtents");
if (extents.size() != m_outD * 2)
throw std::invalid_argument(
"The OutputExtents parameter must have " +
Strings::toString(m_outD * 2) +
" entries "
"(2 for each dimension in the OUTPUT workspace).");
m_minExtents.clear();
m_maxExtents.clear();
for (size_t d = 0; d < m_outD; d++) {
m_minExtents.push_back(extents[d * 2]);
m_maxExtents.push_back(extents[d * 2 + 1]);
}
m_numBins = this->getProperty("OutputBins");
if (m_numBins.size() != m_outD)
throw std::invalid_argument("The OutputBins parameter must have 1 entry "
"for each dimension in the OUTPUT workspace.");
m_NormalizeBasisVectors = this->getProperty("NormalizeBasisVectors");
m_transformScaling.clear();
// Create the dimensions based on the strings from the user
for (size_t i = 0; i < dimChars.size(); i++) {
std::string propName = "BasisVector0";
propName[11] = dimChars[i];
try {
makeBasisVectorFromString(getPropertyValue(propName));
} catch (std::exception &e) {
throw std::invalid_argument("Error parsing the " + propName +
" parameter: " + std::string(e.what()));
}
}
// Number of output binning dimensions found
m_outD = m_binDimensions.size();
if (m_outD == 0)
throw std::runtime_error(
"No output dimensions were found in the MDEventWorkspace. Cannot bin!");
// Get the Translation parameter
std::vector<double> translVector;
try {
translVector = getProperty("Translation");
} catch (std::exception &e) {
throw std::invalid_argument("Error parsing the Translation parameter: " +
std::string(e.what()));
}
// Default to 0,0,0 when not specified
if (translVector.empty())
translVector.resize(m_inWS->getNumDims(), 0);
m_translation = VMD(translVector);
if (m_translation.getNumDims() != m_inWS->getNumDims())
throw std::invalid_argument(
"The number of dimensions in the Translation parameter is "
"not consistent with the number of dimensions in the input workspace.");
// Validate
if (m_outD > m_inWS->getNumDims())
throw std::runtime_error(
"More output dimensions were specified than input dimensions "
"exist in the MDEventWorkspace. Cannot bin!");
if (m_binningScaling.size() != m_outD)
throw std::runtime_error(
"Inconsistent number of entries in scaling vector.");
}
//----------------------------------------------------------------------------------------------
/** Loads the dimensions and create the coordinate transform, using the inputs.
* This is for the general (i.e. non-aligned) case
*/
void SlicingAlgorithm::createGeneralTransform() {
// Process all the input properties
this->processGeneralTransformProperties();
// Number of input dimensions
size_t inD = m_inWS->getNumDims();
// ----- Make the basis vectors orthogonal -------------------------
bool ForceOrthogonal = getProperty("ForceOrthogonal");
if (ForceOrthogonal && m_bases[0].getNumDims() == 3 && m_bases.size() >= 2) {
std::vector<VMD> firstTwo = m_bases;
firstTwo.resize(2, VMD(3));
std::vector<VMD> ortho = VMD::makeVectorsOrthogonal(firstTwo);
// Set the bases back
ortho.resize(m_bases.size(), VMD(3));
m_bases = ortho;
g_log.information() << "Basis vectors forced to be orthogonal: ";
for (size_t i = 0; i < m_bases.size(); i++)
g_log.information() << m_bases[i].toString(",") << "; ";
g_log.information() << std::endl;
}
// Now, convert the original vector to the coordinates of the ORIGNAL ws, if
// any
if (m_originalWS) {
CoordTransform const *toOrig = m_inWS->getTransformToOriginal();
m_translation = toOrig->applyVMD(m_translation);
}
// OK now find the min/max coordinates of the edges in the INPUT workspace
m_inputMinPoint = m_translation;
for (size_t d = 0; d < m_outD; d++)
// Translate from the outCoords=(0,0,0) to outCoords=(min,min,min)
m_inputMinPoint += (m_bases[d] * m_binDimensions[d]->getMinimum());
// std::cout << m_inputMinPoint << " m_inputMinPoint " << std::endl;
// Create the CoordTransformAffine for BINNING with these basis vectors
DataObjects::CoordTransformAffine *ct =
new DataObjects::CoordTransformAffine(inD, m_outD);
// Note: the scaling makes the coordinate correspond to a bin index
// ct->buildOrthogonal(m_inputMinPoint, this->m_bases,
// VMD(this->m_binningScaling));
ct->buildNonOrthogonal(m_inputMinPoint, this->m_bases,
VMD(this->m_binningScaling) / VMD(m_transformScaling));
this->m_transform = ct;
// Transformation original->binned
DataObjects::CoordTransformAffine *ctFrom =
new DataObjects::CoordTransformAffine(inD, m_outD);
// ctFrom->buildOrthogonal(m_translation, this->m_bases,
// VMD(m_transformScaling));
ctFrom->buildNonOrthogonal(m_translation, this->m_bases,
VMD(m_transformScaling) / VMD(m_transformScaling));
m_transformFromOriginal = ctFrom;
// Validate
if (m_transform->getInD() != inD)
throw std::invalid_argument(
"The number of input dimensions in the CoordinateTransform "
"object is not consistent with the number of dimensions in the input "
"workspace.");
if (m_transform->getOutD() != m_outD)
throw std::invalid_argument(
"The number of output dimensions in the CoordinateTransform "
"object is not consistent with the number of dimensions specified in "
"the OutDimX, etc. properties.");
// Now the reverse transformation
m_transformToOriginal = NULL;
if (m_outD == inD) {
// Can't reverse transform if you lost dimensions.
DataObjects::CoordTransformAffine *ctTo =
new DataObjects::CoordTransformAffine(inD, m_outD);
Matrix<coord_t> fromMatrix = ctFrom->getMatrix();
Matrix<coord_t> toMatrix = fromMatrix;
// Invert the affine matrix to get the reverse transformation
toMatrix.Invert();
ctTo->setMatrix(toMatrix);
m_transformToOriginal = ctTo;
}
}
//----------------------------------------------------------------------------------------------
/** Generate a MDHistoDimension_sptr from a comma-sep string (for AlignedDim0,
*etc.)
* Must be called in order X,Y,Z,T.
*
* @param str :: name,minimum,maximum,number_of_bins
*/
void SlicingAlgorithm::makeAlignedDimensionFromString(const std::string &str) {
if (str.empty()) {
throw std::runtime_error(
"Empty string passed to one of the AlignedDim0 parameters.");
} else {
// Strip spaces
std::string input = Strings::strip(str);
if (input.size() < 4)
throw std::invalid_argument(
"Dimensions string is too short to be valid: " + str);
// Find the 3rd comma from the end
size_t n = std::string::npos;
for (size_t i = 0; i < 3; i++) {
n = input.find_last_of(",", n);
if (n == std::string::npos)
throw std::invalid_argument("Wrong number of values (4 are expected) "
"in the dimensions string: " +
str);
if (n == 0)
throw std::invalid_argument("Dimension string starts with a comma: " +
str);
n--;
}
// The name is everything before the 3rd comma from the end
std::string name = input.substr(0, n + 1);
name = Strings::strip(name);
// And split the rest of it
input = input.substr(n + 2);
std::vector<std::string> strs;
boost::split(strs, input, boost::is_any_of(","));
if (strs.size() != 3)
throw std::invalid_argument(
"Wrong number of values (3 are expected) after the name "
"in the dimensions string: " +
str);
// Extract the arguments
coord_t min, max;
int numBins = 0;
Strings::convert(strs[0], min);
Strings::convert(strs[1], max);
Strings::convert(strs[2], numBins);
if (name.size() == 0)
throw std::invalid_argument("Name should not be blank.");
if (min >= max)
throw std::invalid_argument("Min should be > max.");
if (numBins < 1)
throw std::invalid_argument("Number of bins should be >= 1.");
// Find the named axis in the input workspace
size_t dim_index = 0;
try {
dim_index = m_inWS->getDimensionIndexByName(name);
} catch (std::runtime_error &) {
// The dimension was not found by name. Try using the ID
try {
dim_index = m_inWS->getDimensionIndexById(name);
} catch (std::runtime_error &) {
throw std::runtime_error("Dimension " + name +
" was not found in the "
"MDEventWorkspace! Cannot continue.");
}
}
// Copy the dimension name, ID and units
IMDDimension_const_sptr inputDim = m_inWS->getDimension(dim_index);
const auto& frame = inputDim->getMDFrame();
m_binDimensions.push_back(MDHistoDimension_sptr(
new MDHistoDimension(inputDim->getName(), inputDim->getDimensionId(),
frame, min, max, numBins)));
// Add the index from which we're binning to the vector
this->m_dimensionToBinFrom.push_back(dim_index);
}
}
//----------------------------------------------------------------------------------------------
/** Using the parameters, create a coordinate transformation
* for aligned cuts
*/
void SlicingAlgorithm::createAlignedTransform() {
std::string dimChars = this->getDimensionChars();
// Validate inputs
bool previousWasEmpty = false;
size_t numDims = 0;
for (size_t i = 0; i < dimChars.size(); i++) {
std::string propName = "AlignedDim0";
propName[10] = dimChars[i];
std::string prop = Strings::strip(getPropertyValue(propName));
if (!prop.empty())
numDims++;
if (!prop.empty() && previousWasEmpty)
throw std::invalid_argument(
"Please enter the AlignedDim parameters in the order 0,1,2, etc.,"
"without skipping any entries.");
previousWasEmpty = prop.empty();
}
// Number of input dimension
size_t inD = m_inWS->getNumDims();
// Validate
if (numDims == 0)
throw std::runtime_error("No output dimensions specified.");
if (numDims > inD)
throw std::runtime_error(
"More output dimensions were specified than input dimensions "
"exist in the MDEventWorkspace.");
// Create the dimensions based on the strings from the user
for (size_t i = 0; i < numDims; i++) {
std::string propName = "AlignedDim0";
propName[10] = dimChars[i];
makeAlignedDimensionFromString(getPropertyValue(propName));
}
// Number of output binning dimensions found
m_outD = m_binDimensions.size();
// Now we build the coordinate transformation object
m_translation = VMD(inD);
m_bases.clear();
std::vector<coord_t> origin(m_outD), scaling(m_outD);
for (size_t d = 0; d < m_outD; d++) {
origin[d] = m_binDimensions[d]->getMinimum();
scaling[d] = 1.0f / m_binDimensions[d]->getBinWidth();
// Origin in the input
m_translation[m_dimensionToBinFrom[d]] = origin[d];
// Create a unit basis vector that corresponds to this
VMD basis(inD);
basis[m_dimensionToBinFrom[d]] = 1.0;
m_bases.push_back(basis);
}
// Transform for binning
m_transform = new DataObjects::CoordTransformAligned(
m_inWS->getNumDims(), m_outD, m_dimensionToBinFrom, origin, scaling);
// Transformation original->binned. There is no offset or scaling!
std::vector<coord_t> unitScaling(m_outD, 1.0);
std::vector<coord_t> zeroOrigin(m_outD, 0.0);
m_transformFromOriginal = new DataObjects::CoordTransformAligned(
inD, m_outD, m_dimensionToBinFrom, zeroOrigin, unitScaling);
// Now the reverse transformation.
if (m_outD == inD) {
// Make the reverse map = if you're in the output dimension "od", what INPUT
// dimension index is that?
Matrix<coord_t> mat = m_transformFromOriginal->makeAffineMatrix();
mat.Invert();
DataObjects::CoordTransformAffine *tmp =
new DataObjects::CoordTransformAffine(inD, m_outD);
tmp->setMatrix(mat);
m_transformToOriginal = tmp;
} else {
// Changed # of dimensions - can't reverse the transform
m_transformToOriginal = NULL;
g_log.warning("SlicingAlgorithm: Your slice will cause the output "
"workspace to have less dimensions than the input. This will "
"affect your ability to create subsequent slices.");
}
}
//-----------------------------------------------------------------------------------------------
/** Read the algorithm properties and creates the appropriate transforms
* for slicing the MDEventWorkspace.
*
* NOTE: The m_inWS member must be set first.
* If the workspace is based on another, e.g. result from BinMD,
* m_inWS will be modified to be the original workspace and the transformations
* will be altered to match.
*
* The m_transform, m_transformFromOriginal and m_transformToOriginal transforms
*will be set.
*/
void SlicingAlgorithm::createTransform() {
if (!m_inWS)
throw std::runtime_error("SlicingAlgorithm::createTransform(): input "
"MDWorkspace must be set first!");
if (boost::dynamic_pointer_cast<MatrixWorkspace>(m_inWS))
throw std::runtime_error(this->name() +
" cannot be run on a MatrixWorkspace!");
// Is the transformation aligned with axes?
m_axisAligned = getProperty("AxisAligned");
// Refer to the original workspace. Make sure that is possible
if (m_inWS->numOriginalWorkspaces() > 0)
m_originalWS = boost::dynamic_pointer_cast<IMDWorkspace>(
m_inWS->getOriginalWorkspace());
if (m_originalWS) {
if (m_axisAligned)
throw std::runtime_error(
"Cannot perform axis-aligned binning on a MDHistoWorkspace. "
"Please use non-axis aligned binning.");
if (m_originalWS->getNumDims() != m_inWS->getNumDims())
throw std::runtime_error(
"SlicingAlgorithm::createTransform(): Cannot propagate "
"a transformation if the number of dimensions has changed.");
if (!m_inWS->getTransformToOriginal())
throw std::runtime_error(
"SlicingAlgorithm::createTransform(): Cannot propagate "
"a transformation. There is no transformation saved from " +
m_inWS->getName() + " back to " + m_originalWS->getName() + ".");
// Fail if the MDHistoWorkspace was modified by binary operation
DataObjects::MDHistoWorkspace_sptr inHisto =
boost::dynamic_pointer_cast<DataObjects::MDHistoWorkspace>(m_inWS);
if (inHisto) {
if (inHisto->getNumExperimentInfo() > 0) {
const Run &run = inHisto->getExperimentInfo(0)->run();
if (run.hasProperty("mdhisto_was_modified")) {
Property *prop = run.getProperty("mdhisto_was_modified");
if (prop) {
if (prop->value() == "1") {
throw std::runtime_error(
"This MDHistoWorkspace was modified by a binary operation "
"(e.g. Plus, Minus). "
"It is not currently possible to rebin a modified "
"MDHistoWorkspace because that requires returning to the "
"original "
"(unmodified) MDEventWorkspace, and so would give incorrect "
"results. "
"Instead, you can use SliceMD and perform operations on the "
"resulting "
"MDEventWorkspaces, which preserve all events. "
"You can override this check by removing the "
"'mdhisto_was_modified' sample log.");
}
}
}
}
}
g_log.notice() << "Performing " << this->name()
<< " on the original workspace, '" << m_originalWS->getName()
<< "'" << std::endl;
}
// Create the coordinate transformation
m_transform = NULL;
if (m_axisAligned)
this->createAlignedTransform();
else
this->createGeneralTransform();
// Finalize, for binnign MDHistoWorkspace
if (m_originalWS) {
// The intermediate workspace is the MDHistoWorkspace being BINNED
m_intermediateWS = m_inWS;
CoordTransform const *originalToIntermediate =
m_intermediateWS->getTransformFromOriginal();
if (originalToIntermediate &&
(m_originalWS->getNumDims() == m_intermediateWS->getNumDims())) {
try {
// The transform from the INPUT to the INTERMEDIATE ws
// intermediate_coords = [OriginalToIntermediate] * [thisToOriginal] *
// these_coords
Matrix<coord_t> matToOriginal =
m_transformToOriginal->makeAffineMatrix();
Matrix<coord_t> matOriginalToIntermediate =
originalToIntermediate->makeAffineMatrix();
Matrix<coord_t> matToIntermediate =
matOriginalToIntermediate * matToOriginal;
m_transformToIntermediate = new DataObjects::CoordTransformAffine(
m_originalWS->getNumDims(), m_intermediateWS->getNumDims());
m_transformToIntermediate->setMatrix(matToIntermediate);
// And now the reverse
matToIntermediate.Invert();
m_transformFromIntermediate = new DataObjects::CoordTransformAffine(
m_intermediateWS->getNumDims(), m_originalWS->getNumDims());
m_transformFromIntermediate->setMatrix(matToIntermediate);
} catch (std::runtime_error &) {
// Ignore error. Have no transform.
}
}
// Replace the input workspace with the original MDEventWorkspace
// for future binning
m_inWS = m_originalWS;
}
//
// if (m_inWS->hasOriginalWorkspace())
// {
// // A was transformed to B
// // Now we transform B to C
// // So we come up with the A -> C transformation
//
// IMDWorkspace_sptr origWS = m_inWS->getOriginalWorkspace();
// g_log.notice() << "Performing " << this->name() << " on the original
// workspace, '" << origWS->getName() << "'" << std::endl;
//
// if (origWS->getNumDims() != m_inWS->getNumDims())
// throw std::runtime_error("SlicingAlgorithm::createTransform():
// Cannot propagate a transformation if the number of dimensions has
// changed.");
//
// // A->C transformation
// CoordTransform * fromOrig =
// CoordTransformAffine::combineTransformations(
// m_inWS->getTransformFromOriginal(), m_transformFromOriginal );
// // C->A transformation
// CoordTransform * toOrig =
// CoordTransformAffine::combineTransformations( m_transformToOriginal,
// m_inWS->getTransformToOriginal() );
// // A->C binning transformation
// CoordTransform * binningTransform =
// CoordTransformAffine::combineTransformations(
// m_inWS->getTransformFromOriginal(), m_transform );
//
// // Replace the transforms
// delete m_transformFromOriginal;
// delete m_transformToOriginal;
// delete m_transform;
// m_transformFromOriginal = fromOrig;
// m_transformToOriginal = toOrig;
// m_transform = binningTransform;
//
// coord_t in[2] = {0,0};
// coord_t out[2] = {0,0};
// m_transform->apply(in, out);
// std::cout << "0,0 gets binningTransformed to " << VMD(2, out) <<
// std::endl;
// in[0] = 10; in[1] = 10;
// m_transform->apply(in, out);
// std::cout << "10,10 gets binningTransformed to " << VMD(2, out) <<
// std::endl;
//
// // Replace the input workspace
// m_inWS = origWS;
// }
}
//----------------------------------------------------------------------------------------------
/** Create an implicit function for picking boxes, based on the indexes in the
* output MDHistoWorkspace.
* This needs to be in the space of the INPUT MDEventWorkspace.
*
* In the most general case, this function assumes ORTHOGONAL BASIS VECTORS!
* However, the following cases handle non-orthogonal vectors:
* 2 bases in 2D space
* 2 bases in 3D space
* 3 bases in 3D space
* 3 bases in 4D space
*
* @param chunkMin :: the minimum index in each dimension to consider "valid"
*(inclusive).
* NULL to use the entire range.
* @param chunkMax :: the maximum index in each dimension to consider "valid"
*(exclusive)
* NULL to use the entire range.
* @return MDImplicitFunction created
*/
MDImplicitFunction *
SlicingAlgorithm::getGeneralImplicitFunction(const size_t *const chunkMin,
const size_t *const chunkMax) {
size_t nd = m_inWS->getNumDims();
// General implicit function
MDImplicitFunction *func = new MDImplicitFunction;
// First origin = min of each basis vector
VMD o1 = m_translation;
// Second origin = max of each basis vector
VMD o2 = m_translation;
// And this the list of basis vectors. Each vertex is given by o1+bases[i].
std::vector<VMD> bases;
VMD x, y, z, t;
for (size_t d = 0; d < m_bases.size(); d++) {
double xMin = m_binDimensions[d]->getMinimum();
double xMax = m_binDimensions[d]->getMaximum();
// Move the position if you're using a chunk
if (chunkMin != NULL)
xMin = m_binDimensions[d]->getX(chunkMin[d]);
if (chunkMax != NULL)
xMax = m_binDimensions[d]->getX(chunkMax[d]);
// Offset the origin by the position along the basis vector
o1 += (m_bases[d] * xMin);
o2 += (m_bases[d] * xMax);
VMD thisBase = m_bases[d] * (xMax - xMin);
bases.push_back(thisBase);
if (d == 0)
x = thisBase;
if (d == 1)
y = thisBase;
if (d == 2)
z = thisBase;
if (d == 3)
t = thisBase;
}
// Dimensionality of the box
size_t boxDim = bases.size();
// Create a Z vector by doing the cross-product of X and Y
if (boxDim == 2 && nd == 3) {
z = x.cross_prod(y);
z.normalize();
}
// Point that is sure to be inside the volume of interest
VMD insidePoint = (o1 + o2) / 2.0;
VMD normal = bases[0];
if (boxDim == 1) {
// 2 planes defined by 1 basis vector
// Your normal = the x vector
func->addPlane(MDPlane(x, o1));
func->addPlane(MDPlane(x * -1.0, o2));
} else if (boxDim == 2 && nd == 2) {
// 4 planes defined by 2 basis vectors (general to non-orthogonal basis
// vectors)
std::vector<VMD> vectors;
// X plane
vectors.clear();
vectors.push_back(x);
func->addPlane(MDPlane(vectors, o1, insidePoint));
func->addPlane(MDPlane(vectors, o2, insidePoint));
// Y plane
vectors.clear();
vectors.push_back(y);
func->addPlane(MDPlane(vectors, o1, insidePoint));
func->addPlane(MDPlane(vectors, o2, insidePoint));
} else if (boxDim == 2 && nd == 3) {
// 4 planes defined by 2 basis vectors (general to non-orthogonal basis
// vectors)
// The vertical = cross-product of X and Y basis vectors
z = x.cross_prod(y);
std::vector<VMD> vectors;
// X plane
vectors.clear();
vectors.push_back(x);
vectors.push_back(z);
func->addPlane(MDPlane(vectors, o1, insidePoint));
func->addPlane(MDPlane(vectors, o2, insidePoint));
// Y plane
vectors.clear();
vectors.push_back(y);
vectors.push_back(z);
func->addPlane(MDPlane(vectors, o1, insidePoint));
func->addPlane(MDPlane(vectors, o2, insidePoint));
} else if (boxDim == 3 && nd == 3) {
// 6 planes defined by 3 basis vectors (general to non-orthogonal basis
// vectors)
VMD xyNormal = x.cross_prod(y);
// Flip the normal if the input bases were NOT right-handed.
if (!MDPlane(xyNormal, o1).isPointBounded(insidePoint))
xyNormal *= -1.0;
func->addPlane(MDPlane(xyNormal, o1));
func->addPlane(MDPlane(xyNormal * -1.0, o2));
VMD xzNormal = z.cross_prod(x);
if (!MDPlane(xzNormal, o1).isPointBounded(insidePoint))
xzNormal *= -1.0;
func->addPlane(MDPlane(xzNormal, o1));
func->addPlane(MDPlane(xzNormal * -1.0, o2));
VMD yzNormal = y.cross_prod(z);
if (!MDPlane(yzNormal, o1).isPointBounded(insidePoint))
yzNormal *= -1.0;
func->addPlane(MDPlane(yzNormal, o1));
func->addPlane(MDPlane(yzNormal * -1.0, o2));
} else if (boxDim == 3 && nd == 4) {
// 6 planes defined by 3 basis vectors, in 4D world. (General to
// non-orthogonal basis vectors)
// Get a vector (in 4D) that is normal to all 3 other bases.
t = VMD::getNormalVector(bases);
// All the planes will have "t" in their plane.
std::vector<VMD> vectors;
// XY plane
vectors.clear();
vectors.push_back(x);
vectors.push_back(y);
vectors.push_back(t);
func->addPlane(MDPlane(vectors, o1, insidePoint));
func->addPlane(MDPlane(vectors, o2, insidePoint));
// XZ plane
vectors.clear();
vectors.push_back(x);
vectors.push_back(z);
vectors.push_back(t);
func->addPlane(MDPlane(vectors, o1, insidePoint));
func->addPlane(MDPlane(vectors, o2, insidePoint));
// YZ plane
vectors.clear();
vectors.push_back(y);
vectors.push_back(z);
vectors.push_back(t);
func->addPlane(MDPlane(vectors, o1, insidePoint));
func->addPlane(MDPlane(vectors, o2, insidePoint));
} else {
// Last-resort, totally general case
// 2*N planes defined by N basis vectors, in any dimensionality workspace.
// Assumes orthogonality!
for (size_t i = 0; i < bases.size(); i++) {
// For each basis vector, make two planes, perpendicular to it and facing
// inwards
func->addPlane(MDPlane(bases[i], o1));
func->addPlane(MDPlane(bases[i] * -1.0, o2));
}
}
return func;
}
//----------------------------------------------------------------------------------------------
/** Create an implicit function for picking boxes, based on the indexes in the
* output MDHistoWorkspace.
* This needs to be in the space of the INPUT MDEventWorkspace
*
* @param chunkMin :: the minimum index in each dimension to consider "valid"
*(inclusive).
* NULL to use the entire range.
* @param chunkMax :: the maximum index in each dimension to consider "valid"
*(exclusive)
* NULL to use the entire range.
* @return MDImplicitFunction created
*/
MDImplicitFunction *
SlicingAlgorithm::getImplicitFunctionForChunk(const size_t *const chunkMin,
const size_t *const chunkMax) {
size_t nd = m_inWS->getNumDims();
if (m_axisAligned) {
std::vector<coord_t> function_min(
nd, -1e30f); // default to all space if the dimension is not specified
std::vector<coord_t> function_max(
nd, +1e30f); // default to all space if the dimension is not specified
for (size_t bd = 0; bd < m_outD; bd++) {
// Dimension in the MDEventWorkspace
size_t d = m_dimensionToBinFrom[bd];
if (chunkMin)
function_min[d] = m_binDimensions[bd]->getX(chunkMin[bd]);
else
function_min[d] = m_binDimensions[bd]->getX(0);
if (chunkMax)
function_max[d] = m_binDimensions[bd]->getX(chunkMax[bd]);
else
function_max[d] =
m_binDimensions[bd]->getX(m_binDimensions[bd]->getNBins());
}
MDBoxImplicitFunction *function =
new MDBoxImplicitFunction(function_min, function_max);
return function;
} else {
// General implicit function
return getGeneralImplicitFunction(chunkMin, chunkMax);
}
}
#if 0
/**
* Create an MDFrame for the Non-Axis-Aligned case
* 1. Check that Q-based frames are not mixed with non-Q-based frames
* 2. Select Q-based frame if available
* @param units: the units
* @returns the unique pointer
*/
Mantid::Geometry::MDFrame_uptr SlicingAlgorithm::createMDFrameForNonAxisAligned(std::string units) const {
}
#endif
} // namespace Mantid
} // namespace DataObjects