refs #9194. Cleanup, document and speedups.

Code/Mantid/Framework/Crystal/inc/MantidCrystal/ConnectedComponentLabeling.h

@@ -10,6 +10,12 @@
 
 namespace Mantid
 {
+
+namespace API
+{
+  class Progress;
+}
+
 namespace Crystal
 {
   namespace ConnectedComponentMappingTypes
@@ -58,12 +64,13 @@ namespace Crystal
     void startLabelingId(const size_t& id);
 
     /// Execute and return clusters
-    boost::shared_ptr<Mantid::API::IMDHistoWorkspace> execute(Mantid::API::IMDHistoWorkspace_sptr ws, BackgroundStrategy * const strategy) const;
+    boost::shared_ptr<Mantid::API::IMDHistoWorkspace> execute(Mantid::API::IMDHistoWorkspace_sptr ws, 
+      BackgroundStrategy * const strategy, Mantid::API::Progress& progress) const;
 
     /// Execute and return clusters, as well as maps to integrated label values
     boost::shared_ptr<Mantid::API::IMDHistoWorkspace> executeAndIntegrate(
       Mantid::API::IMDHistoWorkspace_sptr ws, BackgroundStrategy * const strategy, ConnectedComponentMappingTypes::LabelIdIntensityMap& labelMap,
-      ConnectedComponentMappingTypes::PositionToLabelIdMap& positionLabelMap) const;
+      ConnectedComponentMappingTypes::PositionToLabelIdMap& positionLabelMap, Mantid::API::Progress& progress) const;
 
     /// Destructor
     virtual ~ConnectedComponentLabeling();
@@ -74,7 +81,8 @@ namespace Crystal
     void calculateDisjointTree(Mantid::API::IMDHistoWorkspace_sptr ws, 
       BackgroundStrategy * const strategy, std::vector<DisjointElement>& neighbourElements,
       ConnectedComponentMappingTypes::LabelIdIntensityMap& labelMap,
-      ConnectedComponentMappingTypes::PositionToLabelIdMap& positionLabelMap) const;
+      ConnectedComponentMappingTypes::PositionToLabelIdMap& positionLabelMap, 
+      Mantid::API::Progress& progress) const;
 
     /// Start labeling index
     size_t m_startId;

Code/Mantid/Framework/Crystal/inc/MantidCrystal/DisjointElement.h

@@ -57,6 +57,17 @@ namespace Crystal
     DisjointElement(const DisjointElement& other);
     /// Assignment operator.
     DisjointElement& operator=(const DisjointElement& other);
+    /// Less than
+    inline bool operator<(const DisjointElement& other) const 
+    {
+      return m_id < other.getId();
+		}
+    /// Greater than
+    inline bool operator>(const DisjointElement& other) const 
+    {
+      return m_id > other.getId();
+    }
+
   private:
     bool hasParent() const;
     int compress();

Code/Mantid/Framework/Crystal/inc/MantidCrystal/PeakBackground.h

@@ -2,15 +2,21 @@
 #define MANTID_CRYSTAL_PEAKBACKGROUND_H_
 
 #include "MantidKernel/System.h"
+#include "MantidKernel/V3D.h"
 #include "MantidAPI/IPeaksWorkspace.h"
 #include "MantidAPI/IMDIterator.h"
 #include "MantidAPI/IMDWorkspace.h"
 #include "MantidCrystal/HardThresholdBackground.h"
+#include <boost/function.hpp>
 
 
 
 namespace Mantid
 {
+namespace API
+{
+  class IPeak;
+}
 namespace Crystal
 {
 
@@ -48,6 +54,8 @@ namespace Crystal
     double m_radiusEstimate;
     /// MD coordinates to use
     Mantid::API::SpecialCoordinateSystem m_mdCoordinates;
+    /// Pointer to member function used for coordinate determination.
+    boost::function<Mantid::Kernel::V3D(const Mantid::API::IPeak*)> m_coordFunction;
 
   public:
     /// Constructor

Code/Mantid/Framework/Crystal/src/ConnectedComponentLabeling.cpp

@@ -1,9 +1,10 @@
 #include "MantidKernel/MultiThreaded.h"
 #include "MantidKernel/V3D.h"
-#include "MantidAPI/IMDHistoWorkspace.h"
 #include "MantidAPI/AlgorithmManager.h"
 #include "MantidAPI/FrameworkManager.h"
+#include "MantidAPI/IMDHistoWorkspace.h"
 #include "MantidAPI/IMDIterator.h"
+#include "MantidAPI/Progress.h"
 #include "MantidCrystal/ConnectedComponentLabeling.h"
 #include "MantidCrystal/BackgroundStrategy.h"
 #include "MantidCrystal/DisjointElement.h"
@@ -12,6 +13,10 @@
 #include <boost/scoped_ptr.hpp>
 #include <stdexcept>
 #include <set>
+#include <algorithm>
+#include <iterator>
+
+
 
 using namespace Mantid::API;
 using namespace Mantid::Kernel;
@@ -50,6 +55,17 @@ namespace Mantid
         }
         return outWS;
       }
