Commit d870d12b authored by Walsh, Michael's avatar Walsh, Michael
Browse files

refactor is passing tests

parent 25590ae0
......@@ -6,68 +6,61 @@
// SPDX - License - Identifier: GPL - 3.0 +
#pragma once
#include "MantidAPI/MatrixWorkspace_fwd.h"
#include "MantidKernel/Logger.h"
namespace Mantid {
namespace Algorithms {
/* This is a simple class originally intended for use solely with FindCenterOfMassPosition2.cpp
*
*
*/
class WorkspaceBoundingBox {
public:
WorkspaceBoundingBox(API::MatrixWorkspace_sptr workspace);
~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){
this.x = x;
this.y = y;
}
WorkspaceBoundingBox(API::MatrixWorkspace_sptr workspace, Kernel::Logger &g_log);
WorkspaceBoundingBox(Kernel::Logger &g_log);
~WorkspaceBoundingBox();
void setCenter(double x, double y) {
this.centerX = x;
this.centerY = y;
}
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 setBounds(double xMin, double xMax, double yMin, double yMax) {
this.xMin = xMin;
this.xMax = xMax;
this.yMin = yMin;
this.yMax = 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 beam_radius, int index, const bool direct_beam);
bool containsPoint(double x, double y);
void normalizePosition(double x, double y);
void updateMinMax(int index);
double calculateDistance();
double calculateRadiusX();
double calculateRadiusY();
double updatePositionAndReturnCount(int index);
int findFirstValidWs(const int numSpec);
bool isValidWs(int index);
bool isOutOfBoundsOfNonDirectBeam(const double beam_radius, int index, const bool direct_beam);
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;
}
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;
/// Logger for this helper
Kernel::Logger &g_log;
};
}
}
\ No newline at end of file
} // namespace Algorithms
} // namespace Mantid
\ No newline at end of file
......@@ -87,97 +87,93 @@ FindCenterOfMassPosition2::sumUsingSpectra(DataObjects::EventWorkspace_const_spt
return inputWS;
}
double initBoundingBox(WorkspaceBoundingBox boundingBox, const double beam_radius, const bool direct_beam) {
double initBoundingBox(WorkspaceBoundingBox &boundingBox, const int numSpec, const double beamRadius,
const bool directBeam) {
double total_count = 0;
for (int i = 0; i < numSpec; i++) {
if(!boundingBox.isValidWs(i))
continue;
boundingBox.updateMinMax(i);
if (boundingBox.isOutOfBoundsOfNonDirectBeam(beam_radius, i, direct_beam))
continue;
else
total_count += boundingBox.updatePositionAndReturnCount(i);
if (!boundingBox.isValidWs(i))
continue;
boundingBox.updateMinMax(i);
if (!boundingBox.isOutOfBoundsOfNonDirectBeam(beamRadius, i, directBeam))
continue;
else
total_count += boundingBox.updatePositionAndReturnCount(i);
}
return total_count;
}
double updateBoundingBox(WorkspaceBoundingBox boundingBox, const double beam_radius, const bool direct_beam) {
double updateBoundingBox(WorkspaceBoundingBox &boundingBox, WorkspaceBoundingBox previousBoundingBox, const int numSpec,
const double beamRadius, const bool directBeam) {
double total_count = 0;
for (int i = 0; i < numSpec; i++) {
if(!boundingBox.isValidWs(i))
continue;
if (!boundingBox.isValidWs(i))
continue;
const V3D position = boundingBox.getWorkspace()->spectrumInfo().position(i);
const V3D position = boundingBox.getWorkspace()->spectrumInfo().position(i);
if (previousBoundingBox.containsPoint(position.X(), position.Y())) {
if (boundingBox.isOutOfBoundsOfNonDirectBeam(beam_radius, i, direct_beam))
continue;
else
total_count += boundingBox.updatePositionAndCount(i);
}
if (previousBoundingBox.containsPoint(position.X(), position.Y())) {
if (!boundingBox.isOutOfBoundsOfNonDirectBeam(beamRadius, i, directBeam))
continue;
else
total_count += boundingBox.updatePositionAndReturnCount(i);
}
}
return total_count;
}
bool equals(double a, double b) {
return fabs(a, b) < std::numeric_limits<double>::lowest;
}
bool equals(double a, double b) { return fabs(a - b) < std::numeric_limits<double>::min(); }
void FindCenterOfMassPosition2::findCenterOfMass(API::MatrixWorkspace_sptr inputWS, double &center_x, double &center_y,
const int numSpec, Progress &progress) {
const double tolerance = getProperty("Tolerance");
const bool direct_beam = getProperty("DirectBeam");
const double beam_radius = getProperty("BeamRadius");
const auto &spectrumInfo = inputWS->spectrumInfo();
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 in:
WorkspaceBoundingBox boundingBox(inputWs);
WorkspaceBoundingBox boundingBox(inputWS, g_log);
boundingBox.setCenter(center_x, center_y);
// Starting values for the bounding box and the center
WorkspaceBoundingBox previousBoundingBox(inputWs);
previousBoundingBox.setBounds(0,0,0,0);
WorkspaceBoundingBox previousBoundingBox(g_log);
previousBoundingBox.setBounds(0., 0., 0., 0.);
// Initialize book-keeping
double distance = -1;
double distance_check = 0;
double total_count = initBoundingBox(boundingBox, beam_radius, direct_beam);
double distanceCheck = 0;
double total_count = initBoundingBox(boundingBox, numSpec, beamRadius, directBeam);
int n_local_minima = 0;
int n_iteration = 0;
// Find center of mass and iterate until we converge
// to within the tolerance specified in meters
bool first_run = true;
while (distance > tolerance || distance < 0) {
// Normalize output to find center of mass position
boundingBox.normalizePosition(total_count, total_count);
// Compute the distance to the previous iteration
distance = boundingBox.calculateDistance();
// Recenter around new mass position
double radius_x = boundingBox.calculateRadiusX()
double radius_x = boundingBox.calculateRadiusX();
double radius_y = boundingBox.calculateRadiusY();
if (!direct_beam && (radius_x <= beam_radius || radius_y <= beam_radius)) {
if (!directBeam && (radius_x <= beamRadius || radius_y <= beamRadius)) {
g_log.error() << "Center of mass falls within the beam center area: "
"stopping here\n";
break;
}
boundingBox.setCenter(boundingBox.getX(), boundingBox.getY());
previousBoundingBox.setBounds(boundingBox.getCenterX() - radius_x,
boundingBox.getCenterX() + radius_x,
boundingBox.getCenterY() - radius_y,
boundingBox.getCenterY() + radius_y);
previousBoundingBox.setBounds(boundingBox.getCenterX() - radius_x, boundingBox.getCenterX() + radius_x,
boundingBox.getCenterY() - radius_y, boundingBox.getCenterY() + radius_y);
// Check to see if we have the same result
// as the previous iteration
if (equals(distance, distance_check)) {
if (equals(distance, distanceCheck)) {
n_local_minima++;
} else {
n_local_minima = 0;
......@@ -196,12 +192,11 @@ void FindCenterOfMassPosition2::findCenterOfMass(API::MatrixWorkspace_sptr input
break;
}
distance_check = distance;
first_run = false;
distanceCheck = distance;
// Count histogram for normalization
boudningBox.setPosition(0,0);
total_count = updateBoundingBox(boundingBox, beam_radius, direct_beam);
boundingBox.setPosition(0, 0);
total_count = updateBoundingBox(boundingBox, previousBoundingBox, numSpec, beamRadius, directBeam);
progress.report("Find Beam Center");
}
......@@ -264,7 +259,7 @@ void FindCenterOfMassPosition2::exec() {
if (inputEventWS) {
inputWS = sumUsingSpectra(inputEventWS, numSpec, progress);
WorkspaceFactory::Instance().initializeFromParent(*inputWSWvl, *inputWS, false);
} else {/diffraction/-/boards?scope=all&assignee_username=wqp
} else {
// Sum up all the wavelength bins
IAlgorithm_sptr childAlg = createChildAlgorithm("Integration");
childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputWSWvl);
......
#include "MantidAlgorithms/WorkspaceBoundingBox.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/SpectrumInfo.h"
WorkspaceBoundingBox(API::MatrixWorkspace_sptr workspace)
:workspace(workspace) {
//maybe calculate first run on init
namespace Mantid {
namespace Algorithms {
WorkspaceBoundingBox::WorkspaceBoundingBox(API::MatrixWorkspace_sptr workspace, Kernel::Logger &g_log)
: workspace(workspace), g_log(g_log) {}
WorkspaceBoundingBox::WorkspaceBoundingBox(Kernel::Logger &g_log) : g_log(g_log) {}
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;
}
bool WorkspaceBoundingBox::isValidWs(int index) {
const auto spectrumInfo = this.workspace->spectrumInfo();
const auto spectrumInfo = this->workspace->spectrumInfo();
if (!spectrumInfo.hasDetectors(index)) {
g_log.warning() << "Workspace index " << i << " has no detector assigned to it - discarding\n";
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.
......@@ -16,7 +40,7 @@ bool WorkspaceBoundingBox::isValidWs(int index) {
return false;
// Get the current spectrum
auto &YIn = inputWS->y(iindex);
auto &YIn = this->workspace->y(index);
// Skip if NaN of inf
if (std::isnan(YIn[m_specID]) || std::isinf(YIn[m_specID]))
return false;
......@@ -24,41 +48,43 @@ bool WorkspaceBoundingBox::isValidWs(int index) {
}
int WorkspaceBoundingBox::findFirstValidWs(const int numSpec) {
const auto spectrumInfo = this.workspace->spectrumInfo();
for(int i = 0; i < numSpec; ++i) {
if(isValidWs(spectrumInfo, i))
const auto spectrumInfo = this->workspace->spectrumInfo();
int i;
for (i = 0; i < numSpec; ++i) {
if (isValidWs(i))
break;
}
return i;
}
double WorkspaceBoundingBox::updatePositionAndReturnCount(int index){
const auto spectrumInfo = this.workspace->spectrumInfo();
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;
this->x += YIn[m_specID] * x;
this->y += YIn[m_specID] * y;
return YIn[m_specID];
}
void WorkspaceBoundingBox::updateMinMax(int index){
const auto spectrumInfo = this.workspace->spectrumInfo();
void WorkspaceBoundingBox::updateMinMax(int index) {
const auto spectrumInfo = this->workspace->spectrumInfo();
double x = spectrumInfo.position(index).X();
double y = spectrumInfo.position(index).Y();
this.xMin = std::min(x, this.xMin);
this.xMax = std::max(x, this.xMax);
this.yMin = std::min(y, this.yMin);
this.yMax = std::max(y, this.yMax);
this->xMin = std::min(x, this->xMin);
this->xMax = std::max(x, this->xMax);
this->yMin = std::min(y, this->yMin);
this->yMax = std::max(y, this->yMax);
}
bool WorkspaceBoundingBox::isOutOfBoundsOfNonDirectBeam(const double beam_radius, int index, const bool direct_beam){
const auto spectrumInfo = this.workspace->spectrumInfo();
bool WorkspaceBoundingBox::isOutOfBoundsOfNonDirectBeam(const double beamRadius, int index, const bool directBeam) {
const auto spectrumInfo = this->workspace->spectrumInfo();
double x = spectrumInfo.position(index).X();
double y = spectrumInfo.position(index).Y();
if (!direct_beam) {
double dx = x - this.centerX;
double dy = y - this.centerY;
if (dx * dx + dy * dy < beam_radius * beam_radius)
if (!directBeam) {
double dx = x - this->centerX;
double dy = y - this->centerY;
if (dx * dx + dy * dy < beamRadius * beamRadius)
return false;
}
return true;
......@@ -67,20 +93,19 @@ bool WorkspaceBoundingBox::isOutOfBoundsOfNonDirectBeam(const double beam_radius
double WorkspaceBoundingBox::calculateDistance() {
return sqrt((centerX - x) * (centerX - x) + (centerY - y) * (centerY - y));
}
double WorkspaceBoundingBox::calculateRadiusX() {
return std::min((x - xMin), (xMax - x));
}
double WorkspaceBoundingBox::calculateRadiusY() {
return std::min((y - yMin), (yMax - y));
double WorkspaceBoundingBox::calculateRadiusX() { return std::min((x - xMin), (xMax - x)); }
double WorkspaceBoundingBox::calculateRadiusY() { return std::min((y - yMin), (yMax - y)); }
void WorkspaceBoundingBox::normalizePosition(double x, double y) {
this->x /= x;
this->y /= y;
}
void WorkspaceBoundingBox::normalizePosition(double x, double y){
this.x /= x;
this.y /= y;
bool WorkspaceBoundingBox::containsPoint(double x, double y) {
if (x > this->xMax || x < this->xMin || y > yMax || y < yMin)
return false;
return true;
}
bool WorkspaceBoundingBox::containsPoint(double x, double y) {
if(x > this.xMax || x < this.xMin || y > yMax, || y < yMin)
return false;
return true;
}
\ No newline at end of file
} // namespace Algorithms
} // namespace Mantid
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment