Merge pull request #32427 from walshmm/sans657_findcenterofmassposition2_cleanup_fix-ornl-next

Sans657 findcenterofmassposition2 cleanup fix ornl next

Framework/Algorithms/CMakeLists.txt

@@ -335,6 +335,7 @@ set(SRC_FILES
     src/WeightingStrategy.cpp
     src/WienerSmooth.cpp
     src/WorkflowAlgorithmRunner.cpp
+    src/WorkspaceBoundingBox.cpp
     src/WorkspaceJoiners.cpp
     src/XDataConverter.cpp
     src/XrayAbsorptionCorrection.cpp)
@@ -684,6 +685,7 @@ set(INC_FILES
     inc/MantidAlgorithms/WeightingStrategy.h
     inc/MantidAlgorithms/WienerSmooth.h
     inc/MantidAlgorithms/WorkflowAlgorithmRunner.h
+    inc/MantidAlgorithms/WorkspaceBoundingBox.h
     inc/MantidAlgorithms/WorkspaceJoiners.h
     inc/MantidAlgorithms/XDataConverter.h
     inc/MantidAlgorithms/XrayAbsorptionCorrection.h)

Framework/Algorithms/inc/MantidAlgorithms/FindCenterOfMassPosition2.h

@@ -11,6 +11,8 @@
 //----------------------------------------------------------------------
 #include "MantidAPI/Algorithm.h"
 #include "MantidAlgorithms/DllConfig.h"
+#include "MantidAlgorithms/WorkspaceBoundingBox.h"
+#include "MantidDataObjects/EventWorkspace_fwd.h"
 
 namespace Mantid {
 namespace Algorithms {
@@ -63,6 +65,14 @@ private:
   void init() override;
   /// Execution code
   void exec() override;
+  /// Helper functions
+  void findCenterOfMass(API::MatrixWorkspace_sptr inputWS, double &centerX, double &centerY, const int numSpec,
+                        API::Progress &progress);
+  API::MatrixWorkspace_sptr sumUsingSpectra(DataObjects::EventWorkspace_const_sptr inputEventWS, const int numSpec,
+                                            API::Progress &progress);
+  void storeOutputWorkspace(double centerX, double centerY);
+  // Iteration cutoff
+  const int m_maxIteration = 200;
 };
 
 } // namespace Algorithms

Framework/Algorithms/inc/MantidAlgorithms/WorkspaceBoundingBox.h

+// Mantid Repository : https://github.com/mantidproject/mantid
+//
+// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
+//   NScD Oak Ridge National Laboratory, European Spallation Source,
+//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
+// SPDX - License - Identifier: GPL - 3.0 +
+#pragma once
+
+#include "MantidAPI/MatrixWorkspace_fwd.h"
+#include "MantidKernel/Logger.h"
+
+namespace Mantid {
+namespace Algorithms {
+namespace {
+// static logger
+Kernel::Logger g_log("WorkspaceBoundingBox");
+} // namespace
+
+/* This is a simple class originally intended for use solely with FindCenterOfMassPosition2.cpp
+ *
+ */
+class WorkspaceBoundingBox {
+public:
+  WorkspaceBoundingBox(API::MatrixWorkspace_sptr workspace);
+  WorkspaceBoundingBox();
+  ~WorkspaceBoundingBox();
+
+  API::MatrixWorkspace_sptr getWorkspace() { return workspace; }
+  double getX() { return x; }
+  double getY() { return y; }
+  double getCenterX() { return centerX; }
+  double getCenterY() { return centerY; }
+  double getXMin() { return xMin; }
+  double getXMax() { return xMax; }
+  double getYMin() { return yMin; }
+  double getYMax() { return yMax; }
+
+  void setPosition(double x, double y);
+  void setCenter(double x, double y);
+  void setBounds(double xMin, double xMax, double yMin, double yMax);
+
+  double calculateDistance();
+  double calculateRadiusX();
+  double calculateRadiusY();
+
+  double updatePositionAndReturnCount(int index);
+  int findFirstValidWs(const int numSpec);
+  bool isValidWs(int index);
+  bool isOutOfBoundsOfNonDirectBeam(const double beamRadius, int index, const bool directBeam);
+  bool containsPoint(double x, double y);
+  void normalizePosition(double x, double y);
+  void updateMinMax(int index);
+
+private:
+  API::MatrixWorkspace_sptr workspace;
+  double x = 0;
+  double y = 0;
+  double centerX;
+  double centerY;
+  double xMin = 0;
+  double xMax = 0;
+  double yMin = 0;
+  double yMax = 0;
+  const int m_specID = 0;
+};
+
+} // namespace Algorithms
+} // namespace Mantid