+
+      template<typename T>
+      T reportEvery(const T& maxReports, const T& maxIterations)
+      {
+        T frequency = maxReports;
+        if (maxIterations >= maxReports)
+        {
+          frequency = maxIterations/maxReports;
+        }
+        return frequency;
+      }
     }
 
     //----------------------------------------------------------------------------------------------
@@ -93,18 +109,21 @@ namespace Mantid
     void ConnectedComponentLabeling::calculateDisjointTree(IMDHistoWorkspace_sptr ws, 
       BackgroundStrategy * const strategy, VecElements& neighbourElements,
       LabelIdIntensityMap& labelMap,
-      PositionToLabelIdMap& positionLabelMap
+      PositionToLabelIdMap& positionLabelMap,
+      Progress& progress
       ) const
     {
 
       VecIndexes allNonBackgroundIndexes;
       allNonBackgroundIndexes.reserve(ws->getNPoints());
 
+      progress.doReport("Pre-processing to filter background out");
+      const size_t nPoints = ws->getNPoints();
+      progress.resetNumSteps(10000, 0.0, 0.25);
       if(m_runMultiThreaded)
       {
-
         std::vector<API::IMDIterator*> iterators = ws->createIterators(getNThreads());
-        const int nthreads = static_cast<int>(iterators.size());
+        const int nthreads = getNThreads();
         std::vector<VecIndexes> manyNonBackgroundIndexes(nthreads);
 
         PARALLEL_FOR_NO_WSP_CHECK()
@@ -118,6 +137,7 @@ namespace Mantid
               if(!strategyCopy->isBackground(iterator))
               {
                 nonBackgroundIndexes.push_back( iterator->getLinearIndex() );
+                progress.report();
               }
             }
             while(iterator->next());
@@ -131,23 +151,36 @@ namespace Mantid
       }
       else
       {
+        progress.resetNumSteps(1, 0.0, 0.5);
         API::IMDIterator *iterator = ws->createIterator(NULL);
         do
         {
           if(!strategy->isBackground(iterator))
           {
             allNonBackgroundIndexes.push_back( iterator->getLinearIndex() );
+            progress.report();
           }
+
         }
         while(iterator->next());
       }
 
       // -------- Perform labeling -----------
+      progress.doReport("Perform connected component labeling");
+
       const size_t maxNeighbours = calculateMaxNeighbours(ws.get());
       IMDIterator* iterator = ws->createIterator(NULL);
       size_t currentLabelCount = m_startId;
