Commit ba992b00 authored by Janik Zikovsky's avatar Janik Zikovsky
Browse files

Refs #3122: Methods for finding the centroid of a sphere. Algorithm started.

parent 4ca22651
......@@ -8,6 +8,7 @@ set ( SRC_FILES
src/IMDBox.cpp
src/MDBin.cpp
src/MDBox.cpp
src/MDCentroidPeaks.cpp
src/MDEWPeakIntegration.cpp
src/MDEventFactory.cpp
src/MDEventWorkspace.cpp
......@@ -38,6 +39,7 @@ set ( INC_FILES
inc/MantidMDEvents/IMDBox.h
inc/MantidMDEvents/MDBin.h
inc/MantidMDEvents/MDBox.h
inc/MantidMDEvents/MDCentroidPeaks.h
inc/MantidMDEvents/MDDimensionExtents.h
inc/MantidMDEvents/MDDimensionStats.h
inc/MantidMDEvents/MDEWPeakIntegration.h
......@@ -61,6 +63,7 @@ set ( TEST_FILES
test/IMDBoxTest.h
test/MDBinTest.h
test/MDBoxTest.h
test/MDCentroidPeaksTest.h
test/MDDimensionStatsTest.h
test/MDEWPeakIntegrationTest.h
test/MDEventFactoryTest.h
......
......@@ -10,7 +10,11 @@
#include "MantidMDEvents/MDEvent.h"
#include <iosfwd>
/// Define to compare performance when tracking the centroid of events when adding (if true) or only in RefreshCache (if false)
/** Define to compare performance when tracking the centroid of events when adding (if true)
* or only in RefreshCache (if false).
* An additional difference is that the centroid saved in the IMDBox object
* will be normalized by the total signal if false.
*/
#define MDBOX_TRACKCENTROID_WHENADDING 1
namespace Mantid
......@@ -89,6 +93,9 @@ namespace MDEvents
/** Sphere (peak) integration */
virtual void integrateSphere(CoordTransform & radiusTransform, const coord_t radiusSquared, signal_t & signal, signal_t & errorSquared) const = 0;
/** Find the centroid around a sphere */
virtual void centroidSphere(CoordTransform & radiusTransform, const coord_t radiusSquared, coord_t * centroid, signal_t & signal) const = 0;
// -------------------------------------------------------------------------------------------
/** Split sub-boxes, if this is possible and neede for this box */
virtual void splitAllIfNeeded(Mantid::Kernel::ThreadScheduler * /*ts*/ = NULL)
......
......@@ -68,6 +68,8 @@ namespace MDEvents
void integrateSphere(CoordTransform & radiusTransform, const coord_t radiusSquared, signal_t & signal, signal_t & errorSquared) const;
void centroidSphere(CoordTransform & radiusTransform, const coord_t radiusSquared, coord_t * centroid, signal_t & signal) const;
void refreshCache(Kernel::ThreadScheduler * /*ts*/ = NULL);
protected:
......
#ifndef MANTID_MDEVENTS_MDCENTROIDPEAKS_H_
#define MANTID_MDEVENTS_MDCENTROIDPEAKS_H_
#include "MantidKernel/System.h"
#include "MantidAPI/Algorithm.h"
namespace Mantid
{
namespace MDEvents
{
/** Find the centroid of single-crystal peaks in a MDEventWorkspace, in order to refine their positions.
*
* @author Janik Zikovsky
* @date 2011-06-01
*/
class DLLExport MDCentroidPeaks : public API::Algorithm
{
public:
MDCentroidPeaks();
~MDCentroidPeaks();
/// Algorithm's name for identification
virtual const std::string name() const { return "MDCentroidPeaks";};
/// Algorithm's version for identification
virtual int version() const { return 1;};
/// Algorithm's category for identification
virtual const std::string category() const { return "General";}
private:
/// Sets documentation strings for this algorithm
virtual void initDocs();
/// Initialise the properties
void init();
/// Run the algorithm
void exec();
};
} // namespace Mantid
} // namespace MDEvents
#endif /* MANTID_MDEVENTS_MDCENTROIDPEAKS_H_ */
......@@ -56,12 +56,15 @@ namespace MDEvents
void integrateSphere(CoordTransform & radiusTransform, const coord_t radiusSquared, signal_t & signal, signal_t & errorSquared) const;
void centroidSphere(CoordTransform & radiusTransform, const coord_t radiusSquared, coord_t * centroid, signal_t & signal) const;
void splitContents(size_t index, Kernel::ThreadScheduler * ts = NULL);
void splitAllIfNeeded(Kernel::ThreadScheduler * ts = NULL);
void refreshCache(Kernel::ThreadScheduler * ts = NULL);
// ======================= Testing/Debugging Methods =================
/** For testing: get 9a copy of) the vector of boxes */
std::vector<IMDBox<MDE, nd>*> getBoxes()
......
......@@ -53,6 +53,9 @@ namespace MDEvents
void integrateSphere(CoordTransform & /*radiusTransform*/, const coord_t /*radiusSquared*/, signal_t & /*signal*/, signal_t & /*errorSquared*/) const
{ throw std::runtime_error("Not implemented."); }
void centroidSphere(CoordTransform & , const coord_t , coord_t * , signal_t & ) const
{ throw std::runtime_error("Not implemented."); }
// --------------------------------------------------------------------------------------------
/// Return which dimension (index) was split
......
......@@ -278,6 +278,36 @@ namespace MDEvents
}
}
/** Find the centroid of all events contained within by doing a weighted average
* of their coordinates.
*
* @param radiusTransform :: nd-to-1 coordinate transformation that converts from these
* dimensions to the distance (squared) from the center of the sphere.
* @param radiusSquared :: radius^2 below which to integrate
* @param[out] centroid :: array of size [nd]; its centroid will be added
* @param[out] signal :: set to the integrated signal
*/
TMDE(
void MDBox)::centroidSphere(CoordTransform & radiusTransform, const coord_t radiusSquared, coord_t * centroid, signal_t & signal) const
{
typename std::vector<MDE>::const_iterator it = data.begin();
typename std::vector<MDE>::const_iterator it_end = data.end();
// For each MDEvent
for (; it != it_end; ++it)
{
coord_t out[nd];
radiusTransform.apply(it->getCenter(), out);
if (out[0] < radiusSquared)
{
double eventSignal = it->getSignal();
signal += eventSignal;
for (size_t d=0; d<nd; d++)
centroid[d] += it->getCenter(d) * eventSignal;
}
}
}
}//namespace MDEvents
......
#include "MantidMDEvents/MDCentroidPeaks.h"
#include "MantidKernel/System.h"
namespace Mantid
{
namespace MDEvents
{
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(MDCentroidPeaks)
using namespace Mantid::Kernel;
using namespace Mantid::API;
//----------------------------------------------------------------------------------------------
/** Constructor
*/
MDCentroidPeaks::MDCentroidPeaks()
{
}
//----------------------------------------------------------------------------------------------
/** Destructor
*/
MDCentroidPeaks::~MDCentroidPeaks()
{
}
//----------------------------------------------------------------------------------------------
/// Sets documentation strings for this algorithm
void MDCentroidPeaks::initDocs()
{
this->setWikiSummary("Find the centroid of single-crystal peaks in a MDEventWorkspace, in order to refine their positions.");
this->setOptionalMessage("Find the centroid of single-crystal peaks in a MDEventWorkspace, in order to refine their positions.");
this->setWikiDescription(
"This algorithm starts with a PeaksWorkspace containing the expected positions of peaks in reciprocal space. "
"It calculates the centroid of the peak by calculating the average of the coordinates of all events within a given "
"radius of the peak, weighted by the weight (signal) of the event."
"\n\n"
"To speed up the calculation, the centroid of the boxes contained within the radius is used (rather than going "
"through all individual events)."
);
}
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void MDCentroidPeaks::init()
{
declareProperty(new WorkspaceProperty<>("InputWorkspace","",Direction::Input), "An input workspace.");
declareProperty(new WorkspaceProperty<>("OutputWorkspace","",Direction::Output), "An output workspace.");
}
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void MDCentroidPeaks::exec()
{
// TODO Auto-generated execute stub
}
} // namespace Mantid
} // namespace MDEvents
......@@ -609,8 +609,6 @@ namespace MDEvents
{
// There is a chance that this part of the box is within integration volume,
// even if no vertex of it is.
IMDBox<MDE, nd> * box = boxes[i];
coord_t boxCenter[nd];
box->getCenter(boxCenter);
......@@ -654,6 +652,43 @@ namespace MDEvents
}
//-----------------------------------------------------------------------------------------------
/** Find the centroid of all events contained within by doing a weighted average
* of their coordinates.
*
* @param radiusTransform :: nd-to-1 coordinate transformation that converts from these
* dimensions to the distance (squared) from the center of the sphere.
* @param radiusSquared :: radius^2 below which to integrate
* @param[out] centroid :: array of size [nd]; its centroid will be added
* @param[out] signal :: set to the integrated signal
*/
TMDE(
void MDGridBox)::centroidSphere(CoordTransform & radiusTransform, const coord_t radiusSquared, coord_t * centroid, signal_t & signal) const
{
for (size_t i=0; i < numBoxes; ++i)
{
// Go through each contained box
IMDBox<MDE, nd> * box = boxes[i];
coord_t boxCenter[nd];
box->getCenter(boxCenter);
// Distance from center to the peak integration center
coord_t out[nd];
radiusTransform.apply(boxCenter, out);
if (out[0] < diagonalSquared*0.72 + radiusSquared)
{
// If the center is closer than the size of the box, then it MIGHT be touching.
// (We multiply by 0.72 (about sqrt(2)) to look for half the diagonal).
// NOTE! Watch out for non-spherical transforms!
// Go down one level to keep centroiding
box->centroidSphere(radiusTransform, radiusSquared, centroid, signal);
}
} // (for each box)
}
}//namespace MDEvents
}//namespace Mantid
......
......@@ -54,6 +54,7 @@ class IMDBoxTester : public IMDBox<MDE,nd>
{}
virtual void integrateSphere(CoordTransform & /*radiusTransform*/, const coord_t /*radiusSquared*/, signal_t & /*signal*/, signal_t & /*errorSquared*/) const {};
virtual void centroidSphere(CoordTransform & /*radiusTransform*/, const coord_t /*radiusSquared*/, coord_t *, signal_t & ) const {};
};
......
......@@ -312,6 +312,52 @@ public:
}
void test_centroidSphere()
{
MDBox<MDEvent<2>,2> b;
MDEvent<2> ev(2.0, 2.0);
ev.setCenter(0, 2.0);
ev.setCenter(1, 3.0);
b.addEvent(ev);
MDEvent<2> ev2(4.0, 4.0);
ev2.setCenter(0, 4.0);
ev2.setCenter(1, 4.0);
b.addEvent(ev2);
// The sphere transformation
bool dimensionsUsed[2] = {true,true};
coord_t center[2] = {0,0};
CoordTransformDistance sphere(2, center, dimensionsUsed);
// Set up the data for the centroid
coord_t centroid[2] = {0,0};
double signal = 0.0;
b.centroidSphere(sphere, 400., centroid, signal);
for (size_t d=0; d<2; d++)
centroid[d] /= signal;
// This should be the weighted centroid
TS_ASSERT_DELTA( signal, 6.000, 0.001);
TS_ASSERT_DELTA( centroid[0], 3.333, 0.001);
TS_ASSERT_DELTA( centroid[1], 3.666, 0.001);
// --- Reset and reduce the radius ------
signal = 0;
for (size_t d=0; d<2; d++)
centroid[d] = 0.0;
b.centroidSphere(sphere, 16., centroid, signal);
for (size_t d=0; d<2; d++)
centroid[d] /= signal;
// Only one event was contained
TS_ASSERT_DELTA( signal, 2.000, 0.001);
TS_ASSERT_DELTA( centroid[0], 2.000, 0.001);
TS_ASSERT_DELTA( centroid[1], 3.000, 0.001);
}
};
......
#ifndef MANTID_MDEVENTS_MDCENTROIDPEAKSTEST_H_
#define MANTID_MDEVENTS_MDCENTROIDPEAKSTEST_H_
#include <cxxtest/TestSuite.h>
#include "MantidKernel/Timer.h"
#include "MantidKernel/System.h"
#include <iostream>
#include <iomanip>
#include "MantidMDEvents/MDCentroidPeaks.h"
using namespace Mantid::MDEvents;
using namespace Mantid::API;
class MDCentroidPeaksTest : public CxxTest::TestSuite
{
public:
void test_Init()
{
MDCentroidPeaks alg;
TS_ASSERT_THROWS_NOTHING( alg.initialize() )
TS_ASSERT( alg.isInitialized() )
}
void xtest_exec()
{
// Name of the output workspace.
std::string outWSName("MDCentroidPeaksTest_OutputWS");
MDCentroidPeaks alg;
TS_ASSERT_THROWS_NOTHING( alg.initialize() )
TS_ASSERT( alg.isInitialized() )
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("REPLACE_PROPERTY_NAME_HERE!!!!", "value") );
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", outWSName) );
TS_ASSERT_THROWS_NOTHING( alg.execute(); );
TS_ASSERT( alg.isExecuted() );
// Retrieve the workspace from data service. TODO: Change to your desired type
Workspace_sptr ws;
TS_ASSERT_THROWS_NOTHING( ws = boost::dynamic_pointer_cast<Workspace>(AnalysisDataService::Instance().retrieve(outWSName)) );
TS_ASSERT(ws);
if (!ws) return;
// TODO: Check the results
// Remove workspace from the data service.
AnalysisDataService::Instance().remove(outWSName);
}
};
#endif /* MANTID_MDEVENTS_MDCENTROIDPEAKSTEST_H_ */
......@@ -848,7 +848,7 @@ public:
do_check_integrateSphere(box, 4.51,4.5, 0.001, 0.0, "Tiny but off the event.");
do_check_integrateSphere(box, 2.0,2.0, 0.49, 0.0, "At a corner but grabbing nothing");
do_check_integrateSphere(box, 4.8,4.5, 0.35, 1.0, "Too small to contain any vertices");
do_check_integrateSphere(box, 5.0,5.0, 1.0, 4.0, "Contains a box completely");
do_check_integrateSphere(box, 5.0,5.0, 1.0, 4.0, "At a corner, containing 4 neighbors");
do_check_integrateSphere(box, 4.5,4.5, 0.9, 1.0, "Contains one box completely");
do_check_integrateSphere(box, 0.5,0.5, 0.9, 1.0, "Contains one box completely, at the edges");
do_check_integrateSphere(box, 9.5,0.5, 0.9, 1.0, "Contains one box completely, at the edges");
......@@ -953,6 +953,76 @@ public:
}
//------------------------------------------------------------------------------------------------
/** For test_integrateSphere
*
* @param box
* @param radius :: radius to integrate
* @param xExpected :: expected centroid
* @param yExpected :: expected centroid
*/
void do_check_centroidSphere(MDGridBox<MDEvent<2>,2> & box, coord_t x, coord_t y,
const coord_t radius,
double numExpected, coord_t xExpected, coord_t yExpected,
std::string message)
{
std::cout << "Centroid of sphere of radius " << radius << " at " << x << "," << y << "------" << message << "--\n";
// The sphere transformation
bool dimensionsUsed[2] = {true,true};
coord_t center[2] = {x,y};
CoordTransformDistance sphere(2, center, dimensionsUsed);
double signal = 0;
coord_t centroid[2] = {0., 0.};
box.centroidSphere(sphere, radius*radius, centroid, signal);
// Normalized
if (signal != 0.0)
{
for (size_t d=0; d<2; d++)
centroid[d] /= signal;
}
TSM_ASSERT_DELTA( message, signal, 1.0*numExpected, 1e-5);
TSM_ASSERT_DELTA( message, centroid[0], xExpected, 1e-5);
TSM_ASSERT_DELTA( message, centroid[1], yExpected, 1e-5);
}
/** Re-used suite of tests */
void test_centroidSphere()
{
// 10x10 sized box
MDGridBox<MDEvent<2>,2> * box_ptr = makeMDGridBox<2>();
feedMDBox<2>(box_ptr, 1);
// Events are at 0.5, 1.5, etc.
MDGridBox<MDEvent<2>,2> & box = *box_ptr;
TS_ASSERT_EQUALS( box.getNPoints(), 10*10);
do_check_centroidSphere(box, 4.5,4.5, 0.5, 1.0, 4.5, 4.5, "Too small to contain any vertices");
do_check_centroidSphere(box, 4.5, 4.5, 0.001, 1.0, 4.5, 4.5, "Tiny but still has an event.");
do_check_centroidSphere(box, 4.51,4.5, 0.001, 0.0, 0.0, 0.0, "Tiny but off the event.");
do_check_centroidSphere(box, 2.0,2.0, 0.49, 0.0, 0.0, 0.0, "At a corner but grabbing nothing");
do_check_centroidSphere(box, 4.8,4.5, 0.35, 1.0, 4.5, 4.5, "Too small to contain any vertices");
do_check_centroidSphere(box, 5.0,5.0, 1.0, 4.0, 5.0, 5.0, "At a corner, containing 4 neighbors");
do_check_centroidSphere(box, 4.5,4.5, 0.9, 1.0, 4.5, 4.5, "Contains one box completely");
do_check_centroidSphere(box, 0.5,0.5, 0.9, 1.0, 0.5, 0.5, "Contains one box completely, at the edges");
do_check_centroidSphere(box, 9.5,0.5, 0.9, 1.0, 9.5, 0.5, "Contains one box completely, at the edges");
do_check_centroidSphere(box, 0.5,9.5, 0.9, 1.0, 0.5, 9.5, "Contains one box completely, at the edges");
do_check_centroidSphere(box, 4.5,9.5, 0.9, 1.0, 4.5, 9.5, "Contains one box completely, at the edges");
do_check_centroidSphere(box, 9.5,9.5, 0.9, 1.0, 9.5, 9.5, "Contains one box completely, at the edges");
do_check_centroidSphere(box, 1.5,1.5, 1.95, 9.0, 1.5, 1.5, "Contains 5 boxes completely, and 4 boxes with a point");
do_check_centroidSphere(box, -1.0,0.5, 1.55, 1.0, 0.5, 0.5, "Off an edge but enough to get an event");
// Now I add an event very near an edge
coord_t center[2] = {0.001, 0.5};
box.addEvent(MDEvent<2>(1.0, 1.0, center));
do_check_integrateSphere(box, -1.0,0.5, 1.01, 1.0, "Off an edge but just barely enough to get an event");
do_check_integrateSphere(box, 0.0,0.5, 0.01, 1.0, "Tiny, but just barely enough to get an event");
}
private:
std::string message;
};
......@@ -968,14 +1038,18 @@ class MDGridBoxTestPerformance : 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 MDGridBoxTestPerformance *createSuite() { return new MDGridBoxTestPerformance(); }
static void destroySuite( MDGridBoxTestPerformance *suite ) { delete suite; }
MDGridBox<MDEvent<3>,3> * box3;
MDGridBox<MDEvent<3>,3> * box3b;
std::vector<MDEvent<3> > events;
void setUp()
MDGridBoxTestPerformance()
{
// Split 5x5x5, 2 deep.
box3 = MDGridBoxTest::makeRecursiveMDGridBox<3>(5,1);
box3b = MDGridBoxTest::makeRecursiveMDGridBox<3>(5,1);
// Make the list of fake events, random dist.
......@@ -998,10 +1072,20 @@ public:
box3b->refreshCache();
}
~MDGridBoxTestPerformance()
{
delete box3b;
}
void setUp()
{
// Make a fresh box.
box3 = MDGridBoxTest::makeRecursiveMDGridBox<3>(5,1);
}
void tearDown()
{
delete box3;
delete box3b;
}
......@@ -1026,70 +1110,118 @@ public:
}
}
/** Smallish sphere in the middle goes partially through lots of boxes */
void test_sphereIntegrate_inTheMiddle()
//-----------------------------------------------------------------------------
/** Do a sphere integration
*
* @param center :: coordinate of the center
* @param radius :: radius
*/
void do_test_sphereIntegrate(coord_t * center, coord_t radius, double expectSignal, double tol)
{
// The sphere transformation
bool dimensionsUsed[3] = {true,true,true};
coord_t center[3] = {2.5, 2.5, 2.5};
CoordTransformDistance sphere(3, center, dimensionsUsed);
// Repeat the integration a lot
double signal, errorSquared;
for (size_t i=0; i < 1000; i++)
{
signal = 0;
errorSquared = 0;
box3b->integrateSphere(sphere, 1.0, signal, errorSquared);
box3b->integrateSphere(sphere, radius*radius, signal, errorSquared);
}
// The expected number of events, given a sphere of radius 1.0
TS_ASSERT_DELTA(signal, (1e6/125)*(4.0*3.14159/3.0), 2000);
TS_ASSERT_DELTA(signal, expectSignal, tol);
TS_ASSERT_DELTA(signal, errorSquared, 1e-3);
}
/** Smallish sphere in the middle goes partially through lots of boxes */
void test_sphereIntegrate_inTheMiddle()
{
coord_t center[3] = {2.5, 2.5, 2.5};
do_test_sphereIntegrate(center, 1.0, (1e6/125)*(4.0*3.14159/3.0), 2000.0);
}
/** Huge sphere containing all within */
void test_sphereIntegrate_inTheMiddle_largeSphere()
{
// The sphere transformation
bool dimensionsUsed[3] = {true,true,true};
coord_t center[3] = {2.5, 2.5, 2.5};
CoordTransformDistance sphere(3, center, dimensionsUsed);
double signal, errorSquared;
for (size_t i=0; i < 1000; i++)
{