use inverse projection matrix for ray generation, add benchmarks for ray generation

.gitignore

 .idea
 .vscode
+cmake-build-*
 cmake-build-debug
 cmake-build-release
 build/

CMakeLists.txt

@@ -119,35 +119,42 @@ if(NOT ELSA_MASTER_PROJECT)
     set(ELSA_TESTING OFF)
 endif(NOT ELSA_MASTER_PROJECT)
 
-if(ELSA_TESTING)
-    message(STATUS "elsa testing is enabled")
+if(ELSA_TESTING OR ELSA_BENCHMARKS)
 
     enable_testing()
     add_subdirectory(thirdparty/Catch2)
 
-    add_custom_target(tests
-        COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure
-        WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
-        COMMENT "Build and run all the tests.")
-
     # add the CMake modules for automatic test discovery
     set(CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/thirdparty/Catch2/contrib" ${CMAKE_MODULE_PATH})
 
-    if(ELSA_COVERAGE)
-        message(STATUS "elsa test coverage is enabled")
+    if(ELSA_TESTING)
+        message(STATUS "elsa testing is enabled")
+
+        add_custom_target(tests
+            COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure
+            WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+            COMMENT "Build and run all the tests.")
+        if(ELSA_COVERAGE)
+            message(STATUS "elsa test coverage is enabled")
+
+            include(CodeCoverage)
+            APPEND_COVERAGE_COMPILER_FLAGS()
+            set(COVERAGE_LCOV_EXCLUDES '${PROJECT_SOURCE_DIR}/thirdparty/*' '/usr/*')
+            SETUP_TARGET_FOR_COVERAGE_LCOV(NAME test_coverage EXECUTABLE ctest)
 
-        include(CodeCoverage)
-        APPEND_COVERAGE_COMPILER_FLAGS()
-        set(COVERAGE_LCOV_EXCLUDES '${PROJECT_SOURCE_DIR}/thirdparty/*' '/usr/*')
-        SETUP_TARGET_FOR_COVERAGE_LCOV(NAME test_coverage EXECUTABLE ctest)
+        else(ELSA_COVERAGE)
+            message(STATUS "elsa test coverage is disabled")
+        endif(ELSA_COVERAGE)
+    endif(ELSA_TESTING)
 
-    else(ELSA_COVERAGE)
-        message(STATUS "elsa test coverage is disabled")
-    endif(ELSA_COVERAGE)
+    if (ELSA_BENCHMARKS)
+        message(STATUS "elsa benchmarks are enabled")
+        add_subdirectory(benchmarks)
+    endif(ELSA_BENCHMARKS)
 
-else(ELSA_TESTING)
+else(ELSA_TESTING OR ELSA_BENCHMARKS)
     message(STATUS "    elsa testing is disabled")
-endif(ELSA_TESTING)
+endif(ELSA_TESTING OR ELSA_BENCHMARKS)
 
 
 # ------------ add code/docs -----------

benchmarks/CMakeLists.txt

+cmake_minimum_required(VERSION 3.9)
+
+# enable ctest and Catch test discovery
+include(CTest)
+include(Catch)
+
+macro(ELSA_BENCHMARK NAME)
+    # Add source to all sources, so we can create a common target to build and run the benchmarks
+    SET(BENCHMARK_SOURCES ${BENCHMARK_SOURCES} bench_${NAME}.cpp)
+
+    # create the test executable
+    add_executable(bench_${NAME} EXCLUDE_FROM_ALL bench_${NAME}.cpp bench_main.cpp)
+    # add catch and the corresponding elsa library
+    target_link_libraries(bench_${NAME} PRIVATE Catch2::Catch2 elsa::all)
+    # enable C++17
+    target_compile_features(bench_${NAME} PUBLIC cxx_std_17)
+    # include helpers
+    target_include_directories(bench_${NAME} PRIVATE ${CMAKE_SOURCE_DIR}/elsa/test_routines/)
+    # discover test for CMake and CTest
+    catch_discover_tests(bench_${NAME} TEST_SPEC "" )
+endmacro(ELSA_BENCHMARK)
+
+ELSA_BENCHMARK(RayGenerationBench)
+ELSA_BENCHMARK(Projectors)
+
+# Add a single executable for all benchmarks, as CTest removes a lot of the output
+add_executable(bench_all EXCLUDE_FROM_ALL bench_main.cpp ${BENCHMARK_SOURCES})
+# add catch and the corresponding elsa library
+target_link_libraries(bench_all PRIVATE Catch2::Catch2 elsa::all)
+# enable C++17
+target_compile_features(bench_all PUBLIC cxx_std_17)
+# include helpers
+target_include_directories(bench_all PRIVATE ${CMAKE_SOURCE_DIR}/elsa/test_routines/)
+
+# Add the custom target to run all the benchmarks
+add_custom_target(benchmarks
+        COMMAND ${CMAKE_CURRENT_BINARY_DIR}/bench_all
+        WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+        USES_TERMINAL
+        COMMENT "Run benchmarks")
+add_dependencies(benchmarks bench_all)
\ No newline at end of file

benchmarks/bench_Projectors.cpp

+/**
+ * \file test_RayGenerationBench.cpp
+ *
+ * \brief Benchmarks for projectors
+ *
+ * \author David Frank
+ */
+
+#define CATCH_CONFIG_ENABLE_BENCHMARKING
+
+#include <catch2/catch.hpp>
+
+#include "Logger.h"
+#include "elsaDefines.h"
+#include "PhantomGenerator.h"
+#include "CircleTrajectoryGenerator.h"
+
+#include "SiddonsMethod.h"
+#include "JosephsMethod.h"
+
+using namespace elsa;
+
+using DataContainer_t = DataContainer<real_t>;
+
+// Some minor fixtures
+template <typename Projector>
+void runProjector2D(index_t coeffsPerDim)
+{
+    // generate 2d phantom
+    IndexVector_t size(2);
+    size << coeffsPerDim, coeffsPerDim;
+    auto phantom = DataContainer_t(DataDescriptor(size));
+    phantom = 0;
+
+    // generate circular trajectory
+    index_t noAngles{180}, arc{360};
+    auto [geometry, sinoDescriptor] = CircleTrajectoryGenerator::createTrajectory(
+        noAngles, phantom.getDataDescriptor(), arc, 20, 20);
+
+    // setup operator for 2d X-ray transform
+    Projector projector(phantom.getDataDescriptor(), *sinoDescriptor, geometry);
+
+    DataContainer_t sinogram(*sinoDescriptor);
+    BENCHMARK("Forward projection")
+    {
+        sinogram = projector.apply(phantom);
+        return sinogram;
+    };
+
+    BENCHMARK("Backward projection") { return projector.applyAdjoint(sinogram); };
+}
+
+template <typename Projector>
+void runProjector3D(index_t coeffsPerDim)
+{
+    // Turn logger off
+    Logger::setLevel(Logger::LogLevel::OFF);
+
+    // generate 2d phantom
+    IndexVector_t size(3);
+    size << coeffsPerDim, coeffsPerDim, coeffsPerDim;
+    auto phantom = DataContainer_t(DataDescriptor(size));
+    phantom = 0;
+
+    // generate circular trajectory
+    index_t noAngles{180}, arc{360};
+    auto [geometry, sinoDescriptor] = CircleTrajectoryGenerator::createTrajectory(
+        noAngles, phantom.getDataDescriptor(), arc, 20, 20);
+
+    // setup operator for 2d X-ray transform
+    Projector projector(phantom.getDataDescriptor(), *sinoDescriptor, geometry);
+
+    DataContainer_t sinogram(*sinoDescriptor);
+    BENCHMARK("Forward projection")
+    {
+        sinogram = projector.apply(phantom);
+        return sinogram;
+    };
+
+    BENCHMARK("Backward projection") { return projector.applyAdjoint(sinogram); };
+}
+
+TEST_CASE("Testing Siddon's projector in 2D")
+{
+    // Turn logger off
+    Logger::setLevel(Logger::LogLevel::OFF);
+
+    using Siddon = SiddonsMethod<real_t>;
+
+    GIVEN("A 8x8 Problem:") { runProjector2D<Siddon>(8); }
+
+    GIVEN("A 16x16 Problem:") { runProjector2D<Siddon>(16); }
+
+    GIVEN("A 32x32 Problem:") { runProjector2D<Siddon>(32); }
+
+    GIVEN("A 64x64 Problem:") { runProjector2D<Siddon>(64); }
+}
+
+TEST_CASE("Testing Siddon's projector in 3D")
+{
+    // Turn logger off
+    Logger::setLevel(Logger::LogLevel::OFF);
+
+    using Siddon = SiddonsMethod<real_t>;
+
+    GIVEN("A 8x8x8 Problem:") { runProjector3D<Siddon>(8); }
+
+    GIVEN("A 16x16x16 Problem:") { runProjector3D<Siddon>(16); }
+
+    GIVEN("A 32x32x32 Problem:") { runProjector3D<Siddon>(32); }
+}
+
+TEST_CASE("Testing Joseph's projector in 2D")
+{
+    // Turn logger off
+    Logger::setLevel(Logger::LogLevel::OFF);
+
+    using Joseph = JosephsMethod<real_t>;
+
+    GIVEN("A 8x8 Problem:") { runProjector2D<Joseph>(8); }
+
+    GIVEN("A 16x16 Problem:") { runProjector2D<Joseph>(16); }
+
+    GIVEN("A 32x32 Problem:") { runProjector2D<Joseph>(32); }
+
+    GIVEN("A 64x64 Problem:") { runProjector2D<Joseph>(64); }
+}
+
+TEST_CASE("Testing Joseph's projector in 3D")
+{
+    // Turn logger off
+    Logger::setLevel(Logger::LogLevel::OFF);
+
+    using Joseph = JosephsMethod<real_t>;
+
+    GIVEN("A 8x8x8 Problem:") { runProjector3D<Joseph>(8); }
+
+    GIVEN("A 16x16x16 Problem:") { runProjector3D<Joseph>(16); }
+
+    GIVEN("A 32x32x32 Problem:") { runProjector3D<Joseph>(32); }
+}
\ No newline at end of file

benchmarks/bench_RayGenerationBench.cpp

+/**
+ * \file test_RayGenerationBench.cpp
+ *
+ * \brief Benchmarks for ray generation
+ *
+ * \author David Frank
+ */
+
+#define CATCH_CONFIG_ENABLE_BENCHMARKING
+
+#include <catch2/catch.hpp>
+#include "Geometry.h"
+#include <string>
+#include <cstdlib>
+
+using namespace elsa;
+static const index_t dimension = 2;
+
+void iterate2D(const Geometry& geo)
+{
+    for (real_t detPixel : std::initializer_list<real_t>{0.5, 2.5, 4.5}) {
+        RealVector_t pixel(1);
+        pixel << detPixel;
+        BENCHMARK("Ray for detector at pixel " + std::to_string(detPixel))
+        {
+            return geo.computeRayTo(pixel);
+        };
+    }
+}
+
+void iterate3D(const Geometry& geo)
+{
+    for (real_t detPixel1 : std::initializer_list<real_t>{0.5, 2.5, 4.5}) {
+        for (real_t detPixel2 : std::initializer_list<real_t>{0.5, 2.5, 4.5}) {
+            RealVector_t pixel(2);
+            pixel << detPixel1, detPixel2;
+            BENCHMARK("Ray for detector at pixel " + std::to_string(detPixel1) + "/"
+                      + std::to_string(detPixel2))
+            {
+                return geo.computeRayTo(pixel);
+            };
+        }
+    }
+}
+
+TEST_CASE("Ray generation for 2D")
+{
+    IndexVector_t volCoeff(2);
+    volCoeff << 5, 5;
+    DataDescriptor ddVol(volCoeff);
+
+    IndexVector_t detCoeff(1);
+    detCoeff << 5;
+    DataDescriptor ddDet(detCoeff);
+
+    real_t s2c = 10;
+    real_t c2d = 4;
+
+    GIVEN("Geometry without offset and rotation")
+    {
+        Geometry g(s2c, c2d, 0, ddVol, ddDet);
+
+        // test outer + central pixels
+        iterate2D(g);
+    }
+
+    GIVEN("Geometry with offset but no rotation")
+    {
+        real_t offset = 2;
+        Geometry g(s2c, c2d, 0, ddVol, ddDet, offset);
+
+        // test outer + central pixels
+        iterate2D(g);
+    }
+
+    GIVEN("Geometry at 90°, but no offset")
+    {
+        real_t angle = pi_t / 2; // 90 degrees
+        Geometry g(s2c, c2d, angle, ddVol, ddDet);
+
+        // test outer + central pixels
+        iterate2D(g);
+    }
+
+    GIVEN("Geometry at 45° with offset")
+    {
+        real_t angle = pi_t / 4; // 45 degrees
+        real_t cx = -1;
+        real_t cy = 2;
+        Geometry g(s2c, c2d, angle, ddVol, ddDet, 0, cx, cy);
+
+        // test outer + central pixels
+        iterate2D(g);
+    }
+}
+
+TEST_CASE("Ray generation for 3D")
+{
+    IndexVector_t volCoeff(3);
+    volCoeff << 5, 5, 5;
+    DataDescriptor ddVol(volCoeff);
+
+    IndexVector_t detCoeff(2);
+    detCoeff << 5, 5;
+    DataDescriptor ddDet(detCoeff);
+
+    real_t s2c = 10;
+    real_t c2d = 4;
+
+    GIVEN("Geometry without offset and rotation")
+    {
+        Geometry g(s2c, c2d, ddVol, ddDet, 0);
+
+        // test outer + central pixels
+        iterate3D(g);
+    }
+
+    GIVEN("Geometry with offset but no rotation")
+    {
+        real_t px = -1;
+        real_t py = 3;
+        Geometry g(s2c, c2d, ddVol, ddDet, 0, 0, 0, px, py);
+
+        // test outer + central pixels
+        iterate3D(g);
+    }
+
+    GIVEN("Geometry at 90°, but no offset")
+    {
+        real_t angle = pi_t / 2;
+        Geometry g(s2c, c2d, ddVol, ddDet, angle);
+
+        // test outer + central pixels
+        iterate3D(g);
+    }
+
+    GIVEN("Geometry at 45°/22.5 with offset")
+    {
+        real_t angle1 = pi_t / 4;
+        real_t angle2 = pi_t / 2;
+        RealVector_t offset(3);
+        offset << 1, -2, -1;
+        Geometry g(s2c, c2d, ddVol, ddDet, angle1, angle2, 0, 0, 0, offset[0], offset[1],
+                   offset[2]);
+
+        // test outer + central pixels
+        iterate3D(g);
+    }
+
+    GIVEN("Geometry at 45/22.5/12.25 with offset")
+    {
+        real_t angle1 = pi_t / 4;
+        real_t angle2 = pi_t / 2;
+        real_t angle3 = pi_t / 8;
+
+        RealMatrix_t rot1(3, 3);
+        rot1 << std::cos(angle1), 0, std::sin(angle1), 0, 1, 0, -std::sin(angle1), 0,
+            std::cos(angle1);
+
+        RealMatrix_t rot2(3, 3);
+        rot2 << std::cos(angle2), -std::sin(angle2), 0, std::sin(angle2), std::cos(angle2), 0, 0, 0,
+            1;
+
+        RealMatrix_t rot3(3, 3);
+        rot3 << std::cos(angle3), 0, std::sin(angle3), 0, 1, 0, -std::sin(angle3), 0,
+            std::cos(angle3);
+
+        RealMatrix_t R = rot1 * rot2 * rot3;
+
+        Geometry g(s2c, c2d, ddVol, ddDet, R);
+
+        // test outer + central pixels
+        iterate3D(g);
+    }
+}
\ No newline at end of file

benchmarks/bench_main.cpp

+#define CATCH_CONFIG_RUNNER
+#define CATCH_CONFIG_ENABLE_BENCHMARKING
+#include <catch2/catch.hpp>
+
+int main(int argc, char* argv[])
+{
+    Catch::Session session; // There must be exactly one instance
+
+    // writing to session.configData() here sets defaults
+    // this is the preferred way to set them
+
+    // Get binary name
+    std::string executable_path = argv[0];
+
+    // Find name or executable (without filesystem as GCC 7.4 doesn't support it)
+    // find last of / or \ (depending on linux or windows systems)
+    auto found = executable_path.find_last_of("/\\");
+    std::string filename = executable_path.substr(found + 1);
+
+    // set reporter and filename to match binary
+    session.configData().reporterName = "console";
+    // session.configData().outputFilename = filename + ".xml";
+
+    int returnCode = session.applyCommandLine(argc, argv);
+    if (returnCode != 0) // Indicates a command line error
+        return returnCode;
+
+    // writing to session.configData() or session.Config() here
+    // overrides command line args
+    // only do this if you know you need to
+
+    int numFailed = session.run();
+
+    // numFailed is clamped to 255 as some unices only use the lower 8 bits.
+    // This clamping has already been applied, so just return it here
+    // You can also do any post run clean-up here
+    return numFailed;
+}

elsa/CMakeLists.txt

@@ -13,7 +13,7 @@ macro(ELSA_TEST NAME)
     # create the test executable
     add_executable(test_${NAME} EXCLUDE_FROM_ALL test_${NAME}.cpp test_main.cpp)
     # add catch and the corresponding elsa library
-    target_link_libraries(test_${NAME} PRIVATE Catch2::Catch2 ${ELSA_MODULE_TARGET_NAME})
+    target_link_libraries(test_${NAME} PRIVATE Catch2::Catch2 ${ELSA_MODULE_TARGET_NAME} elsa::test_routines)
     # enable C++17
     target_compile_features(test_${NAME} PUBLIC cxx_std_17)
     # include helpers
@@ -50,6 +50,7 @@ if(ELSA_BUILD_CUDA_PROJECTORS)
     add_subdirectory(projectors_cuda)
 endif(ELSA_BUILD_CUDA_PROJECTORS)
 add_subdirectory(generators)
+add_subdirectory(test_routines)
 
 
 # library to build and add all registered components of the library

elsa/projectors/Geometry.cpp

@@ -148,14 +148,9 @@ namespace elsa
         RealVector_t homP(_objectDimension);
         homP << p, 1;
 
-        // solve for ray direction
-        RealVector_t rd = (_P.block(0, 0, _objectDimension, _objectDimension))
-                              .colPivHouseholderQr()
-                              .solve(homP)
-                              .normalized();
-        rd.normalize();
-
-        return std::make_pair(ro, rd);
+        // multiplication of inverse projection matrix and homogeneous detector coordinate
+        auto rd = (_Pinv * homP).normalized();
+        return std::make_pair(ro, rd.head(_objectDimension));
     }
 
     const RealMatrix_t& Geometry::getProjectionMatrix() const { return _P; }

