WorkspaceBoundingBox.cpp

Go to the documentation of this file.

#include "MantidAlgorithms/WorkspaceBoundingBox.h"
#include "MantidAPI/MatrixWorkspace.h"
 
namespace Mantid::Algorithms {
 
namespace {
constexpr int HISTOGRAM_INDEX{0};
}
 
WorkspaceBoundingBox::WorkspaceBoundingBox(const API::MatrixWorkspace_const_sptr &workspace,
                                           const double integrationRadius, const double beamRadius,
                                           const bool ignoreDirectBeam, const double cenX, const double cenY)
    : m_workspace(workspace), m_integrationRadiusSq(integrationRadius * integrationRadius),
      m_beamRadiusSq(beamRadius * beamRadius), m_ignoreDirectBeam(ignoreDirectBeam) {
  if (m_workspace->y(0).size() != 1)
    throw std::runtime_error("This object only works with integrated workspaces");
 
  m_spectrumInfo = &workspace->spectrumInfo();
  m_numSpectra = m_workspace->getNumberHistograms();
 
  this->setCenterPrev(cenX, cenY);
  this->initOverallRangeAndFindFirstCenter();
}
 
WorkspaceBoundingBox::~WorkspaceBoundingBox() = default;
 
Kernel::V3D WorkspaceBoundingBox::position(const std::size_t index) const {
  if (!m_spectrumInfo)
    throw std::runtime_error("SpectrumInfo object is not initialized");
 
  return m_spectrumInfo->position(index);
}
 
// find the min/max for x/y coords in set of valid spectra, update position of bounding box
void WorkspaceBoundingBox::initOverallRangeAndFindFirstCenter() {
  // reset temporary place for the calculation
  this->resetIntermediatePosition();
  double totalCount = 0;
 
  for (std::size_t i = 0; i < m_numSpectra; i++) {
    if (!isValidIndex(i))
      continue;
 
    updateMinMax(i);
 
    if (includeInIntegration(i))
      totalCount += updatePositionAndReturnCount(i);
  }
  this->normalizePosition(totalCount);
}
 
// In subsequent iterations check if spectra fits in the normalized bounding box(generated by previous iterations)
// If so, update position
void WorkspaceBoundingBox::findNewCenterPosition() {
  // reset temporary place for the calculation
  this->resetIntermediatePosition();
  double totalCount = 0;
 
  const auto &spectrumInfo = m_workspace->spectrumInfo();
  for (std::size_t i = 0; i < m_numSpectra; i++) {
    if (!isValidIndex(i))
      continue;
 
    const V3D pos = spectrumInfo.position(i);
 
    if (this->symmetricRegionContainsPoint(pos.X(), pos.Y()) && includeInIntegration(pos)) {
      totalCount += updatePositionAndReturnCount(i);
    }
  }
  this->normalizePosition(totalCount);
}
 
double WorkspaceBoundingBox::countsValue(const std::size_t index) const {
  return m_workspace->y(index)[HISTOGRAM_INDEX];
}
 
void WorkspaceBoundingBox::resetIntermediatePosition() {
  this->m_centerXPosCurr = 0.;
  this->m_centerYPosCurr = 0.;
}
 
void WorkspaceBoundingBox::setCenterPrev(const double x, const double y) {
  this->m_centerXPosPrev = x;
  this->m_centerYPosPrev = y;
}
 
void WorkspaceBoundingBox::setBounds(const double xMin, const double xMax, const double yMin, const double yMax) {
  this->m_xBoxMin = xMin;
  this->m_xBoxMax = xMax;
  this->m_yBoxMin = yMin;
  this->m_yBoxMax = yMax;
}
 
bool WorkspaceBoundingBox::isValidIndex(const std::size_t index) const {
  if (!m_spectrumInfo)
    throw std::runtime_error("SpectrumInfo object is not initialized");
  if (!m_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 (this->m_spectrumInfo->isMonitor(index) || this->m_spectrumInfo->isMasked(index))
    return false;
 
  // Get the current spectrum
  const auto YIn = this->countsValue(index);
  // Skip if NaN of inf
  if (std::isnan(YIn) || std::isinf(YIn))
    return false;
  return true;
}
 
double WorkspaceBoundingBox::updatePositionAndReturnCount(const std::size_t index) {
  const auto counts = this->countsValue(index);
  const auto &pos = this->position(index);
 
  this->m_centerXPosCurr += counts * pos.X();
  this->m_centerYPosCurr += counts * pos.Y();
 
  return counts;
}
 
void WorkspaceBoundingBox::updateMinMax(const std::size_t index) {
  const auto &pos = this->position(index);
  const double x = pos.X();
  const double y = pos.Y();
 
  this->m_xPosMin = std::min(x, this->m_xPosMin);
  this->m_xPosMax = std::max(x, this->m_xPosMax);
  this->m_yPosMin = std::min(y, this->m_yPosMin);
  this->m_yPosMax = std::max(y, this->m_yPosMax);
}
 
bool WorkspaceBoundingBox::includeInIntegration(const std::size_t index) {
  return this->includeInIntegration(this->position(index));
}
 
bool WorkspaceBoundingBox::includeInIntegration(const Kernel::V3D &position) {
  const double dx = position.X() - this->m_centerXPosPrev;
  const double dy = position.Y() - this->m_centerYPosPrev;
  const double distanceFromCenterSq = dx * dx + dy * dy;
 
  if (m_ignoreDirectBeam) {
    if (distanceFromCenterSq < m_beamRadiusSq)
      return false;
  }
 
  return bool(distanceFromCenterSq < m_integrationRadiusSq);
}
 
double WorkspaceBoundingBox::distanceFromPrevious() const {
  const auto xExtent = (m_centerXPosPrev - m_centerXPosCurr);
  const auto yExtent = (m_centerYPosPrev - m_centerYPosCurr);
  return sqrt(xExtent * xExtent + yExtent * yExtent);
}
 
bool WorkspaceBoundingBox::centerOfMassWithinBeamCenter() {
  if (m_ignoreDirectBeam) {
    const double radiusX = this->calculateRadiusX();
    const double radiusY = this->calculateRadiusY();
    if (radiusX * radiusX <= m_beamRadiusSq || radiusY * radiusY <= m_beamRadiusSq) {
      return true;
    }
  }
  return false;
}
 
void WorkspaceBoundingBox::prepareCenterCalculation() {
  this->setCenterPrev(m_centerXPosCurr, m_centerYPosCurr);
 
  const double radiusX = this->calculateRadiusX();
  const double radiusY = this->calculateRadiusY();
  this->setBounds(m_centerXPosCurr - radiusX, m_centerXPosCurr + radiusX, m_centerYPosCurr - radiusY,
                  m_centerYPosCurr + radiusY);
}
 
double WorkspaceBoundingBox::calculateRadiusX() const {
  return std::min((m_centerXPosCurr - m_xPosMin), (m_xPosMax - m_centerXPosCurr));
}
 
double WorkspaceBoundingBox::calculateRadiusY() const {
  return std::min((m_centerYPosCurr - m_yPosMin), (m_yPosMax - m_centerYPosCurr));
}
 
void WorkspaceBoundingBox::normalizePosition(const double totalCounts) {
  this->m_centerXPosCurr /= std::fabs(totalCounts);
  this->m_centerYPosCurr /= std::fabs(totalCounts);
}
 
bool WorkspaceBoundingBox::symmetricRegionContainsPoint(double x, double y) {
  return (x <= this->m_xBoxMax && x >= this->m_xBoxMin && y <= m_yBoxMax && y >= m_yBoxMin);
}
 
} // namespace Mantid::Algorithms