diff --git a/Code/Mantid/Framework/Kernel/inc/MantidKernel/Utils.h b/Code/Mantid/Framework/Kernel/inc/MantidKernel/Utils.h
index 10d89f4b94a5b7839e1caf51088d592a54b9910c..e53ef353320531563f25afcfef32b993c83153e8 100644
--- a/Code/Mantid/Framework/Kernel/inc/MantidKernel/Utils.h
+++ b/Code/Mantid/Framework/Kernel/inc/MantidKernel/Utils.h
@@ -286,14 +286,14 @@ inline void getIndicesFromLinearIndex(const size_t linear_index,
* @param subject_indices
* @param num_bins
* @param index_max
- * @param range
+ * @param widths : width in pixels per dimension
* @return True if the are neighbours, otherwise false.
*/
inline bool isNeighbourOfSubject(const size_t ndims,
const size_t neighbour_linear_index,
const size_t *subject_indices,
const size_t *num_bins,
- const size_t *index_max, const size_t range = 1) {
+ const size_t *index_max, const std::vector<int>& widths) {
std::vector<size_t> neighbour_indices(ndims);
Utils::NestedForLoop::GetIndicesFromLinearIndex(ndims, neighbour_linear_index,
num_bins, index_max,
@@ -303,7 +303,7 @@ inline bool isNeighbourOfSubject(const size_t ndims,
long double diff =
std::abs(static_cast<long double>(subject_indices[ind]) -
static_cast<long double>(neighbour_indices[ind]));
- if (diff > range) {
+ if (diff > widths[ind]/2) {
return false;
}
}
diff --git a/Code/Mantid/Framework/MDAlgorithms/src/SmoothMD.cpp b/Code/Mantid/Framework/MDAlgorithms/src/SmoothMD.cpp
index 8d7fcbc8123fa68abe576e181cec2a481d40b865..5412eff169a7d857a95bbc3165a3adea942b0c7d 100644
--- a/Code/Mantid/Framework/MDAlgorithms/src/SmoothMD.cpp
+++ b/Code/Mantid/Framework/MDAlgorithms/src/SmoothMD.cpp
@@ -31,11 +31,13 @@ using namespace Mantid::MDEvents;
typedef std::vector<int> WidthVector;
// Typedef for an optional md histo workspace
-typedef boost::optional<IMDHistoWorkspace_const_sptr> OptionalIMDHistoWorkspace_const_sptr;
+typedef boost::optional<IMDHistoWorkspace_const_sptr>
+ OptionalIMDHistoWorkspace_const_sptr;
// Typedef for a smoothing function
typedef boost::function<IMDHistoWorkspace_sptr(
- IMDHistoWorkspace_const_sptr, const WidthVector &, OptionalIMDHistoWorkspace_const_sptr)> SmoothFunction;
+ IMDHistoWorkspace_const_sptr, const WidthVector &,
+ OptionalIMDHistoWorkspace_const_sptr)> SmoothFunction;
// Typedef for a smoothing function map keyed by name.
typedef std::map<std::string, SmoothFunction> SmoothFunctionMap;
@@ -57,12 +59,13 @@ std::vector<std::string> functions() {
* Maps a function name to a smoothing function
* @return function map
*/
-SmoothFunctionMap makeFunctionMap(Mantid::MDAlgorithms::SmoothMD * instance) {
+SmoothFunctionMap makeFunctionMap(Mantid::MDAlgorithms::SmoothMD *instance) {
SmoothFunctionMap map;
- map.insert(std::make_pair("Hat", boost::bind( &Mantid::MDAlgorithms::SmoothMD::hatSmooth, instance, _1, _2, _3 )));
+ map.insert(std::make_pair(
+ "Hat", boost::bind(&Mantid::MDAlgorithms::SmoothMD::hatSmooth, instance,
+ _1, _2, _3)));
return map;
}
-
}
namespace Mantid {
@@ -105,15 +108,18 @@ const std::string SmoothMD::summary() const {
* @param weightingWS : Weighting workspace (optional)
* @return Smoothed MDHistoWorkspace
*/
-IMDHistoWorkspace_sptr SmoothMD::hatSmooth(IMDHistoWorkspace_const_sptr toSmooth,
- const WidthVector &widthVector, OptionalIMDHistoWorkspace_const_sptr weightingWS) {
+IMDHistoWorkspace_sptr
+SmoothMD::hatSmooth(IMDHistoWorkspace_const_sptr toSmooth,
+ const WidthVector &widthVector,
+ OptionalIMDHistoWorkspace_const_sptr weightingWS) {
const bool useWeights = weightingWS.is_initialized();
uint64_t nPoints = toSmooth->getNPoints();
- Progress progress(this, 0, 1, size_t( double(nPoints) * 1.1 ) );
+ Progress progress(this, 0, 1, size_t(double(nPoints) * 1.1));
// Create the output workspace.
IMDHistoWorkspace_sptr outWS = toSmooth->clone();
- progress.reportIncrement(size_t(double(nPoints) * 0.1)); // Report ~10% progress
+ progress.reportIncrement(
+ size_t(double(nPoints) * 0.1)); // Report ~10% progress
const int nThreads = Mantid::API::FrameworkManager::Instance()
.getNumOMPThreads(); // NThreads to Request
@@ -129,51 +135,48 @@ IMDHistoWorkspace_sptr SmoothMD::hatSmooth(IMDHistoWorkspace_const_sptr toSmooth
do {
// Gets all vertex-touching neighbours
- size_t iteratorIndex = iterator->getLinearIndex();
+ size_t iteratorIndex = iterator->getLinearIndex();
- if(useWeights) {
+ if (useWeights) {
- // Check that we could measuer here.
- if ( (*weightingWS)->getSignalAt(iteratorIndex) == 0 ) {
+ // Check that we could measuer here.
+ if ((*weightingWS)->getSignalAt(iteratorIndex) == 0) {
- outWS->setSignalAt(iteratorIndex, std::numeric_limits<double>::quiet_NaN());
+ outWS->setSignalAt(iteratorIndex,
+ std::numeric_limits<double>::quiet_NaN());
- outWS->setErrorSquaredAt(iteratorIndex, std::numeric_limits<double>::quiet_NaN());
+ outWS->setErrorSquaredAt(iteratorIndex,
+ std::numeric_limits<double>::quiet_NaN());
- continue; // Skip we couldn't measure here.
- }
+ continue; // Skip we couldn't measure here.
+ }
}
std::vector<size_t> neighbourIndexes =
- iterator->findNeighbourIndexesByWidth(
- widthVector.front()); // TODO we should use the whole width vector
- // not just the first element.
+ iterator->findNeighbourIndexesByWidth(widthVector);
+
size_t nNeighbours = neighbourIndexes.size();
double sumSignal = iterator->getSignal();
double sumSqError = iterator->getError();
for (size_t i = 0; i < neighbourIndexes.size(); ++i) {
- if(useWeights) {
- if ( (*weightingWS)->getSignalAt(neighbourIndexes[i]) == 0 )
- {
- // Nothing measured here. We cannot use that neighbouring point.
- nNeighbours -= 1;
- continue;
- }
+ if (useWeights) {
+ if ((*weightingWS)->getSignalAt(neighbourIndexes[i]) == 0) {
+ // Nothing measured here. We cannot use that neighbouring point.
+ nNeighbours -= 1;
+ continue;
+ }
}
sumSignal += toSmooth->getSignalAt(neighbourIndexes[i]);
double error = toSmooth->getErrorAt(neighbourIndexes[i]);
sumSqError += (error * error);
-
}
// Calculate the mean
- outWS->setSignalAt(iteratorIndex,
- sumSignal / double(nNeighbours + 1));
+ outWS->setSignalAt(iteratorIndex, sumSignal / double(nNeighbours + 1));
// Calculate the sample variance
outWS->setErrorSquaredAt(iteratorIndex,
sumSqError / double(nNeighbours + 1));
-
progress.report();
} while (iterator->next());
@@ -184,7 +187,6 @@ IMDHistoWorkspace_sptr SmoothMD::hatSmooth(IMDHistoWorkspace_const_sptr toSmooth
return outWS;
}
-
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
@@ -219,7 +221,8 @@ void SmoothMD::init() {
docBuffer.str());
declareProperty(new WorkspaceProperty<API::IMDHistoWorkspace>(
- "InputNormalizationWorkspace", "", Direction::Input, PropertyMode::Optional),
+ "InputNormalizationWorkspace", "", Direction::Input,
+ PropertyMode::Optional),
"Multidimensional weighting workspace. Optional.");
declareProperty(new WorkspaceProperty<API::IMDHistoWorkspace>(
@@ -236,15 +239,21 @@ void SmoothMD::exec() {
IMDHistoWorkspace_sptr toSmooth = this->getProperty("InputWorkspace");
// Get the input weighting workspace
- IMDHistoWorkspace_sptr weightingWS = this->getProperty("InputNormalizationWorkspace");
+ IMDHistoWorkspace_sptr weightingWS =
+ this->getProperty("InputNormalizationWorkspace");
OptionalIMDHistoWorkspace_const_sptr optionalWeightingWS;
- if(weightingWS)
- {
- optionalWeightingWS = weightingWS;
+ if (weightingWS) {
+ optionalWeightingWS = weightingWS;
}
// Get the width vector
- const std::vector<int> widthVector = this->getProperty("WidthVector");
+ std::vector<int> widthVector = this->getProperty("WidthVector");
+ if (widthVector.size() == 1) {
+ // Pad the width vector out to the right size if only one entry has been
+ // provided.
+ widthVector =
+ std::vector<int>(toSmooth->getNumDims(), widthVector.front());
+ }
// Find the choosen smooth operation
const std::string smoothFunctionName = this->getProperty("Function");
@@ -264,55 +273,67 @@ std::map<std::string, std::string> SmoothMD::validateInputs() {
std::map<std::string, std::string> product;
+ IMDHistoWorkspace_sptr toSmoothWs = this->getProperty("InputWorkspace");
+
// Check the width vector
const std::string widthVectorPropertyName = "WidthVector";
std::vector<int> widthVector = this->getProperty(widthVectorPropertyName);
- for (auto it = widthVector.begin(); it != widthVector.end(); ++it) {
- const int widthEntry = *it;
- if (widthEntry % 2 == 0) {
- std::stringstream message;
- message << widthVectorPropertyName
- << " entries must be odd numbers. Bad entry is " << widthEntry;
- product.insert(std::make_pair(widthVectorPropertyName, message.str()));
+
+ if (widthVector.size() != 1 && widthVector.size() != toSmoothWs->getNumDims()) {
+ product.insert(std::make_pair(widthVectorPropertyName,
+ widthVectorPropertyName +
+ " can either have one entry or needs to "
+ "have entries for each dimension of the "
+ "InputWorkspace."));
+ } else {
+ for (auto it = widthVector.begin(); it != widthVector.end(); ++it) {
+ const int widthEntry = *it;
+ if (widthEntry % 2 == 0) {
+ std::stringstream message;
+ message << widthVectorPropertyName
+ << " entries must be odd numbers. Bad entry is " << widthEntry;
+ product.insert(std::make_pair(widthVectorPropertyName, message.str()));
+ }
}
}
// Check the dimensionality of the normalization workspace
- const std::string normalisationWorkspacePropertyName = "InputNormalizationWorkspace";
- IMDHistoWorkspace_sptr toSmoothWs = this->getProperty("InputWorkspace");
- IMDHistoWorkspace_sptr normWs = this->getProperty(normalisationWorkspacePropertyName);
- if(normWs)
- {
- const size_t nDimsNorm = normWs->getNumDims();
- const size_t nDimsSmooth = toSmoothWs->getNumDims();
- if(nDimsNorm != nDimsSmooth) {
+ const std::string normalisationWorkspacePropertyName =
+ "InputNormalizationWorkspace";
+
+ IMDHistoWorkspace_sptr normWs =
+ this->getProperty(normalisationWorkspacePropertyName);
+ if (normWs) {
+ const size_t nDimsNorm = normWs->getNumDims();
+ const size_t nDimsSmooth = toSmoothWs->getNumDims();
+ if (nDimsNorm != nDimsSmooth) {
+ std::stringstream message;
+ message << normalisationWorkspacePropertyName
+ << " has a different number of dimensions than InputWorkspace. "
+ "Shapes of inputs must be the same. Cannot continue "
+ "smoothing.";
+ product.insert(
+ std::make_pair(normalisationWorkspacePropertyName, message.str()));
+ } else {
+ // Loop over dimensions and check nbins.
+ for (size_t i = 0; i < nDimsNorm; ++i) {
+ const size_t nBinsNorm = normWs->getDimension(i)->getNBins();
+ const size_t nBinsSmooth = toSmoothWs->getDimension(i)->getNBins();
+ if (nBinsNorm != nBinsSmooth) {
std::stringstream message;
message << normalisationWorkspacePropertyName
- << " has a different number of dimensions than InputWorkspace. Shapes of inputs must be the same. Cannot continue smoothing.";
- product.insert(std::make_pair(normalisationWorkspacePropertyName, message.str()));
- }
- else
- {
- // Loop over dimensions and check nbins.
- for(size_t i = 0; i < nDimsNorm; ++i)
- {
- const size_t nBinsNorm = normWs->getDimension(i)->getNBins();
- const size_t nBinsSmooth = toSmoothWs->getDimension(i)->getNBins();
- if(nBinsNorm != nBinsSmooth)
- {
- std::stringstream message;
- message << normalisationWorkspacePropertyName
- << ". Number of bins from dimension with index "<< i << " do not match. "
- << nBinsSmooth << " expected. Got " << nBinsNorm
- << ". Shapes of inputs must be the same. Cannot continue smoothing.";
- product.insert(std::make_pair(normalisationWorkspacePropertyName, message.str()));
- break;
- }
- }
+ << ". Number of bins from dimension with index " << i
+ << " do not match. " << nBinsSmooth << " expected. Got "
+ << nBinsNorm << ". Shapes of inputs must be the same. Cannot "
+ "continue smoothing.";
+ product.insert(std::make_pair(normalisationWorkspacePropertyName,
+ message.str()));
+ break;
+ }
}
+ }
}
-
return product;
}
diff --git a/Code/Mantid/Framework/MDAlgorithms/test/SmoothMDTest.h b/Code/Mantid/Framework/MDAlgorithms/test/SmoothMDTest.h
index dd00e34212beb1784f171bbf0606f961a7421fbc..1920aed9aa6e530c223a01cf488d7b52dba3a52f 100644
--- a/Code/Mantid/Framework/MDAlgorithms/test/SmoothMDTest.h
+++ b/Code/Mantid/Framework/MDAlgorithms/test/SmoothMDTest.h
@@ -75,6 +75,24 @@ public:
alg.execute(), std::runtime_error &);
}
+ void test_width_vector_must_not_be_arbitrary_size() {
+ auto toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(
+ 1 /*signal*/, 2 /*numDims*/, 3 /*numBins in each dimension*/);
+
+ std::vector<int> badWidths(11, 3); // odd number value = 3, but size of 11 has no meaning
+
+ SmoothMD alg;
+ alg.setChild(true);
+ alg.initialize();
+ alg.setPropertyValue("OutputWorkspace", "dummy");
+ alg.setProperty("InputWorkspace", toSmooth);
+ alg.setProperty(
+ "WidthVector",
+ badWidths); // Width vector is the wrong size
+ TSM_ASSERT_THROWS("Size of with vector is wrong should throw.", alg.execute(),
+ std::runtime_error &);
+ }
+
void test_simple_smooth_hat_function() {
auto toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(
1 /*signal*/, 2 /*numDims*/, 3 /*numBins in each dimension*/);
@@ -207,6 +225,59 @@ public:
TS_ASSERT_EQUALS(z, out->getSignalAt(12));
}
+
+ void test_smooth_hat_function_mixed_widths() {
+
+ auto toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(
+ 2 /*signal*/, 2 /*numDims*/, 5 /*numBins in each dimension*/);
+ toSmooth->setSignalAt(2, 1.0);
+ toSmooth->setSignalAt(22, 1.0);
+
+ /*
+ 2D MDHistoWorkspace Input
+
+ 2 - 2 - 1 - 2 - 2
+ 2 - 2 - 2 - 2 - 2
+ 2 - 2 - 2 - 2 - 2
+ 2 - 2 - 2 - 2 - 2
+ 2 - 2 - 1 - 2 - 2
+
+ */
+
+ SmoothMD alg;
+ alg.setChild(true);
+ alg.initialize();
+ std::vector<int> widthVector;
+ widthVector.push_back(3); // 3 = width in zeroth dimension
+ widthVector.push_back(5); // 5 = width in first dimension
+ alg.setProperty("WidthVector", widthVector);
+ alg.setProperty("InputWorkspace", toSmooth);
+ alg.setPropertyValue("OutputWorkspace", "dummy");
+ alg.execute();
+ IMDHistoWorkspace_sptr out = alg.getProperty("OutputWorkspace");
+
+ /*
+ Smoothing region for centre point
+
+ 2 - |2 - 1 - 2| - 2
+ 2 - |2 - 2 - 2| - 2
+ 2 - |2 -|2|- 2| - 2
+ 2 - |2 - 2 - 2| - 2
+ 2 - |2 - 1 - 2| - 2
+
+ 3 by 5 with.
+
+ average for centre point should be [ (3 * 5) - 2 ] * 2 / (3 * 5) = 28 / 15
+
+ */
+
+ // Check vertexes
+ double expectedSmoothedValue = 28.0 / 15;
+
+ TS_ASSERT_EQUALS(expectedSmoothedValue, out->getSignalAt(12));
+ }
+
+
void test_dimensional_check_of_weight_ws() {
MDHistoWorkspace_sptr a =
diff --git a/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspaceIterator.h b/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspaceIterator.h
index af649ba80233cf8053c8e1f078d28634dde785f6..5e88c423bc9abef3d551a2ec6e14a3f95637fb13 100644
--- a/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspaceIterator.h
+++ b/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspaceIterator.h
@@ -7,12 +7,13 @@
#include "MantidGeometry/MDGeometry/MDImplicitFunction.h"
#include "MantidMDEvents/SkippingPolicy.h"
#include <map>
+#include <vector>
namespace Mantid {
namespace MDEvents {
// Typdef for a map for mapping width of neighbours (key) to permutations needed in the calcualtion.
-typedef std::map<size_t, std::vector<int64_t> > PermutationsMap;
+typedef std::map<std::vector<int>, std::vector<int64_t> > PermutationsMap;
/** An implementation of IMDIterator that iterates through
a MDHistoWorkspace. It treats the bin in the workspace as
@@ -113,7 +114,9 @@ public:
std::vector<size_t> findNeighbourIndexesFaceTouching() const;
- std::vector<size_t> findNeighbourIndexesByWidth(const int width) const;
+ std::vector<size_t> findNeighbourIndexesByWidth(const int& width) const;
+
+ std::vector<size_t> findNeighbourIndexesByWidth(const std::vector<int>& widths) const;
virtual bool isWithinBounds(size_t index) const;
@@ -166,7 +169,7 @@ protected:
SkippingPolicy_scptr m_skippingPolicy;
/// Create or fetch permutations relating to a given neighbour width.
- std::vector<int64_t> createPermutations(const int width) const;
+ std::vector<int64_t> createPermutations(const std::vector<int>& widths) const;
};
} // namespace MDEvents
diff --git a/Code/Mantid/Framework/MDEvents/src/MDHistoWorkspaceIterator.cpp b/Code/Mantid/Framework/MDEvents/src/MDHistoWorkspaceIterator.cpp
index d1510c15ba339c55ef6f90d7c19a79d320b60a2d..a7fdba24e77080849e78677a9f8d68e2d519cf81 100644
--- a/Code/Mantid/Framework/MDEvents/src/MDHistoWorkspaceIterator.cpp
+++ b/Code/Mantid/Framework/MDEvents/src/MDHistoWorkspaceIterator.cpp
@@ -5,6 +5,8 @@
#include "MantidGeometry/MDGeometry/IMDDimension.h"
#include <algorithm>
#include <utility>
+#include <functional>
+#include <numeric>
using namespace Mantid::Kernel;
using namespace Mantid::API;
@@ -418,6 +420,7 @@ MDHistoWorkspaceIterator::findNeighbourIndexesFaceTouching() const {
m_indexMax, m_index);
std::vector<size_t> neighbourIndexes; // Accumulate neighbour indexes.
+ std::vector<int> widths(m_nd, 3); // Face touching width is always 3 in each dimension
for (size_t i = 0; i < m_permutationsFaceTouching.size(); ++i) {
if (m_permutationsFaceTouching[i] == 0) {
continue;
@@ -426,7 +429,7 @@ MDHistoWorkspaceIterator::findNeighbourIndexesFaceTouching() const {
size_t neighbour_index = m_pos + m_permutationsFaceTouching[i];
if (neighbour_index < m_ws->getNPoints() &&
Utils::isNeighbourOfSubject(m_nd, neighbour_index, m_index,
- m_indexMaker, m_indexMax)) {
+ m_indexMaker, m_indexMax, widths)) {
neighbourIndexes.push_back(neighbour_index);
}
}
@@ -449,30 +452,36 @@ bool MDHistoWorkspaceIterator::isWithinBounds(size_t index) const {
*the iterator is moved and the method is called,
* we can cache the results, and re-use them as the only factors are the and the
*dimensionality, the width (n-neighbours).
- *
- * @return
+ * @param widths : vector of integer widths.
+ * @return index permutations
*/
std::vector<int64_t>
-MDHistoWorkspaceIterator::createPermutations(const int width) const {
+MDHistoWorkspaceIterator::createPermutations(const std::vector<int>& widths) const {
// look-up
- PermutationsMap::iterator it = m_permutationsVertexTouchingMap.find(width);
+ PermutationsMap::iterator it = m_permutationsVertexTouchingMap.find(widths);
if (it == m_permutationsVertexTouchingMap.end()) {
- if (width % 2 == 0) {
+ if (widths[0] % 2 == 0) {
throw std::invalid_argument("MDHistoWorkspaceIterator::"
"findNeighbourIndexesByWidth, width must "
"always be an even number");
}
+ if (widths.size() != m_nd) {
+ throw std::invalid_argument("MDHistoWorkspaceIterator::"
+ "findNeighbourIndexesByWidth, size of widths must be the same as the number of dimensions.");
+ }
int64_t offset = 1;
// Size of block will be width ^ nd
std::vector<int64_t> permutationsVertexTouching;
- permutationsVertexTouching.reserve(integerPower(width, m_nd));
+ // Calculate maximum permutations size.
+ int product = std::accumulate(widths.begin(), widths.end(), 1, std::multiplies<int>());
+ permutationsVertexTouching.reserve(product);
- int centreIndex = (width) / 2; // Deliberately truncate to get centre index
+ int centreIndex = widths[0] / 2; // Deliberately truncate to get centre index
- for (int i = 0; i < width; ++i) {
+ for (int i = 0; i < widths[0]; ++i) {
// for width = 3 : -1, 0, 1
// for width = 5 : -2, -1, 0, 1, 2
permutationsVertexTouching.push_back(centreIndex - i);
@@ -485,7 +494,7 @@ MDHistoWorkspaceIterator::createPermutations(const int width) const {
offset * static_cast<int64_t>(m_ws->getDimension(j - 1)->getNBins());
size_t nEntries = permutationsVertexTouching.size();
- for (int k = 1; k <= width / 2; ++k) {
+ for (int k = 1; k <= widths[j] / 2; ++k) {
for (size_t m = 0; m < nEntries; m++) {
permutationsVertexTouching.push_back((offset * k) +
permutationsVertexTouching[m]);
@@ -495,23 +504,38 @@ MDHistoWorkspaceIterator::createPermutations(const int width) const {
}
}
- m_permutationsVertexTouchingMap.insert(std::make_pair(width, permutationsVertexTouching));
+ m_permutationsVertexTouchingMap.insert(std::make_pair(widths, permutationsVertexTouching));
}
// In either case, get the result.
- return m_permutationsVertexTouchingMap[width];
+ return m_permutationsVertexTouchingMap[widths];
}
/**
* Find vertex-touching neighbours.
+ *
+ * Expands out the width to make a n-dimensional vector filled with the requested width.
+ *
* @param width : Odd number of pixels for all dimensions.
* @return collection of indexes.
*/
std::vector<size_t>
-MDHistoWorkspaceIterator::findNeighbourIndexesByWidth(const int width) const {
+MDHistoWorkspaceIterator::findNeighbourIndexesByWidth(const int & width) const {
+
+ return this->findNeighbourIndexesByWidth(std::vector<int>(m_nd, width));
+}
+
+
+/**
+ * Find vertex-touching neighbours.
+ * @param widths : Vector containing odd number of pixels per dimension. Entries match dimensions of iterator.
+ * @return collection of indexes.
+ */
+std::vector<size_t>
+MDHistoWorkspaceIterator::findNeighbourIndexesByWidth(const std::vector<int>& widths) const {
// Find existing or create required index permutations.
- std::vector<int64_t> permutationsVertexTouching = createPermutations(width);
+ std::vector<int64_t> permutationsVertexTouching = createPermutations(widths);
Utils::NestedForLoop::GetIndicesFromLinearIndex(m_nd, m_pos, m_indexMaker,
m_indexMax, m_index);
@@ -527,7 +551,7 @@ MDHistoWorkspaceIterator::findNeighbourIndexesByWidth(const int width) const {
size_t neighbour_index = m_pos + permutationsVertexTouching[i];
if (neighbour_index < m_ws->getNPoints() &&
Utils::isNeighbourOfSubject(m_nd, neighbour_index, m_index,
- m_indexMaker, m_indexMax, width / 2)) {
+ m_indexMaker, m_indexMax, widths) ) {
neighbourIndexes.push_back(neighbour_index);
}
}
diff --git a/Code/Mantid/Framework/MDEvents/test/MDHistoWorkspaceIteratorTest.h b/Code/Mantid/Framework/MDEvents/test/MDHistoWorkspaceIteratorTest.h
index 50a722742a840ac2f69914a755a8aa90818017e1..2f5bb66efc31c7b79817d11465b1b8807fd03395 100644
--- a/Code/Mantid/Framework/MDEvents/test/MDHistoWorkspaceIteratorTest.h
+++ b/Code/Mantid/Framework/MDEvents/test/MDHistoWorkspaceIteratorTest.h
@@ -800,16 +800,93 @@ public:
neighbourIndexes = it->findNeighbourIndexesByWidth(width);
TS_ASSERT_EQUALS(8, neighbourIndexes.size());
// Is on an edge
- TSM_ASSERT( "Neighbour at index 0 is 5", doesContainIndex(neighbourIndexes, 5));
- TSM_ASSERT( "Neighbour at index 0 is 6", doesContainIndex(neighbourIndexes, 6));
- TSM_ASSERT( "Neighbour at index 0 is 7", doesContainIndex(neighbourIndexes, 7));
- TSM_ASSERT( "Neighbour at index 0 is 9", doesContainIndex(neighbourIndexes, 9));
- TSM_ASSERT( "Neighbour at index 0 is 10", doesContainIndex(neighbourIndexes, 10));
- TSM_ASSERT( "Neighbour at index 0 is 11", doesContainIndex(neighbourIndexes, 11));
- TSM_ASSERT( "Neighbour at index 0 is 13", doesContainIndex(neighbourIndexes, 13));
- TSM_ASSERT( "Neighbour at index 0 is 14", doesContainIndex(neighbourIndexes, 14));
+ TSM_ASSERT( "Neighbour at index is 5", doesContainIndex(neighbourIndexes, 5));
+ TSM_ASSERT( "Neighbour at index is 6", doesContainIndex(neighbourIndexes, 6));
+ TSM_ASSERT( "Neighbour at index is 7", doesContainIndex(neighbourIndexes, 7));
+ TSM_ASSERT( "Neighbour at index is 9", doesContainIndex(neighbourIndexes, 9));
+ TSM_ASSERT( "Neighbour at index is 10", doesContainIndex(neighbourIndexes, 10));
+ TSM_ASSERT( "Neighbour at index is 11", doesContainIndex(neighbourIndexes, 11));
+ TSM_ASSERT( "Neighbour at index is 13", doesContainIndex(neighbourIndexes, 13));
+ TSM_ASSERT( "Neighbour at index is 14", doesContainIndex(neighbourIndexes, 14));
}
+ void test_neighbours_2d_vertex_touching_by_width_vector()
+ {
+ const size_t nd = 2;
+ std::vector<int> widthVector;
+ widthVector.push_back(5);
+ widthVector.push_back(3);
+
+ MDHistoWorkspace_sptr ws = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.0, nd, 4);
+ /*
+ 2D MDHistoWorkspace
+
+ 0 - 1 - 2 - 3
+ 4 - 5 - 6 - 7
+ 8 - 9 -10 -11
+ 12-13 -14 -15
+ */
+ MDHistoWorkspaceIterator * it = new MDHistoWorkspaceIterator(ws);
+
+ // At initial position
+ /*
+ |0| - 1 - 2 - 3
+ 4 - 5 - 6 - 7
+ 8 - 9 -10 -11
+ 12-13 -14 -15
+ */
+ std::vector<size_t> neighbourIndexes = it->findNeighbourIndexesByWidth(widthVector);
+ TS_ASSERT_EQUALS(5, neighbourIndexes.size());
+ // Is on an edge
+ TSM_ASSERT( "Neighbour at index is 1", doesContainIndex(neighbourIndexes, 1));
+ TSM_ASSERT( "Neighbour at index is 2", doesContainIndex(neighbourIndexes, 2));
+ TSM_ASSERT( "Neighbour at index is 4", doesContainIndex(neighbourIndexes, 4));
+ TSM_ASSERT( "Neighbour at index is 5", doesContainIndex(neighbourIndexes, 5));
+ TSM_ASSERT( "Neighbour at index is 6", doesContainIndex(neighbourIndexes, 6));
+
+
+
+ // At centreish position
+ /*
+ 0 - 1 - 2 - 3
+ 4 - |5| - 6 - 7
+ 8 - 9 -10 -11
+ 12-13 -14 -15
+ */
+ it->jumpTo(5);
+ neighbourIndexes = it->findNeighbourIndexesByWidth(widthVector);
+ TS_ASSERT_EQUALS(11, neighbourIndexes.size());
+ // Is on an edge
+ for(int i = 0; i < 12; ++i)
+ {
+ if(i == 5)
+ {
+ continue; // skip over the current index of the iterator.
+ }
+ std::stringstream buffer;
+ buffer << "Neighbour at index 5 should include " << i;
+ TSM_ASSERT( buffer.str(), doesContainIndex(neighbourIndexes, i));
+ }
+
+ // At end position
+ /*
+ 0 - 1 - 2 - 3
+ 4 - 5 - 6 - 7
+ 8 - 9 -10 -11
+ 12-13 -14 -|15|
+ */
+ it->jumpTo(15);
+ neighbourIndexes = it->findNeighbourIndexesByWidth(widthVector);
+ TS_ASSERT_EQUALS(5, neighbourIndexes.size());
+ // Is on an edge
+ TSM_ASSERT( "Neighbour at index is 9", doesContainIndex(neighbourIndexes, 9));
+ TSM_ASSERT( "Neighbour at index is 10", doesContainIndex(neighbourIndexes, 10));
+ TSM_ASSERT( "Neighbour at index is 11", doesContainIndex(neighbourIndexes, 11));
+ TSM_ASSERT( "Neighbour at index is 13", doesContainIndex(neighbourIndexes, 13));
+ TSM_ASSERT( "Neighbour at index is 14", doesContainIndex(neighbourIndexes, 14));
+ }
+
+
void test_neighbours_3d_vertex_touching_width()
{
const size_t nd = 3;