diff --git a/Code/Mantid/Framework/Crystal/inc/MantidCrystal/ConnectedComponentLabeling.h b/Code/Mantid/Framework/Crystal/inc/MantidCrystal/ConnectedComponentLabeling.h index 0c5141545f4e..f9c7b6634a9c 100644 --- a/Code/Mantid/Framework/Crystal/inc/MantidCrystal/ConnectedComponentLabeling.h +++ b/Code/Mantid/Framework/Crystal/inc/MantidCrystal/ConnectedComponentLabeling.h @@ -37,6 +37,7 @@ namespace Crystal { public: ConnectedComponentLabeling(); + size_t getStartLabelId() const; void startLabelingId(const size_t& id); boost::shared_ptr execute(Mantid::API::IMDHistoWorkspace_sptr ws, BackgroundStrategy * const strategy) const; virtual ~ConnectedComponentLabeling(); diff --git a/Code/Mantid/Framework/Crystal/src/ConnectedComponentLabeling.cpp b/Code/Mantid/Framework/Crystal/src/ConnectedComponentLabeling.cpp index caf78d5a73ee..5698d6f48a56 100644 --- a/Code/Mantid/Framework/Crystal/src/ConnectedComponentLabeling.cpp +++ b/Code/Mantid/Framework/Crystal/src/ConnectedComponentLabeling.cpp @@ -38,6 +38,14 @@ namespace Mantid m_startId = id; } + /** + @return: The start label id. + */ + size_t ConnectedComponentLabeling::getStartLabelId() const + { + return m_startId; + } + //---------------------------------------------------------------------------------------------- /** Destructor */ diff --git a/Code/Mantid/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp b/Code/Mantid/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp index 49268f60df4e..38b85cbf7c7d 100644 --- a/Code/Mantid/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp +++ b/Code/Mantid/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp @@ -14,6 +14,9 @@ Uses connected component analysis to integrate peaks in an PeaksWorkspace over a #include "MantidKernel/Utils.h" #include "MantidDataObjects/PeaksWorkspace.h" #include +#include +#include +#include using namespace Mantid::API; using namespace Mantid::Kernel; @@ -21,6 +24,24 @@ using namespace Mantid::DataObjects; namespace { + typedef boost::tuple SignalErrorSQPair; + typedef std::map LabelIdIntensityMap; + typedef std::map PositionToLabelIdMap; + + class IsNearPeak +{ +public: + IsNearPeak(const V3D& coordinates, const double& thresholdDistance ) : m_coordinates(coordinates), m_thresholdDistance(thresholdDistance) + {} + bool operator()( const PositionToLabelIdMap::value_type& v ) const + { + return v.first.distance(m_coordinates) < m_thresholdDistance; + } +private: + const V3D m_coordinates; + const double m_thresholdDistance; +}; + class PeakBackground : public Mantid::Crystal::HardThresholdBackground { private: @@ -39,6 +60,9 @@ namespace { if(!HardThresholdBackground::isBackground(iterator) ) { + const VMD& center = iterator->getCenter(); + V3D temp(center[0], center[1], center[2]); // This assumes dims 1, 2, and 3 in the workspace correspond to positions. + for(size_t i = 0; i < m_peaksWS->getNumberPeaks(); ++i) { V3D coords; @@ -54,12 +78,11 @@ namespace { coords= m_peaksWS->getPeak(i).getHKL(); } - const VMD& center = iterator->getCenter(); - V3D temp(center[0], center[1], center[2]); if(coords.distance(temp) < m_radiusEstimate) { return false; } + } } @@ -125,18 +148,10 @@ namespace Mantid void IntegratePeaksUsingClusters::init() { declareProperty(new WorkspaceProperty("InputWorkspace","",Direction::Input), "Input md workspace."); - std::vector propOptions; - propOptions.push_back("NoNormalization"); - propOptions.push_back("NumberOfEventsNormalization"); - propOptions.push_back("VolumeNormalization"); - declareProperty("Normalization", "NoNormalization",boost::make_shared(propOptions), - "Normalization to use." - ); - - declareProperty(new WorkspaceProperty("PeaksWorkspace","", Direction::Input),"A PeaksWorkspace containing the peaks to integrate."); - declareProperty(new PropertyWithValue("RadiusEstimate", 0.1, Direction::Input), "Estimate of Peak Radius. This is optional and should not affect the ouput of the integration, but it will remove noise from the resulting MD cluster image."); - declareProperty(new PropertyWithValue("Threshold", 0, Direction::Input), "Threshold"); - //declareProperty(new WorkspaceProperty("OutputWorkspace","",Direction::Output), "An output integrated peaks workspace."); + declareProperty(new WorkspaceProperty("PeaksWorkspace","", Direction::Input),"A PeaksWorkspace containing the peaks to integrate."); + declareProperty(new PropertyWithValue("RadiusEstimate", 0.1, Direction::Input), "Estimate of Peak Radius. Points beyond this radius will not be considered, so caution towards the larger end."); + declareProperty(new PropertyWithValue("Threshold", 0, Direction::Input), "Threshold signal above which to consider peaks"); + declareProperty(new WorkspaceProperty("OutputWorkspace","",Direction::Output), "An output integrated peaks workspace."); declareProperty(new WorkspaceProperty("OutputWorkspaceMD","",Direction::Output), "MDHistoWorkspace containing the clusters."); } @@ -146,31 +161,87 @@ namespace Mantid void IntegratePeaksUsingClusters::exec() { IMDHistoWorkspace_sptr mdWS = getProperty("InputWorkspace"); - PeaksWorkspace_sptr peaksWS = getProperty("PeaksWorkspace"); + PeaksWorkspace_sptr inPeakWS = getProperty("PeaksWorkspace"); + Mantid::DataObjects::PeaksWorkspace_sptr peakWS = getProperty("OutputWorkspace"); + if (peakWS != inPeakWS) + peakWS = inPeakWS->clone(); const SpecialCoordinateSystem mdCoordinates = mdWS->getSpecialCoordinateSystem(); - - - std::string strNormalization = getPropertyValue("Normalization"); - MDNormalization normalization = NoNormalization; - if(strNormalization == "NumberOfEventsNormalization") - { - normalization = NumEventsNormalization; - } - else if(strNormalization == "VolumeNormalization") - { - normalization = VolumeNormalization; - } + // TODO. check special coordinates have been set. If unknown throw. const double threshold = getProperty("Threshold"); const double radiusEstimate = getProperty("RadiusEstimate"); - PeakBackground background(peaksWS, radiusEstimate, threshold, normalization, mdCoordinates); + PeakBackground background(peakWS, radiusEstimate, threshold, NoNormalization, mdCoordinates); //HardThresholdBackground background(threshold, normalization); ConnectedComponentLabeling analysis; IMDHistoWorkspace_sptr clusters = analysis.execute(mdWS, &background); - //setProperty("OutputWorkspace", peaksWS); + /* + Note that the following may be acheived better inside the clustering utility at the same time as the + cluster workspace is populated. + + Accumulate intesity values for each peak cluster and key by label_id + */ + LabelIdIntensityMap labelMap; + PositionToLabelIdMap positionMap; + + for(size_t i = 0; i < clusters->getNPoints(); ++i) + { + const size_t& label_id = static_cast(clusters->getSignalAt(i)); + + const double& signal = mdWS->getSignalAt(i); + double errorSQ = mdWS->getErrorAt(i); + errorSQ *=errorSQ; + if(label_id >= analysis.getStartLabelId()) + { + if(labelMap.find(label_id) != labelMap.end()) + { + SignalErrorSQPair current = labelMap[label_id]; + labelMap[label_id] = SignalErrorSQPair(current.get<0>() + signal, current.get<1>() + errorSQ); + } + else + { + labelMap[label_id] = SignalErrorSQPair(signal, errorSQ); + + const VMD& center = mdWS->getCenter(i); + V3D temp(center[0], center[1], center[2]); + positionMap[temp] = label_id; //Record charcteristic position of the cluster. + } + } + } + + + for(size_t i =0; i < peakWS->getNumberPeaks(); ++i) + { + IPeak& peak = peakWS->getPeak(i); + V3D coords; + if(mdCoordinates==QLab) + { + coords= peakWS->getPeak(i).getQLabFrame(); + } + else if(mdCoordinates==QSample) + { + coords= peakWS->getPeak(i).getQSampleFrame(); + } + else if(mdCoordinates==Mantid::API::HKL) + { + coords= peakWS->getPeak(i).getHKL(); + } + + /* Now find the label corresponding to these coordinates. Use the characteristic coordinates of the coordinates + recorded earlier to do this. A better implemention would be a direct lookup. + */ + IsNearPeak nearPeak(coords, radiusEstimate); + auto iterator = std::find_if(positionMap.begin(), positionMap.end(), nearPeak); + if(iterator != positionMap.end()) + { + peak.setIntensity(labelMap[ iterator->second ].get<0>()); + peak.setSigmaIntensity(labelMap[ iterator->second ].get<1>()); + } + } + + setProperty("OutputWorkspace", peakWS); setProperty("OutputWorkspaceMD", clusters); } diff --git a/Code/Mantid/Framework/Crystal/test/ConnectedComponentLabelingTest.h b/Code/Mantid/Framework/Crystal/test/ConnectedComponentLabelingTest.h index 390d573d7b56..7591ebef52e1 100644 --- a/Code/Mantid/Framework/Crystal/test/ConnectedComponentLabelingTest.h +++ b/Code/Mantid/Framework/Crystal/test/ConnectedComponentLabelingTest.h @@ -80,6 +80,20 @@ class ConnectedComponentLabelingTest: public CxxTest::TestSuite FrameworkManager::Instance(); } + void test_default_start_label_id() + { + ConnectedComponentLabeling ccl; + TSM_ASSERT_EQUALS("Start Label Id should be 1 by default", 1, ccl.getStartLabelId()); + } + + void test_set_get_start_label_id() + { + ConnectedComponentLabeling ccl; + const size_t startLabelId = 10; + ccl.startLabelingId(startLabelId); + TS_ASSERT_EQUALS(startLabelId, ccl.getStartLabelId()) + } + void test_1d_one_node() { IMDHistoWorkspace_sptr inWS = MDEventsTestHelper::makeFakeMDHistoWorkspace(1, 1, 1); // Single node. Simpliest possible test case diff --git a/Code/Mantid/Framework/Crystal/test/IntegratePeaksUsingClustersTest.h b/Code/Mantid/Framework/Crystal/test/IntegratePeaksUsingClustersTest.h index 7febbc1ec977..c9fc569d88bd 100644 --- a/Code/Mantid/Framework/Crystal/test/IntegratePeaksUsingClustersTest.h +++ b/Code/Mantid/Framework/Crystal/test/IntegratePeaksUsingClustersTest.h @@ -4,8 +4,24 @@ #include #include "MantidCrystal/IntegratePeaksUsingClusters.h" +#include "MantidTestHelpers/MDEventsTestHelper.h" +#include "MantidTestHelpers/WorkspaceCreationHelper.h" +#include "MantidTestHelpers/ComponentCreationHelper.h" +#include "MantidAPI/FrameworkManager.h" +#include "MantidAPI/AlgorithmManager.h" +#include "MantidAPI/Workspace.h" +#include "MantidKernel/V3D.h" +#include "MantidDataObjects/PeaksWorkspace.h" -using Mantid::Crystal::IntegratePeaksUsingClusters; +#include +#include + +using namespace Mantid::Crystal; +using namespace Mantid::API; +using namespace Mantid::Kernel; +using namespace Mantid::MDEvents; +using namespace Mantid::Geometry; +using namespace Mantid::DataObjects; class IntegratePeaksUsingClustersTest : public CxxTest::TestSuite { @@ -15,6 +31,10 @@ class IntegratePeaksUsingClustersTest : public CxxTest::TestSuite static IntegratePeaksUsingClustersTest *createSuite() { return new IntegratePeaksUsingClustersTest(); } static void destroySuite( IntegratePeaksUsingClustersTest *suite ) { delete suite; } + IntegratePeaksUsingClustersTest() + { + FrameworkManager::Instance(); + } void test_Init() { @@ -22,38 +42,126 @@ class IntegratePeaksUsingClustersTest : public CxxTest::TestSuite TS_ASSERT_THROWS_NOTHING( alg.initialize() ) TS_ASSERT( alg.isInitialized() ) } - - void test_exec() + + void test_peaks_workspace_mandatory() { - // Name of the output workspace. - std::string outWSName("IntegratePeaksUsingClustersTest_OutputWS"); - - IntegratePeaksUsingClusters 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 = AnalysisDataService::Instance().retrieveWS(outWSName) ); - TS_ASSERT(ws); - if (!ws) return; - - // TODO: Check the results - - // Remove workspace from the data service. - AnalysisDataService::Instance().remove(outWSName); + IMDHistoWorkspace_sptr mdws = MDEventsTestHelper::makeFakeMDHistoWorkspace(1,1); + + IntegratePeaksUsingClusters alg; + alg.setRethrows(true); + alg.initialize(); + alg.setProperty("InputWorkspace", mdws); + alg.setPropertyValue("OutputWorkspaceMD", "out_md"); + alg.setPropertyValue("OutputWorkspace", "out_peaks"); + TSM_ASSERT_THROWS("PeaksWorkspace required", alg.execute(), std::runtime_error&); } - - void test_Something() + + void test_input_md_workspace_mandatory() + { + auto peaksws = WorkspaceCreationHelper::createPeaksWorkspace(); + + IntegratePeaksUsingClusters alg; + alg.setRethrows(true); + alg.initialize(); + alg.setProperty("PeaksWorkspace", peaksws); + alg.setPropertyValue("OutputWorkspaceMD", "out_md"); + alg.setPropertyValue("OutputWorkspace", "out_peaks"); + TSM_ASSERT_THROWS("InputWorkspace required", alg.execute(), std::runtime_error&); + } + + void add_fake_md_peak(Workspace_sptr mdws, const size_t& nEvents, const double& h, const double& k, const double& l, const double& radius) + { + auto fakeMDEventDataAlg = AlgorithmManager::Instance().create("FakeMDEventData"); + fakeMDEventDataAlg->setChild(true); + fakeMDEventDataAlg->initialize(); + fakeMDEventDataAlg->setProperty("InputWorkspace", mdws); + std::stringstream peakstream; + peakstream << nEvents << ", " << h << ", " << k << ", " << l << ", " << radius; + fakeMDEventDataAlg->setPropertyValue("PeakParams", peakstream.str()); + fakeMDEventDataAlg->execute(); + } + + boost::tuple + make_peak_and_md_ws(const std::vector& hklValues, + const double& min, const double& max, const double& peakRadius=1, + const size_t nEventsInPeak=1000, const size_t& nBins=100) { - TSM_ASSERT( "You forgot to write a test!", 0); + Instrument_sptr inst = ComponentCreationHelper::createTestInstrumentRectangular(1, 100, 0.05); + + auto mdworkspaceAlg = AlgorithmManager::Instance().create("CreateMDWorkspace"); + mdworkspaceAlg->setChild(true); + mdworkspaceAlg->initialize(); + mdworkspaceAlg->setProperty("Dimensions", 3); + std::vector extents = boost::assign::list_of(min)(max)(min)(max)(min)(max).convert_to_container >(); + mdworkspaceAlg->setProperty("Extents", extents); + mdworkspaceAlg->setPropertyValue("Names", "H,K,L"); + mdworkspaceAlg->setPropertyValue("Units", "-,-,-"); + mdworkspaceAlg->setPropertyValue("OutputWorkspace", "IntegratePeaksMDTest_MDEWS"); + mdworkspaceAlg->execute(); + Workspace_sptr mdws = mdworkspaceAlg->getProperty("OutputWorkspace"); + + // --- Make a fake PeaksWorkspace --- + PeaksWorkspace_sptr peakWS(new PeaksWorkspace()); + peakWS->setInstrument(inst); + + for(size_t i = 0; iaddPeak(peak); + + add_fake_md_peak(mdws, nEventsInPeak, h, k, l, peakRadius); + } + + auto binMDAlg = AlgorithmManager::Instance().create("BinMD"); + binMDAlg->setChild(true); + binMDAlg->initialize(); + binMDAlg->setProperty("InputWorkspace", mdws); + binMDAlg->setPropertyValue("OutputWorkspace", "output_ws"); + binMDAlg->setProperty("AxisAligned", true); + + std::stringstream dimensionstring; + dimensionstring << "," << min << ", " << max << "," << nBins ; + + binMDAlg->setPropertyValue("AlignedDim0", "H" + dimensionstring.str()); + binMDAlg->setPropertyValue("AlignedDim1", "K" + dimensionstring.str()); + binMDAlg->setPropertyValue("AlignedDim2", "L" + dimensionstring.str()); + binMDAlg->execute(); + + Workspace_sptr temp = binMDAlg->getProperty("OutputWorkspace"); + IMDHistoWorkspace_sptr outMDWS = boost::dynamic_pointer_cast(temp); + return boost::tuple(outMDWS, peakWS); } + void test_integrateSinglePeak() + { + std::vector hklValues; + hklValues.push_back(V3D(1,1,1)); + auto inputWorkspaces = make_peak_and_md_ws(hklValues, -10, 10); + + auto mdWS = inputWorkspaces.get<0>(); + auto peaksWS = inputWorkspaces.get<1>(); + + IntegratePeaksUsingClusters alg; + alg.initialize(); + alg.setChild(true); + alg.setProperty("InputWorkspace", mdWS); + alg.setProperty("PeaksWorkspace", peaksWS); + alg.setProperty("Threshold", 1000); + alg.setProperty("RadiusEstimate", 1.1); + alg.setPropertyValue("OutputWorkspace", "out_ws"); + alg.setPropertyValue("OutputWorkspaceMD", "out_ws_md"); + alg.execute(); + IPeaksWorkspace_sptr out_peaks = alg.getProperty("OutputWorkspace"); + IMDHistoWorkspace_sptr out_clusters = alg.getProperty("OutputWorkspaceMD"); + } + };