Skip to content
Snippets Groups Projects
Commit ca696540 authored by Owen Arnold's avatar Owen Arnold
Browse files

refs #11393. Normalization workspace option.

parent c66bd0c0
No related merge requests found
...@@ -4,6 +4,7 @@ ...@@ -4,6 +4,7 @@
#include "MantidKernel/System.h" #include "MantidKernel/System.h"
#include "MantidAPI/Algorithm.h" #include "MantidAPI/Algorithm.h"
#include <boost/shared_ptr.hpp> #include <boost/shared_ptr.hpp>
#include <boost/optional.hpp>
namespace Mantid namespace Mantid
{ {
...@@ -49,7 +50,7 @@ namespace MDAlgorithms ...@@ -49,7 +50,7 @@ namespace MDAlgorithms
std::map<std::string, std::string> validateInputs(); std::map<std::string, std::string> validateInputs();
boost::shared_ptr<Mantid::API::IMDHistoWorkspace> hatSmooth(boost::shared_ptr<const Mantid::API::IMDHistoWorkspace> toSmooth, boost::shared_ptr<Mantid::API::IMDHistoWorkspace> hatSmooth(boost::shared_ptr<const Mantid::API::IMDHistoWorkspace> toSmooth,
const std::vector<int> &widthVector); const std::vector<int> &widthVector, boost::optional<boost::shared_ptr<const Mantid::API::IMDHistoWorkspace> > weighting);
private: private:
......
...@@ -17,6 +17,7 @@ ...@@ -17,6 +17,7 @@
#include <string> #include <string>
#include <sstream> #include <sstream>
#include <utility> #include <utility>
#include <limits>
#include <boost/function.hpp> #include <boost/function.hpp>
#include <boost/bind.hpp> #include <boost/bind.hpp>
#include <boost/scoped_ptr.hpp> #include <boost/scoped_ptr.hpp>
...@@ -26,9 +27,17 @@ using namespace Mantid::Kernel; ...@@ -26,9 +27,17 @@ using namespace Mantid::Kernel;
using namespace Mantid::API; using namespace Mantid::API;
using namespace Mantid::MDEvents; using namespace Mantid::MDEvents;
// Typedef for with vector
typedef std::vector<int> WidthVector; 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( typedef boost::function<IMDHistoWorkspace_sptr(
IMDHistoWorkspace_const_sptr, const WidthVector &)> 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; typedef std::map<std::string, SmoothFunction> SmoothFunctionMap;
namespace { namespace {
...@@ -50,7 +59,7 @@ std::vector<std::string> functions() { ...@@ -50,7 +59,7 @@ std::vector<std::string> functions() {
*/ */
SmoothFunctionMap makeFunctionMap(Mantid::MDAlgorithms::SmoothMD * instance) { SmoothFunctionMap makeFunctionMap(Mantid::MDAlgorithms::SmoothMD * instance) {
SmoothFunctionMap map; SmoothFunctionMap map;
map.insert(std::make_pair("Hat", boost::bind( &Mantid::MDAlgorithms::SmoothMD::hatSmooth, instance, _1, _2 ))); map.insert(std::make_pair("Hat", boost::bind( &Mantid::MDAlgorithms::SmoothMD::hatSmooth, instance, _1, _2, _3 )));
return map; return map;
} }
...@@ -93,11 +102,13 @@ const std::string SmoothMD::summary() const { ...@@ -93,11 +102,13 @@ const std::string SmoothMD::summary() const {
* width. * width.
* @param toSmooth : Workspace to smooth * @param toSmooth : Workspace to smooth
* @param widthVector : Width vector * @param widthVector : Width vector
* @param weightingWS : Weighting workspace (optional)
* @return Smoothed MDHistoWorkspace * @return Smoothed MDHistoWorkspace
*/ */
IMDHistoWorkspace_sptr SmoothMD::hatSmooth(IMDHistoWorkspace_const_sptr toSmooth, IMDHistoWorkspace_sptr SmoothMD::hatSmooth(IMDHistoWorkspace_const_sptr toSmooth,
const WidthVector &widthVector) { const WidthVector &widthVector, OptionalIMDHistoWorkspace_const_sptr weightingWS) {
const bool useWeights = weightingWS.is_initialized();
uint64_t nPoints = toSmooth->getNPoints(); 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. // Create the output workspace.
...@@ -118,31 +129,53 @@ IMDHistoWorkspace_sptr SmoothMD::hatSmooth(IMDHistoWorkspace_const_sptr toSmooth ...@@ -118,31 +129,53 @@ IMDHistoWorkspace_sptr SmoothMD::hatSmooth(IMDHistoWorkspace_const_sptr toSmooth
do { do {
// Gets all vertex-touching neighbours // 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 = std::vector<size_t> neighbourIndexes =
iterator->findNeighbourIndexesByWidth( iterator->findNeighbourIndexesByWidth(
widthVector.front()); // TODO we should use the whole width vector widthVector.front()); // TODO we should use the whole width vector
// not just the first element. // not just the first element.
const size_t nNeighbours = neighbourIndexes.size(); size_t nNeighbours = neighbourIndexes.size();
double sumSignal = iterator->getSignal(); double sumSignal = iterator->getSignal();
double sumSqError = iterator->getError(); double sumSqError = iterator->getError();
for (size_t i = 0; i < neighbourIndexes.size(); ++i) { 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]); sumSignal += toSmooth->getSignalAt(neighbourIndexes[i]);
double error = toSmooth->getErrorAt(neighbourIndexes[i]); double error = toSmooth->getErrorAt(neighbourIndexes[i]);
sumSqError += (error * error); sumSqError += (error * error);
} }
// Calculate the mean // Calculate the mean
outWS->setSignalAt(iterator->getLinearIndex(), outWS->setSignalAt(iteratorIndex,
sumSignal / double(nNeighbours + 1)); sumSignal / double(nNeighbours + 1));
// Calculate the sample variance // Calculate the sample variance
outWS->setErrorSquaredAt(iterator->getLinearIndex(), outWS->setErrorSquaredAt(iteratorIndex,
sumSqError / double(nNeighbours + 1)); sumSqError / double(nNeighbours + 1));
progress.report(); progress.report();
} while (iterator->next()); } while (iterator->next());
PARALLEL_END_INTERUPT_REGION PARALLEL_END_INTERUPT_REGION
} }
...@@ -185,6 +218,10 @@ void SmoothMD::init() { ...@@ -185,6 +218,10 @@ void SmoothMD::init() {
Direction::Input), Direction::Input),
docBuffer.str()); docBuffer.str());
declareProperty(new WorkspaceProperty<API::IMDHistoWorkspace>(
"InputNormalizationWorkspace", "", Direction::Input, PropertyMode::Optional),
"Multidimensional weighting workspace. Optional.");
declareProperty(new WorkspaceProperty<API::IMDHistoWorkspace>( declareProperty(new WorkspaceProperty<API::IMDHistoWorkspace>(
"OutputWorkspace", "", Direction::Output), "OutputWorkspace", "", Direction::Output),
"An output smoothed MDHistoWorkspace."); "An output smoothed MDHistoWorkspace.");
...@@ -197,6 +234,15 @@ void SmoothMD::exec() { ...@@ -197,6 +234,15 @@ void SmoothMD::exec() {
// Get the input workspace to smooth // Get the input workspace to smooth
IMDHistoWorkspace_sptr toSmooth = this->getProperty("InputWorkspace"); 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 // Get the width vector
const std::vector<int> widthVector = this->getProperty("WidthVector"); const std::vector<int> widthVector = this->getProperty("WidthVector");
...@@ -205,7 +251,7 @@ void SmoothMD::exec() { ...@@ -205,7 +251,7 @@ void SmoothMD::exec() {
SmoothFunctionMap functionMap = makeFunctionMap(this); SmoothFunctionMap functionMap = makeFunctionMap(this);
SmoothFunction smoothFunction = functionMap[smoothFunctionName]; SmoothFunction smoothFunction = functionMap[smoothFunctionName];
// invoke the smoothing operation // invoke the smoothing operation
auto smoothed = smoothFunction(toSmooth, widthVector); auto smoothed = smoothFunction(toSmooth, widthVector, optionalWeightingWS);
setProperty("OutputWorkspace", smoothed); setProperty("OutputWorkspace", smoothed);
} }
...@@ -218,6 +264,7 @@ std::map<std::string, std::string> SmoothMD::validateInputs() { ...@@ -218,6 +264,7 @@ std::map<std::string, std::string> SmoothMD::validateInputs() {
std::map<std::string, std::string> product; std::map<std::string, std::string> product;
// Check the width vector
const std::string widthVectorPropertyName = "WidthVector"; const std::string widthVectorPropertyName = "WidthVector";
std::vector<int> widthVector = this->getProperty(widthVectorPropertyName); std::vector<int> widthVector = this->getProperty(widthVectorPropertyName);
for (auto it = widthVector.begin(); it != widthVector.end(); ++it) { for (auto it = widthVector.begin(); it != widthVector.end(); ++it) {
...@@ -229,6 +276,43 @@ std::map<std::string, std::string> SmoothMD::validateInputs() { ...@@ -229,6 +276,43 @@ std::map<std::string, std::string> SmoothMD::validateInputs() {
product.insert(std::make_pair(widthVectorPropertyName, message.str())); 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) {
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; return product;
} }
......
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
#include "MantidTestHelpers/MDEventsTestHelper.h" #include "MantidTestHelpers/MDEventsTestHelper.h"
#include "MantidAPI/IMDHistoWorkspace.h" #include "MantidAPI/IMDHistoWorkspace.h"
#include <vector> #include <vector>
#include <boost/math/special_functions/fpclassify.hpp>
using Mantid::MDAlgorithms::SmoothMD; using Mantid::MDAlgorithms::SmoothMD;
using namespace Mantid::API; using namespace Mantid::API;
...@@ -26,196 +26,299 @@ public: ...@@ -26,196 +26,299 @@ public:
} }
void test_function_is_of_right_type() { void test_function_is_of_right_type() {
SmoothMD alg; SmoothMD alg;
alg.initialize(); alg.initialize();
TSM_ASSERT_THROWS("Function can only be of known types for SmoothMD", alg.setProperty("Function", "magic_function"), std::invalid_argument&); 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() void test_reject_negative_width_vector_entry() {
{ SmoothMD alg;
SmoothMD alg; alg.initialize();
alg.initialize(); TSM_ASSERT_THROWS("N-pixels contains zero",
TSM_ASSERT_THROWS("N-pixels contains zero", alg.setProperty("WidthVector", std::vector<int>(1, 0)), std::invalid_argument&); alg.setProperty("WidthVector", std::vector<int>(1, 0)),
std::invalid_argument &);
} }
void test_mandatory_width_vector_entry() void test_mandatory_width_vector_entry() {
{ SmoothMD alg;
SmoothMD alg; alg.initialize();
alg.initialize(); TSM_ASSERT_THROWS("Empty WidthVector",
TSM_ASSERT_THROWS("Empty WidthVector", alg.setProperty("WidthVector", std::vector<int>()), std::invalid_argument&); alg.setProperty("WidthVector", std::vector<int>()),
std::invalid_argument &);
} }
void test_width_entry_must_be_odd() void test_width_entry_must_be_odd() {
{ auto toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(
auto toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(1 /*signal*/, 1 /*numDims*/, 4 /*numBins in each dimension*/); 1 /*signal*/, 1 /*numDims*/, 4 /*numBins in each dimension*/);
SmoothMD alg; SmoothMD alg;
alg.setChild(true); alg.setChild(true);
alg.initialize(); alg.initialize();
alg.setPropertyValue("OutputWorkspace", "dummy"); alg.setPropertyValue("OutputWorkspace", "dummy");
alg.setProperty("InputWorkspace", toSmooth); alg.setProperty("InputWorkspace", toSmooth);
alg.setProperty("WidthVector", std::vector<int>(1, 4)); // Width vector contains even number == 4 alg.setProperty(
TSM_ASSERT_THROWS("One bad entry. Should throw.", alg.execute(), std::runtime_error&); "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_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
*/
std::vector<int> widthVector; SmoothMD alg;
widthVector.push_back(3); // OK alg.setChild(true);
widthVector.push_back(5); // OK alg.initialize();
widthVector.push_back(2); // Not OK 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");
alg.setProperty("WidthVector", widthVector); // Width vector contains even number /*
TSM_ASSERT_THROWS("Some good entries, but should still throw", alg.execute(), std::runtime_error&); 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_simple_smooth_hat_function() void test_smooth_hat_function_3_pix_width() {
{ auto toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(
auto toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(1 /*signal*/, 2 /*numDims*/, 3 /*numBins in each dimension*/); 1 /*signal*/, 2 /*numDims*/, 3 /*numBins in each dimension*/);
toSmooth->setSignalAt(4, 2.0);
/* /*
2D MDHistoWorkspace Input 2D MDHistoWorkspace Input
1 - 1 - 1 1 - 1 - 1
1 - 1 - 1 1 - 2 - 1
1 - 1 - 1 1 - 1 - 1
*/ */
SmoothMD alg; SmoothMD alg;
alg.setChild(true); alg.setChild(true);
alg.initialize(); alg.initialize();
std::vector<int> widthVector(1, 3); std::vector<int> widthVector(1, 3);
alg.setProperty("WidthVector", widthVector); alg.setProperty("WidthVector", widthVector);
alg.setProperty("InputWorkspace", toSmooth); alg.setProperty("InputWorkspace", toSmooth);
alg.setPropertyValue("OutputWorkspace", "dummy"); alg.setPropertyValue("OutputWorkspace", "dummy");
alg.execute(); alg.execute();
IMDHistoWorkspace_sptr out = alg.getProperty("OutputWorkspace"); IMDHistoWorkspace_sptr out = alg.getProperty("OutputWorkspace");
/* /*
2D MDHistoWorkspace Expected 2D MDHistoWorkspace Expected
1 - 1 - 1 5/4 - 7/6 - 5/4
1 - 1 - 1 7/6 - 10/9 - 7/6
1 - 1 - 1 5/4 - 7/6 - 5/4
*/ */
for(size_t i = 0; i < out->getNPoints(); ++i)
{ TS_ASSERT_EQUALS(5.0 / 4, out->getSignalAt(0));
TS_ASSERT_EQUALS(1, out->getSignalAt(i)); TS_ASSERT_EQUALS(7.0 / 6, out->getSignalAt(1));
TS_ASSERT_EQUALS(1, out->getErrorAt(i)); 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_3_pix_width() void test_dimensional_check_of_weight_ws() {
{
auto toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(1 /*signal*/, 2 /*numDims*/, 3 /*numBins in each dimension*/);
toSmooth->setSignalAt(4, 2.0);
/* MDHistoWorkspace_sptr a =
2D MDHistoWorkspace Input MDEventsTestHelper::makeFakeMDHistoWorkspace(1.0 /*signal value*/, 1 /*dimensionality*/,
9);
1 - 1 - 1 MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(
1 - 2 - 1 1.0 /*signal value*/, 2 /*dimensionality*/, 9); // one dimension larger
1 - 1 - 1
*/
SmoothMD alg; SmoothMD alg;
alg.setChild(true); alg.setChild(true);
alg.initialize(); alg.initialize();
std::vector<int> widthVector(1, 3); std::vector<int> widthVector(1, 3); // Smooth with width == 3
alg.setProperty("WidthVector", widthVector); alg.setProperty("WidthVector", widthVector);
alg.setProperty("InputWorkspace", toSmooth); alg.setProperty("InputWorkspace", a);
alg.setProperty("InputNormalizationWorkspace", b);
alg.setPropertyValue("OutputWorkspace", "dummy"); alg.setPropertyValue("OutputWorkspace", "dummy");
alg.execute();
IMDHistoWorkspace_sptr out = alg.getProperty("OutputWorkspace");
/*
2D MDHistoWorkspace Expected
5/4 - 7/6 - 5/4 TSM_ASSERT_THROWS("Input unsmoothed and input Normalisation workspaces must have the same dimensionality",
7/6 - 10/9 - 7/6 alg.execute(), std::runtime_error &);
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() void test_shape_check_of_weight_ws() {
{
auto toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(1 /*signal*/, 2 /*numDims*/, 5 /*numBins in each dimension*/);
toSmooth->setSignalAt(12, 4.0);
/* const size_t nd = 1;
2D MDHistoWorkspace Input
1 - 1 - 1 - 1 - 1 MDHistoWorkspace_sptr a =
1 - 1 - 1 - 1 - 1 MDEventsTestHelper::makeFakeMDHistoWorkspace(1.0 /*signal value*/, nd,
1 - 1 - 4 - 1 - 1 10);
1 - 1 - 1 - 1 - 1
1 - 1 - 1 - 1 - 1
*/ MDHistoWorkspace_sptr b = MDEventsTestHelper::makeFakeMDHistoWorkspace(
1.0 /*signal value*/, nd, 10 + 1); // one bin longer
SmoothMD alg; SmoothMD alg;
alg.setChild(true); alg.setChild(true);
alg.initialize(); alg.initialize();
std::vector<int> widthVector(1, 5); // Smooth with width == 5 std::vector<int> widthVector(1, 3); // Smooth with width == 3
alg.setProperty("WidthVector", widthVector); alg.setProperty("WidthVector", widthVector);
alg.setProperty("InputWorkspace", toSmooth); alg.setProperty("InputWorkspace", a);
alg.setProperty("InputNormalizationWorkspace", b);
alg.setPropertyValue("OutputWorkspace", "dummy"); alg.setPropertyValue("OutputWorkspace", "dummy");
alg.execute();
IMDHistoWorkspace_sptr out = alg.getProperty("OutputWorkspace"); TSM_ASSERT_THROWS("Input unsmoothed and input Normalisation workspaces must have the same shape",
alg.execute(), std::runtime_error &);
/*
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_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 { class SmoothMDTestPerformance : public CxxTest::TestSuite {
private: private:
IMDHistoWorkspace_sptr m_toSmooth;
IMDHistoWorkspace_sptr m_toSmooth;
public: public:
// This pair of boilerplate methods prevent the suite being created statically // This pair of boilerplate methods prevent the suite being created statically
// This means the constructor isn't called when running other tests // This means the constructor isn't called when running other tests
static SmoothMDTestPerformance *createSuite() { return new SmoothMDTestPerformance(); } static SmoothMDTestPerformance *createSuite() {
return new SmoothMDTestPerformance();
}
static void destroySuite(SmoothMDTestPerformance *suite) { delete suite; } static void destroySuite(SmoothMDTestPerformance *suite) { delete suite; }
SmoothMDTestPerformance() SmoothMDTestPerformance() {
{ m_toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(
m_toSmooth = MDEventsTestHelper::makeFakeMDHistoWorkspace(1 /*signal*/, 2 /*numDims*/, 500 /*numBins in each dimension*/); 1 /*signal*/, 2 /*numDims*/, 500 /*numBins in each dimension*/);
} }
void test_execute_hat_function() { void test_execute_hat_function() {
...@@ -230,8 +333,20 @@ public: ...@@ -230,8 +333,20 @@ public:
IMDHistoWorkspace_sptr out = alg.getProperty("OutputWorkspace"); IMDHistoWorkspace_sptr out = alg.getProperty("OutputWorkspace");
TS_ASSERT(out); 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_ */ #endif /* MANTID_MDALGORITHMS_SMOOTHMDTEST_H_ */
...@@ -12,6 +12,8 @@ Description ...@@ -12,6 +12,8 @@ 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**. 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.
.. figure:: /images/PreSmooth.png .. figure:: /images/PreSmooth.png
:alt: PreSmooth.png :alt: PreSmooth.png
:width: 400px :width: 400px
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment