diff --git a/Code/Mantid/Framework/API/inc/MantidAPI/IMDHistoWorkspace.h b/Code/Mantid/Framework/API/inc/MantidAPI/IMDHistoWorkspace.h
index e58fa0a91b9c87bc7f9bf5890d0cb114e8e4e71d..5da50570e195a3d5e866a35b85f2470a02a5ecc6 100644
--- a/Code/Mantid/Framework/API/inc/MantidAPI/IMDHistoWorkspace.h
+++ b/Code/Mantid/Framework/API/inc/MantidAPI/IMDHistoWorkspace.h
@@ -78,6 +78,8 @@ namespace API
virtual double & operator [](const size_t & index) = 0;
virtual void setCoordinateSystem(const Mantid::API::SpecialCoordinateSystem coordinateSystem) = 0;
+
+ boost::shared_ptr<IMDHistoWorkspace> cloneShape() const = 0;
protected:
virtual const std::string toString() const;
};
diff --git a/Code/Mantid/Framework/Crystal/inc/MantidCrystal/ConnectedComponentLabeling.h b/Code/Mantid/Framework/Crystal/inc/MantidCrystal/ConnectedComponentLabeling.h
index c483331e55e52cb0c3c68e0906c991a79e6116c1..39d11cf2bbb8457e1dfb6953ba081fdedf7a7fa2 100644
--- a/Code/Mantid/Framework/Crystal/inc/MantidCrystal/ConnectedComponentLabeling.h
+++ b/Code/Mantid/Framework/Crystal/inc/MantidCrystal/ConnectedComponentLabeling.h
@@ -11,6 +11,10 @@
namespace Mantid
{
+ namespace Kernel
+ {
+ class Logger;
+ }
namespace API
{
@@ -101,6 +105,8 @@ namespace Mantid
/// Run multithreaded
const boost::optional<int> m_nThreads;
+ Mantid::Kernel::Logger& m_logger;
+
};
} // namespace Crystal
diff --git a/Code/Mantid/Framework/Crystal/inc/MantidCrystal/IntegratePeaksUsingClusters.h b/Code/Mantid/Framework/Crystal/inc/MantidCrystal/IntegratePeaksUsingClusters.h
index 87f970edd3c91ce532f6ef93f9c344f5fbe555ba..139fd5637ace28ac8eb13bb61a6b537bad67ce62 100644
--- a/Code/Mantid/Framework/Crystal/inc/MantidCrystal/IntegratePeaksUsingClusters.h
+++ b/Code/Mantid/Framework/Crystal/inc/MantidCrystal/IntegratePeaksUsingClusters.h
@@ -3,6 +3,7 @@
#include "MantidKernel/System.h"
#include "MantidAPI/Algorithm.h"
+#include "MantidAPI/IMDWorkspace.h"
namespace Mantid
{
@@ -42,6 +43,7 @@ namespace Crystal
virtual const std::string category() const;
private:
+ Mantid::API::MDNormalization getNormalization();
virtual void initDocs();
void init();
void exec();
@@ -53,4 +55,4 @@ namespace Crystal
} // namespace Crystal
} // namespace Mantid
-#endif /* MANTID_CRYSTAL_INTEGRATEPEAKSUSINGCLUSTERS_H_ */
\ No newline at end of file
+#endif /* MANTID_CRYSTAL_INTEGRATEPEAKSUSINGCLUSTERS_H_ */
diff --git a/Code/Mantid/Framework/Crystal/src/ConnectedComponentLabeling.cpp b/Code/Mantid/Framework/Crystal/src/ConnectedComponentLabeling.cpp
index 7d3bd25e0548ebc3bf0e302d375e576fa94bcf28..abca8d2386923e5b2b21b74778324b904df0901d 100644
--- a/Code/Mantid/Framework/Crystal/src/ConnectedComponentLabeling.cpp
+++ b/Code/Mantid/Framework/Crystal/src/ConnectedComponentLabeling.cpp
@@ -1,4 +1,5 @@
#include "MantidKernel/MultiThreaded.h"
+#include "MantidKernel/Logger.h"
#include "MantidKernel/V3D.h"
#include "MantidAPI/AlgorithmManager.h"
#include "MantidAPI/FrameworkManager.h"
@@ -45,6 +46,27 @@ namespace Mantid
return maxNeighbours;
}
+ /**
+ * Get the maximum number of clusters possible in an MD space. For use as an offset.
+ * @param ws : Workspace to interogate
+ * @param nIterators : number of equally sized iterators working on the workspace
+ * @return Maximum number of clusters.
+ */
+ size_t calculateMaxClusters(IMDHistoWorkspace const * const ws, size_t nIterators)
+ {
+ size_t maxClusters = 1;
+ for(size_t i = 0; i < ws->getNumDims(); ++i)
+ {
+ maxClusters *= ws->getDimension(i)->getNBins()/2;
+ }
+ maxClusters /= nIterators;
+ if(maxClusters == 0)
+ {
+ maxClusters = ws->getNPoints();
+ }
+ return maxClusters;
+ }
+
/**
* Helper non-member to clone the input workspace
* @param inWS: To clone
@@ -113,6 +135,18 @@ namespace Mantid
typedef std::set<DisjointPair> SetDisjointPair;
typedef std::vector<EdgeIndexPair> VecEdgeIndexPair;
+ /**
+ * Free function performing the CCL implementation over a range defined by the iterator.
+ *
+ * @param iterator : Iterator giving access the the image
+ * @param strategy : Strategy for identifying background
+ * @param neighbourElements : Grid of DisjointElements the same size as the image
+ * @param progress : Progress object to update
+ * @param maxNeighbours : Maximum number of neighbours each element may have. Determined by dimensionality.
+ * @param startLabelId : Start label index to increment
+ * @param edgeIndexVec : Vector of edge index pairs. To identify elements across iterator boundaries to resolve later.
+ * @return
+ */
size_t doConnectedComponentLabeling(
IMDIterator* iterator,
BackgroundStrategy * const strategy,
@@ -123,6 +157,7 @@ namespace Mantid
)
{
size_t currentLabelCount = startLabelId;
+ strategy->configureIterator(iterator); // Set up such things as desired Normalization.
do
{
if(!strategy->isBackground(iterator))
@@ -205,7 +240,7 @@ namespace Mantid
* @param nThreads : Optional argument of number of threads to use.
*/
ConnectedComponentLabeling::ConnectedComponentLabeling(const size_t& startId, const boost::optional<int> nThreads)
- : m_startId(startId), m_nThreads(nThreads)
+ : m_startId(startId), m_nThreads(nThreads), m_logger(Kernel::Logger::get("ConnectedComponentLabeling"))
{
if(m_nThreads.is_initialized() && m_nThreads.get() < 0)
{
@@ -274,30 +309,34 @@ namespace Mantid
VecElements neighbourElements(ws->getNPoints());
const size_t maxNeighbours = calculateMaxNeighbours(ws.get());
+
progress.doReport("Identifying clusters");
size_t frequency = reportEvery<size_t>(10000, ws->getNPoints());
- progress.resetNumSteps(frequency, 0.0, 0.5);
+ progress.resetNumSteps(frequency, 0.0, 0.8);
// For each process maintains pair of index from within process bounds to index outside process bounds
const int nThreadsToUse = getNThreads();
+
if(nThreadsToUse > 1)
{
std::vector<API::IMDIterator*> iterators = ws->createIterators(nThreadsToUse);
- const int nthreads = getNThreads();
- std::vector<VecEdgeIndexPair> parallelEdgeVec(nthreads);
+
+ const size_t maxClustersPossible = calculateMaxClusters(ws.get(), nThreadsToUse);
+
+ std::vector<VecEdgeIndexPair> parallelEdgeVec(nThreadsToUse);
- std::vector<std::map<size_t, boost::shared_ptr<Cluster> > > parallelClusterMapVec(nthreads);
+ std::vector<std::map<size_t, boost::shared_ptr<Cluster> > > parallelClusterMapVec(nThreadsToUse);
// ------------- Stage One. Local CCL in parallel.
+ m_logger.debug("Parallel solve local CCL");
PARALLEL_FOR_NO_WSP_CHECK()
- for(int i = 0; i < nthreads; ++i)
+ for(int i = 0; i < nThreadsToUse; ++i)
{
API::IMDIterator* iterator = iterators[i];
- iterator->setNormalization(NoNormalization);
boost::scoped_ptr<BackgroundStrategy> strategy(baseStrategy->clone()); // local strategy
VecEdgeIndexPair& edgeVec = parallelEdgeVec[i]; // local edge indexes
- const size_t startLabel = m_startId + (i * ws->getNPoints()); // Ensure that label ids are totally unique within each parallel unit.
+ const size_t startLabel = m_startId + (i * maxClustersPossible); // Ensure that label ids are totally unique within each parallel unit.
const size_t endLabel = doConnectedComponentLabeling(iterator, strategy.get(), neighbourElements, progress, maxNeighbours, startLabel, edgeVec);
// Create clusters from labels.
@@ -332,6 +371,7 @@ namespace Mantid
}
// Percolate minimum label across boundaries for indexes where there is ambiguity.
+ m_logger.debug("Percolate minimum label across boundaries");
std::vector<boost::shared_ptr<Cluster> > incompleteClusterVec;
incompleteClusterVec.reserve(1000);
std::set<size_t> usedLabels;
@@ -369,7 +409,8 @@ namespace Mantid
// ------------- Stage 3 In parallel, process each incomplete cluster.
progress.doReport("Merging clusters across processors");
- progress.resetNumSteps(incompleteClusterVec.size(), 0.5, 0.75);
+ m_logger.debug("Merging clusters across processors");
+ progress.resetNumSteps(incompleteClusterVec.size(), 0.8, 0.9);
PARALLEL_FOR_NO_WSP_CHECK()
for(int i = 0; i < static_cast<int>(incompleteClusterVec.size()); ++i)
{
@@ -377,13 +418,16 @@ namespace Mantid
cluster->toUniformMinimum(neighbourElements);
progress.report();
}
-
- // Now combine clusters and add the resolved clusters to the clustermap.
+ m_logger.debug("Remove duplicates");
+ m_logger.debug() << incompleteClusterVec.size() << " clusters to reconstruct" << std::endl;
+ // Now combine clusters and add the resolved clusters to the clusterMap.
+ SetIds usedOwningClusterIds;
for(size_t i = 0; i < incompleteClusterVec.size(); ++i)
{
const size_t label = incompleteClusterVec[i]->getLabel();
- if(!does_contain_key(clusterMap, label))
+ if(usedOwningClusterIds.find(label) == usedOwningClusterIds.end())
{
+ usedOwningClusterIds.insert(label);
clusterMap.insert(std::make_pair(label, incompleteClusterVec[i]));
}
else
@@ -409,7 +453,6 @@ namespace Mantid
// Associate the member DisjointElements with a cluster. Involves looping back over iterator.
iterator->jumpTo(0); // Reset
- iterator->setNormalization(NoNormalization); // TODO: check that this is a valid assumption.
std::set<size_t> labelIds;
do
{
diff --git a/Code/Mantid/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp b/Code/Mantid/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp
index 09b701642ee91d4c8d490880eb1b9d0cf2738b23..0936f5ff036b3c86b1113e3194038fd4d489f25a 100644
--- a/Code/Mantid/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp
+++ b/Code/Mantid/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp
@@ -140,6 +140,16 @@ namespace Mantid
declareProperty(
new PropertyWithValue<double>("Threshold", 0, compositeValidator, Direction::Input),
"Threshold signal above which to consider peaks");
+
+ std::vector<std::string> normalizations(3);
+ normalizations[0] = "NoNormalization";
+ normalizations[1] = "VolumeNormalization";
+ normalizations[2] = "NumEventsNormalization";
+
+ declareProperty("Normalization",normalizations[1],
+ Kernel::IValidator_sptr(new Kernel::ListValidator<std::string>(normalizations)),
+ "Normalization to use with Threshold. Defaults to VolumeNormalization to account for different binning.");
+
declareProperty(new WorkspaceProperty<IPeaksWorkspace>("OutputWorkspace", "", Direction::Output),
"An output integrated peaks workspace.");
declareProperty(
@@ -147,6 +157,30 @@ namespace Mantid
"MDHistoWorkspace containing the labeled clusters used by the algorithm.");
}
+
+ /**
+ * Get the normalization. For use with iterators + background strategies.
+ * @return Chosen normalization
+ */
+ MDNormalization IntegratePeaksUsingClusters::getNormalization()
+ {
+ std::string normProp = getPropertyValue("Normalization");
+ Mantid::API::MDNormalization normalization;
+ if (normProp == "NoNormalization")
+ {
+ normalization = NoNormalization;
+ }
+ else if (normProp == "VolumeNormalization")
+ {
+ normalization = VolumeNormalization;
+ }
+ else
+ {
+ normalization = NumEventsNormalization;
+ }
+ return normalization;
+ }
+
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
@@ -178,12 +212,12 @@ namespace Mantid
const double threshold = getProperty("Threshold");
// Make a background strategy for the CCL analysis to use.
- HardThresholdBackground backgroundStrategy(threshold,NoNormalization);
+ HardThresholdBackground backgroundStrategy(threshold, this->getNormalization());
// CCL. Multi-processor version.
ConnectedComponentLabeling analysis;
Progress progress(this, 0, 1, 1);
- // Peform CCL.
+ // Perform CCL.
ClusterTuple clusters = analysis.executeAndFetchClusters(mdWS, &backgroundStrategy, progress);
// Extract the clusters
ConnectedComponentMappingTypes::ClusterMap& clusterMap = clusters.get<1>();
@@ -195,14 +229,15 @@ namespace Mantid
std::map<size_t, size_t> labelsTakenByPeaks;
progress.doReport("Performing Peak Integration");
- progress.resetNumSteps(peakWS->getNumberPeaks(), 0.75, 1);
+ g_log.information("Starting Integration");
+ progress.resetNumSteps(peakWS->getNumberPeaks(), 0.9, 1);
PARALLEL_FOR1(peakWS)
for(int i = 0; i < peakWS->getNumberPeaks(); ++i)
{
PARALLEL_START_INTERUPT_REGION
IPeak& peak = peakWS->getPeak(i);
const V3D& peakCenterInMDFrame = peakTransform->transformPeak(peak);
- const Mantid::signal_t signalValue = outHistoWS->getSignalAtVMD(peakCenterInMDFrame, NoNormalization);
+ const Mantid::signal_t signalValue = outHistoWS->getSignalAtVMD(peakCenterInMDFrame, NoNormalization); // No normalization when extracting label ids!
if(boost::math::isnan(signalValue))
{
g_log.warning() << "Warning: image for integration is off edge of detector for peak " << i << std::endl;
diff --git a/Code/Mantid/Framework/Crystal/test/ConnectedComponentLabelingTest.h b/Code/Mantid/Framework/Crystal/test/ConnectedComponentLabelingTest.h
index 516972f75e36b9f5001cf3697894d1626c98844b..d3ee6c2a55f215b1bf5e9757f5fcae7321d6c806 100644
--- a/Code/Mantid/Framework/Crystal/test/ConnectedComponentLabelingTest.h
+++ b/Code/Mantid/Framework/Crystal/test/ConnectedComponentLabelingTest.h
@@ -106,7 +106,7 @@ public:
MockBackgroundStrategy mockStrategy;
EXPECT_CALL(mockStrategy, isBackground(_)).Times(static_cast<int>(inWS->getNPoints())*2).WillRepeatedly(Return(false));// A filter that passes everything.
-
+ EXPECT_CALL(mockStrategy, configureIterator(_)).Times(1);
size_t labelingId = 1;
int multiThreaded = 1;
ConnectedComponentLabeling ccl(labelingId, multiThreaded);
@@ -128,7 +128,7 @@ public:
MockBackgroundStrategy mockStrategy;
EXPECT_CALL(mockStrategy, isBackground(_)).Times(static_cast<int>(inWS->getNPoints())*2).WillRepeatedly(Return(false));// A filter that passes everything.
-
+ EXPECT_CALL(mockStrategy, configureIterator(_)).Times(1);
size_t labelingId = 2;
int multiThreaded = 1;
ConnectedComponentLabeling ccl(labelingId, multiThreaded);
@@ -151,7 +151,7 @@ public:
IMDHistoWorkspace_sptr inWS = MDEventsTestHelper::makeFakeMDHistoWorkspace(1, 1, 6); // Makes a 1 by 6 md ws with identical signal values.
MockBackgroundStrategy mockStrategy;
-
+ EXPECT_CALL(mockStrategy, configureIterator(_)).Times(1);
/*
* We use the is background strategy to set up two disconected blocks for us.
* */
@@ -192,7 +192,7 @@ public:
IMDHistoWorkspace_sptr inWS = MDEventsTestHelper::makeFakeMDHistoWorkspace(1, 1, 5); // Makes a 1 by 5 md ws with identical signal values.
MockBackgroundStrategy mockStrategy;
-
+ EXPECT_CALL(mockStrategy, configureIterator(_)).Times(1);
/*
* We use the is background strategy to set up three disconected blocks for us.
* */
@@ -230,6 +230,7 @@ public:
IMDHistoWorkspace_sptr inWS = MDEventsTestHelper::makeFakeMDHistoWorkspace(1, 2, 4); // Makes a 4 by 4 grid.
MockBackgroundStrategy mockStrategy;
+ EXPECT_CALL(mockStrategy, configureIterator(_)).Times(1);
EXPECT_CALL(mockStrategy, isBackground(_)).WillRepeatedly(Return(false));// Nothing is treated as background
size_t labelingId = 1;
@@ -249,6 +250,7 @@ public:
IMDHistoWorkspace_sptr inWS = MDEventsTestHelper::makeFakeMDHistoWorkspace(1, 2, 3); // Makes a 3 by 3 grid.
MockBackgroundStrategy mockStrategy;
+ EXPECT_CALL(mockStrategy, configureIterator(_)).Times(1);
/*
* We treat alternate cells as background, which actually should result in a single object. Think of a chequered flag.
@@ -292,6 +294,7 @@ public:
IMDHistoWorkspace_sptr inWS = MDEventsTestHelper::makeFakeMDHistoWorkspace(1, 3, 3); // Makes a 3 by 3 by 3 grid. All populated with a single value.
MockBackgroundStrategy mockStrategy;
+ EXPECT_CALL(mockStrategy, configureIterator(_)).Times(1);
/*
* We treat alternate cells as background, which actually should result in a single object. Think of a chequered flag.
diff --git a/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspace.h b/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspace.h
index 2406ab0d3e25b01f26d40d6548ade6245969f238..d0cd3edb8679160e7ddc3a0b19f96f87ceba0d8c 100644
--- a/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspace.h
+++ b/Code/Mantid/Framework/MDEvents/inc/MantidMDEvents/MDHistoWorkspace.h
@@ -400,6 +400,10 @@ namespace MDEvents
/// Get the size of an element in the HistoWorkspace.
static size_t sizeOfElement();
+
+ /// Virtual constructor.
+ boost::shared_ptr<IMDHistoWorkspace> cloneShape() const;
+
private:
void initVertexesArray();
diff --git a/Code/Mantid/Framework/MDEvents/src/MDHistoWorkspace.cpp b/Code/Mantid/Framework/MDEvents/src/MDHistoWorkspace.cpp
index 8fe513c57d6cf62e49ac856032027c1f84709aab..669866466d960a47b8034087c6932a127831dae4 100644
--- a/Code/Mantid/Framework/MDEvents/src/MDHistoWorkspace.cpp
+++ b/Code/Mantid/Framework/MDEvents/src/MDHistoWorkspace.cpp
@@ -1308,6 +1308,11 @@ namespace MDEvents
return (3 * sizeof(signal_t) ) + sizeof(bool);
}
+ boost::shared_ptr<IMDHistoWorkspace> MDHistoWorkspace::clone() const;
+ {
+ return boost::make_
+ }
+
} // namespace Mantid
} // namespace MDEvents