elsa/projectors/tests/test_Geometry.cpp

@@ -129,8 +129,9 @@ SCENARIO("Testing 2D geometries")
                     auto c = g.getCameraCenter();
                     REQUIRE((ro - c).sum() == Approx(0));
 
-                    real_t factor =
-                        (std::abs(rd[0]) > 0) ? ((pixel[0] - ro[0] - px) / rd[0]) : (s2c + c2d);
+                    real_t factor = (std::abs(rd[0]) >= Approx(0.001))
+                                        ? ((pixel[0] - ro[0] - px) / rd[0])
+                                        : (s2c + c2d);
                     real_t detCoordY = ro[1] + factor * rd[1];
                     REQUIRE(detCoordY == Approx(ddVol.getLocationOfOrigin()[1] + c2d));
                 }

elsa/test_routines/CMakeLists.txt

 cmake_minimum_required(VERSION 3.9)
 
+# set the name of the module
+set(ELSA_MODULE_NAME test_routines)
+set(ELSA_MODULE_TARGET_NAME elsa_test_routines)
+set(ELSA_MODULE_EXPORT_TARGET elsa_${ELSA_MODULE_NAME}Targets)
+
+
+# list all the headers of the module
+set(MODULE_HEADERS
+        testHelpers.h)
+
+# list all the code files of the module
+set(MODULE_SOURCES
+        )
+
+
+# build the module library
+add_library(${ELSA_MODULE_TARGET_NAME} ${MODULE_HEADERS} ${MODULE_SOURCES})
+add_library(elsa::${ELSA_MODULE_NAME} ALIAS ${ELSA_MODULE_TARGET_NAME})
+
+target_link_libraries(${ELSA_MODULE_TARGET_NAME} PUBLIC elsa_core elsa_logging)
+
+target_include_directories(${ELSA_MODULE_TARGET_NAME}
+        PUBLIC
+        $<INSTALL_INTERFACE:include/elsa/${ELSA_MODULE_NAME}>
+        $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+        )
+
+# require C++17
+target_compile_features(${ELSA_MODULE_TARGET_NAME} PUBLIC cxx_std_17)
+# set -fPIC
+set_target_properties(${ELSA_MODULE_TARGET_NAME} PROPERTIES POSITION_INDEPENDENT_CODE ON)
+
+# This is only used privately, so don't bother exporting it and installing it
\ No newline at end of file

elsa/test_routines/testHelpers.h

@@ -5,67 +5,107 @@
 #include <random>
 #include "elsaDefines.h"
 
-/**
- * \brief comparing two number types for approximate equality for complex and regular number
- *
- * \tparam T - arithmetic data type
- * \return true if same number
- *
- * Use example in test case: REQUIRE(checkSameNumbers(a, b));
- * The CHECK(...) assertion in the function ensures that the values are reported when the test fails
- */
-template <typename T>
-bool checkSameNumbers(T left, T right, int epsilonFactor = 1)
+#include <iomanip>
+#include <limits>
+#include <cassert>
+
+namespace elsa
 {
-    using numericalBaseType = elsa::GetFloatingPointType_t<T>;
+    /**
+     * \brief Epsilon value for our test suit
+     */
+    static constexpr elsa::real_t epsilon = static_cast<elsa::real_t>(0.0001);
+
+    /**
+     * \brief comparing two number types for approximate equality for complex and regular number
+     *
+     * \tparam T - arithmetic data type
+     * \return true if same number
+     *
+     * Use example in test case: REQUIRE(checkSameNumbers(a, b));
+     * The CHECK(...) assertion in the function ensures that the values are reported when the test
+     * fails
+     */
+    template <typename T>
+    bool checkSameNumbers(T left, T right, int epsilonFactor = 1)
+    {
+        using numericalBaseType = elsa::GetFloatingPointType_t<T>;
 
-    numericalBaseType eps = std::numeric_limits<numericalBaseType>::epsilon()
-                            * static_cast<numericalBaseType>(epsilonFactor)
-                            * static_cast<numericalBaseType>(100);
+        numericalBaseType eps = std::numeric_limits<numericalBaseType>::epsilon()
+                                * static_cast<numericalBaseType>(epsilonFactor)
+                                * static_cast<numericalBaseType>(100);
 
-    if constexpr (std::is_same_v<T,
-                                 std::complex<float>> || std::is_same_v<T, std::complex<double>>) {
-        CHECK(Approx(left.real()).epsilon(eps) == right.real());
-        CHECK(Approx(left.imag()).epsilon(eps) == right.imag());
-        return Approx(left.real()).epsilon(eps) == right.real()
-               && Approx(left.imag()).epsilon(eps) == right.imag();
-    } else {
-        CHECK(Approx(left).epsilon(eps) == right);
-        return Approx(left).epsilon(eps) == right;
+        if constexpr (std::is_same_v<
+                          T, std::complex<float>> || std::is_same_v<T, std::complex<double>>) {
+            CHECK(Approx(left.real()).epsilon(eps) == right.real());
+            CHECK(Approx(left.imag()).epsilon(eps) == right.imag());
+            return Approx(left.real()).epsilon(eps) == right.real()
+                   && Approx(left.imag()).epsilon(eps) == right.imag();
+        } else {
+            CHECK(Approx(left).epsilon(eps) == right);
+            return Approx(left).epsilon(eps) == right;
+        }
     }
-}
 
-/**
- * \brief Generates a random Eigen matrix for different data_t types with integer values limited to
- * a certain range