-      for(size_t ii = 0; ii < allNonBackgroundIndexes.size(); ++ii)
+      const size_t nIndexesToProcess= allNonBackgroundIndexes.size();
+      const size_t maxReports = 100;
+      const size_t frequency = reportEvery(maxReports, nIndexesToProcess);
+      progress.resetNumSteps(100, 0.25, 0.5);
+      for(size_t ii = 0; ii < nIndexesToProcess ; ++ii)
       {
+        if(ii % frequency == 0)
+        {
+          progress.doReport();
+        }
         size_t& currentIndex = allNonBackgroundIndexes[ii];
         iterator->jumpTo(currentIndex);
 
@@ -184,17 +217,10 @@ namespace Mantid
         else
         {
           // Choose the lowest neighbour index as the parent.
-          size_t parentIndex = nonEmptyNeighbourIndexes[0];
-          for (size_t i = 1; i < nonEmptyNeighbourIndexes.size(); ++i)
-          {
-            size_t neighIndex = nonEmptyNeighbourIndexes[i];
-            if (neighbourElements[neighIndex].getId() < neighbourElements[parentIndex].getId())
-            {
-              parentIndex = i;
-            }
-          }
+          const VecElements::iterator& minIt = std::min_element(neighbourElements.begin(), neighbourElements.end());
+          const size_t& parentIndex = std::distance(neighbourElements.begin(), minIt);
           // Get the chosen parent
-          DisjointElement& parentElement = neighbourElements[parentIndex];
+          DisjointElement& parentElement = *minIt;
           // Make this element a copy of the parent
           neighbourElements[currentIndex] = parentElement;
           // Union remainder parents with the chosen parent
@@ -213,51 +239,60 @@ namespace Mantid
     }
 
     boost::shared_ptr<Mantid::API::IMDHistoWorkspace> ConnectedComponentLabeling::execute(
-      IMDHistoWorkspace_sptr ws, BackgroundStrategy * const strategy) const
+      IMDHistoWorkspace_sptr ws, BackgroundStrategy * const strategy, Progress& progress) const
     {
       VecElements neighbourElements(ws->getNPoints());
 
       // Perform the bulk of the connected component analysis, but don't collapse the elements yet.
       LabelIdIntensityMap labelMap; // This will not get used.
       PositionToLabelIdMap positionLabelMap; // This will not get used.
-      calculateDisjointTree(ws, strategy, neighbourElements, labelMap, positionLabelMap);
+      calculateDisjointTree(ws, strategy, neighbourElements, labelMap, positionLabelMap, progress);
 
       // Create the output workspace from the input workspace
       IMDHistoWorkspace_sptr outWS = cloneInputWorkspace(ws);
 
+      progress.doReport("Generating cluster image");
+      const int nIndexesToProcess = static_cast<int>(neighbourElements.size());
+      progress.resetNumSteps(nIndexesToProcess, 0.5, 0.75);
       // Set each pixel to the root of each disjointed element.
       PARALLEL_FOR_NO_WSP_CHECK()
-      for (int i = 0; i < static_cast<int>(neighbourElements.size()); ++i)
+      for (int i = 0; i < nIndexesToProcess; ++i)
       {
-        //std::cout << "Element\t" << i << " Id: \t" << neighbourElements[i].getId() << " This location:\t"<< &neighbourElements[i] << " Root location:\t" << neighbourElements[i].getParent() << " Root Id:\t" <<  neighbourElements[i].getRoot() << std::endl;
         if(!neighbourElements[i].isEmpty())
         {
           outWS->setSignalAt(i, neighbourElements[i].getRoot());
+          progress.doReport();
         }
         else
         {
           outWS->setSignalAt(i, 0);
         }
         outWS->setErrorSquaredAt(i, 0);
+
       }
 
       return outWS;
     }
 
     boost::shared_ptr<Mantid::API::IMDHistoWorkspace> ConnectedComponentLabeling::executeAndIntegrate(
       IMDHistoWorkspace_sptr ws, BackgroundStrategy * const strategy, LabelIdIntensityMap& labelMap,
-      PositionToLabelIdMap& positionLabelMap) const
+      PositionToLabelIdMap& positionLabelMap, Progress& progress) const
     {
       VecElements neighbourElements(ws->getNPoints());
 
       // Perform the bulk of the connected component analysis, but don't collapse the elements yet.
-      calculateDisjointTree(ws, strategy, neighbourElements, labelMap, positionLabelMap);
+      calculateDisjointTree(ws, strategy, neighbourElements, labelMap, positionLabelMap, progress);
 
       // Create the output workspace from the input workspace
       IMDHistoWorkspace_sptr outWS = cloneInputWorkspace(ws);
 
+      progress.doReport("Integrating clusters and generating cluster image");
+      const size_t nIterations = neighbourElements.size();
+      const size_t maxReports = 100;
+      const size_t frequency = reportEvery(maxReports, nIterations);
+      progress.resetNumSteps(maxReports, 0.5, 0.75);
       // Set each pixel to the root of each disjointed element.
-      for (size_t i = 0; i < neighbourElements.size(); ++i)
+      for (size_t i = 0; i < nIterations; ++i)
       {
         if(!neighbourElements[i].isEmpty())
         {
@@ -266,7 +301,7 @@ namespace Mantid
           errorSQ *=errorSQ; // Error squared at index
           const size_t& labelId = neighbourElements[i].getRoot();
           // Set the output cluster workspace signal value
-          outWS->setSignalAt(i, labelId);
+          outWS->setSignalAt(i, static_cast<Mantid::signal_t>(labelId));
 
           SignalErrorSQPair current = labelMap[labelId];
           // Sum labels. This is integration!
@@ -279,6 +314,10 @@ namespace Mantid
           outWS->setSignalAt(i, 0);
         }
         outWS->setErrorSquaredAt(i, 0);
+        if(i % frequency == 0)
+        {
+          progress.doReport();
+        }
       }
 
       return outWS;

Code/Mantid/Framework/Crystal/src/IntegratePeaksUsingClusters.cpp

@@ -1,7 +1,15 @@
 /*WIKI*
-Handles integration of arbitary single crystal peaks shapes.
+Integrates arbitary shaped single crystal peaks defined on an [[MDHistoWorkspace]] using connected component analysis to determine
+regions of interest around each peak of the [[PeaksWorkspace]]. The output is an integrated [[PeaksWorkspace]] as well as an image 
+containing the labels assigned to each cluster for diagnostic and visualisation purposes.
 
-Uses connected component analysis to integrate peaks in an PeaksWorkspace over an MDHistoWorkspace of data.
+A threshold for the Peak should be defined below which, parts of the image are treated as background. In addition, a radius estimate
+is required to dispose of those clusters which are not to do with peaks, and also to associate clusters in the image with a peak center.
+You can view the radius estimate as a radius cut-off.
+
+This algorithm uses an imaging technique, and it is therefore important that the MDHistoWorkspace you are using is binned to a sufficient
+resolution via [[BinMD]]. You can overlay the intergrated peaks workspace in the [[MantidPlot:_SliceViewer#Viewing_Peaks_Workspaces|Slice Viewer]] over
+the generated Cluster Labeled OutputWorkspaceMD to see what the interation region used for each peak amounts to. 
 *WIKI*/
 
 #include "MantidCrystal/IntegratePeaksUsingClusters.h"
@@ -11,6 +19,7 @@ Uses connected component analysis to integrate peaks in an PeaksWorkspace over a
 #include "MantidAPI/WorkspaceProperty.h"
 #include "MantidAPI/IMDIterator.h"
 #include "MantidAPI/AlgorithmManager.h"
+#include "MantidAPI/Progress.h"
 #include "MantidKernel/MultiThreaded.h"
 #include "MantidKernel/CompositeValidator.h"
 #include "MantidKernel/MandatoryValidator.h"
@@ -110,7 +119,7 @@ namespace Mantid
 
       declareProperty(new PropertyWithValue<double>("Threshold", 0, positiveValidator->clone(), Direction::Input), "Threshold signal above which to consider peaks");
       declareProperty(new WorkspaceProperty<IPeaksWorkspace>("OutputWorkspace","",Direction::Output), "An output integrated peaks workspace."); 
-      declareProperty(new WorkspaceProperty<IMDHistoWorkspace>("OutputWorkspaceMD","",Direction::Output), "MDHistoWorkspace containing the clusters.");
+      declareProperty(new WorkspaceProperty<IMDHistoWorkspace>("OutputWorkspaceMD","",Direction::Output), "MDHistoWorkspace containing the labeled clusters used by the algorithm.");
     }
 
     //----------------------------------------------------------------------------------------------
