diff --git a/Framework/Algorithms/CMakeLists.txt b/Framework/Algorithms/CMakeLists.txt
index 3a653bd17cbad0a5692b857deddd8863c572e721..217daf7a5c9e08a4125e223ed79229622f92002f 100644
--- a/Framework/Algorithms/CMakeLists.txt
+++ b/Framework/Algorithms/CMakeLists.txt
@@ -142,6 +142,7 @@ set ( SRC_FILES
src/He3TubeEfficiency.cpp
src/IQTransform.cpp
src/IdentifyNoisyDetectors.cpp
+ src/ImggTomographicReconstruction.cpp
src/IntegrateByComponent.cpp
src/Integration.cpp
src/InterpolatingRebin.cpp
@@ -268,6 +269,7 @@ set ( SRC_FILES
src/TimeAtSampleStrategyDirect.cpp
src/TimeAtSampleStrategyElastic.cpp
src/TimeAtSampleStrategyIndirect.cpp
+ src/Tomography/FBPTomopy.cpp
src/Transpose.cpp
src/UnGroupWorkspace.cpp
src/UnaryOperation.cpp
@@ -283,6 +285,12 @@ set ( SRC_FILES
src/XDataConverter.cpp
)
+# C as they originally come from a different project (tomopy)
+set ( C_SRC_FILES
+ src/Tomography/tomopy/fbp.c
+ src/Tomography/tomopy/utils.c
+)
+
set ( INC_FILES
inc/MantidAlgorithms/AbsorptionCorrection.h
inc/MantidAlgorithms/AddLogDerivative.h
@@ -428,6 +436,7 @@ set ( INC_FILES
inc/MantidAlgorithms/He3TubeEfficiency.h
inc/MantidAlgorithms/IQTransform.h
inc/MantidAlgorithms/IdentifyNoisyDetectors.h
+ inc/MantidAlgorithms/ImggTomographicReconstruction.h
inc/MantidAlgorithms/IntegrateByComponent.h
inc/MantidAlgorithms/Integration.h
inc/MantidAlgorithms/InterpolatingRebin.h
@@ -557,6 +566,9 @@ set ( INC_FILES
inc/MantidAlgorithms/TimeAtSampleStrategyDirect.h
inc/MantidAlgorithms/TimeAtSampleStrategyElastic.h
inc/MantidAlgorithms/TimeAtSampleStrategyIndirect.h
+ inc/MantidAlgorithms/Tomography/FBPTomopy.h
+ inc/MantidAlgorithms/Tomography/tomopy/fbp.h
+ inc/MantidAlgorithms/Tomography/tomopy/utils.h
inc/MantidAlgorithms/Transpose.h
inc/MantidAlgorithms/UnGroupWorkspace.h
inc/MantidAlgorithms/UnaryOperation.h
@@ -581,7 +593,7 @@ set(SRC_UNITY_IGNORE_FILES src/AlignDetectors.cpp
if(UNITY_BUILD)
include(UnityBuild)
- enable_unity_build(Algorithms SRC_FILES SRC_UNITY_IGNORE_FILES 10)
+ enable_unity_build(Algorithms SRC_FILES C_SRC_FILES SRC_UNITY_IGNORE_FILES 10)
endif(UNITY_BUILD)
set ( TEST_FILES
@@ -722,6 +734,7 @@ set ( TEST_FILES
He3TubeEfficiencyTest.h
IQTransformTest.h
IdentifyNoisyDetectorsTest.h
+ ImggTomographicReconstructionTest.h
IntegrateByComponentTest.h
IntegrationTest.h
InterpolatingRebinTest.h
@@ -832,6 +845,7 @@ set ( TEST_FILES
TimeAtSampleStrategyDirectTest.h
TimeAtSampleStrategyElasticTest.h
TimeAtSampleStrategyIndirectTest.h
+ Tomography/FBPTomopyTest.h
TransposeTest.h
UnGroupWorkspaceTest.h
UnaryOperationTest.h
@@ -848,7 +862,7 @@ set ( TEST_FILES
set ( TEST_PY_FILES NormaliseToUnityTest.py )
if (COVERALLS)
- foreach( loop_var ${SRC_FILES} ${INC_FILES})
+ foreach( loop_var ${SRC_FILES} ${C_SRC_FILES} ${INC_FILES})
set_property(GLOBAL APPEND PROPERTY COVERAGE_SRCS "${CMAKE_CURRENT_SOURCE_DIR}/${loop_var}")
endforeach(loop_var)
endif()
@@ -856,7 +870,7 @@ endif()
# Add a precompiled header where they are supported
enable_precompiled_headers ( inc/MantidAlgorithms/PrecompiledHeader.h SRC_FILES )
# Add the target for this directory
-add_library ( Algorithms ${SRC_FILES} ${INC_FILES})
+add_library ( Algorithms ${SRC_FILES} ${C_SRC_FILES} ${INC_FILES})
# Set the name of the generated library
set_target_properties ( Algorithms PROPERTIES OUTPUT_NAME MantidAlgorithms
COMPILE_DEFINITIONS "IN_MANTID_ALGORITHMS" )
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/ImggTomographicReconstruction.h b/Framework/Algorithms/inc/MantidAlgorithms/ImggTomographicReconstruction.h
new file mode 100644
index 0000000000000000000000000000000000000000..6edb67009036ed4738f71087988a8b14ea9cbe79
--- /dev/null
+++ b/Framework/Algorithms/inc/MantidAlgorithms/ImggTomographicReconstruction.h
@@ -0,0 +1,81 @@
+#ifndef MANTID_ALGORITHMS_IMGGTOMOGRAPHICRECONSTRUCTION_H_
+#define MANTID_ALGORITHMS_IMGGTOMOGRAPHICRECONSTRUCTION_H_
+
+#include "MantidAlgorithms/DllConfig.h"
+#include "MantidAPI/Algorithm.h"
+#include "MantidAPI/WorkspaceGroup_fwd.h"
+
+namespace Mantid {
+namespace Algorithms {
+
+/**
+ ImggTomographicReconstruction: reconstruct volumes from a sequence
+ of 2D projections using tomographic reconstruction methods.
+
+ Copyright © 2016 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
+ National Laboratory & European Spallation Source
+
+ 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 MANTID_ALGORITHMS_DLL ImggTomographicReconstruction final
+ : public API::Algorithm {
+public:
+ const std::string name() const override final;
+
+ int version() const override final;
+
+ const std::string category() const override final;
+
+ const std::string summary() const override final;
+
+private:
+ void init() override final;
+
+ void exec() override final;
+
+ bool processGroups() override final;
+
+ std::map<std::string, std::string> validateInputs() override;
+
+ std::unique_ptr<std::vector<float>>
+ prepareProjectionAngles(API::WorkspaceGroup_const_sptr wks, double minAngle,
+ double maxAngle) const;
+
+ std::unique_ptr<std::vector<float>>
+ prepareInputData(size_t totalSize, API::WorkspaceGroup_const_sptr wsg);
+
+ std::unique_ptr<std::vector<float>> prepareDataVol(size_t totalSize);
+
+ std::unique_ptr<std::vector<float>> prepareCenters(int cor, size_t totalSize);
+
+ size_t xSizeProjections(API::WorkspaceGroup_const_sptr wks) const;
+
+ size_t pSizeProjections(API::WorkspaceGroup_const_sptr wks) const;
+
+ size_t ySizeProjections(API::WorkspaceGroup_const_sptr wks) const;
+
+ API::WorkspaceGroup_sptr buildOutputWks(const std::vector<float> &dataVol,
+ size_t xsize, size_t ysize,
+ size_t sliceSize);
+};
+
+} // namespace Algorithms
+} // namespace Mantid
+
+#endif /* MANTID_ALGORITHMS_IMGGTOMOGRAPHICRECONSTRUCTION_H_ */
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/Tomography/FBPTomopy.h b/Framework/Algorithms/inc/MantidAlgorithms/Tomography/FBPTomopy.h
new file mode 100644
index 0000000000000000000000000000000000000000..0d295f82a82b2a2e9195e343c9b8031950720be0
--- /dev/null
+++ b/Framework/Algorithms/inc/MantidAlgorithms/Tomography/FBPTomopy.h
@@ -0,0 +1,45 @@
+#ifndef MANTID_ALGORITHMS_TOMOGRAPHY_FBPTOMOPY_H_
+#define MANTID_ALGORITHMS_TOMOGRAPHY_FBPTOMOPY_H_
+
+#include "MantidAlgorithms/DllConfig.h"
+
+namespace Mantid {
+namespace Algorithms {
+namespace Tomography {
+
+/**
+ Wrapper for C function(s) that implement tomopy reconstruction
+ methods.
+
+ Copyright © 2016 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
+ National Laboratory & European Spallation Source
+
+ 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>
+*/
+
+/// FBP - Filtered Back-Projection routine as implemented in tomopy
+void MANTID_ALGORITHMS_DLL FBPTomopy(const float *data, int dy, int dt, int dx,
+ const float *center, const float *theta,
+ float *recon, int ngridx, int ngridy);
+
+} // namespace Tomography
+} // namespace Algorithms
+} // namespace Mantid
+
+#endif /* MANTID_ALGORITHMS_TOMOGRAPHY_FBPTOMOPY_H_ */
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/Tomography/tomopy/fbp.h b/Framework/Algorithms/inc/MantidAlgorithms/Tomography/tomopy/fbp.h
new file mode 100644
index 0000000000000000000000000000000000000000..ade127406daee9f2e30484d34ec32badcf192fbc
--- /dev/null
+++ b/Framework/Algorithms/inc/MantidAlgorithms/Tomography/tomopy/fbp.h
@@ -0,0 +1,34 @@
+/**
+ fbp function from Tomopy - FBP reconstruction method
+
+ Copyright © 2016 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
+ National Laboratory & European Spallation Source
+
+ 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>
+*/
+#ifndef MANTID_ALGORITHMS_TOMOGRAPHY_TOMOPY_FBP_H_
+#define MANTID_ALGORITHMS_TOMOGRAPHY_TOMOPY_FBP_H_
+
+extern "C" {
+void fbp(const float *data, int dy, int dt, int dx, const float *center,
+ const float *theta, float *recon, int ngridx, int ngridy,
+ const char * /*fname*/, const float * /*filter_par*/);
+}
+
+#endif /* MANTID_ALGORITHMS_TOMOGRAPHY_TOMOPY_FBP_H_ */
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/Tomography/tomopy/utils.h b/Framework/Algorithms/inc/MantidAlgorithms/Tomography/tomopy/utils.h
new file mode 100644
index 0000000000000000000000000000000000000000..607160a3be2b50c29d7f6644b2c601fa160cd1ee
--- /dev/null
+++ b/Framework/Algorithms/inc/MantidAlgorithms/Tomography/tomopy/utils.h
@@ -0,0 +1,103 @@
+/**
+ Common functions for different reconstruction methods in Tomopy
+
+ Copyright © 2016 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
+ National Laboratory & European Spallation Source
+
+ 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>
+*/
+/**
+ * This file is substantially modified but originally based on the
+ * file utils.h from tomopy (available from
+ * https:github.com/tomopy/tomopy/blob/master/src/utils.h) which is:
+ */
+/*
+ Copyright (c) 2015, UChicago Argonne, LLC. All rights reserved.
+
+ Copyright 2015. UChicago Argonne, LLC. This software was produced
+ under U.S. Government contract DE-AC02-06CH11357 for Argonne National
+ Laboratory (ANL), which is operated by UChicago Argonne, LLC for the
+ U.S. Department of Energy. The U.S. Government has rights to use,
+ reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR
+ UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR
+ ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is
+ modified to produce derivative works, such modified software should
+ be clearly marked, so as not to confuse it with the version available
+ from ANL.
+
+ Additionally, redistribution and use in source and binary forms, with
+ or without modification, are permitted provided that the following
+ conditions are met:
+
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+ * Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in
+ the documentation and/or other materials provided with the
+ distribution.
+
+ * Neither the name of UChicago Argonne, LLC, Argonne National
+ Laboratory, ANL, the U.S. Government, nor the names of its
+ contributors may be used to endorse or promote products derived
+ from this software without specific prior written permission.
+
+ THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS
+ "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago
+ Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ POSSIBILITY OF SUCH DAMAGE.
+*/
+#ifndef MANTID_ALGORITHMS_TOMOGRAPHY_TOMOPY_UTILS_H_
+#define MANTID_ALGORITHMS_TOMOGRAPHY_TOMOPY_UTILS_H_
+
+void preprocessing(int ngridx, int ngridy, int dz, float center, float *mov,
+ float *gridx, float *gridy);
+
+int calc_quadrant(float theta_p);
+
+void calc_coords(int ngridx, int ngridy, float xi, float yi, float sin_p,
+ float cos_p, const float *gridx, const float *gridy,
+ float *coordx, float *coordy);
+
+void trim_coords(int ngridx, int ngridy, const float *coordx,
+ const float *coordy, const float *gridx, const float *gridy,
+ int *asize, float *ax, float *ay, int *bsize, float *bx,
+ float *by);
+
+void sort_intersections(int ind_condition, int asize, const float *ax,
+ const float *ay, int bsize, const float *bx,
+ const float *by, int *csize, float *coorx,
+ float *coory);
+
+void calc_dist(int ngridx, int ngridy, int csize, const float *coorx,
+ const float *coory, int *indi, float *dist);
+
+void calc_simdata(int s, int p, int d, int ngridx, int ngridy, int dt, int dx,
+ int csize, const int *indi, const float *dist,
+ const float *model, float *simdata);
+
+#endif /* MANTID_ALGORITHMS_TOMOGRAPHY_TOMOPY_UTILS_H_ */
diff --git a/Framework/Algorithms/src/ImggTomographicReconstruction.cpp b/Framework/Algorithms/src/ImggTomographicReconstruction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e6c8ac35556d2f529e5455c5d297677b7d4f6cb7
--- /dev/null
+++ b/Framework/Algorithms/src/ImggTomographicReconstruction.cpp
@@ -0,0 +1,412 @@
+#include "MantidAlgorithms/ImggTomographicReconstruction.h"
+#include "MantidAlgorithms/Tomography/FBPTomopy.h"
+#include "MantidAPI/AlgorithmManager.h"
+#include "MantidAPI/CommonBinsValidator.h"
+#include "MantidAPI/FileProperty.h"
+#include "MantidAPI/HistogramValidator.h"
+#include "MantidAPI/MatrixWorkspace.h"
+#include "MantidAPI/WorkspaceUnitValidator.h"
+#include "MantidAPI/WorkspaceFactory.h"
+#include "MantidAPI/WorkspaceGroup.h"
+#include "MantidDataObjects/Workspace2D.h"
+#include "MantidKernel/BoundedValidator.h"
+#include "MantidKernel/CompositeValidator.h"
+#include "MantidKernel/ListValidator.h"
+#include "MantidKernel/make_unique.h"
+
+namespace Mantid {
+namespace Algorithms {
+
+// Register the algorithm into the AlgorithmFactory
+DECLARE_ALGORITHM(ImggTomographicReconstruction)
+
+namespace {
+// Just to have explicit double => float casts in std::copy/transform
+struct DoubleToFloatStd {
+ float operator()(const double &dblValue) const {
+ return static_cast<float>(dblValue);
+ }
+};
+}
+
+//----------------------------------------------------------------------------------------------
+
+/// Algorithms name for identification. @see Algorithm::name
+const std::string ImggTomographicReconstruction::name() const {
+ return "ImggTomographicReconstruction";
+}
+
+/// Algorithm's version for identification. @see Algorithm::version
+int ImggTomographicReconstruction::version() const { return 1; }
+
+/// Algorithm's category for identification. @see Algorithm::category
+const std::string ImggTomographicReconstruction::category() const {
+ return "Diffraction\\Imaging;Diffraction\\Tomography";
+}
+
+/// Algorithm's summary for use in the GUI and help. @see Algorithm::summary
+const std::string ImggTomographicReconstruction::summary() const {
+ return "Reconstruct a 3D volume from 2D imaging projection data";
+}
+
+namespace {
+const std::string PROP_INPUT_WS = "InputWorkspace";
+const std::string PROP_METHOD = "Method";
+const std::string PROP_OUTPUT_WS = "OutputWorkspace";
+const std::string PROP_COR = "CenterOfRotation";
+const std::string PROP_RELAXATION_PARAM = "RelaxationParameter";
+const std::string PROP_MAX_CORES = "MaximumCores";
+const std::string PROP_MIN_PROJ_ANGLE = "MinProjectionAngle";
+const std::string PROP_MAX_PROJ_ANGLE = "MaxProjectionAngle";
+}
+
+//----------------------------------------------------------------------------------------------
+/** Initialize the algorithm's properties.
+ */
+void ImggTomographicReconstruction::init() {
+
+ declareProperty(
+ Kernel::make_unique<API::WorkspaceProperty<API::Workspace>>(
+ PROP_INPUT_WS, "", Kernel::Direction::Input),
+ "Group of workspace holding images (with one spectrum per pixel row).");
+
+ std::vector<std::string> methods{"FBP (tomopy)"};
+ declareProperty(
+ PROP_METHOD, methods.front(),
+ boost::make_shared<Kernel::ListValidator<std::string>>(methods),
+ "Reconstruction method", Kernel::Direction::Input);
+
+ declareProperty(
+ Kernel::make_unique<API::WorkspaceProperty<API::WorkspaceGroup>>(
+ PROP_OUTPUT_WS, "", Kernel::Direction::Output),
+ "Output reconstructed volume, as a group of workspaces where "
+ "each workspace holds one slice of the volume.");
+
+ auto zeroOrPosInt = boost::make_shared<Kernel::BoundedValidator<int>>();
+ zeroOrPosInt->setLower(-1);
+ declareProperty(PROP_COR, -1, zeroOrPosInt,
+ "Center of rotation for the reconstruction (in pixels).");
+
+ auto zeroOrPosDbl = boost::make_shared<Kernel::BoundedValidator<double>>();
+ zeroOrPosDbl->setLower(0.0);
+ declareProperty(PROP_RELAXATION_PARAM, 0.5, zeroOrPosDbl,
+ "Relaxation parameter for the reconstruction method.");
+
+ zeroOrPosInt->setLower(0);
+ declareProperty(PROP_MAX_CORES, 0, zeroOrPosInt,
+ "Maximum number of cores to use for parallel runs. Leave "
+ "empty to use all available cores.");
+
+ declareProperty(PROP_MIN_PROJ_ANGLE, 0.0, zeroOrPosDbl,
+ "Minimum projection angle.");
+
+ declareProperty(PROP_MAX_PROJ_ANGLE, 180.0, zeroOrPosDbl,
+ "Maximum projection angle (assuming a uniform angle increase "
+ "from first to last projection.");
+}
+
+std::map<std::string, std::string>
+ImggTomographicReconstruction::validateInputs() {
+ std::map<std::string, std::string> result;
+
+ API::Workspace_const_sptr inWks = getProperty(PROP_INPUT_WS);
+ API::WorkspaceGroup_const_sptr inGrp =
+ boost::dynamic_pointer_cast<const API::WorkspaceGroup>(inWks);
+ if (!inGrp) {
+ result[PROP_INPUT_WS] = "The current version of this algorithm only "
+ "supports input workspaces of type WorkspaceGroup";
+ } else {
+ if (inGrp->size() < 2) {
+ result[PROP_INPUT_WS] = "The input workspace must have at least two "
+ "workspaces (projection images)";
+ }
+
+ auto first = inGrp->getItem(0);
+ auto fwks = boost::dynamic_pointer_cast<API::MatrixWorkspace>(first);
+ if (!fwks) {
+ result[PROP_INPUT_WS] =
+ "Unable to get a matrix workspace from the first "
+ "item of the input workspace group " +
+ first->getTitle() +
+ ". It must contain workspaces of type MatrixWorkspace";
+ } else {
+ int cor = getProperty(PROP_COR);
+ size_t bsize = fwks->blocksize();
+ if (cor < 0 || cor >= static_cast<int>(bsize)) {
+ result[PROP_COR] =
+ "The center of rotation must be between 0 and the "
+ "number of columns in the input projection images (0 to " +
+ std::to_string(bsize - 1) + ")";
+ }
+ }
+ // Not validating requirements on all input workspaces here (there could be
+ // many)
+ }
+
+ double minAngle = getProperty(PROP_MIN_PROJ_ANGLE);
+ double maxAngle = getProperty(PROP_MAX_PROJ_ANGLE);
+
+ if (minAngle >= maxAngle) {
+ result[PROP_MIN_PROJ_ANGLE] = PROP_MIN_PROJ_ANGLE +
+ " cannot be equal to or lower than " +
+ PROP_MAX_PROJ_ANGLE;
+ }
+ return result;
+}
+
+//----------------------------------------------------------------------------------------------
+/** Execute the algorithm.
+ */
+void ImggTomographicReconstruction::exec() {
+ throw std::runtime_error(
+ "This algorithm cannot be executed with a single workspace as input");
+}
+
+bool ImggTomographicReconstruction::processGroups() {
+ API::Workspace_const_sptr inWks = getProperty("InputWorkspace");
+
+ API::WorkspaceGroup_const_sptr wks =
+ boost::dynamic_pointer_cast<const API::WorkspaceGroup>(inWks);
+
+ if (!wks) {
+ g_log.error(
+ "Could not retrieve the input workspace as a workspace group: ");
+ return false;
+ }
+
+ // TODO: apply validators here on every input image/workspace
+ for (size_t idx = 0; idx < wks->size(); idx++) {
+ auto item = wks->getItem(0);
+ auto mWS = boost::dynamic_pointer_cast<API::MatrixWorkspace>(item);
+ if (!mWS) {
+ throw std::runtime_error("Unable to get a matrix workspace from the "
+ "element of te workspace group with title " +
+ item->getTitle());
+ }
+
+ auto wsValidator = boost::make_shared<Kernel::CompositeValidator>();
+ wsValidator->add<API::CommonBinsValidator>();
+ wsValidator->add<API::HistogramValidator>();
+ // Probably we won't need this:
+ // wsValidator->add<API::WorkspaceUnitValidator>("Label");
+ const std::string validation = wsValidator->isValid(mWS);
+ if (validation != "") {
+ throw std::runtime_error(
+ "Validation of input image / matrix workspace failed: " + validation);
+ }
+ }
+
+ double minAngle = getProperty(PROP_MIN_PROJ_ANGLE);
+ double maxAngle = getProperty(PROP_MAX_PROJ_ANGLE);
+ auto angles = prepareProjectionAngles(wks, minAngle, maxAngle);
+
+ // these values are expected as 'int' by the tompy routines
+ const int ysize = static_cast<int>(ySizeProjections(wks));
+ const int projSize = static_cast<int>(angles->size());
+ const int xsize = static_cast<int>(xSizeProjections(wks));
+ // total size of input data in voxels
+ size_t totalInSize = ysize * projSize * xsize;
+ // total size of the reconstructed volume
+ size_t totalReconSize = ysize * ysize * xsize;
+
+ auto inVol = prepareInputData(totalInSize, wks);
+ auto reconVol = prepareDataVol(totalReconSize);
+
+ int cor = getProperty(PROP_COR);
+ auto centers = prepareCenters(cor, ysize);
+
+ Mantid::Algorithms::Tomography::FBPTomopy(
+ inVol->data(), ysize, projSize, xsize, centers->data(), angles->data(),
+ reconVol->data(), xsize, ysize);
+
+ size_t expectedVox = ysize * ysize * xsize;
+ if (reconVol->size() != expectedVox) {
+ std::stringstream stream;
+ stream << std::string("The reconstructed volume data block does not "
+ "have the expected dimensions. It has ")
+ << reconVol->size() << " voxels, whereas I was expecting: " << ysize
+ << " slices by " << ysize << " rows by " << xsize
+ << " columns = " << expectedVox << " voxels in total";
+ throw std::runtime_error(stream.str());
+ }
+
+ const auto outputGrp = buildOutputWks(*reconVol, xsize, ysize, ysize);
+ setProperty(PROP_OUTPUT_WS, outputGrp);
+
+ g_log.notice() << "Finished reconstruction of volume from workspace "
+ << wks->getTitle() << " with " << projSize
+ << " input projections, " << ysize << " rows by " << xsize
+ << " columns." << std::endl;
+
+ // This is an ugly workaround for the issue that processGroups()
+ // does not fully finish the algorithm:
+ // https://github.com/mantidproject/mantid/issues/11361
+ // A similar workaround is used in CompareWorkspaces and others
+ // Store output workspace in AnalysisDataService and finish execution
+ if (!isChild())
+ this->store();
+ setExecuted(true);
+ notificationCenter().postNotification(
+ new FinishedNotification(this, this->isExecuted()));
+
+ return true;
+}
+
+std::unique_ptr<std::vector<float>>
+ImggTomographicReconstruction::prepareProjectionAngles(
+ API::WorkspaceGroup_const_sptr wks, double minAngle,
+ double maxAngle) const {
+
+ auto angles = Kernel::make_unique<std::vector<float>>(wks->size());
+
+ auto vec = *angles;
+ double factor = (maxAngle - minAngle);
+ for (size_t idx = 0; idx < wks->size(); ++idx) {
+ vec[idx] = static_cast<float>(minAngle +
+ factor * static_cast<double>(idx) /
+ static_cast<double>(idx - 1));
+ }
+
+ return angles;
+}
+
+std::unique_ptr<std::vector<float>>
+ImggTomographicReconstruction::prepareInputData(
+ size_t totalSize, API::WorkspaceGroup_const_sptr wksg) {
+ auto data = prepareDataVol(totalSize);
+
+ if (!wksg || 0 == wksg->size())
+ return data;
+
+ auto first = wksg->getItem(0);
+ auto fwks = boost::dynamic_pointer_cast<API::MatrixWorkspace>(first);
+ if (!fwks) {
+ throw std::runtime_error(
+ "Unable to get a matrix workspace from the first "
+ "item of the input workspace group " +
+ first->getTitle() +
+ ". It must contain workspaces of type MatrixWorkspace");
+ }
+
+ size_t ysize = fwks->getNumberHistograms();
+ size_t xsize = fwks->blocksize();
+ const size_t oneSliceSize = xsize * ysize;
+
+ PARALLEL_FOR_NO_WSP_CHECK()
+ for (int slice = 0; slice < static_cast<int>(wksg->size()); ++slice) {
+ size_t startSlice = slice * oneSliceSize;
+ for (size_t row = 0; row < ysize; ++row) {
+ const Mantid::API::ISpectrum *specRow = fwks->getSpectrum(row);
+ const auto &dataY = specRow->readY();
+ size_t startRow = startSlice + row * ysize;
+ // MSVC will produce C4244 warnings in <xutility> (double=>float
+ // converstion)
+ // std::copy(dataY.begin(), dataY.end(), data->begin() + startRow);
+ std::transform(dataY.begin(), dataY.end(), data->begin() + startRow,
+ DoubleToFloatStd());
+ }
+ }
+
+ return data;
+}
+
+std::unique_ptr<std::vector<float>>
+ImggTomographicReconstruction::prepareDataVol(size_t totalSize) {
+ return Kernel::make_unique<std::vector<float>>(totalSize);
+}
+
+std::unique_ptr<std::vector<float>>
+ImggTomographicReconstruction::prepareCenters(int cor, size_t totalSize) {
+ auto centers = Kernel::make_unique<std::vector<float>>(totalSize);
+
+ for (auto &cnt : *centers) {
+ cnt = static_cast<float>(cor);
+ }
+
+ return centers;
+}
+
+size_t ImggTomographicReconstruction::xSizeProjections(
+ API::WorkspaceGroup_const_sptr wks) const {
+ auto first = wks->getItem(0);
+
+ auto wksItem = boost::dynamic_pointer_cast<API::MatrixWorkspace>(first);
+ if (!wksItem) {
+ throw std::runtime_error("Unable to get a matrix workspace from the first "
+ "item of the workspace group " +
+ first->getTitle());
+ }
+
+ return wksItem->getNumberHistograms();
+}
+
+size_t ImggTomographicReconstruction::pSizeProjections(
+ API::WorkspaceGroup_const_sptr wks) const {
+ return wks->size();
+}
+
+size_t ImggTomographicReconstruction::ySizeProjections(
+ API::WorkspaceGroup_const_sptr wks) const {
+ auto first = wks->getItem(0);
+
+ auto wksItem = boost::dynamic_pointer_cast<API::MatrixWorkspace>(first);
+ if (!wksItem) {
+ throw std::runtime_error("Unable to get a matrix workspace from the first "
+ "item of the workspace group " +
+ first->getTitle());
+ }
+
+ return wksItem->blocksize();
+}
+
+/**
+ * Transfer data from a numpy/tomopy style data volume to a workspace
+ * group of MatrixWorkspaces. This will typically have ysize slices
+ * (images) of dimensions ysize rows by xsize columns.
+
+ *
+ * @param dataVol a 3D data volume that may have been generated as
+ * output from a reconstruction
+ * @param xsize image x dimension or number of columns (bins)
+ * @param ysize image y dimension or number of rows (spectra)
+ * @param sliceSize number of slices
+ *
+ * @return a workspace group with reconstruction slices
+ */
+API::WorkspaceGroup_sptr
+ImggTomographicReconstruction::buildOutputWks(const std::vector<float> &dataVol,
+ size_t xsize, size_t ysize,
+ size_t sliceSize) {
+
+ // auto wsGroup = boost::make_shared<API::WorkspaceGroup>();
+ auto wsGroup = API::WorkspaceGroup_sptr(new API::WorkspaceGroup());
+ wsGroup->setTitle("Reconstructed volume from imaging projection data");
+
+ const size_t oneSliceSize = xsize * ysize;
+ PARALLEL_FOR_NO_WSP_CHECK()
+ for (int slice = 0; slice < static_cast<int>(sliceSize); ++slice) {
+ // individual slices as Workspace2D/MatrixWorkspace
+ DataObjects::Workspace2D_sptr sliceWS =
+ boost::dynamic_pointer_cast<DataObjects::Workspace2D>(
+ API::WorkspaceFactory::Instance().create("Workspace2D", ysize,
+ xsize + 1, xsize));
+ size_t startSlice = slice * oneSliceSize;
+ for (size_t row = 0; row < ysize; ++row) {
+ Mantid::API::ISpectrum *specRow = sliceWS->getSpectrum(row);
+ auto &dataX = specRow->dataX();
+ std::fill(dataX.begin(), dataX.end(), static_cast<double>(row));
+
+ size_t startRow = startSlice + row * ysize;
+ size_t endRow = startRow + xsize;
+ auto &dataY = specRow->dataY();
+ std::transform(dataVol.begin() + startRow, dataVol.begin() + endRow,
+ dataY.begin(), DoubleToFloatStd());
+ }
+ wsGroup->addWorkspace(sliceWS);
+ }
+
+ return wsGroup;
+}
+
+} // namespace Algorithms
+} // namespace Mantid
diff --git a/Framework/Algorithms/src/Tomography/FBPTomopy.cpp b/Framework/Algorithms/src/Tomography/FBPTomopy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..807ece8f2922621c7f88bcf1332627d10141e76c
--- /dev/null
+++ b/Framework/Algorithms/src/Tomography/FBPTomopy.cpp
@@ -0,0 +1,15 @@
+#include "MantidAlgorithms/Tomography/FBPTomopy.h"
+#include "MantidAlgorithms/Tomography/tomopy/fbp.h"
+
+namespace Mantid {
+namespace Algorithms {
+namespace Tomography {
+
+void FBPTomopy(const float *data, int dy, int dt, int dx, const float *center,
+ const float *theta, float *recon, int ngridx, int ngridy) {
+ fbp(data, dy, dt, dx, center, theta, recon, ngridx, ngridy, nullptr, nullptr);
+}
+
+} // namespace Tomography
+} // namespace Algorithms
+} // namespace Mantid
diff --git a/Framework/Algorithms/src/Tomography/tomopy/fbp.c b/Framework/Algorithms/src/Tomography/tomopy/fbp.c
new file mode 100644
index 0000000000000000000000000000000000000000..b0451a11807ab74b2091d415604b5a2c13189552
--- /dev/null
+++ b/Framework/Algorithms/src/Tomography/tomopy/fbp.c
@@ -0,0 +1,160 @@
+/**
+ This file is originally based on the implementation of the FBP (Filtered
+ Back-Projection) method of tomopy, available from
+ https:github.com/tomopy/tomopy/blob/master/src/fbp.c, which is:
+ */
+/*
+ Copyright (c) 2015, UChicago Argonne, LLC. All rights reserved.
+
+ Copyright 2015. UChicago Argonne, LLC. This software was produced
+ under U.S. Government contract DE-AC02-06CH11357 for Argonne National
+ Laboratory (ANL), which is operated by UChicago Argonne, LLC for the
+ U.S. Department of Energy. The U.S. Government has rights to use,
+ reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR
+ UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR
+ ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is
+ modified to produce derivative works, such modified software should
+ be clearly marked, so as not to confuse it with the version available
+ from ANL.
+
+ Additionally, redistribution and use in source and binary forms, with
+ or without modification, are permitted provided that the following
+ conditions are met:
+
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+ * Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in
+ the documentation and/or other materials provided with the
+ distribution.
+
+ * Neither the name of UChicago Argonne, LLC, Argonne National
+ Laboratory, ANL, the U.S. Government, nor the names of its
+ contributors may be used to endorse or promote products derived
+ from this software without specific prior written permission.
+
+ THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS
+ "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago
+ Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#include "MantidAlgorithms/Tomography/tomopy/utils.h"
+
+#include <assert.h>
+#include <math.h>
+#include <stdlib.h>
+
+void fbp(const float *data, int dy, int dt, int dx, const float *center,
+ const float *theta, float *recon, int ngridx, int ngridy,
+ const char *fname, const float *filter_par) {
+
+ int s, p, d, n;
+ int quadrant;
+ float theta_p, sin_p, cos_p;
+ float mov, xi, yi;
+ int asize, bsize, csize;
+ int ind_data, ind_recon;
+
+ float *gridx, *gridy, *coordx, *coordy, *ax, *ay, *bx, *by, *coorx, *coory,
+ *dist;
+ int *indi;
+
+ gridx = (float *)malloc((size_t)(ngridx + 1) * sizeof(float));
+ gridy = (float *)malloc((size_t)(ngridy + 1) * sizeof(float));
+ coordx = (float *)malloc((size_t)(ngridy + 1) * sizeof(float));
+ coordy = (float *)malloc((size_t)(ngridx + 1) * sizeof(float));
+ ax = (float *)malloc((size_t)(ngridx + ngridy) * sizeof(float));
+ ay = (float *)malloc((size_t)(ngridx + ngridy) * sizeof(float));
+ bx = (float *)malloc((size_t)(ngridx + ngridy) * sizeof(float));
+ by = (float *)malloc((size_t)(ngridx + ngridy) * sizeof(float));
+ coorx = (float *)malloc((size_t)(ngridx + ngridy) * sizeof(float));
+ coory = (float *)malloc((size_t)(ngridx + ngridy) * sizeof(float));
+ dist = (float *)malloc((size_t)(ngridx + ngridy) * sizeof(float));
+ indi = (int *)malloc((size_t)(ngridx + ngridy) * sizeof(int));
+
+ /* unused arguments: */
+ (void)fname;
+ (void)filter_par;
+
+ assert(coordx != NULL && coordy != NULL && ax != NULL && ay != NULL &&
+ by != NULL && bx != NULL && coorx != NULL && coory != NULL &&
+ dist != NULL && indi != NULL);
+
+ /* For each slice */
+ for (s = 0; s < dy; s++) {
+ preprocessing(ngridx, ngridy, dx, center[s], &mov, gridx,
+ gridy); /* Outputs: mov, gridx, gridy */
+
+ /* For each projection angle */
+ for (p = 0; p < dt; p++) {
+ /*
+ Calculate the sin and cos values
+ of the projection angle and find
+ at which quadrant on the cartesian grid.
+ */
+ theta_p = (float)fmod(theta[p], 2 * M_PI);
+ quadrant = calc_quadrant(theta_p);
+ sin_p = sinf(theta_p);
+ cos_p = cosf(theta_p);
+
+ /* For each detector pixel */
+ for (d = 0; d < dx; d++) {
+ /* Calculate coordinates */
+ xi = (float)(-ngridx - ngridy);
+ yi = (float)(1 - dx) / 2.0f + (float)d + mov;
+ calc_coords(ngridx, ngridy, xi, yi, sin_p, cos_p, gridx, gridy, coordx,
+ coordy);
+
+ /* Merge the (coordx, gridy) and (gridx, coordy) */
+ trim_coords(ngridx, ngridy, coordx, coordy, gridx, gridy, &asize, ax,
+ ay, &bsize, bx, by);
+
+ /*
+ Sort the array of intersection points (ax, ay) and
+ (bx, by). The new sorted intersection points are
+ stored in (coorx, coory). Total number of points
+ are csize.
+ */
+ sort_intersections(quadrant, asize, ax, ay, bsize, bx, by, &csize,
+ coorx, coory);
+
+ /*
+ Calculate the distances (dist) between the
+ intersection points (coorx, coory). Find the
+ indices of the pixels on the reconstruction grid.
+ */
+ calc_dist(ngridx, ngridy, csize, coorx, coory, indi, dist);
+
+ /* Update */
+ ind_recon = s * ngridx * ngridy;
+ ind_data = d + p * dx + s * dt * dx;
+ for (n = 0; n < csize - 1; n++) {
+ recon[indi[n] + ind_recon] += data[ind_data] * dist[n];
+ }
+ }
+ }
+ }
+
+ free(gridx);
+ free(gridy);
+ free(coordx);
+ free(coordy);
+ free(ax);
+ free(ay);
+ free(bx);
+ free(by);
+ free(coorx);
+ free(coory);
+ free(dist);
+ free(indi);
+}
diff --git a/Framework/Algorithms/src/Tomography/tomopy/utils.c b/Framework/Algorithms/src/Tomography/tomopy/utils.c
new file mode 100644
index 0000000000000000000000000000000000000000..6b5248d5bac5f037c5203170a183ab5218b05ed0
--- /dev/null
+++ b/Framework/Algorithms/src/Tomography/tomopy/utils.c
@@ -0,0 +1,207 @@
+/**
+ This file is originally based on the implementation of the FBP (Filtered
+ Back-Projection) method of tomopy, available from
+ https:github.com/tomopy/tomopy/blob/master/src/fbp.c, which is:
+ */
+/*
+ Copyright (c) 2015, UChicago Argonne, LLC. All rights reserved.
+
+ Copyright 2015. UChicago Argonne, LLC. This software was produced
+ under U.S. Government contract DE-AC02-06CH11357 for Argonne National
+ Laboratory (ANL), which is operated by UChicago Argonne, LLC for the
+ U.S. Department of Energy. The U.S. Government has rights to use,
+ reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR
+ UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR
+ ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is
+ modified to produce derivative works, such modified software should
+ be clearly marked, so as not to confuse it with the version available
+ from ANL.
+
+ Additionally, redistribution and use in source and binary forms, with
+ or without modification, are permitted provided that the following
+ conditions are met:
+
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+ * Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in
+ the documentation and/or other materials provided with the
+ distribution.
+
+ * Neither the name of UChicago Argonne, LLC, Argonne National
+ Laboratory, ANL, the U.S. Government, nor the names of its
+ contributors may be used to endorse or promote products derived
+ from this software without specific prior written permission.
+
+ THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS
+ "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago
+ Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#include "MantidAlgorithms/Tomography/tomopy/utils.h"
+
+#include <math.h>
+
+void preprocessing(int ry, int rz, int num_pixels, float center, float *mov,
+ float *gridx, float *gridy) {
+ int i;
+
+ for (i = 0; i <= ry; i++) {
+ gridx[i] = ((float)-ry / 2.f + (float)i);
+ }
+
+ for (i = 0; i <= rz; i++) {
+ gridy[i] = ((float)-rz / 2.f + (float)i);
+ }
+
+ *mov = ((float)num_pixels - 1) / 2.0f - center;
+ if (*mov - floor(*mov) < 0.01f) {
+ *mov += 0.01f;
+ }
+ *mov += 0.5f;
+}
+
+int calc_quadrant(float theta_p) {
+ int quadrant;
+ if ((theta_p >= 0 && theta_p < M_PI / 2) ||
+ (theta_p >= M_PI && theta_p < 3 * M_PI / 2) ||
+ (theta_p >= -M_PI && theta_p < -M_PI / 2) ||
+ (theta_p >= -2 * M_PI && theta_p < -3 * M_PI / 2)) {
+ quadrant = 1;
+ } else {
+ quadrant = 0;
+ }
+ return quadrant;
+}
+
+void calc_coords(int ry, int rz, float xi, float yi, float sin_p, float cos_p,
+ const float *gridx, const float *gridy, float *coordx,
+ float *coordy) {
+ float srcx, srcy, detx, dety;
+ float slope, islope;
+ int n;
+
+ srcx = xi * cos_p - yi * sin_p;
+ srcy = xi * sin_p + yi * cos_p;
+ detx = -xi * cos_p - yi * sin_p;
+ dety = -xi * sin_p + yi * cos_p;
+
+ slope = (srcy - dety) / (srcx - detx);
+ islope = 1 / slope;
+ for (n = 0; n <= ry; n++) {
+ coordy[n] = slope * (gridx[n] - srcx) + srcy;
+ }
+ for (n = 0; n <= rz; n++) {
+ coordx[n] = islope * (gridy[n] - srcy) + srcx;
+ }
+}
+
+void trim_coords(int ry, int rz, const float *coordx, const float *coordy,
+ const float *gridx, const float *gridy, int *asize, float *ax,
+ float *ay, int *bsize, float *bx, float *by) {
+ int n;
+
+ *asize = 0;
+ *bsize = 0;
+ for (n = 0; n <= rz; n++) {
+ if (coordx[n] >= gridx[0] + 1e-2) {
+ if (coordx[n] <= gridx[ry] - 1e-2) {
+ ax[*asize] = coordx[n];
+ ay[*asize] = gridy[n];
+ (*asize)++;
+ }
+ }
+ }
+ for (n = 0; n <= ry; n++) {
+ if (coordy[n] >= gridy[0] + 1e-2) {
+ if (coordy[n] <= gridy[rz] - 1e-2) {
+ bx[*bsize] = gridx[n];
+ by[*bsize] = coordy[n];
+ (*bsize)++;
+ }
+ }
+ }
+}
+
+void sort_intersections(int ind_condition, int asize, const float *ax,
+ const float *ay, int bsize, const float *bx,
+ const float *by, int *csize, float *coorx,
+ float *coory) {
+ int i = 0, j = 0, k = 0;
+ int a_ind;
+ while (i < asize && j < bsize) {
+ a_ind = (ind_condition) ? i : (asize - 1 - i);
+ if (ax[a_ind] < bx[j]) {
+ coorx[k] = ax[a_ind];
+ coory[k] = ay[a_ind];
+ i++;
+ k++;
+ } else {
+ coorx[k] = bx[j];
+ coory[k] = by[j];
+ j++;
+ k++;
+ }
+ }
+ while (i < asize) {
+ a_ind = (ind_condition) ? i : (asize - 1 - i);
+ coorx[k] = ax[a_ind];
+ coory[k] = ay[a_ind];
+ i++;
+ k++;
+ }
+ while (j < bsize) {
+ coorx[k] = bx[j];
+ coory[k] = by[j];
+ j++;
+ k++;
+ }
+ *csize = asize + bsize;
+}
+
+void calc_dist(int ry, int rz, int csize, const float *coorx,
+ const float *coory, int *indi, float *dist) {
+ int n;
+
+ for (n = 0; n < csize - 1; n++) {
+ int i1, i2;
+ float x1, x2;
+ float diffx, diffy, midx, midy;
+ int indx, indy;
+
+ diffx = coorx[n + 1] - coorx[n];
+ diffy = coory[n + 1] - coory[n];
+ dist[n] = (float)sqrt(diffx * diffx + diffy * diffy);
+ midx = (coorx[n + 1] + coorx[n]) / 2;
+ midy = (coory[n + 1] + coory[n]) / 2;
+ x1 = midx + (float)ry / 2.f;
+ x2 = midy + (float)rz / 2.f;
+ i1 = (int)(midx + ry / 2.);
+ i2 = (int)(midy + rz / 2.);
+ indx = i1 - (i1 > x1);
+ indy = i2 - (i2 > x2);
+ indi[n] = indy + (indx * rz);
+ }
+}
+
+void calc_simdata(int s, int p, int d, int ry, int rz, int dt, int dx,
+ int csize, const int *indi, const float *dist,
+ const float *model, float *simdata) {
+ int n;
+
+ int index_model = s * ry * rz;
+ int index_data = d + p * dx + s * dt * dx;
+ for (n = 0; n < csize - 1; n++) {
+ simdata[index_data] += model[indi[n] + index_model] * dist[n];
+ }
+}
diff --git a/Framework/Algorithms/test/ImggTomographicReconstructionTest.h b/Framework/Algorithms/test/ImggTomographicReconstructionTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..a6d196ddc4b9ea0f997c9b2df42ba3088a1536a6
--- /dev/null
+++ b/Framework/Algorithms/test/ImggTomographicReconstructionTest.h
@@ -0,0 +1,286 @@
+#ifndef MANTID_ALGORITHMS_IMGGTOMOGRAPHICRECONSTRUCTIONTEST_H_
+#define MANTID_ALGORITHMS_IMGGTOMOGRAPHICRECONSTRUCTIONTEST_H_
+
+#include <cxxtest/TestSuite.h>
+
+#include "MantidAlgorithms/ImggTomographicReconstruction.h"
+#include "MantidAPI/AlgorithmManager.h"
+#include "MantidKernel/Exception.h"
+
+#include "MantidTestHelpers/WorkspaceCreationHelper.h"
+
+#include <Poco/File.h>
+
+using Mantid::Algorithms::ImggTomographicReconstruction;
+using namespace Mantid;
+
+class ImggTomographicReconstructionTest : 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 ImggTomographicReconstructionTest *createSuite() {
+ return new ImggTomographicReconstructionTest();
+ }
+
+ static void destroySuite(ImggTomographicReconstructionTest *suite) {
+ delete suite;
+ }
+
+ void test_init() {
+ ImggTomographicReconstruction alg;
+ TS_ASSERT_THROWS_NOTHING(alg.initialize());
+ TS_ASSERT(alg.isInitialized());
+
+ double relax;
+ TS_ASSERT_THROWS_NOTHING(relax = alg.getProperty("RelaxationParameter"));
+ TS_ASSERT_EQUALS(relax, 0.5);
+ }
+
+ void test_errors_options() {
+ auto alg = API::AlgorithmManager::Instance().create(
+ "ImggTomographicReconstruction");
+
+ TS_ASSERT_THROWS_NOTHING(
+ alg->setPropertyValue("OutputWorkspace", "_unused_for_child"));
+
+ TS_ASSERT_THROWS(
+ alg->setPropertyValue("BitDepth", "this_is_wrong_you_must_fail"),
+ std::runtime_error);
+ }
+
+ void test_exec_fails_inexistent_workspace() {
+ ImggTomographicReconstruction alg;
+ TS_ASSERT_THROWS_NOTHING(alg.initialize());
+ TS_ASSERT(alg.isInitialized());
+ TS_ASSERT_THROWS(
+ alg.setPropertyValue("InputWorkspace", "inexistent_workspace_fails"),
+ std::invalid_argument);
+
+ TS_ASSERT_THROWS(alg.execute(), std::runtime_error);
+ TS_ASSERT(!alg.isExecuted());
+ }
+
+ void test_exec_fails_wrong_workspace() {
+ API::MatrixWorkspace_sptr a =
+ WorkspaceCreationHelper::CreateWorkspaceSingleValue(3);
+
+ ImggTomographicReconstruction alg;
+ TS_ASSERT_THROWS_NOTHING(alg.initialize());
+ TS_ASSERT(alg.isInitialized());
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("InputWorkspace", a));
+
+ TS_ASSERT_THROWS(alg.execute(), std::runtime_error);
+ TS_ASSERT(!alg.isExecuted());
+
+ API::MatrixWorkspace_sptr wsSingle =
+ WorkspaceCreationHelper::Create2DWorkspace(10, 10);
+
+ ImggTomographicReconstruction algTwo;
+ TS_ASSERT_THROWS_NOTHING(algTwo.initialize());
+ TS_ASSERT(algTwo.isInitialized());
+ TS_ASSERT_THROWS_NOTHING(algTwo.setProperty("InputWorkspace", wsSingle));
+
+ TS_ASSERT_THROWS(algTwo.execute(), std::runtime_error);
+ TS_ASSERT(!algTwo.isExecuted());
+ }
+
+ void test_exec_fails_single_proj() {
+ const std::string projectionsGrpName("only_one_projection");
+ API::WorkspaceGroup_sptr projectionsGrp =
+ WorkspaceCreationHelper::CreateWorkspaceGroup(1, 4, 4,
+ projectionsGrpName);
+
+ ImggTomographicReconstruction alg;
+ TS_ASSERT_THROWS_NOTHING(alg.initialize());
+ TS_ASSERT(alg.isInitialized());
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("InputWorkspace", projectionsGrp));
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("MinProjectionAngle", 0.0));
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("MaxProjectionAngle", 260.0));
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("RelaxationParameter", 1.25));
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("CenterOfRotation", 4));
+
+ TS_ASSERT_THROWS(alg.execute(), std::runtime_error);
+ TS_ASSERT(!alg.isExecuted());
+ }
+
+ void test_exec_fails_wrong_center() {
+ const std::string projectionsGrpName("only_two_small_projections");
+ API::WorkspaceGroup_sptr projectionsGrp =
+ WorkspaceCreationHelper::CreateWorkspaceGroup(2, 4, 4,
+ projectionsGrpName);
+
+ ImggTomographicReconstruction alg;
+ TS_ASSERT_THROWS_NOTHING(alg.initialize());
+ TS_ASSERT(alg.isInitialized());
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("InputWorkspace", projectionsGrp));
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("MinProjectionAngle", 0.0));
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("MaxProjectionAngle", 260.0));
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("RelaxationParameter", 1.25));
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("CenterOfRotation", 10000));
+
+ // should throw because of the wrong center / outside of image dimensions
+ TS_ASSERT_THROWS(alg.execute(), std::runtime_error);
+ TS_ASSERT(!alg.isExecuted());
+ }
+
+ void test_exec_runs() {
+ const std::string projectionsGrpName("only_four_proj");
+ int ysize = 16;
+ int xsize = 16;
+ int numProj = 4;
+ API::WorkspaceGroup_sptr projectionsGrp =
+ WorkspaceCreationHelper::CreateWorkspaceGroup(numProj, ysize, xsize,
+ projectionsGrpName);
+
+ for (size_t proj = 0; proj < static_cast<size_t>(proj); ++proj) {
+ API::Workspace_sptr ws;
+ TS_ASSERT_THROWS_NOTHING(ws = projectionsGrp->getItem(proj));
+ API::MatrixWorkspace_sptr projWS;
+ TS_ASSERT_THROWS_NOTHING(
+ projWS = boost::dynamic_pointer_cast<API::MatrixWorkspace>(ws));
+ MantidVec &dataY = projWS->dataY(7);
+ TS_ASSERT_THROWS_NOTHING(std::fill(dataY.begin(), dataY.end(), 5000.0));
+ for (size_t idx = 0; idx < static_cast<size_t>(ysize); ++idx) {
+ projWS->dataY(ysize)[ysize - idx - 1] = 987.6;
+ }
+ }
+
+ ImggTomographicReconstruction alg;
+ // getProperty with workspaceGroup isn't working with
+ // Algorithm::processGroups. So we need retrieveWS<>() below:
+ const std::string reconName = "recon_1";
+ auto recon = runWithValidCenter(alg, projectionsGrp, 7, reconName);
+ // Fix when this works well
+ TSM_ASSERT("Expected that getProperty would return nullptr with "
+ "WorkspaceGroup when using processGroups()",
+ !recon);
+
+ TS_ASSERT_THROWS_NOTHING(
+ recon = API::AnalysisDataService::Instance()
+ .retrieveWS<API::WorkspaceGroup>(reconName));
+ TSM_ASSERT("The reconstruction/result workspace is not valid", recon);
+
+ size_t grpSize = recon->size();
+ TSM_ASSERT_EQUALS(
+ "The number of items in the output/reconstruction workspace is wrong",
+ grpSize, ysize);
+
+ const size_t pix1x = 5;
+ const size_t pix1y = 14;
+ const double referencePix1 = 8.0;
+ const size_t pix2x = 7;
+ const size_t pix2y = 8;
+ const double referencePix2 = 8.0;
+ for (size_t idx = 0; idx < grpSize; ++idx) {
+ auto wks = recon->getItem(idx);
+ TSM_ASSERT(
+ "The output workspace group should have valid slice workspaces", wks);
+ auto sliceWS = boost::dynamic_pointer_cast<API::MatrixWorkspace>(wks);
+ TSM_ASSERT("The slice workspaces should be of type MatrixWorkspace",
+ sliceWS);
+
+ TSM_ASSERT_EQUALS("Unexpected number of columns in output slices",
+ sliceWS->blocksize(), xsize);
+ TSM_ASSERT_EQUALS("Unexpected number of rows in output slices",
+ sliceWS->getNumberHistograms(), ysize);
+ TSM_ASSERT_DELTA("Unexpected value in output pixel",
+ sliceWS->readY(pix1y)[pix1x], referencePix1, 1e-4);
+ TSM_ASSERT_DELTA("Unexpected value in output pixel",
+ sliceWS->readY(pix2y)[pix2x], referencePix2, 1e-4);
+ }
+ }
+
+ void test_exec_runs0s() {
+ const std::string projectionsGrpName("a_couple_0_images");
+ int ysize = 8;
+ int xsize = 8;
+ int numProj = 2;
+ API::WorkspaceGroup_sptr projectionsGrp =
+ WorkspaceCreationHelper::CreateWorkspaceGroup(numProj, ysize, xsize,
+ projectionsGrpName);
+
+ for (size_t proj = 0; proj < static_cast<size_t>(proj); ++proj) {
+ API::Workspace_sptr ws;
+ TS_ASSERT_THROWS_NOTHING(ws = projectionsGrp->getItem(proj));
+ API::MatrixWorkspace_sptr projWS;
+ TS_ASSERT_THROWS_NOTHING(
+ projWS = boost::dynamic_pointer_cast<API::MatrixWorkspace>(ws));
+ for (size_t row = 0; row < static_cast<size_t>(ysize); ++row) {
+ MantidVec &dataY = projWS->dataY(row);
+ TS_ASSERT_THROWS_NOTHING(std::fill(dataY.begin(), dataY.end(), 0.0));
+ }
+ }
+
+ ImggTomographicReconstruction alg;
+ const std::string reconName = "recon_0";
+ auto recon = runWithValidCenter(alg, projectionsGrp, 7, reconName);
+ // Fix when this works well
+ TSM_ASSERT("Expected that getProperty would return nullptr with "
+ "WorkspaceGroup when using processGroups()",
+ !recon);
+
+ TS_ASSERT_THROWS_NOTHING(
+ recon = API::AnalysisDataService::Instance()
+ .retrieveWS<API::WorkspaceGroup>(reconName));
+ TSM_ASSERT("The reconstruction/result workspace is not valid", recon);
+
+ size_t grpSize = recon->size();
+ TSM_ASSERT_EQUALS(
+ "The number of items in the output/reconstruction workspace is wrong",
+ grpSize, ysize);
+
+ const size_t pix1x = 5;
+ const size_t pix1y = 2;
+ const double referencePix1 = 0.0;
+ const size_t pix2x = 7;
+ const size_t pix2y = 7;
+ const double referencePix2 = 0.0;
+ for (size_t idx = 0; idx < grpSize; ++idx) {
+ auto wks = recon->getItem(idx);
+ TSM_ASSERT(
+ "The output workspace group should have valid slice workspaces", wks);
+ auto sliceWS = boost::dynamic_pointer_cast<API::MatrixWorkspace>(wks);
+ TSM_ASSERT("The slice workspaces should be of type MatrixWorkspace",
+ sliceWS);
+
+ TSM_ASSERT_EQUALS("Unexpected number of columns in output slices",
+ sliceWS->blocksize(), xsize);
+ TSM_ASSERT_EQUALS("Unexpected number of rows in output slices",
+ sliceWS->getNumberHistograms(), ysize);
+ TSM_ASSERT_DELTA("Unexpected value in output pixel",
+ sliceWS->readY(pix1y)[pix1x], referencePix1, 1e-4);
+ TSM_ASSERT_DELTA("Unexpected value in output pixel",
+ sliceWS->readY(pix2y)[pix2x], referencePix2, 1e-4);
+ }
+ }
+
+private:
+ API::WorkspaceGroup_sptr runWithValidCenter(API::Algorithm &alg,
+ API::WorkspaceGroup_sptr wksg,
+ int center,
+ const std::string &outName) {
+ TS_ASSERT_THROWS_NOTHING(alg.initialize());
+ TS_ASSERT(alg.isInitialized());
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("InputWorkspace", wksg));
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("MinProjectionAngle", 0.0));
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("MaxProjectionAngle", 180.0));
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("RelaxationParameter", 1.25));
+ TS_ASSERT_THROWS_NOTHING(alg.setProperty("CenterOfRotation", center));
+
+ TS_ASSERT_THROWS_NOTHING(alg.setPropertyValue("OutputWorkspace", outName));
+ TSM_ASSERT_THROWS_NOTHING(
+ "execute() throwed for an algorithm with a "
+ "supposedly correct center parameter for which no "
+ "failure was expected",
+ alg.execute());
+ TSM_ASSERT("The algorithm execution didn't finish successfully when no "
+ "issues where expected ",
+ alg.isExecuted());
+
+ API::WorkspaceGroup_sptr result = alg.getProperty("OutputWorkspace");
+ return result;
+ }
+};
+
+#endif /* MANTID_ALGORITHMS_IMGGTOMOGRAPHICRECONSTRUCTIONTEST_H_ */
diff --git a/Framework/Algorithms/test/Tomography/FBPTomopyTest.h b/Framework/Algorithms/test/Tomography/FBPTomopyTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..091844e61345c8e001ad15d24440fb212fcd8982
--- /dev/null
+++ b/Framework/Algorithms/test/Tomography/FBPTomopyTest.h
@@ -0,0 +1,114 @@
+#ifndef MANTID_ALGORITHMS_FBPTOMOPYTEST_H_
+#define MANTID_ALGORITHMS_FBPTOMOPYTEST_H_
+
+#include <cxxtest/TestSuite.h>
+
+#include "MantidAlgorithms/Tomography/FBPTomopy.h"
+#include "MantidAPI/AlgorithmManager.h"
+#include "MantidKernel/Exception.h"
+
+#include "MantidTestHelpers/WorkspaceCreationHelper.h"
+
+#include <numeric>
+
+using namespace Mantid::Algorithms::Tomography;
+
+class FBPTomopyTest : 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 FBPTomopyTest *createSuite() { return new FBPTomopyTest(); }
+
+ static void destroySuite(FBPTomopyTest *suite) { delete suite; }
+
+ void test_null() {
+ FBPTomopy(nullptr, 0, 0, 0, nullptr, nullptr, nullptr, 0, 0);
+ }
+
+ void test_small_buffer_flat() {
+ const size_t numProj = 3;
+ const size_t ysize = 8;
+ const size_t xsize = 8;
+ const size_t projSize = numProj * ysize * xsize;
+ std::array<float, projSize> projImages;
+ std::fill(projImages.begin(), projImages.end(), 33.0f);
+
+ const size_t reconSize = ysize * ysize * xsize;
+ std::array<float, reconSize> reconVol;
+ std::fill(reconVol.begin(), reconVol.end(), 0.0f);
+
+ const size_t numAngles = numProj;
+ std::array<float, numAngles> angles{{0.0, 90.0, 180.0}};
+
+ const size_t numCenters = ysize;
+ std::array<float, numCenters> centers{{4, 3, 3, 4, 4, 4, 4, 3}};
+
+ FBPTomopy(projImages.data(), ysize, numProj, xsize, centers.data(),
+ angles.data(), reconVol.data(), xsize, ysize);
+
+ TS_ASSERT_DELTA(reconVol.front(), 66.7519, 1);
+ TS_ASSERT_DELTA(reconVol[1], 67.9431, 1);
+ TS_ASSERT_DELTA(reconVol[2], 90.9802, 1);
+ TS_ASSERT_DELTA(reconVol[3], 111.1035, 1);
+ TS_ASSERT_DELTA(reconVol[4], 88.6847, 1);
+ TS_ASSERT_DELTA(reconVol[20], 102.3709, 1);
+ TS_ASSERT_DELTA(reconVol[50], 97.4173, 1);
+ TS_ASSERT_DELTA(reconVol[100], 97.4172, 1);
+ TS_ASSERT_DELTA(reconVol[150], 99.8273, 1);
+ TS_ASSERT_DELTA(reconVol[projSize - 5], 82.4248, 1);
+ TS_ASSERT_DELTA(reconVol[projSize - 4], 74.1906, 1);
+ TS_ASSERT_DELTA(reconVol[projSize - 3], 60.5043, 1);
+ TS_ASSERT_DELTA(reconVol[projSize - 2], 67.9431, 1);
+ TS_ASSERT_DELTA(reconVol.back(), 66.7519, 1);
+ }
+
+ void test_buffer_idx() {
+ const size_t numProj = 8;
+ const size_t ysize = 16;
+ const size_t xsize = 16;
+ const size_t projSize = numProj * ysize * xsize;
+ std::array<float, projSize> projImages;
+ std::iota(projImages.begin(), projImages.end(), 0.0f);
+
+ // inconsistent/stressing values
+ std::fill(projImages.begin() + 300, projImages.begin() + 400, 333.0f);
+ std::fill(projImages.begin() + 600, projImages.begin() + 850, 999.0f);
+ std::fill(projImages.begin() + 990, projImages.begin() + 1100, 1000.0f);
+ std::fill(projImages.begin() + 1500, projImages.begin() + 1700, -444.0f);
+ std::fill(projImages.begin() + 1900, projImages.begin() + 2000, 765.0f);
+
+ const size_t reconSize = ysize * ysize * xsize;
+ std::array<float, reconSize> reconVol;
+ std::fill(reconVol.begin(), reconVol.end(), 0.0f);
+
+ const size_t numAngles = numProj;
+ std::array<float, numAngles> angles{
+ {0.0, 45.0, 90.0, 135.0, 180.0, 225.0, 270.0, 315.0}};
+
+ const size_t numCenters = ysize;
+ std::array<float, numCenters> centers;
+ std::fill(centers.begin(), centers.end(), 7.5f);
+
+ FBPTomopy(projImages.data(), ysize, numProj, xsize, centers.data(),
+ angles.data(), reconVol.data(), xsize, ysize);
+
+ TS_ASSERT_DELTA(reconVol.front(), 241.6610, 1);
+ TS_ASSERT_DELTA(reconVol[1], 286.6727, 1);
+ TS_ASSERT_DELTA(reconVol[2], 392.6853, 1);
+ TS_ASSERT_DELTA(reconVol[3], 350.0282, 1);
+ TS_ASSERT_DELTA(reconVol[4], 429.8395, 1);
+ TS_ASSERT_DELTA(reconVol[200], 516.5272, 1);
+ TS_ASSERT_DELTA(reconVol[500], 1202.6435, 1);
+ TS_ASSERT_DELTA(reconVol[1000], 3604.5090, 1);
+ TS_ASSERT_DELTA(reconVol[1500], 8092.9765, 1);
+ TS_ASSERT_DELTA(reconVol[1900], 7941.9135, 1);
+ TS_ASSERT_DELTA(reconVol[projSize - 5], 6829.1318, 1);
+ TS_ASSERT_DELTA(reconVol[projSize - 4], 5662.6342, 1);
+ TS_ASSERT_DELTA(reconVol[projSize - 3], 5747.1845, 1);
+ TS_ASSERT_DELTA(reconVol[projSize - 2], 6034.2163, 1);
+ TS_ASSERT_DELTA(reconVol.back(), 3736.5483, 1);
+ }
+};
+
+#endif /* MANTID_ALGORITHM_FBPTOMOPYTEST_H_ */
diff --git a/docs/source/algorithms/ImggTomographicReconstruction-v1.rst b/docs/source/algorithms/ImggTomographicReconstruction-v1.rst
new file mode 100644
index 0000000000000000000000000000000000000000..709b0f076e925ef5a3d3f2e776e70e12860a4d8b
--- /dev/null
+++ b/docs/source/algorithms/ImggTomographicReconstruction-v1.rst
@@ -0,0 +1,125 @@
+
+.. algorithm::
+
+.. summary::
+
+.. alias::
+
+.. properties::
+
+Description
+-----------
+
+.. warning:: This is an early, experimental version of the algorithm.
+
+The input and output workspaces are workspace groups where every
+element is an image workspace. In the input workspace every image is a
+2D projection from a different angle whereas in the output workspace
+every image is a slice of a reconstructed 3D volume. The input
+workspace must have one image workspace per projection from a
+tomography imaging experiment. The output workspace will have one
+image workspace for every slice of the output reconstructed volume.
+
+The following method is supported: FBP (following the TomoPy
+implementation [TomoPy2014]).
+
+The implementation of TomoPy methods are based on the TomoPy source
+code available from https://github.com/tomopy/tomopy/, which is:
+
+::
+
+ Copyright 2015. UChicago Argonne, LLC. This software was produced
+ under U.S. Government contract DE-AC02-06CH11357 for Argonne National
+ Laboratory (ANL), which is operated by UChicago Argonne, LLC for the
+ U.S. Department of Energy. The U.S. Government has rights to use,
+ reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR
+ UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR
+ ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is
+ modified to produce derivative works, such modified software should
+ be clearly marked, so as not to confuse it with the version available
+ from ANL.
+
+ Additionally, redistribution and use in source and binary forms, with
+ or without modification, are permitted provided that the following
+ conditions are met:
+
+ * Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+ * Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in
+ the documentation and/or other materials provided with the
+ distribution.
+
+ * Neither the name of UChicago Argonne, LLC, Argonne National
+ Laboratory, ANL, the U.S. Government, nor the names of its
+ contributors may be used to endorse or promote products derived
+ from this software without specific prior written permission.
+
+ THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS
+ "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago
+ Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ POSSIBILITY OF SUCH DAMAGE.
+
+References:
+
+.. [TomoPy2014] Gursoy D, De Carlo F, Xiao X,
+ Jacobsen C. (2014). TomoPy: a framework for the analysis of
+ synchrotron tomographic data. J. Synchrotron Rad. 21. 1188-1193
+ doi:10.1107/S1600577514013939
+
+
+Usage
+-----
+
+**Example - ReconstructProjections**
+
+.. testcode:: ReconstructProjections
+
+ # Note: you would load FITS images like this:
+ # wsg = LoadFITS(Filename='FITS_small_01.fits', LoadAsRectImg=1, OutputWorkspace='projections')
+ wsg_name = 'projections'
+ # Produce 16 projections with 32x32 pixels
+ projections = []
+ unit_label=UnitFactory.create('Label')
+ unit_label.setLabel('width','cm')
+ for proj in range(0, 8):
+ wks_name = 'wks_proj_' + str(proj)
+ wks = CreateSampleWorkspace(NumBanks=32, BankPixelWidth=1, XMin=0, XMax=32, BinWidth=1, OutputWorkspace=wks_name)
+ wks.getAxis(0).setUnit('Label')
+ projections.append(wks)
+ wsg_proj = GroupWorkspaces(projections, OutputWorkspace=wsg_name)
+ wsg_reconstructed = ImggTomographicReconstruction(InputWorkspace=wsg_proj, CenterOfRotation=15)
+ rows = wsg_reconstructed.getItem(0).getNumberHistograms()
+ columns = wsg_reconstructed.getItem(0).blocksize()
+ print ("The output reconstruction has {0} slices of {1} by {2} pixels".
+ format(wsg_reconstructed.size(), rows, columns))
+ slice_idx = 2
+ coord_x = 8
+ coord_y = 15
+ print ("Value of pixel at coordinate ({0},{1}) in slice {2}: {3:.1f}".
+ format(coord_x, coord_y, slice_idx,
+ wsg_reconstructed.getItem(2).readY(coord_y)[coord_x]))
+
+.. testcleanup:: ReconstructProjections
+
+ DeleteWorkspace(wsg_name)
+
+Output:
+
+.. testoutput:: ReconstructProjections
+
+ The output reconstruction has 32 slices of 32 by 32 pixels
+ Value of pixel at coordinate (8,15) in slice 2: 2.4
+
+.. categories::
+
+.. sourcelink::
diff --git a/docs/source/release/v3.7.0/diffraction.rst b/docs/source/release/v3.7.0/diffraction.rst
index be06a2a910fadf08cb1c2da07a8e5654f30edcee..c6187f92f77795633cdcdcbee37a6ae83b56c0eb 100644
--- a/docs/source/release/v3.7.0/diffraction.rst
+++ b/docs/source/release/v3.7.0/diffraction.rst
@@ -97,9 +97,16 @@ Imaging
- The new algorithm :ref:`ImggAggregateWavelengths <algm-ImggAggregateWavelengths>`
aggregates stacks of images from wavelength dependent data.
+- The algorithm `ImggTomographicReconstruction
+ <algm-ImggTomographicReconstruction>` has been introduced. This is a
+ first experimental version that implements the Filtered
+ Back-Projection (FBP) reconstruction method using the FBP
+ implementation of the `TomoPy package
+ <http://www.aps.anl.gov/tomopy/>`_.
- Images loaded as Mantid workspaces can now be saved into FITS files
using the algorithm :ref:`SaveFITS <algm-SaveFITS>`.
+
Improvements in the tomographic reconstruction graphical user interface:
- New capabilities added when visualizing stacks of images:
diff --git a/docs/source/release/v3.7.0/framework.rst b/docs/source/release/v3.7.0/framework.rst
index 7a7f122b47bef5f6f64d4e76328cc227a4a44ed2..846bbf2bb2e7787b3333074035f923cf02981310 100644
--- a/docs/source/release/v3.7.0/framework.rst
+++ b/docs/source/release/v3.7.0/framework.rst
@@ -28,9 +28,7 @@ New
:ref:`GetDetOffsetsMultiPeaks <algm-GetDetOffsetsMultiPeaks>`, :ref:`CalibrateRectangularDetectors <algm-CalibrateRectangularDetectors>`, *et al* and
minimizes the difference between the *DIFC* of the instrument and
calibration by moving and rotating instrument components.
-
- :ref:`CorrectTOF <algm-CorrectTOF>` applies to the time-of-flight correction which considers the specified elastic peak position.
-
- :ref:`EnggFitDIFCFromPeaks <algm-AlignComponents>` fits GSAS calibration
parameters (DIFA, DIFC, TZERO) from peaks fitted using
:ref:`EnggFitPeaks <algm-EnggFitPeaks>`.
@@ -42,6 +40,9 @@ New
- :ref:`ImggAggregateWavelengths <algm-ImggAggregateWavelengths>`
aggregates stacks of images from wavelength dependent imaging
into one or more output bands.
+- :ref:`ImggTomographicReconstruction
+ <algm-ImggTomographicReconstruction>` implements a method for 3D
+ tomographic reconstruction from projection images.
- :ref:`SaveFITS <algm-SaveFITS>` saves images in FITS format.
Renamed