diff --git a/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/MDHistoWorkspaceIterator.h b/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/MDHistoWorkspaceIterator.h
index 4eae0c223115210a887c3023af9a217476b9237f..5a914cbb28f10b974934fbdf3e2352fc8f2c0e46 100644
--- a/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/MDHistoWorkspaceIterator.h
+++ b/Code/Mantid/Framework/DataObjects/inc/MantidDataObjects/MDHistoWorkspaceIterator.h
@@ -6,10 +6,15 @@
#include "MantidDataObjects/MDHistoWorkspace.h"
#include "MantidGeometry/MDGeometry/MDImplicitFunction.h"
#include "MantidDataObjects/SkippingPolicy.h"
+#include <map>
+#include <vector>
namespace Mantid {
namespace DataObjects {
+// Typdef for a map for mapping width of neighbours (key) to permutations needed in the calcualtion.
+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
a box containing a single "event", at the center of each bin
@@ -109,8 +114,14 @@ public:
std::vector<size_t> findNeighbourIndexesFaceTouching() 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;
+ size_t permutationCacheSize() const;
+
protected:
/// The MDHistoWorkspace being iterated.
const MDHistoWorkspace *m_ws;
@@ -148,12 +159,17 @@ protected:
/// Array to find indices from linear indices
size_t *m_indexMaker;
- /// Neighbour finding permutations.
- mutable std::vector<int64_t> m_permutationsVertexTouching;
+ /// Neigbour finding permutations for face touching neighbours (3 by 3 width).
mutable std::vector<int64_t> m_permutationsFaceTouching;
+ /// Neighbour finding permutations map for vertex touching. Keyed via the width (n-pixels) of neighbours required.
+ mutable PermutationsMap m_permutationsVertexTouchingMap;
+
/// Skipping policy.
SkippingPolicy_scptr m_skippingPolicy;
+
+ /// Create or fetch permutations relating to a given neighbour width.
+ std::vector<int64_t> createPermutations(const std::vector<int>& widths) const;
};
} // namespace DataObjects
diff --git a/Code/Mantid/Framework/DataObjects/src/MDHistoWorkspaceIterator.cpp b/Code/Mantid/Framework/DataObjects/src/MDHistoWorkspaceIterator.cpp
index 0213d77f9be76d9fd52887e0033a19ad5d34d5ed..c30405614890f9203f1efb2565a8265682b7cd2d 100644
--- a/Code/Mantid/Framework/DataObjects/src/MDHistoWorkspaceIterator.cpp
+++ b/Code/Mantid/Framework/DataObjects/src/MDHistoWorkspaceIterator.cpp
@@ -3,6 +3,10 @@
#include "MantidKernel/VMD.h"
#include "MantidKernel/Utils.h"
#include "MantidGeometry/MDGeometry/IMDDimension.h"
+#include <algorithm>
+#include <utility>
+#include <functional>
+#include <numeric>
using namespace Mantid::Kernel;
using namespace Mantid::API;
@@ -10,15 +14,7 @@ using Mantid::Geometry::IMDDimension_const_sptr;
namespace Mantid {
namespace DataObjects {
-namespace {
-size_t integerPower(const size_t base, const size_t pow) {
- size_t result = 1;
- for (size_t i = 0; i < pow; ++i) {
- result *= base;
- }
- return result;
-}
-}
+
//----------------------------------------------------------------------------------------------
/** Constructor
@@ -159,39 +155,24 @@ MDHistoWorkspaceIterator::init(const MDHistoWorkspace *workspace,
next();
}
- // --- Calculate index permutations for neighbour finding vertex touching ---
- auto temp = std::vector<int64_t>(integerPower(3, m_nd), 0);
- m_permutationsVertexTouching.swap(temp);
-
// --- Calculate index permutations for neighbour finding face touching ---
- temp = std::vector<int64_t>(2 * m_nd);
+ auto temp = std::vector<int64_t>(2 * m_nd);
m_permutationsFaceTouching.swap(temp);
int64_t offset = 1;
- m_permutationsVertexTouching[0] = 0;
- m_permutationsVertexTouching[1] = 1;
- m_permutationsVertexTouching[2] = -1;
m_permutationsFaceTouching[0] = -1;
m_permutationsFaceTouching[1] = 1;
// Figure out what possible indexes deltas to generate indexes that are next
// to the current one.
- size_t nVertexTouchingPermutations = 3;
for (size_t j = 1; j < m_nd; ++j) {
offset =
offset * static_cast<int64_t>(m_ws->getDimension(j - 1)->getNBins());
- size_t counter = nVertexTouchingPermutations;
- for (size_t k = 0; k < nVertexTouchingPermutations; k += 1, counter += 2) {
- int64_t newVariant = m_permutationsVertexTouching[k] + offset;
- m_permutationsVertexTouching[counter] = newVariant;
- m_permutationsVertexTouching[counter + 1] = (-1 * newVariant);
- }
m_permutationsFaceTouching[j * 2] = offset;
m_permutationsFaceTouching[(j * 2) + 1] = -offset;
- nVertexTouchingPermutations *= 3;
}
}
@@ -412,24 +393,7 @@ size_t MDHistoWorkspaceIterator::getLinearIndex() const { return m_pos; }
*/
std::vector<size_t> MDHistoWorkspaceIterator::findNeighbourIndexes() const {
- Utils::NestedForLoop::GetIndicesFromLinearIndex(m_nd, m_pos, m_indexMaker,
- m_indexMax, m_index);
-
- // Filter out indexes that are are not actually neighbours.
- std::vector<size_t> neighbourIndexes; // Accumulate neighbour indexes.
- for (size_t i = 0; i < m_permutationsVertexTouching.size(); ++i) {
- if (m_permutationsVertexTouching[i] == 0) {
- continue;
- }
-
- size_t neighbour_index = m_pos + m_permutationsVertexTouching[i];
- if (neighbour_index < m_ws->getNPoints() &&
- Utils::isNeighbourOfSubject(m_nd, neighbour_index, m_index,
- m_indexMaker, m_indexMax)) {
- neighbourIndexes.push_back(neighbour_index);
- }
- }
- return neighbourIndexes;
+ return this->findNeighbourIndexesByWidth(3 /*immediate neighbours only*/);
}
/**
@@ -448,6 +412,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;
@@ -456,7 +421,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);
}
}
@@ -472,5 +437,133 @@ bool MDHistoWorkspaceIterator::isWithinBounds(size_t index) const {
return index >= m_begin && index < m_max;
}
+/**
+ * This is to create the permutations needed to operate find neighbours in the
+ *vertex-touching schenarios
+ * Rather than having to fabricate what the possible permutations are each time
+ *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).
+ * @param widths : vector of integer widths.
+ * @return index permutations
+ */
+std::vector<int64_t>
+MDHistoWorkspaceIterator::createPermutations(const std::vector<int>& widths) const {
+ // look-up
+ PermutationsMap::iterator it = m_permutationsVertexTouchingMap.find(widths);
+ if (it == m_permutationsVertexTouchingMap.end()) {
+
+ 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;
+ // Calculate maximum permutations size.
+ int product = std::accumulate(widths.begin(), widths.end(), 1, std::multiplies<int>());
+ permutationsVertexTouching.reserve(product);
+
+ int centreIndex = widths[0] / 2; // Deliberately truncate to get centre index
+
+ 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);
+ }
+
+ // Figure out what possible indexes deltas to generate indexes that are next
+ // to the current one.
+ for (size_t j = 1; j < m_nd; ++j) {
+ offset =
+ offset * static_cast<int64_t>(m_ws->getDimension(j - 1)->getNBins());
+
+ size_t nEntries = permutationsVertexTouching.size();
+ for (int k = 1; k <= widths[j] / 2; ++k) {
+ for (size_t m = 0; m < nEntries; m++) {
+ permutationsVertexTouching.push_back((offset * k) +
+ permutationsVertexTouching[m]);
+ permutationsVertexTouching.push_back((offset * k * (-1)) +
+ permutationsVertexTouching[m]);
+ }
+ }
+ }
+
+ m_permutationsVertexTouchingMap.insert(std::make_pair(widths, permutationsVertexTouching));
+ }
+
+ // In either case, get the result.
+ 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 {
+
+ 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(widths);
+
+ Utils::NestedForLoop::GetIndicesFromLinearIndex(m_nd, m_pos, m_indexMaker,
+ m_indexMax, m_index);
+
+ // Filter out indexes that are are not actually neighbours.
+ // Accumulate neighbour indexes.
+ std::vector<size_t> neighbourIndexes;
+ for (size_t i = 0; i < permutationsVertexTouching.size(); ++i) {
+ if (permutationsVertexTouching[i] == 0) {
+ continue;
+ }
+
+ 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, widths) ) {
+ neighbourIndexes.push_back(neighbour_index);
+ }
+ }
+
+ // Remove duplicates
+ std::sort(neighbourIndexes.begin(), neighbourIndexes.end());
+ neighbourIndexes.erase(
+ std::unique(neighbourIndexes.begin(), neighbourIndexes.end()),
+ neighbourIndexes.end());
+ return neighbourIndexes;
+}
+
+/**
+ *
+ * @return The size of the permutation cache.
+ */
+size_t MDHistoWorkspaceIterator::permutationCacheSize() const
+{
+ return m_permutationsVertexTouchingMap.size();
+}
+
} // namespace Mantid
} // namespace DataObjects
diff --git a/Code/Mantid/Framework/DataObjects/test/MDHistoWorkspaceIteratorTest.h b/Code/Mantid/Framework/DataObjects/test/MDHistoWorkspaceIteratorTest.h
index d8949c9a26ad8f969bbb9a509add03679c30c678..e39b618a81543c6159385c992ad147b8fcfead09 100644
--- a/Code/Mantid/Framework/DataObjects/test/MDHistoWorkspaceIteratorTest.h
+++ b/Code/Mantid/Framework/DataObjects/test/MDHistoWorkspaceIteratorTest.h
@@ -15,6 +15,7 @@
#include "MantidGeometry/MDGeometry/MDPlane.h"
#include <boost/assign/list_of.hpp>
#include <boost/function.hpp>
+#include <boost/bind.hpp>
using namespace Mantid;
using namespace Mantid::DataObjects;
@@ -279,7 +280,9 @@ public:
histoIt->getLinearIndex());
}
- bool doesContainIndex(const std::vector<size_t>& container, const size_t element)
+ //template<typename ContainerType, typename ElementType>
+ template<class ContainerType>
+ bool doesContainIndex(const ContainerType& container, const typename ContainerType::value_type element)
{
return std::find(container.begin(), container.end(), element) != container.end();
}
@@ -298,6 +301,7 @@ public:
TS_ASSERT(!iterator.isWithinBounds(end));
}
+
void do_test_neighbours_1d(
boost::function<std::vector<size_t>(MDHistoWorkspaceIterator*)> findNeighbourMemberFunction)
{
@@ -321,14 +325,14 @@ public:
std::vector<size_t> neighbourIndexes = findNeighbourMemberFunction(it);
TS_ASSERT_EQUALS(1, neighbourIndexes.size());
// should be on edge
- TSM_ASSERT( "Neighbour at index 0 is 1", doesContainIndex(neighbourIndexes, 1));
+ TSM_ASSERT( "Neighbour at index 0 is 1", doesContainIndex(neighbourIndexes, size_t(1)));
// Go to intermediate position
/*
0 - 1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9
- ^
- |
- */
+ ^
+ |
+ */
it->next();
neighbourIndexes = findNeighbourMemberFunction(it);
TS_ASSERT_EQUALS(2, neighbourIndexes.size());
@@ -339,9 +343,9 @@ public:
// Go to last position
/*
0 - 1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9
- ^
- |
- */
+ ^
+ |
+ */
it->jumpTo(9);
neighbourIndexes = findNeighbourMemberFunction(it);
TSM_ASSERT( "Neighbour at index 9 is 8", doesContainIndex(neighbourIndexes, 8));
@@ -671,6 +675,294 @@ public:
}
+ void test_neighbours_1d_with_width()
+ {
+
+ // This is the width to use
+ const int width = 5;
+
+ const size_t nd = 1;
+ MDHistoWorkspace_sptr ws = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.0, nd, 10);
+ /*
+ 1D MDHistoWorkspace
+
+ 0 - 1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9
+
+ */
+
+ MDHistoWorkspaceIterator * it = new MDHistoWorkspaceIterator(ws);
+
+ // At first position
+ /*
+ 0 - 1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9
+ ^
+ |
+ */
+
+ std::vector<size_t> neighbourIndexes = it->findNeighbourIndexesByWidth(width);
+ TS_ASSERT_EQUALS(2, neighbourIndexes.size());
+ // should be on edge
+ TSM_ASSERT( "Neighbours at index 0 includes 1", doesContainIndex(neighbourIndexes, 1));
+ TSM_ASSERT( "Neighbours at index 0 includes 2", doesContainIndex(neighbourIndexes, 1));
+
+ // Go to intermediate position
+ /*
+ 0 - 1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9
+ ^
+ |
+ */
+ it->next();
+ neighbourIndexes = it->findNeighbourIndexesByWidth(width);
+ TS_ASSERT_EQUALS(3, neighbourIndexes.size());
+ // should be on edge
+ TSM_ASSERT( "Neighbours at index 1 includes 0", doesContainIndex(neighbourIndexes, 0));
+ TSM_ASSERT( "Neighbours at index 1 includes 2", doesContainIndex(neighbourIndexes, 2));
+ TSM_ASSERT( "Neighbours at index 1 includes 3", doesContainIndex(neighbourIndexes, 3));
+
+ // Go to last position
+ /*
+ 0 - 1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9
+ ^
+ |
+ */
+ it->jumpTo(9);
+ neighbourIndexes = it->findNeighbourIndexesByWidth(width);
+ TS_ASSERT_EQUALS(2, neighbourIndexes.size());
+ TSM_ASSERT( "Neighbours at index 9 includes 8", doesContainIndex(neighbourIndexes, 8));
+ TSM_ASSERT( "Neighbours at index 9 includes 7", doesContainIndex(neighbourIndexes, 7));
+ }
+
+ void test_neighbours_2d_vertex_touching_by_width()
+ {
+ const size_t nd = 2;
+ const int width = 5;
+ 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(width);
+ TS_ASSERT_EQUALS(8, neighbourIndexes.size());
+ // Is on an edge
+ TSM_ASSERT( "Neighbour at index 0 is 1", doesContainIndex(neighbourIndexes, 1));
+ TSM_ASSERT( "Neighbour at index 0 is 2", doesContainIndex(neighbourIndexes, 2));
+ TSM_ASSERT( "Neighbour at index 0 is 4", doesContainIndex(neighbourIndexes, 4));
+ 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 8", doesContainIndex(neighbourIndexes, 8));
+ TSM_ASSERT( "Neighbour at index 0 is 9", doesContainIndex(neighbourIndexes, 9));
+ TSM_ASSERT( "Neighbour at index 0 is 10", doesContainIndex(neighbourIndexes, 10));
+
+
+ // 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(width);
+ TS_ASSERT_EQUALS(15, neighbourIndexes.size());
+ // Is on an edge
+ for(int i = 0; i < 16; ++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(width);
+ TS_ASSERT_EQUALS(8, neighbourIndexes.size());
+ // Is on an edge
+ 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;
+ const int width = 5;
+ MDHistoWorkspace_sptr ws = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.0, nd, 4);
+ /*
+ 3D MDHistoWorkspace
+
+ [[[ 0 1 2 3]
+ [ 4 5 6 7]
+ [ 8 9 10 11]
+ [12 13 14 15]]
+
+ [[16 17 18 19]
+ [20 21 22 23]
+ [24 25 26 27]
+ [28 29 30 31]]
+
+ [[32 33 34 35]
+ [36 37 38 39]
+ [40 41 42 43]
+ [44 45 46 47]]
+
+ [[48 49 50 51]
+ [52 53 54 55]
+ [56 57 58 59]
+ [60 61 62 63]]]
+ */
+
+ MDHistoWorkspaceIterator * it = new MDHistoWorkspaceIterator(ws);
+
+ // Start at Index = 0
+ std::vector<size_t> neighbourIndexes = it->findNeighbourIndexesByWidth(width);
+ TS_ASSERT_EQUALS(26, neighbourIndexes.size());
+ // Is on an edge
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 1));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 2));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 4));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 5));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 6));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 8));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 9));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 10));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 16));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 17));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 18));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 20));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 21));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 22));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 24));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 25));
+ TS_ASSERT(doesContainIndex(neighbourIndexes, 26));
+ }
+
+ void test_cache()
+ {
+ const size_t nd = 1;
+ MDHistoWorkspace_sptr ws = MDEventsTestHelper::makeFakeMDHistoWorkspace(1.0, nd, 10);
+ /*
+ 1D MDHistoWorkspace
+
+ 0 - 1 - 2 - 3 - 4 - 5 - 6 - 7 - 8 - 9
+
+ */
+
+ MDHistoWorkspaceIterator * it = new MDHistoWorkspaceIterator(ws);
+ TSM_ASSERT_EQUALS("Empty cache expected", 0, it->permutationCacheSize());
+ it->findNeighbourIndexesByWidth(3);
+ TSM_ASSERT_EQUALS("One cache item expected", 1, it->permutationCacheSize());
+ it->findNeighbourIndexesByWidth(3);
+ TSM_ASSERT_EQUALS("One cache item expected", 1, it->permutationCacheSize()); // Same item, no change to cache
+ it->findNeighbourIndexesByWidth(5);
+ TSM_ASSERT_EQUALS("Two cache entries expected", 2, it->permutationCacheSize());
+ }
+
+
};
class MDHistoWorkspaceIteratorTestPerformance: public CxxTest::TestSuite
@@ -799,6 +1091,15 @@ public:
} while (iterator.next());
}
+ void test_findNeighboursByWidth()
+ {
+ MDHistoWorkspaceIterator iterator(small_ws, new SkipNothing());
+ do
+ {
+ iterator.findNeighbourIndexesByWidth(5);
+ } while (iterator.next());
+ }
+
};
#endif /* MANTID_DATAOBJECTS_MDHISTOWORKSPACEITERATORTEST_H_ */
diff --git a/Code/Mantid/Framework/Kernel/inc/MantidKernel/Utils.h b/Code/Mantid/Framework/Kernel/inc/MantidKernel/Utils.h
index 449f16ed831385307d7204cae96ad528b3b1ac6e..e53ef353320531563f25afcfef32b993c83153e8 100644
--- a/Code/Mantid/Framework/Kernel/inc/MantidKernel/Utils.h
+++ b/Code/Mantid/Framework/Kernel/inc/MantidKernel/Utils.h
@@ -286,13 +286,14 @@ inline void getIndicesFromLinearIndex(const size_t linear_index,
* @param subject_indices
* @param num_bins
* @param index_max
+ * @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 *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,
@@ -302,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 > 1) {
+ if (diff > widths[ind]/2) {
return false;
}
}
diff --git a/Code/Mantid/Framework/MDAlgorithms/CMakeLists.txt b/Code/Mantid/Framework/MDAlgorithms/CMakeLists.txt
index 3e92b53e508838bb7f92d2ea1a2b9f2f41c90df6..1adef4fb0bb0a5e237dcf7fd76138dff4742de8a 100644
--- a/Code/Mantid/Framework/MDAlgorithms/CMakeLists.txt
+++ b/Code/Mantid/Framework/MDAlgorithms/CMakeLists.txt
@@ -104,6 +104,7 @@ set ( SRC_FILES
src/SetMDUsingMask.cpp
src/SliceMD.cpp
src/SlicingAlgorithm.cpp
+ src/SmoothMD.cpp
src/ThresholdMD.cpp
src/TransformMD.cpp
src/UnaryOperationMD.cpp
@@ -219,6 +220,7 @@ set ( INC_FILES
inc/MantidMDAlgorithms/SetMDUsingMask.h
inc/MantidMDAlgorithms/SliceMD.h
inc/MantidMDAlgorithms/SlicingAlgorithm.h
+ inc/MantidMDAlgorithms/SmoothMD.h
inc/MantidMDAlgorithms/ThresholdMD.h
inc/MantidMDAlgorithms/TransformMD.h
inc/MantidMDAlgorithms/UnaryOperationMD.h
@@ -317,6 +319,7 @@ set ( TEST_FILES
SimulateResolutionConvolvedModelTest.h
SliceMDTest.h
SlicingAlgorithmTest.h
+ SmoothMDTest.h
Strontium122Test.h
ThresholdMDTest.h
TobyFitBMatrixTest.h
diff --git a/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/SmoothMD.h b/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/SmoothMD.h
new file mode 100644
index 0000000000000000000000000000000000000000..61c447729a39e417ab6799c08af66a17f952a1cf
--- /dev/null
+++ b/Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/SmoothMD.h
@@ -0,0 +1,66 @@
+#ifndef MANTID_MDALGORITHMS_SMOOTHMD_H_
+#define MANTID_MDALGORITHMS_SMOOTHMD_H_
+
+#include "MantidKernel/System.h"
+#include "MantidAPI/Algorithm.h"
+#include <boost/shared_ptr.hpp>
+#include <boost/optional.hpp>
+
+namespace Mantid
+{
+ namespace API
+ {
+ class IMDHistoWorkspace;
+ }
+namespace MDAlgorithms
+{
+
+ /** SmoothMD : Algorithm for smoothing MDHistoWorkspaces
+
+ Copyright © 2015 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge National Laboratory & European Spallation Source
+
+ This file is part of Mantid.
+
+ Mantid is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ Mantid is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+ File change history is stored at: <https://github.com/mantidproject/mantid>
+ Code Documentation is available at: <http://doxygen.mantidproject.org>
+ */
+ class DLLExport SmoothMD : public API::Algorithm
+ {
+ public:
+ SmoothMD();
+ virtual ~SmoothMD();
+
+ virtual const std::string name() const;
+ virtual int version() const;
+ virtual const std::string category() const;
+ virtual const std::string summary() const;
+ std::map<std::string, std::string> validateInputs();
+
+ boost::shared_ptr<Mantid::API::IMDHistoWorkspace> hatSmooth(boost::shared_ptr<const Mantid::API::IMDHistoWorkspace> toSmooth,
+ const std::vector<int> &widthVector, boost::optional<boost::shared_ptr<const Mantid::API::IMDHistoWorkspace> > weighting);
+
+ private:
+
+ void init();
+ void exec();
+
+ };
+
+
+} // namespace MDAlgorithms
+} // namespace Mantid
+
+#endif /* MANTID_MDALGORITHMS_SMOOTHMD_H_ */
diff --git a/Code/Mantid/Framework/MDAlgorithms/src/SmoothMD.cpp b/Code/Mantid/Framework/MDAlgorithms/src/SmoothMD.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8e46f38e9691b290259c7e4b658b438b5117d0b1
--- /dev/null
+++ b/Code/Mantid/Framework/MDAlgorithms/src/SmoothMD.cpp
@@ -0,0 +1,341 @@
+#include "MantidMDAlgorithms/SmoothMD.h"
+#include "MantidAPI/FrameworkManager.h"
+#include "MantidAPI/IMDHistoWorkspace.h"
+#include "MantidAPI/IMDIterator.h"
+#include "MantidAPI/Progress.h"
+#include "MantidKernel/ArrayProperty.h"
+#include "MantidKernel/ArrayBoundedValidator.h"
+#include "MantidKernel/CompositeValidator.h"
+#include "MantidKernel/ListValidator.h"
+#include "MantidKernel/MandatoryValidator.h"
+#include "MantidKernel/PropertyWithValue.h"
+#include "MantidKernel/MultiThreaded.h"
+#include "MantidDataObjects/MDHistoWorkspaceIterator.h"
+#include <boost/make_shared.hpp>
+#include <vector>
+#include <map>
+#include <string>
+#include <sstream>
+#include <utility>
+#include <limits>
+#include <boost/function.hpp>
+#include <boost/bind.hpp>
+#include <boost/scoped_ptr.hpp>
+#include <boost/tuple/tuple.hpp>
+
+using namespace Mantid::Kernel;
+using namespace Mantid::API;
+using namespace Mantid::DataObjects;
+
+// Typedef for with vector
+typedef std::vector<int> WidthVector;
+
+// Typedef for an optional md histo workspace
+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;
+
+// Typedef for a smoothing function map keyed by name.
+typedef std::map<std::string, SmoothFunction> SmoothFunctionMap;
+
+namespace {
+
+/**
+ * @brief functions
+ * @return Allowed smoothing functions
+ */
+std::vector<std::string> functions() {
+ std::vector<std::string> propOptions;
+ propOptions.push_back("Hat");
+ // propOptions.push_back("Gaussian");
+ return propOptions;
+}
+
+/**
+ * Maps a function name to a smoothing function
+ * @return function map
+ */
+SmoothFunctionMap makeFunctionMap(Mantid::MDAlgorithms::SmoothMD *instance) {
+ SmoothFunctionMap map;
+ map.insert(std::make_pair(
+ "Hat", boost::bind(&Mantid::MDAlgorithms::SmoothMD::hatSmooth, instance,
+ _1, _2, _3)));
+ return map;
+}
+}
+
+namespace Mantid {
+namespace MDAlgorithms {
+
+// Register the algorithm into the AlgorithmFactory
+DECLARE_ALGORITHM(SmoothMD)
+
+//----------------------------------------------------------------------------------------------
+/** Constructor
+ */
+SmoothMD::SmoothMD() {}
+
+//----------------------------------------------------------------------------------------------
+/** Destructor
+ */
+SmoothMD::~SmoothMD() {}
+
+//----------------------------------------------------------------------------------------------
+
+/// Algorithms name for identification. @see Algorithm::name
+const std::string SmoothMD::name() const { return "SmoothMD"; }
+
+/// Algorithm's version for identification. @see Algorithm::version
+int SmoothMD::version() const { return 1; }
+
+/// Algorithm's category for identification. @see Algorithm::category
+const std::string SmoothMD::category() const { return "MDAlgorithms"; }
+
+/// Algorithm's summary for use in the GUI and help. @see Algorithm::summary
+const std::string SmoothMD::summary() const {
+ return "Smooth an MDHistoWorkspace according to a weight function";
+}
+
+/**
+ * Hat function smoothing. All weights even. Hat function boundaries beyond
+ * width.
+ * @param toSmooth : Workspace to smooth
+ * @param widthVector : Width vector
+ * @param weightingWS : Weighting workspace (optional)
+ * @return Smoothed MDHistoWorkspace
+ */
+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));
+ // Create the output workspace.
+ IMDHistoWorkspace_sptr outWS = toSmooth->clone();
+ progress.reportIncrement(
+ size_t(double(nPoints) * 0.1)); // Report ~10% progress
+
+ const int nThreads = Mantid::API::FrameworkManager::Instance()
+ .getNumOMPThreads(); // NThreads to Request
+
+ auto iterators = toSmooth->createIterators(nThreads, NULL);
+
+ PARALLEL_FOR_NO_WSP_CHECK()
+ for (int it = 0; it < int(iterators.size()); ++it) {
+
+ PARALLEL_START_INTERUPT_REGION
+ boost::scoped_ptr<MDHistoWorkspaceIterator> iterator(
+ dynamic_cast<MDHistoWorkspaceIterator *>(iterators[it]));
+
+ do {
+ // Gets all vertex-touching neighbours
+ size_t iteratorIndex = iterator->getLinearIndex();
+
+ if (useWeights) {
+
+ // Check that we could measuer here.
+ if ((*weightingWS)->getSignalAt(iteratorIndex) == 0) {
+
+ outWS->setSignalAt(iteratorIndex,
+ std::numeric_limits<double>::quiet_NaN());
+
+ outWS->setErrorSquaredAt(iteratorIndex,
+ std::numeric_limits<double>::quiet_NaN());
+
+ continue; // Skip we couldn't measure here.
+ }
+ }
+
+ std::vector<size_t> neighbourIndexes =
+ 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;
+ }
+ }
+ sumSignal += toSmooth->getSignalAt(neighbourIndexes[i]);
+ double error = toSmooth->getErrorAt(neighbourIndexes[i]);
+ sumSqError += (error * error);
+ }
+
+ // Calculate the mean
+ outWS->setSignalAt(iteratorIndex, sumSignal / double(nNeighbours + 1));
+ // Calculate the sample variance
+ outWS->setErrorSquaredAt(iteratorIndex,
+ sumSqError / double(nNeighbours + 1));
+
+ progress.report();
+
+ } while (iterator->next());
+ PARALLEL_END_INTERUPT_REGION
+ }
+ PARALLEL_CHECK_INTERUPT_REGION
+
+ return outWS;
+}
+
+//----------------------------------------------------------------------------------------------
+/** Initialize the algorithm's properties.
+ */
+void SmoothMD::init() {
+ declareProperty(new WorkspaceProperty<API::IMDHistoWorkspace>(
+ "InputWorkspace", "", Direction::Input),
+ "An input MDHistoWorkspace to smooth.");
+
+ auto widthVectorValidator = boost::make_shared<CompositeValidator>();
+ auto boundedValidator =
+ boost::make_shared<ArrayBoundedValidator<int>>(1, 100);
+ widthVectorValidator->add(boundedValidator);
+ widthVectorValidator->add(
+ boost::make_shared<MandatoryValidator<std::vector<int>>>());
+
+ declareProperty(new ArrayProperty<int>("WidthVector", widthVectorValidator,
+ Direction::Input),
+ "Width vector. Either specify the width in n-pixels for each "
+ "dimension, or provide a single entry (n-pixels) for all "
+ "dimensions.");
+
+ const auto allFunctionTypes = functions();
+ const std::string first = allFunctionTypes.front();
+
+ std::stringstream docBuffer;
+ docBuffer << "Smoothing function. Defaults to " << first;
+ declareProperty(
+ new PropertyWithValue<std::string>(
+ "Function", first,
+ boost::make_shared<ListValidator<std::string>>(allFunctionTypes),
+ Direction::Input),
+ docBuffer.str());
+
+ declareProperty(new WorkspaceProperty<API::IMDHistoWorkspace>(
+ "InputNormalizationWorkspace", "", Direction::Input,
+ PropertyMode::Optional),
+ "Multidimensional weighting workspace. Optional.");
+
+ declareProperty(new WorkspaceProperty<API::IMDHistoWorkspace>(
+ "OutputWorkspace", "", Direction::Output),
+ "An output smoothed MDHistoWorkspace.");
+}
+
+//----------------------------------------------------------------------------------------------
+/** Execute the algorithm.
+ */
+void SmoothMD::exec() {
+
+ // Get the input workspace to smooth
+ IMDHistoWorkspace_sptr toSmooth = this->getProperty("InputWorkspace");
+
+ // Get the input weighting workspace
+ IMDHistoWorkspace_sptr weightingWS =
+ this->getProperty("InputNormalizationWorkspace");
+ OptionalIMDHistoWorkspace_const_sptr optionalWeightingWS;
+ if (weightingWS) {
+ optionalWeightingWS = weightingWS;
+ }
+
+ // Get the width vector
+ 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");
+ SmoothFunctionMap functionMap = makeFunctionMap(this);
+ SmoothFunction smoothFunction = functionMap[smoothFunctionName];
+ // invoke the smoothing operation
+ auto smoothed = smoothFunction(toSmooth, widthVector, optionalWeightingWS);
+
+ setProperty("OutputWorkspace", smoothed);
+}
+
+/**
+ * validateInputs
+ * @return map of property names to errors.
+ */
+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);
+
+ 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 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
+ << ". 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;
+}
+
+} // namespace MDAlgorithms
+} // namespace Mantid
diff --git a/Code/Mantid/Framework/MDAlgorithms/test/SmoothMDTest.h b/Code/Mantid/Framework/MDAlgorithms/test/SmoothMDTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..db55e4e826616d2197f001639a12268b5d1e36f2
--- /dev/null
+++ b/Code/Mantid/Framework/MDAlgorithms/test/SmoothMDTest.h
@@ -0,0 +1,423 @@
+#ifndef MANTID_MDALGORITHMS_SMOOTHMDTEST_H_
+#define MANTID_MDALGORITHMS_SMOOTHMDTEST_H_
+
+#include <cxxtest/TestSuite.h>
+#include "MantidMDAlgorithms/SmoothMD.h"
+#include "MantidDataObjects/MDHistoWorkspace.h"
+#include "MantidTestHelpers/MDEventsTestHelper.h"
+#include "MantidAPI/IMDHistoWorkspace.h"
+#include <vector>
+#include <boost/math/special_functions/fpclassify.hpp>
+
+using Mantid::MDAlgorithms::SmoothMD;
+using namespace Mantid::API;
+using namespace Mantid::DataObjects;
+
+class SmoothMDTest : public CxxTest::TestSuite {
+public:
+ // This pair of boilerplate methods prevent the suite being created statically
+ // This means the constructor isn't called when running other tests
+ static SmoothMDTest *createSuite() { return new SmoothMDTest(); }
+ static void destroySuite(SmoothMDTest *suite) { delete suite; }
+
+ void test_Init() {
+ SmoothMD alg;
+ TS_ASSERT_THROWS_NOTHING(alg.initialize())
+ }
+
+ void test_function_is_of_right_type() {
+ SmoothMD alg;
+ alg.initialize();
+ TSM_ASSERT_THROWS("Function can only be of known types for SmoothMD",
+ alg.setProperty("Function", "magic_function"),
+ std::invalid_argument &);
+ }
+
+ void test_reject_negative_width_vector_entry() {
+ SmoothMD alg;
+ alg.initialize();
+ TSM_ASSERT_THROWS("N-pixels contains zero",
+ alg.setProperty("WidthVector", std::vector<int>(1, 0)),
+ std::invalid_argument &);
+ }
+
+ void test_mandatory_width_vector_entry() {
+ SmoothMD alg;
+ alg.initialize();
+ TSM_ASSERT_THROWS("Empty WidthVector",
+ alg.setProperty("WidthVector", std::vector<int>()),
+ std::invalid_argument &);
+ }
+
+ void test_width_entry_must_be_odd() {
+ auto toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(
+ 1 /*signal*/, 1 /*numDims*/, 4 /*numBins in each dimension*/);
+
+ SmoothMD alg;
+ alg.setChild(true);
+ alg.initialize();
+ alg.setPropertyValue("OutputWorkspace", "dummy");
+ alg.setProperty("InputWorkspace", toSmooth);
+ alg.setProperty(
+ "WidthVector",
+ std::vector<int>(1, 4)); // Width vector contains even number == 4
+ TSM_ASSERT_THROWS("One bad entry. Should throw.", alg.execute(),
+ std::runtime_error &);
+
+ std::vector<int> widthVector;
+ widthVector.push_back(3); // OK
+ widthVector.push_back(5); // OK
+ widthVector.push_back(2); // Not OK
+
+ alg.setProperty("WidthVector",
+ widthVector); // Width vector contains even number
+ TSM_ASSERT_THROWS("Some good entries, but should still throw",
+ 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*/);
+
+ /*
+ 2D MDHistoWorkspace Input
+
+ 1 - 1 - 1
+ 1 - 1 - 1
+ 1 - 1 - 1
+ */
+
+ SmoothMD alg;
+ alg.setChild(true);
+ alg.initialize();
+ std::vector<int> widthVector(1, 3);
+ alg.setProperty("WidthVector", widthVector);
+ alg.setProperty("InputWorkspace", toSmooth);
+ alg.setPropertyValue("OutputWorkspace", "dummy");
+ alg.execute();
+ IMDHistoWorkspace_sptr out = alg.getProperty("OutputWorkspace");
+
+ /*
+ 2D MDHistoWorkspace Expected
+
+ 1 - 1 - 1
+ 1 - 1 - 1
+ 1 - 1 - 1
+ */
+ for (size_t i = 0; i < out->getNPoints(); ++i) {
+ TS_ASSERT_EQUALS(1, out->getSignalAt(i));
+ TS_ASSERT_EQUALS(1, out->getErrorAt(i));
+ }
+ }
+
+ void test_smooth_hat_function_3_pix_width() {
+ auto toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(
+ 1 /*signal*/, 2 /*numDims*/, 3 /*numBins in each dimension*/);
+ toSmooth->setSignalAt(4, 2.0);
+
+ /*
+ 2D MDHistoWorkspace Input
+
+ 1 - 1 - 1
+ 1 - 2 - 1
+ 1 - 1 - 1
+ */
+
+ SmoothMD alg;
+ alg.setChild(true);
+ alg.initialize();
+ std::vector<int> widthVector(1, 3);
+ alg.setProperty("WidthVector", widthVector);
+ alg.setProperty("InputWorkspace", toSmooth);
+ alg.setPropertyValue("OutputWorkspace", "dummy");
+ alg.execute();
+ IMDHistoWorkspace_sptr out = alg.getProperty("OutputWorkspace");
+
+ /*
+ 2D MDHistoWorkspace Expected
+
+ 5/4 - 7/6 - 5/4
+ 7/6 - 10/9 - 7/6
+ 5/4 - 7/6 - 5/4
+ */
+
+ TS_ASSERT_EQUALS(5.0 / 4, out->getSignalAt(0));
+ TS_ASSERT_EQUALS(7.0 / 6, out->getSignalAt(1));
+ TS_ASSERT_EQUALS(10.0 / 9, out->getSignalAt(4));
+ }
+
+ void test_smooth_hat_function_5_pix_width() {
+ auto toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(
+ 1 /*signal*/, 2 /*numDims*/, 5 /*numBins in each dimension*/);
+ toSmooth->setSignalAt(12, 4.0);
+
+ /*
+ 2D MDHistoWorkspace Input
+
+ 1 - 1 - 1 - 1 - 1
+ 1 - 1 - 1 - 1 - 1
+ 1 - 1 - 4 - 1 - 1
+ 1 - 1 - 1 - 1 - 1
+ 1 - 1 - 1 - 1 - 1
+
+ */
+
+ SmoothMD alg;
+ alg.setChild(true);
+ alg.initialize();
+ std::vector<int> widthVector(1, 5); // Smooth with width == 5
+ alg.setProperty("WidthVector", widthVector);
+ alg.setProperty("InputWorkspace", toSmooth);
+ alg.setPropertyValue("OutputWorkspace", "dummy");
+ alg.execute();
+ IMDHistoWorkspace_sptr out = alg.getProperty("OutputWorkspace");
+
+ /*
+ 2D MDHistoWorkspace Expected
+
+ key:
+ x = 12/9
+ y = 18/15
+ z = 28/25
+ ` = ignore
+
+ x - ` - y - ` - x
+ ` - ` - ` - ` - `
+ y - ` - z - ` - y
+ ` - ` - ` - ` - `
+ x - ` - y - ` - x
+ */
+
+ // Check vertexes
+ double x = 12.0 / 9;
+ TS_ASSERT_EQUALS(x, out->getSignalAt(0));
+ TS_ASSERT_EQUALS(x, out->getSignalAt(4));
+ TS_ASSERT_EQUALS(x, out->getSignalAt(20));
+ TS_ASSERT_EQUALS(x, out->getSignalAt(24));
+
+ // Check edges
+ double y = 18.0 / 15;
+ TS_ASSERT_EQUALS(y, out->getSignalAt(2));
+ TS_ASSERT_EQUALS(y, out->getSignalAt(10));
+ TS_ASSERT_EQUALS(y, out->getSignalAt(14));
+ TS_ASSERT_EQUALS(y, out->getSignalAt(22));
+
+ // Check centre
+ double z = 28.0 / 25;
+ 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 =
+ MDEventsTestHelper::makeFakeMDHistoWorkspace(1.0 /*signal value*/, 1 /*dimensionality*/,
+ 9);
+
+ MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(
+ 1.0 /*signal value*/, 2 /*dimensionality*/, 9); // one dimension larger
+
+ SmoothMD alg;
+ alg.setChild(true);
+ alg.initialize();
+ std::vector<int> widthVector(1, 3); // Smooth with width == 3
+ alg.setProperty("WidthVector", widthVector);
+ alg.setProperty("InputWorkspace", a);
+ alg.setProperty("InputNormalizationWorkspace", b);
+ alg.setPropertyValue("OutputWorkspace", "dummy");
+
+ TSM_ASSERT_THROWS("Input unsmoothed and input Normalisation workspaces must have the same dimensionality",
+ alg.execute(), std::runtime_error &);
+
+ }
+
+ void test_shape_check_of_weight_ws() {
+
+ const size_t nd = 1;
+
+ MDHistoWorkspace_sptr a =
+ MDEventsTestHelper::makeFakeMDHistoWorkspace(1.0 /*signal value*/, nd,
+ 10);
+
+ MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(
+ 1.0 /*signal value*/, nd, 10 + 1); // one bin longer
+
+ SmoothMD alg;
+ alg.setChild(true);
+ alg.initialize();
+ std::vector<int> widthVector(1, 3); // Smooth with width == 3
+ alg.setProperty("WidthVector", widthVector);
+ alg.setProperty("InputWorkspace", a);
+ alg.setProperty("InputNormalizationWorkspace", b);
+ alg.setPropertyValue("OutputWorkspace", "dummy");
+
+ TSM_ASSERT_THROWS("Input unsmoothed and input Normalisation workspaces must have the same shape",
+ alg.execute(), std::runtime_error &);
+
+ }
+
+ void test_smooth_with_normalization_guidance() {
+
+ const size_t nd = 1;
+ MDHistoWorkspace_sptr toSmooth =
+ MDEventsTestHelper::makeFakeMDHistoWorkspace(2.0 /*signal value*/, nd,
+ 10);
+ toSmooth->setSignalAt(7, 3);
+
+ MDHistoWorkspace_sptr normWs = MDEventsTestHelper::makeFakeMDHistoWorkspace(
+ 1.0 /*signal value*/, nd, 10);
+ normWs->setSignalAt(9, 0);
+
+ /*
+ 1D MDHistoWorkspace for normalization
+
+ 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 0
+
+ 1D MDHistoWorkspace for smoothing
+
+ 2 - 2 - 2 - 2 - 2 - 2 - 2 - 3 - 2 - 2
+ */
+
+ SmoothMD alg;
+ alg.setChild(true);
+ alg.initialize();
+ std::vector<int> widthVector(1, 3); // Smooth with width == 3
+ alg.setProperty("WidthVector", widthVector);
+ alg.setProperty("InputWorkspace", toSmooth);
+ alg.setProperty("InputNormalizationWorkspace", normWs);
+ alg.setPropertyValue("OutputWorkspace", "dummy");
+ alg.execute();
+ IMDHistoWorkspace_sptr out = alg.getProperty("OutputWorkspace");
+
+ TSM_ASSERT_EQUALS("Second index should have a smoothed using 2 "
+ "neighbours nothing ignored",
+ (toSmooth->getSignalAt(0) + toSmooth->getSignalAt(1) + toSmooth->getSignalAt(2)) / 3,
+ out->getSignalAt(1));
+
+
+ TSM_ASSERT_EQUALS("Second to last index should have a smoothed using 1 "
+ "neighbour only neighour at 9 should be ignored",
+ (toSmooth->getSignalAt(8) + toSmooth->getSignalAt(7)) / 2,
+ out->getSignalAt(8));
+
+ TSM_ASSERT("Last index should have a smoothed Value of NaN",
+ boost::math::isnan(out->getSignalAt(9)));
+ }
+};
+
+class SmoothMDTestPerformance : public CxxTest::TestSuite {
+private:
+ IMDHistoWorkspace_sptr m_toSmooth;
+
+public:
+ // This pair of boilerplate methods prevent the suite being created statically
+ // This means the constructor isn't called when running other tests
+ static SmoothMDTestPerformance *createSuite() {
+ return new SmoothMDTestPerformance();
+ }
+ static void destroySuite(SmoothMDTestPerformance *suite) { delete suite; }
+
+ SmoothMDTestPerformance() {
+ m_toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(
+ 1 /*signal*/, 2 /*numDims*/, 500 /*numBins in each dimension*/);
+ }
+
+ void test_execute_hat_function() {
+ SmoothMD alg;
+ alg.setChild(true);
+ alg.initialize();
+ std::vector<int> widthVector(1, 5); // Smooth with width == 5
+ alg.setProperty("WidthVector", widthVector);
+ alg.setProperty("InputWorkspace", m_toSmooth);
+ alg.setPropertyValue("OutputWorkspace", "dummy");
+ alg.execute();
+ IMDHistoWorkspace_sptr out = alg.getProperty("OutputWorkspace");
+ TS_ASSERT(out);
+ }
+
+ void test_execute_hat_function_with_normalisation() {
+ SmoothMD alg;
+ alg.setChild(true);
+ alg.initialize();
+ std::vector<int> widthVector(1, 3); // Smooth with width == 3
+ alg.setProperty("WidthVector", widthVector);
+ alg.setProperty("InputWorkspace", m_toSmooth);
+ alg.setProperty("InputNormalizationWorkspace", m_toSmooth);
+ alg.setPropertyValue("OutputWorkspace", "dummy");
+ alg.execute();
+ IMDHistoWorkspace_sptr out = alg.getProperty("OutputWorkspace");
+ TS_ASSERT(out);
+ }
+};
+
+#endif /* MANTID_MDALGORITHMS_SMOOTHMDTEST_H_ */
diff --git a/Code/Mantid/docs/source/algorithms/SmoothMD-v1.rst b/Code/Mantid/docs/source/algorithms/SmoothMD-v1.rst
new file mode 100644
index 0000000000000000000000000000000000000000..57ac9b1775402e29a0e426e15184f763389176f3
--- /dev/null
+++ b/Code/Mantid/docs/source/algorithms/SmoothMD-v1.rst
@@ -0,0 +1,57 @@
+
+.. algorithm::
+
+.. summary::
+
+.. alias::
+
+.. properties::
+
+Description
+-----------
+
+Provides smoothing of :ref:`MDHistoWorkspace <MDHistoWorkspace>` in n-dimensions. The WidthVector relates to the number of pixels to include in the width for each dimension. *WidthVector* **must contain entries that are odd numbers**.
+
+A *InputNormalizationWorkspace* may optionally be provided. Such workspaces must have exactly the same shape as the *InputWorkspace*. Where the signal values from this workspace are zero, the corresponding smoothed value will be NaN. Any un-smoothed values from the *InputWorkspace* corresponding to zero in the *InputNormalizationWorkspace* will be ignored during neighbour calculations, so effectively omitted from the smoothing altogether.
+Note that the NormalizationWorkspace is not changed, and needs to be smoothed as well, using the same parameters and *InputNormalizationWorkspace* as the original data.
+
+.. figure:: /images/PreSmooth.png
+ :alt: PreSmooth.png
+ :width: 400px
+ :align: center
+
+ No smoothing
+
+.. figure:: /images/Smoothed.png
+ :alt: PreSmooth.png
+ :width: 400px
+ :align: center
+
+ Smooth with WidthVector=5
+
+
+Usage
+-----
+
+**Example - SmoothMD**
+
+.. testcode:: SmoothMDExample
+
+ ws = CreateMDWorkspace(Dimensions=2, Extents=[-10,10,-10,10], Names='A,B', Units='U,U')
+ FakeMDEventData(InputWorkspace=ws, PeakParams='100000,-5,0,1')
+ FakeMDEventData(InputWorkspace=ws, PeakParams='100000,5,0,1')
+ histogram = BinMD(InputWorkspace=ws, AlignedDim0='A,-10,10,50', AlignedDim1='B,-10,10,50', OutputExtents='-10,10,-10,10,-10,10', OutputBins='10,10,10')
+ # plotSlice(histogram)
+ smoothed = SmoothMD(InputWorkspace=histogram, WidthVector=5, Function='Hat')
+ # plotSlice(smoothed)
+
+ print 'Smoothed has %i points' % smoothed.getNPoints()
+
+Output:
+
+.. testoutput:: SmoothMDExample
+
+ Smoothed has 2500 points
+
+.. categories::
+
diff --git a/Code/Mantid/docs/source/images/PreSmooth.png b/Code/Mantid/docs/source/images/PreSmooth.png
new file mode 100644
index 0000000000000000000000000000000000000000..fe335f515ed1bea97acd4138906631466c6bfed5
Binary files /dev/null and b/Code/Mantid/docs/source/images/PreSmooth.png differ
diff --git a/Code/Mantid/docs/source/images/Smoothed.png b/Code/Mantid/docs/source/images/Smoothed.png
new file mode 100644
index 0000000000000000000000000000000000000000..c819c6891cdb7b50ef78c64de5e05219712fed0b
Binary files /dev/null and b/Code/Mantid/docs/source/images/Smoothed.png differ