@@ -141,17 +150,21 @@ namespace Mantid
 
       const double threshold = getProperty("Threshold");
       const double radiusEstimate = getProperty("RadiusEstimate");
-      PeakBackground background(peakWS, radiusEstimate, threshold, NoNormalization, mdCoordinates);
-      //HardThresholdBackground background(threshold, normalization);
+      PeakBackground backgroundStrategy(peakWS, radiusEstimate, threshold, NoNormalization, mdCoordinates);
 
       ConnectedComponentLabeling analysis; 
       LabelIdIntensityMap labelMap;
       PositionToLabelIdMap positionMap;
-      IMDHistoWorkspace_sptr clusters = analysis.executeAndIntegrate(mdWS, &background, labelMap, positionMap);
+
+      Progress progress(this, 0, 1, 1);
+      IMDHistoWorkspace_sptr clusters = analysis.executeAndIntegrate(mdWS, &backgroundStrategy, labelMap, positionMap, progress);
 
       // Link integrated values up with peaks.
+      const size_t nPeaks = peakWS->getNumberPeaks();
+      progress.resetNumSteps(nPeaks, 0, 1);
+      progress.doReport("Writing out PeaksWorkspace");
       PARALLEL_FOR1(peakWS)
-      for(int i =0; i < peakWS->getNumberPeaks(); ++i)
+      for(int i =0; i < nPeaks; ++i)
       {
         PARALLEL_START_INTERUPT_REGION
         IPeak& peak = peakWS->getPeak(i);
@@ -179,6 +192,7 @@ namespace Mantid
           peak.setIntensity(labelMap[ iterator->second ].get<0>());
           peak.setSigmaIntensity(labelMap[ iterator->second ].get<1>());
         }
+        progress.report();
         PARALLEL_END_INTERUPT_REGION
       }
       PARALLEL_CHECK_INTERUPT_REGION