Framework/Algorithms/src/FindCenterOfMassPosition2.cpp

@@ -32,7 +32,9 @@ using namespace Geometry;
 using namespace DataObjects;
 
 void FindCenterOfMassPosition2::init() {
-  auto wsValidator = std::make_shared<CompositeValidator>();
+  const auto wsValidator = std::make_shared<CompositeValidator>();
+  const auto positiveDouble = std::make_shared<BoundedValidator<double>>();
+
   wsValidator->add<HistogramValidator>();
   declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, wsValidator));
   declareProperty("Output", "",
@@ -51,7 +53,6 @@ void FindCenterOfMassPosition2::init() {
                   "center of mass "
                   "of the scattering data will be computed by excluding the beam area.");
 
-  auto positiveDouble = std::make_shared<BoundedValidator<double>>();
   positiveDouble->setLower(0);
   declareProperty("BeamRadius", 0.0155, positiveDouble,
                   "Radius of the beam area, in meters, used the exclude the "
@@ -59,186 +60,176 @@ void FindCenterOfMassPosition2::init() {
                   "the center of mass of the scattering pattern.");
 }
 
-void FindCenterOfMassPosition2::exec() {
-  MatrixWorkspace_sptr inputWSWvl = getProperty("InputWorkspace");
-  MatrixWorkspace_sptr inputWS;
+/** Integrate events in the inputEventWS to determine thier sum and error values
+ *
+ *  @param inputEventWS  :: event workspace to ingegrate
+ *  @param numSpec  :: number of events in the workspace to iterate through
+ *  @param progress :: object for reporting progress of the operation
+ *  @return workspace containing the calculated x, y, and e data over numSpec
+ */
+API::MatrixWorkspace_sptr
+FindCenterOfMassPosition2::sumUsingSpectra(DataObjects::EventWorkspace_const_sptr inputEventWS, const int numSpec,
+                                           Progress &progress) {
+  std::vector<double> yValues(numSpec);
+  std::vector<double> eValues(numSpec);
+
+  PARALLEL_FOR_NO_WSP_CHECK()
+  for (int i = 0; i < numSpec; i++) {
+    double sum_i(0), err_i(0);
+    progress.report("Integrating events");
+    const EventList &el = inputEventWS->getSpectrum(i);
+    el.integrate(0, 0, true, sum_i, err_i);
+    yValues[i] = sum_i;
+    eValues[i] = err_i;
+  }
 
-  // Option to exclude beam area
-  bool direct_beam = getProperty("DirectBeam");
+  IAlgorithm_sptr algo = createChildAlgorithm("CreateWorkspace", 0.7, 1.0);
+  algo->setProperty<std::vector<double>>("DataX", std::vector<double>(2, 0.0));
+  algo->setProperty<std::vector<double>>("DataY", yValues);
+  algo->setProperty<std::vector<double>>("DataE", eValues);
+  algo->setProperty<int>("NSpec", numSpec);
+  algo->execute();
 
-  // TODO: Need an input for the X bin to use, assume 0 for now
-  int specID = 0;
-  // Initial center location
-  double center_x = getProperty("CenterX");
-  double center_y = getProperty("CenterY");
-  const double tolerance = getProperty("Tolerance");
-  // Iteration cutoff
-  int max_iteration = 200;
-  // Radius of the beam area, in pixels
-  double beam_radius = getProperty("BeamRadius");
+  MatrixWorkspace_sptr inputWS = algo->getProperty("OutputWorkspace");
+  return inputWS;
+}
 
-  // Get the number of monitors. We assume that all monitors are stored in the
-  // first spectra
-  const auto numSpec = static_cast<int>(inputWSWvl->getNumberHistograms());
+// find the min/max for x/y coords in set of valid spectra, update position of bounding box
+double initBoundingBox(WorkspaceBoundingBox &boundingBox, const int numSpec, const double beamRadius,
+                       const bool directBeam) {
+  double totalCount = 0;
+  for (int i = 0; i < numSpec; i++) {
+    if (!boundingBox.isValidWs(i))
+      continue;
 
-  // Set up the progress reporting object
-  Progress progress(this, 0.0, 1.0, max_iteration);
+    boundingBox.updateMinMax(i);
 
-  EventWorkspace_const_sptr inputEventWS = std::dynamic_pointer_cast<const EventWorkspace>(inputWSWvl);
-  if (inputEventWS) {
-    std::vector<double> y_values(numSpec);
-    std::vector<double> e_values(numSpec);
-
-    PARALLEL_FOR_NO_WSP_CHECK()
-    for (int i = 0; i < numSpec; i++) {
-      double sum_i(0), err_i(0);
-      progress.report("Integrating events");
-      const EventList &el = inputEventWS->getSpectrum(i);
-      el.integrate(0, 0, true, sum_i, err_i);
-      y_values[i] = sum_i;
-      e_values[i] = err_i;
-    }
+    if (!boundingBox.isOutOfBoundsOfNonDirectBeam(beamRadius, i, directBeam))
+      continue;
+    else
+      totalCount += boundingBox.updatePositionAndReturnCount(i);
+  }
+  return totalCount;
+}
 
-    IAlgorithm_sptr algo = createChildAlgorithm("CreateWorkspace", 0.7, 1.0);
-    algo->setProperty<std::vector<double>>("DataX", std::vector<double>(2, 0.0));
-    algo->setProperty<std::vector<double>>("DataY", y_values);
-    algo->setProperty<std::vector<double>>("DataE", e_values);
-    algo->setProperty<int>("NSpec", numSpec);
-    algo->execute();
+// In subsequent iterations check if spectra fits in the normalized bounding box(generated by previous iterations)
+// If so, update position
+double updateBoundingBox(WorkspaceBoundingBox &boundingBox, WorkspaceBoundingBox previousBoundingBox, const int numSpec,
+                         const double beamRadius, const bool directBeam) {
+  double totalCount = 0;
+  for (int i = 0; i < numSpec; i++) {
+    if (!boundingBox.isValidWs(i))
+      continue;
 
-    inputWS = algo->getProperty("OutputWorkspace");
-    WorkspaceFactory::Instance().initializeFromParent(*inputWSWvl, *inputWS, false);
-  } else {
-    // Sum up all the wavelength bins
-    IAlgorithm_sptr childAlg = createChildAlgorithm("Integration");
-    childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputWSWvl);
-    childAlg->executeAsChildAlg();
-    inputWS = childAlg->getProperty("OutputWorkspace");
+    const V3D position = boundingBox.getWorkspace()->spectrumInfo().position(i);
+
+    if (previousBoundingBox.containsPoint(position.X(), position.Y())) {
+      if (!boundingBox.isOutOfBoundsOfNonDirectBeam(beamRadius, i, directBeam))
+        continue;
+      else
+        totalCount += boundingBox.updatePositionAndReturnCount(i);
+    }
   }
+  return totalCount;
+}
+
+bool equals(double a, double b) { return fabs(a - b) < std::numeric_limits<double>::min(); }
+
+/** Iterates through spectrum in the input workspace finding the center of mass until
+ *  we converge to within the tolerance specified in meters
+ *
+ *  @param inputWS  :: workspace to find the center of mass of
+ *  @param centerX  :: save the real center of mass x coord here
+ *  @param centerY  :: save the real center of mass y coord here
+ *  @param numSpec  :: number of spectrum in the workspace to iterate through
+ *  @param progress :: object for reporting progress of the operation
+ */
+void FindCenterOfMassPosition2::findCenterOfMass(API::MatrixWorkspace_sptr inputWS, double &centerX, double &centerY,
+                                                 const int numSpec, Progress &progress) {
+  const double tolerance = getProperty("Tolerance");
+  const bool directBeam = getProperty("DirectBeam");
+  const double beamRadius = getProperty("BeamRadius");
 
   // Define box around center of mass so that only pixels in an area
   // _centered_ on the latest center position are considered. At each
   // iteration we will recompute the bounding box, and we will make
-  // it as large as possible. The largest box is defined as:
-  double xmin0 = 0;
-  double xmax0 = 0;
-  double ymin0 = 0;
-  double ymax0 = 0;
+  // it as large as possible. The largest box is defined in:
+  WorkspaceBoundingBox boundingBox(inputWS);
+  boundingBox.setCenter(centerX, centerY);
 
   // Starting values for the bounding box and the center
-  double xmin = xmin0;
-  double xmax = xmax0;
-  double ymin = ymin0;
-  double ymax = ymax0;
+  WorkspaceBoundingBox previousBoundingBox;
+  previousBoundingBox.setBounds(0., 0., 0., 0.);
 
   // Initialize book-keeping
   double distance = -1;
-  double distance_check = 0;
-  int n_local_minima = 0;
-  int n_iteration = 0;
+  double distanceCheck = 0;
+  double totalCount = initBoundingBox(boundingBox, numSpec, beamRadius, directBeam);
+
+  int totalLocalMinima = 0;
+  int totalIterations = 0;
 
   // Find center of mass and iterate until we converge
-  // to within a quarter of a pixel
-  bool first_run = true;
+  // to within the tolerance specified in meters
   while (distance > tolerance || distance < 0) {
-    // Count histogram for normalization
-    double total_count = 0;
-    double position_x = 0;
-    double position_y = 0;
-
-    const auto &spectrumInfo = inputWS->spectrumInfo();
-    for (int i = 0; i < numSpec; i++) {
-      if (!spectrumInfo.hasDetectors(i)) {
-        g_log.warning() << "Workspace index " << i << " has no detector assigned to it - discarding\n";
-        continue;
-      }
-      // Skip if we have a monitor or if the detector is masked.
-      if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i))
-        continue;
-
-      // Get the current spectrum
-      auto &YIn = inputWS->y(i);
-      // Skip if NaN of inf
-      if (std::isnan(YIn[specID]) || std::isinf(YIn[specID]))
-        continue;
-      const V3D pos = spectrumInfo.position(i);
-      double x = pos.X();
-      double y = pos.Y();
-
-      if (first_run) {
-        xmin0 = std::min(x, xmin0);
-        xmax0 = std::max(x, xmax0);
-        ymin0 = std::min(y, ymin0);
-        ymax0 = std::max(y, ymax0);
-      }
-
-      if (first_run || (x >= xmin && x <= xmax && y >= ymin && y <= ymax)) {
-        if (!direct_beam) {
-          double dx = x - center_x;
-          double dy = y - center_y;
-          if (dx * dx + dy * dy < beam_radius * beam_radius)
-            continue;
-        }
-        position_x += YIn[specID] * x;
-        position_y += YIn[specID] * y;
-        total_count += YIn[specID];
-      }
-    }
-
     // Normalize output to find center of mass position
-    position_x /= total_count;
-    position_y /= total_count;
-
+    boundingBox.normalizePosition(totalCount, totalCount);
     // Compute the distance to the previous iteration
-    distance =
-        sqrt((center_x - position_x) * (center_x - position_x) + (center_y - position_y) * (center_y - position_y));
-
-    // Modify the bounding box around the detector region used to
-    // compute the center of mass so that it is centered around
-    // the new center of mass position.
-    double radius_x = std::min((position_x - xmin0), (xmax0 - position_x));
-    double radius_y = std::min((position_y - ymin0), (ymax0 - position_y));
+    distance = boundingBox.calculateDistance();
+    // Recenter around new mass position
+    double radiusX = boundingBox.calculateRadiusX();
+    double radiusY = boundingBox.calculateRadiusY();
 
-    if (!direct_beam && (radius_x <= beam_radius || radius_y <= beam_radius)) {
+    if (!directBeam && (radiusX <= beamRadius || radiusY <= beamRadius)) {
       g_log.error() << "Center of mass falls within the beam center area: "
                        "stopping here\n";
       break;
     }
 
-    center_x = position_x;
-    center_y = position_y;
+    boundingBox.setCenter(boundingBox.getX(), boundingBox.getY());
+    previousBoundingBox.setBounds(boundingBox.getCenterX() - radiusX, boundingBox.getCenterX() + radiusX,
+                                  boundingBox.getCenterY() - radiusY, boundingBox.getCenterY() + radiusY);
 
-    xmin = center_x - radius_x;
-    xmax = center_x + radius_x;
-    ymin = center_y - radius_y;
-    ymax = center_y + radius_y;
-
-    // Sanity check to avoid getting stuck
-    if (distance == distance_check) {
-      n_local_minima++;
+    // Check to see if we have the same result
+    // as the previous iteration
+    if (equals(distance, distanceCheck)) {
+      totalLocalMinima++;
     } else {
-      n_local_minima = 0;
+      totalLocalMinima = 0;
     }
 
     // Quit if we found the exact same distance five times in a row.
-    if (n_local_minima > 5) {
+    if (totalLocalMinima > 5) {
       g_log.warning() << "Found the same or equivalent center of mass locations "
                          "more than 5 times in a row: stopping here\n";
       break;
     }
 
     // Quit if we haven't converged after the maximum number of iterations.
-    if (++n_iteration > max_iteration) {
-      g_log.warning() << "More than " << max_iteration << " iteration to find beam center: stopping here\n";
+    if (++totalIterations > m_maxIteration) {
+      g_log.warning() << "More than " << m_maxIteration << " iteration to find beam center: stopping here\n";
       break;
     }
 
-    distance_check = distance;
-    first_run = false;
+    distanceCheck = distance;
+
+    // Count histogram for normalization
+    boundingBox.setPosition(0, 0);
+    totalCount = updateBoundingBox(boundingBox, previousBoundingBox, numSpec, beamRadius, directBeam);
 
     progress.report("Find Beam Center");
   }
+  centerX = boundingBox.getCenterX();
+  centerY = boundingBox.getCenterY();
+}
 
+/** Package the algorithm outputs one of two ways depending on whether or
+ *  not it was given an input EventWorkspace to start with
+ *
+ *  @param centerX  :: center of mass x coord to package
+ *  @param centerY  :: center of mass y coord to package
+ */
+void FindCenterOfMassPosition2::storeOutputWorkspace(double centerX, double centerY) {
   std::string output = getProperty("Output");
 
   // If an output workspace name was given, create a TableWorkspace with the
@@ -257,9 +248,9 @@ void FindCenterOfMassPosition2::exec() {
     m_result->addColumn("double", "Value");
 
     Mantid::API::TableRow row = m_result->appendRow();
-    row << "X (m)" << center_x;
+    row << "X (m)" << centerX;
     row = m_result->appendRow();
-    row << "Y (m)" << center_y;
+    row << "Y (m)" << centerY;
 
     setProperty("OutputWorkspace", m_result);
   } else {
@@ -268,12 +259,41 @@ void FindCenterOfMassPosition2::exec() {
       declareProperty(std::make_unique<ArrayProperty<double>>("CenterOfMass", std::make_shared<NullValidator>(),
                                                               Direction::Output));
     std::vector<double> center_of_mass;
-    center_of_mass.emplace_back(center_x);
-    center_of_mass.emplace_back(center_y);
+    center_of_mass.emplace_back(centerX);
+    center_of_mass.emplace_back(centerY);
     setProperty("CenterOfMass", center_of_mass);
   }
 
-  g_log.information() << "Center of Mass found at x=" << center_x << " y=" << center_y << '\n';
+  g_log.information() << "Center of Mass found at x=" << centerX << " y=" << centerY << '\n';
+}
+
+void FindCenterOfMassPosition2::exec() {
+  const MatrixWorkspace_sptr inputWSWvl = getProperty("InputWorkspace");
+  double centerX = getProperty("CenterX");
+  double centerY = getProperty("CenterY");
+
+  MatrixWorkspace_sptr inputWS;
+
+  // Get the number of monitors. We assume that all monitors are stored in the
+  // first spectra
+  const auto numSpec = static_cast<int>(inputWSWvl->getNumberHistograms());
+  // Set up the progress reporting object
+  Progress progress(this, 0.0, 1.0, m_maxIteration);
+  EventWorkspace_const_sptr inputEventWS = std::dynamic_pointer_cast<const EventWorkspace>(inputWSWvl);
+
+  if (inputEventWS) {
+    inputWS = sumUsingSpectra(inputEventWS, numSpec, progress);
+    WorkspaceFactory::Instance().initializeFromParent(*inputWSWvl, *inputWS, false);
+  } else {
+    // Sum up all the wavelength bins
+    IAlgorithm_sptr childAlg = createChildAlgorithm("Integration");
+    childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputWSWvl);
+    childAlg->executeAsChildAlg();
+    inputWS = childAlg->getProperty("OutputWorkspace");
+  }
+
+  findCenterOfMass(inputWS, centerX, centerY, numSpec, progress);
+  storeOutputWorkspace(centerX, centerY);
 }
 
 } // namespace Algorithms

Framework/Algorithms/src/WorkspaceBoundingBox.cpp

+#include "MantidAlgorithms/WorkspaceBoundingBox.h"
+#include "MantidAPI/MatrixWorkspace.h"
+#include "MantidAPI/SpectrumInfo.h"
+
+namespace Mantid {
+namespace Algorithms {
+
+WorkspaceBoundingBox::WorkspaceBoundingBox(API::MatrixWorkspace_sptr workspace) : workspace(workspace) {}
+
+WorkspaceBoundingBox::WorkspaceBoundingBox() {}
+
+WorkspaceBoundingBox::~WorkspaceBoundingBox() {}
+
+void WorkspaceBoundingBox::setPosition(double x, double y) {
+  this->x = x;
+  this->y = y;
+}
+
+void WorkspaceBoundingBox::setCenter(double x, double y) {
+  this->centerX = x;
+  this->centerY = y;
+}
+
+void WorkspaceBoundingBox::setBounds(double xMin, double xMax, double yMin, double yMax) {
+  this->xMin = xMin;
+  this->xMax = xMax;
+  this->yMin = yMin;
+  this->yMax = yMax;
+}
+
+/** Performs checks on the spectrum located at index to determine if
+ *  it is acceptable to be operated on
+ *
+ *  @param index :: index of spectrum data
+ *  @return true/false if its valid
+ */
+bool WorkspaceBoundingBox::isValidWs(int index) {
+  const auto spectrumInfo = this->workspace->spectrumInfo();
+  if (!spectrumInfo.hasDetectors(index)) {
+    g_log.warning() << "Workspace index " << index << " has no detector assigned to it - discarding\n";
+    return false;
+  }
+  // Skip if we have a monitor or if the detector is masked.
+  if (spectrumInfo.isMonitor(index) || spectrumInfo.isMasked(index))
+    return false;
+
+  // Get the current spectrum
+  auto &YIn = this->workspace->y(index);
+  // Skip if NaN of inf
+  if (std::isnan(YIn[m_specID]) || std::isinf(YIn[m_specID]))
+    return false;
+  return true;
+}
+
+/** Searches for the first valid spectrum info in member variable `workspace`
+ *
+ *  @param numSpec :: the number of spectrum in the workspace to search through
+ *  @return index of first valid spectrum
+ */
+int WorkspaceBoundingBox::findFirstValidWs(const int numSpec) {
+  const auto spectrumInfo = this->workspace->spectrumInfo();
+  int i;
+  for (i = 0; i < numSpec; ++i) {
+    if (isValidWs(i))
+      break;
+  }
+  return i;
+}
+
+/** Sets member variables x/y to new x/y based on
+ *  spectrum info and historgram data at the given index
+ *
+ *  @param index :: index of spectrum data
+ *  @return number of points of histogram data at index
+ */
+double WorkspaceBoundingBox::updatePositionAndReturnCount(int index) {
+  const auto spectrumInfo = this->workspace->spectrumInfo();
+  auto &YIn = this->workspace->y(index);
+  double x = spectrumInfo.position(index).X();
+  double y = spectrumInfo.position(index).Y();
+  this->x += YIn[m_specID] * x;
+  this->y += YIn[m_specID] * y;
+  return YIn[m_specID];
+}
+
+/** Compare current mins and maxs to the coordinates of the spectrum at index
+ *  expnd mins and maxs to include this spectrum
+ *