refs #9194. First implementation of algorithm.

Lots of work still to do yet. But on some basic tests generates expected component images (MD workspaces). One thing that will need to be solved yet, is what strategy for background detection to use. Also speedups will be required across the board.
mantidproject · Mar 27, 2014 · aebfb6e · aebfb6e
1 parent 2611b7d
commit aebfb6e
Show file tree

Hide file tree

Showing 5 changed files with 302 additions and 15 deletions.
diff --git a/Code/Mantid/Framework/Crystal/CMakeLists.txt b/Code/Mantid/Framework/Crystal/CMakeLists.txt
@@ -21,6 +21,7 @@ set ( SRC_FILES
 	src/IndexPeaks.cpp
 	src/IndexSXPeaks.cpp
 	src/IntegratePeakTimeSlices.cpp
+	src/IntegratePeaksUsingClusters.cpp
 	src/LatticeErrors.cpp
 	src/LoadHKL.cpp
 	src/LoadIsawPeaks.cpp
@@ -83,6 +84,7 @@ set ( INC_FILES
 	inc/MantidCrystal/IndexPeaks.h
 	inc/MantidCrystal/IndexSXPeaks.h
 	inc/MantidCrystal/IntegratePeakTimeSlices.h
+	inc/MantidCrystal/IntegratePeaksUsingClusters.h
 	inc/MantidCrystal/LatticeErrors.h
 	inc/MantidCrystal/LoadHKL.h
 	inc/MantidCrystal/LoadIsawPeaks.h
@@ -142,9 +144,10 @@ set ( TEST_FILES
 	IndexPeaksTest.h
 	IndexSXPeaksTest.h
 	IntegratePeakTimeSlicesTest.h
+	IntegratePeaksUsingClustersTest.h
 	LoadHKLTest.h
 	LoadIsawPeaksTest.h
-        LoadIsawSpectrumTest.h
+	LoadIsawSpectrumTest.h
 	LoadIsawUBTest.h
 	MaskPeaksWorkspaceTest.h
 	NormaliseVanadiumTest.h

diff --git a/Code/Mantid/Framework/Crystal/inc/MantidCrystal/IntegratePeaksUsingClusters.h b/Code/Mantid/Framework/Crystal/inc/MantidCrystal/IntegratePeaksUsingClusters.h
@@ -0,0 +1,56 @@
+#ifndef MANTID_CRYSTAL_INTEGRATEPEAKSUSINGCLUSTERS_H_
+#define MANTID_CRYSTAL_INTEGRATEPEAKSUSINGCLUSTERS_H_
+
+#include "MantidKernel/System.h"
+#include "MantidAPI/Algorithm.h"
+
+namespace Mantid
+{
+namespace Crystal
+{
+
+  /** IntegratePeaksUsingClusters : Uses clustering to integrate peaks.
+    
+    Copyright &copy; 2014 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
+
+    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 IntegratePeaksUsingClusters  : public API::Algorithm
+  {
+  public:
+    IntegratePeaksUsingClusters();
+    virtual ~IntegratePeaksUsingClusters();
+
+    virtual const std::string name() const;
+    virtual int version() const;
+    virtual const std::string category() const;
+
+  private:
+    virtual void initDocs();
+    void init();
+    void exec();
+
+
+  };
+
+
+} // namespace Crystal
+} // namespace Mantid
+
+#endif  /* MANTID_CRYSTAL_INTEGRATEPEAKSUSINGCLUSTERS_H_ */
diff --git a/Code/Mantid/Framework/Crystal/src/DisjointElement.cpp b/Code/Mantid/Framework/Crystal/src/DisjointElement.cpp
@@ -37,11 +37,7 @@ namespace Mantid
     DisjointElement::DisjointElement(const DisjointElement& other) :
         m_parent(other.m_parent), m_rank(other.m_rank), m_id(other.m_id)
     {
-      if (m_rank > 0)
-      {
-        throw std::logic_error(
-            "This is a parent node. Children cannot be copied, leading to possible inconsistencies.");
-      }
+
       // Don't point to copy object as parent if copy object is it's own parent.
       if (other.m_parent == &other)
       {
@@ -59,11 +55,6 @@ namespace Mantid
       if (this != &other)
       {
 
-        if (other.m_rank > 0)
-        {
-          throw std::logic_error(
-              "This is a parent node. Children cannot be copied, leading to possible inconsistencies.");
-        }
 
         m_parent = other.m_parent;
         m_rank = other.m_rank;
@@ -183,10 +174,7 @@ namespace Mantid
     void DisjointElement::unionWith(DisjointElement* other)
     {
 
-      if (this->getId() == other->getId())
-      {
-        throw std::logic_error("Trying to union two elements with the same Id");
-      }
+
 
       if (other->getRoot() != this->getRoot()) // Check sets do not already have the same root before continuing
       {

diff --git a/Code/Mantid/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp b/Code/Mantid/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp
@@ -0,0 +1,180 @@
+/*WIKI*
+Handles integration of arbitary single crystal peaks shapes.
+
+Uses connected component analysis to integrate peaks in an PeaksWorkspace over an MDHistoWorkspace of data.
+*WIKI*/
+
+#include "MantidCrystal/IntegratePeaksUsingClusters.h"
+#include "MantidCrystal/ConnectedComponentLabeling.h"
+#include "MantidCrystal/HardThresholdBackground.h"
+#include "MantidAPI/IMDHistoWorkspace.h"
+#include "MantidAPI/WorkspaceProperty.h"
+#include "MantidAPI/IMDIterator.h"
+#include "MantidKernel/ListValidator.h"
+#include "MantidKernel/Utils.h"
+#include "MantidDataObjects/PeaksWorkspace.h"
+#include <boost/make_shared.hpp>
+
+using namespace Mantid::API;
+using namespace Mantid::Kernel;
+using namespace Mantid::DataObjects;
+
+namespace
+{
+  class PeakBackground : public Mantid::Crystal::HardThresholdBackground
+  {
+  private:
+
+    IPeaksWorkspace_const_sptr m_peaksWS;
+    const double m_radiusEstimate;
+    const SpecialCoordinateSystem m_mdCoordinates;
+
+  public:
+    PeakBackground(IPeaksWorkspace_const_sptr peaksWS, const double& radiusEstimate, const double& thresholdSignal, const Mantid::API::MDNormalization normalisation, const SpecialCoordinateSystem coordinates) 
+      : Mantid::Crystal::HardThresholdBackground(thresholdSignal, normalisation), m_peaksWS(peaksWS), m_radiusEstimate(radiusEstimate), m_mdCoordinates(coordinates)
+    {
+    }
+
+    virtual bool isBackground(Mantid::API::IMDIterator* iterator) const
+    {
+      if(!HardThresholdBackground::isBackground(iterator) )
+      {
+        for(size_t i = 0; i < m_peaksWS->getNumberPeaks(); ++i)
+        {
+          V3D coords;
+          if(m_mdCoordinates==QLab)
+          {
+            coords= m_peaksWS->getPeak(i).getQLabFrame();
+          }
+          else if(m_mdCoordinates==QSample)
+          {
+            coords= m_peaksWS->getPeak(i).getQSampleFrame();
+          }
+          else if(m_mdCoordinates==Mantid::API::HKL)
+          {
+            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;
+          }
+        }
+
+      }
+      return true;
+    }
+
+    void configureIterator(Mantid::API::IMDIterator* const iterator) const
+    {
+    }
+
+    virtual ~PeakBackground()
+    {
+    }
+  };
+}
+
+namespace Mantid
+{
+  namespace Crystal
+  {
+
+    // Register the algorithm into the AlgorithmFactory
+    DECLARE_ALGORITHM(IntegratePeaksUsingClusters)
+
+
+
+    //----------------------------------------------------------------------------------------------
+    /** Constructor
+    */
+    IntegratePeaksUsingClusters::IntegratePeaksUsingClusters()
+    {
+    }
+
+    //----------------------------------------------------------------------------------------------
+    /** Destructor
+    */
+    IntegratePeaksUsingClusters::~IntegratePeaksUsingClusters()
+    {
+    }
+
+
+    //----------------------------------------------------------------------------------------------
+    /// Algorithm's name for identification. @see Algorithm::name
+    const std::string IntegratePeaksUsingClusters::name() const { return "IntegratePeaksUsingClusters";};
+
+    /// Algorithm's version for identification. @see Algorithm::version
+    int IntegratePeaksUsingClusters::version() const { return 1;};
+
+    /// Algorithm's category for identification. @see Algorithm::category
+    const std::string IntegratePeaksUsingClusters::category() const { return "MDAlgorithms";}
+
+    //----------------------------------------------------------------------------------------------
+    /// Sets documentation strings for this algorithm
+    void IntegratePeaksUsingClusters::initDocs()
+    {
+      this->setWikiSummary("Integrate single crystal peaks using connected component analysis");
+      this->setOptionalMessage(this->getWikiSummary());
+    }
+
+    //----------------------------------------------------------------------------------------------
+    /** Initialize the algorithm's properties.
+    */
+    void IntegratePeaksUsingClusters::init()
+    {
+      declareProperty(new WorkspaceProperty<IMDHistoWorkspace>("InputWorkspace","",Direction::Input), "Input md workspace.");
+      std::vector<std::string> propOptions;
+      propOptions.push_back("NoNormalization");
+      propOptions.push_back("NumberOfEventsNormalization");
+      propOptions.push_back("VolumeNormalization");
+      declareProperty("Normalization", "NoNormalization",boost::make_shared<StringListValidator>(propOptions),
+        "Normalization to use."
+        );
+
+      declareProperty(new WorkspaceProperty<PeaksWorkspace>("PeaksWorkspace","", Direction::Input),"A PeaksWorkspace containing the peaks to integrate.");
+      declareProperty(new PropertyWithValue<double>("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<double>("Threshold", 0, Direction::Input), "Threshold");
+      //declareProperty(new WorkspaceProperty<PeaksWorkspace>("OutputWorkspace","",Direction::Output), "An output integrated peaks workspace."); 
+      declareProperty(new WorkspaceProperty<IMDHistoWorkspace>("OutputWorkspaceMD","",Direction::Output), "MDHistoWorkspace containing the clusters.");
+    }
+
+    //----------------------------------------------------------------------------------------------
+    /** Execute the algorithm.
+    */
+    void IntegratePeaksUsingClusters::exec()
+    {
+      IMDHistoWorkspace_sptr mdWS = getProperty("InputWorkspace");
+      PeaksWorkspace_sptr peaksWS = getProperty("PeaksWorkspace");
+
+      const SpecialCoordinateSystem mdCoordinates = mdWS->getSpecialCoordinateSystem();
+
+
+      std::string strNormalization = getPropertyValue("Normalization");
+      MDNormalization normalization = NoNormalization;
+      if(strNormalization == "NumberOfEventsNormalization")
+      {
+        normalization = NumEventsNormalization;
+      }
+      else if(strNormalization == "VolumeNormalization")
+      {
+        normalization = VolumeNormalization;
+      }
+
+      const double threshold = getProperty("Threshold");
+      const double radiusEstimate = getProperty("RadiusEstimate");
+      PeakBackground background(peaksWS, radiusEstimate, threshold, normalization, mdCoordinates);
+      //HardThresholdBackground background(threshold, normalization);
+
+      ConnectedComponentLabeling analysis;
+      IMDHistoWorkspace_sptr clusters = analysis.execute(mdWS, &background);
+
+      //setProperty("OutputWorkspace", peaksWS);
+      setProperty("OutputWorkspaceMD", clusters);
+    }
+
+
+
+  } // namespace Crystal
+} // namespace Mantid
diff --git a/Code/Mantid/Framework/Crystal/test/IntegratePeaksUsingClustersTest.h b/Code/Mantid/Framework/Crystal/test/IntegratePeaksUsingClustersTest.h
@@ -0,0 +1,60 @@
+#ifndef MANTID_CRYSTAL_INTEGRATEPEAKSUSINGCLUSTERSTEST_H_
+#define MANTID_CRYSTAL_INTEGRATEPEAKSUSINGCLUSTERSTEST_H_
+
+#include <cxxtest/TestSuite.h>
+
+#include "MantidCrystal/IntegratePeaksUsingClusters.h"
+
+using Mantid::Crystal::IntegratePeaksUsingClusters;
+
+class IntegratePeaksUsingClustersTest : 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 IntegratePeaksUsingClustersTest *createSuite() { return new IntegratePeaksUsingClustersTest(); }
+  static void destroySuite( IntegratePeaksUsingClustersTest *suite ) { delete suite; }
+
+
+  void test_Init()
+  {
+    IntegratePeaksUsingClusters alg;
+    TS_ASSERT_THROWS_NOTHING( alg.initialize() )
+    TS_ASSERT( alg.isInitialized() )
+  }
+
+  void test_exec()
+  {
+    // 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<Workspace>(outWSName) );
+    TS_ASSERT(ws);
+    if (!ws) return;
+
+    // TODO: Check the results
+
+    // Remove workspace from the data service.
+    AnalysisDataService::Instance().remove(outWSName);
+  }
+
+  void test_Something()
+  {
+    TSM_ASSERT( "You forgot to write a test!", 0);
+  }
+
+
+};
+
+
+#endif /* MANTID_CRYSTAL_INTEGRATEPEAKSUSINGCLUSTERSTEST_H_ */