LCOV - code coverage report
Current view: top level - projectors_cuda/tests - test_JosephsMethodCUDA.cpp (source / functions) Hit Total Coverage
Test: test_coverage.info.cleaned Lines: 1099 1126 97.6 %
Date: 2021-12-08 03:24:48 Functions: 25 37 67.6 %

          Line data    Source code
       1             : #include "doctest/doctest.h"
       2             : 
       3             : #include "JosephsMethodCUDA.h"
       4             : #include "Geometry.h"
       5             : #include "VolumeDescriptor.h"
       6             : #include "Logger.h"
       7             : #include "PlanarDetectorDescriptor.h"
       8             : #include "testHelpers.h"
       9             : 
      10             : #include <array>
      11             : 
      12             : using namespace elsa;
      13             : using namespace elsa::geometry;
      14             : using namespace doctest;
      15             : 
      16             : // TODO(dfrank): remove this and replace with checkApproxEq
      17             : using doctest::Approx;
      18             : 
      19             : /*
      20             :  * this function declaration can be used in conjunction with decltype to deduce the
      21             :  * template parameter of a templated class at compile time
      22             :  *
      23             :  * the template parameter must be a typename
      24             :  */
      25             : template <template <typename> typename T, typename data_t>
      26             : constexpr data_t return_data_t(const T<data_t>&);
      27             : 
      28          34 : TEST_CASE_TEMPLATE("JosephsMethodCUDA: Calls to functions of super class", TestType,
      29             :                    JosephsMethodCUDA<float>, JosephsMethodCUDA<double>)
      30             : {
      31             :     // Turn logger of
      32           4 :     Logger::setLevel(Logger::LogLevel::OFF);
      33             : 
      34             :     using data_t = decltype(return_data_t(std::declval<TestType>()));
      35           8 :     GIVEN("A projector")
      36             :     {
      37           8 :         IndexVector_t volumeDims(2), sinoDims(2);
      38           4 :         const index_t volSize = 10;
      39           4 :         const index_t detectorSize = 10;
      40           4 :         const index_t numImgs = 10;
      41           4 :         volumeDims << volSize, volSize;
      42           4 :         sinoDims << detectorSize, numImgs;
      43           8 :         VolumeDescriptor volumeDescriptor(volumeDims);
      44           8 :         DataContainer<data_t> volume(volumeDescriptor);
      45           4 :         volume = 0;
      46             : 
      47           4 :         auto stc = SourceToCenterOfRotation{20 * volSize};
      48           4 :         auto ctr = CenterOfRotationToDetector{volSize};
      49             : 
      50           8 :         std::vector<Geometry> geom;
      51          44 :         for (index_t i = 0; i < numImgs; i++) {
      52          40 :             real_t angle = static_cast<real_t>(i) * 2 * pi_t / 10;
      53          40 :             geom.emplace_back(stc, ctr, Radian{angle}, VolumeData2D{Size2D{volumeDims}},
      54             :                               SinogramData2D{Size2D{sinoDims}});
      55             :         }
      56             : 
      57           8 :         PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
      58           8 :         DataContainer<data_t> sino(sinoDescriptor);
      59           4 :         sino = 0;
      60             : 
      61           8 :         TestType fast(volumeDescriptor, sinoDescriptor);
      62           8 :         TestType slow(volumeDescriptor, sinoDescriptor, false);
      63             : 
      64           8 :         WHEN("Projector is cloned")
      65             :         {
      66           8 :             auto fastClone = fast.clone();
      67           8 :             auto slowClone = slow.clone();
      68           8 :             auto sinoClone = DataContainer<data_t>(sinoDescriptor);
      69           4 :             sinoClone = 0;
      70           8 :             auto volumeClone = DataContainer<data_t>(volumeDescriptor);
      71           4 :             volumeClone = 0;
      72             : 
      73           8 :             THEN("Results do not change (may still be slightly different due to summation being "
      74             :                  "performed in a different order)")
      75             :             {
      76           4 :                 fast.apply(volume, sino);
      77           4 :                 fastClone->apply(volume, sinoClone);
      78           4 :                 REQUIRE_UNARY(isCwiseApprox(sino, sinoClone));
      79             : 
      80           4 :                 slowClone->apply(volume, sinoClone);
      81           4 :                 REQUIRE_UNARY(isCwiseApprox(sino, sinoClone));
      82             : 
      83           4 :                 fast.applyAdjoint(sino, volume);
      84           4 :                 fastClone->applyAdjoint(sino, volumeClone);
      85           4 :                 REQUIRE_UNARY(isCwiseApprox(volume, volumeClone));
      86             : 
      87           4 :                 slow.applyAdjoint(sino, volume);
      88           4 :                 slowClone->applyAdjoint(sino, volumeClone);
      89           4 :                 REQUIRE_UNARY(isCwiseApprox(volume, volumeClone));
      90             :             }
      91             :         }
      92             :     }
      93           4 : }
      94             : 
      95         110 : TEST_CASE_TEMPLATE("JosephsMethodCUDA: 2D setup with a single ray", TestType,
      96             :                    JosephsMethodCUDA<float>, JosephsMethodCUDA<double>)
      97             : {
      98             :     // Turn logger of
      99          80 :     Logger::setLevel(Logger::LogLevel::OFF);
     100             : 
     101             :     using data_t = decltype(return_data_t(std::declval<TestType>()));
     102             : 
     103         128 :     GIVEN("A 5x5 volume and detector of size 1 with 1 angle")
     104             :     {
     105             :         // Domain setup
     106          48 :         const index_t volSize = 5;
     107             : 
     108          96 :         IndexVector_t volumeDims(2);
     109          48 :         volumeDims << volSize, volSize;
     110             : 
     111          96 :         VolumeDescriptor volumeDescriptor(volumeDims);
     112          96 :         DataContainer<data_t> volume(volumeDescriptor);
     113             : 
     114             :         // range setup
     115          48 :         const index_t detectorSize = 1;
     116          48 :         const index_t numImgs = 1;
     117             : 
     118          96 :         IndexVector_t sinoDims(2);
     119          48 :         sinoDims << detectorSize, numImgs;
     120             : 
     121             :         // Setup geometry
     122          48 :         auto stc = SourceToCenterOfRotation{20 * volSize};
     123          48 :         auto ctr = CenterOfRotationToDetector{volSize};
     124         144 :         auto volData = VolumeData2D{Size2D{volumeDims}};
     125         144 :         auto sinoData = SinogramData2D{Size2D{sinoDims}};
     126             : 
     127          96 :         std::vector<Geometry> geom;
     128             : 
     129          56 :         GIVEN("A basic geometry setup")
     130             :         {
     131           8 :             geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData));
     132             : 
     133          16 :             PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
     134          16 :             DataContainer<data_t> sino(sinoDescriptor);
     135             : 
     136          16 :             TestType slow(volumeDescriptor, sinoDescriptor, false);
     137          16 :             TestType fast(volumeDescriptor, sinoDescriptor);
     138             : 
     139          12 :             WHEN("Sinogram conatainer is not zero initialized and we project through an empty "
     140             :                  "volume")
     141             :             {
     142           4 :                 volume = 0;
     143           4 :                 sino = 1;
     144             : 
     145           8 :                 THEN("Result is zero")
     146             :                 {
     147           4 :                     fast.apply(volume, sino);
     148           4 :                     DataContainer<data_t> zero(sinoDescriptor);
     149           4 :                     zero = 0;
     150           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     151             : 
     152           4 :                     sino = 1;
     153           4 :                     slow.apply(volume, sino);
     154           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     155             :                 }
     156             :             }
     157             : 
     158          12 :             WHEN("Volume container is not zero initialized and we backproject from an empty "
     159             :                  "sinogram")
     160             :             {
     161           4 :                 sino = 0;
     162           4 :                 volume = 1;
     163             : 
     164           8 :                 THEN("Result is zero")
     165             :                 {
     166           4 :                     fast.applyAdjoint(sino, volume);
     167           4 :                     DataContainer<data_t> zero(volumeDescriptor);
     168           4 :                     zero = 0;
     169           4 :                     REQUIRE_UNARY(isCwiseApprox(volume, zero));
     170             : 
     171           4 :                     volume = 1;
     172           4 :                     slow.applyAdjoint(sino, volume);
     173           4 :                     REQUIRE_UNARY(isCwiseApprox(volume, zero));
     174             :                 }
     175             :             }
     176             :         }
     177             : 
     178          64 :         GIVEN("Rays not intersecting the bounding box are present")
     179             :         {
     180          20 :             WHEN("Tracing along a y-axis-aligned ray with a negative x-coordinate of origin")
     181             :             {
     182           4 :                 geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData),
     183             :                                   PrincipalPointOffset{}, RotationOffset2D{volSize, 0});
     184             : 
     185           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
     186           8 :                 DataContainer<data_t> sino(sinoDescriptor);
     187             : 
     188           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
     189           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
     190             : 
     191           8 :                 THEN("Result of forward projection is zero")
     192             :                 {
     193           4 :                     fast.apply(volume, sino);
     194           8 :                     DataContainer<data_t> zero(sinoDescriptor);
     195           4 :                     zero = 0;
     196           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     197             : 
     198           4 :                     slow.apply(volume, sino);
     199           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     200             : 
     201           8 :                     AND_THEN("Result of backprojection is zero")
     202             :                     {
     203           4 :                         fast.applyAdjoint(sino, volume);
     204           4 :                         DataContainer<data_t> zero(volumeDescriptor);
     205           4 :                         zero = 0;
     206           4 :                         REQUIRE_UNARY(isCwiseApprox(volume, zero));
     207             : 
     208           4 :                         slow.applyAdjoint(sino, volume);
     209           4 :                         REQUIRE_UNARY(isCwiseApprox(volume, zero));
     210             :                     }
     211             :                 }
     212             :             }
     213             : 
     214          20 :             WHEN("Tracing along a y-axis-aligned ray with a x-coordinate of origin beyond the "
     215             :                  "bounding "
     216             :                  "box")
     217             :             {
     218           4 :                 geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData),
     219             :                                   PrincipalPointOffset{}, RotationOffset2D{-volSize, 0});
     220             : 
     221           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
     222           8 :                 DataContainer<data_t> sino(sinoDescriptor);
     223             : 
     224           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
     225           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
     226             : 
     227           8 :                 THEN("Result of forward projection is zero")
     228             :                 {
     229           8 :                     DataContainer<data_t> zero(sinoDescriptor);
     230           4 :                     zero = 0;
     231             : 
     232           4 :                     fast.apply(volume, sino);
     233           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     234             : 
     235           4 :                     slow.apply(volume, sino);
     236           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     237             : 
     238           8 :                     AND_THEN("Result of backprojection is zero")
     239             :                     {
     240           4 :                         DataContainer<data_t> zero(volumeDescriptor);
     241           4 :                         zero = 0;
     242             : 
     243           4 :                         fast.applyAdjoint(sino, volume);
     244           4 :                         REQUIRE_UNARY(isCwiseApprox(volume, zero));
     245             : 
     246           4 :                         slow.applyAdjoint(sino, volume);
     247           4 :                         REQUIRE_UNARY(isCwiseApprox(volume, zero));
     248             :                     }
     249             :                 }
     250             :             }
     251             : 
     252          20 :             WHEN("Tracing along a x-axis-aligned ray with a negative y-coordinate of origin")
     253             :             {
     254           8 :                 geom.emplace_back(stc, ctr, Radian{pi_t / 2}, std::move(volData),
     255           4 :                                   std::move(sinoData), PrincipalPointOffset{},
     256             :                                   RotationOffset2D{0, volSize});
     257             : 
     258           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
     259           8 :                 DataContainer<data_t> sino(sinoDescriptor);
     260             : 
     261           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
     262           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
     263             : 
     264           8 :                 THEN("Result of forward projection is zero")
     265             :                 {
     266           8 :                     DataContainer<data_t> zero(sinoDescriptor);
     267           4 :                     zero = 0;
     268             : 
     269           4 :                     fast.apply(volume, sino);
     270           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     271             : 
     272           4 :                     slow.apply(volume, sino);
     273           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     274             : 
     275           8 :                     AND_THEN("Result of backprojection is zero")
     276             :                     {
     277           4 :                         DataContainer<data_t> zero(volumeDescriptor);
     278           4 :                         zero = 0;
     279             : 
     280           4 :                         fast.applyAdjoint(sino, volume);
     281           4 :                         REQUIRE_UNARY(isCwiseApprox(volume, zero));
     282             : 
     283           4 :                         slow.applyAdjoint(sino, volume);
     284           4 :                         REQUIRE_UNARY(isCwiseApprox(volume, zero));
     285             :                     }
     286             :                 }
     287             :             }
     288             : 
     289          20 :             WHEN("Tracing along a x-axis-aligned ray with a y-coordinate of origin beyond the "
     290             :                  "bounding "
     291             :                  "box")
     292             :             {
     293           8 :                 geom.emplace_back(stc, ctr, Radian{pi_t / 2}, std::move(volData),
     294           4 :                                   std::move(sinoData), PrincipalPointOffset{},
     295             :                                   RotationOffset2D{0, -volSize});
     296             : 
     297           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
     298           8 :                 DataContainer<data_t> sino(sinoDescriptor);
     299             : 
     300           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
     301           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
     302             : 
     303           8 :                 THEN("Result of forward projection is zero")
     304             :                 {
     305           8 :                     DataContainer<data_t> zero(sinoDescriptor);
     306           4 :                     zero = 0;
     307             : 
     308           4 :                     fast.apply(volume, sino);
     309           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     310             : 
     311           4 :                     slow.apply(volume, sino);
     312           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     313             : 
     314           8 :                     AND_THEN("Result of backprojection is zero")
     315             :                     {
     316           4 :                         DataContainer<data_t> zero(volumeDescriptor);
     317           4 :                         zero = 0;
     318             : 
     319           4 :                         fast.applyAdjoint(sino, volume);
     320           4 :                         REQUIRE_UNARY(isCwiseApprox(volume, zero));
     321             : 
     322           4 :                         slow.applyAdjoint(sino, volume);
     323           4 :                         REQUIRE_UNARY(isCwiseApprox(volume, zero));
     324             :                     }
     325             :                 }
     326             :             }
     327             :         }
     328             : 
     329             :         // Expected results
     330         288 :         Eigen::Matrix<data_t, Eigen::Dynamic, 1> backProjections[2];
     331          48 :         backProjections[0].resize(volSize * volSize);
     332          48 :         backProjections[1].resize(volSize * volSize);
     333             : 
     334          48 :         constexpr index_t numCases = 4;
     335          48 :         const std::array<real_t, numCases> angles = {0., 90., 180., 270.};
     336             : 
     337             :         // Testing for axis-aligned rays
     338          56 :         GIVEN("Given rays along the axes")
     339             :         {
     340             : 
     341          16 :             backProjections[0] << 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0,
     342          16 :                 1, 0, 0;
     343          16 :             backProjections[1] << 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0,
     344          16 :                 0, 0, 0;
     345             : 
     346          40 :             for (std::size_t i = 0; i < numCases; i++) {
     347          36 :                 WHEN("An axis-aligned ray with a fixed angle passes through the center of a pixel")
     348             :                 {
     349           8 :                     INFO("An axis-aligned ray with an angle of ", angles[i],
     350             :                          " radians passes through the center of a pixel");
     351           4 :                     geom.emplace_back(stc, ctr, Degree{angles[i]}, std::move(volData),
     352           4 :                                       std::move(sinoData));
     353             : 
     354           8 :                     PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
     355           8 :                     DataContainer<data_t> sino(sinoDescriptor);
     356             : 
     357           8 :                     TestType fast(volumeDescriptor, sinoDescriptor);
     358           8 :                     TestType slow(volumeDescriptor, sinoDescriptor, false);
     359             : 
     360           8 :                     THEN("The result of projecting through a pixel is exactly the pixel value")
     361             :                     {
     362          24 :                         for (index_t j = 0; j < volSize; j++) {
     363          20 :                             volume = 0;
     364          20 :                             if (i % 2 == 0)
     365          20 :                                 volume(volSize / 2, j) = 1;
     366             :                             else
     367           0 :                                 volume(j, volSize / 2) = 1;
     368             : 
     369          20 :                             fast.apply(volume, sino);
     370             :                             // Using doubles significantly increases interpolation accuracy
     371             :                             // For example: when using floats, points very near the border (less
     372             :                             // than ~1/500th of a pixel away from the border) are rounded to
     373             :                             // actually lie on the border. This then yields more accurate results
     374             :                             // when using floats in some of the axis-aligned test cases, despite the
     375             :                             // lower interpolation accuracy.
     376             :                             //
     377             :                             //  => different requirements for floats and doubles, looser
     378             :                             // requirements for doubles
     379          20 :                             REQUIRE_EQ(sino[0], Approx(1.0));
     380             : 
     381          20 :                             slow.apply(volume, sino);
     382          20 :                             REQUIRE_EQ(sino[0], Approx(1.0));
     383             :                         }
     384             : 
     385           8 :                         AND_THEN(
     386             :                             "The backprojection sets the values of all hit pixels to the detector "
     387             :                             "value")
     388             :                         {
     389           4 :                             fast.applyAdjoint(sino, volume);
     390           4 :                             REQUIRE_UNARY(isCwiseApprox(
     391             :                                 volume,
     392             :                                 DataContainer<data_t>(volumeDescriptor, backProjections[i % 2])));
     393             : 
     394           4 :                             slow.applyAdjoint(sino, volume);
     395           4 :                             REQUIRE_UNARY(isCwiseApprox(
     396             :                                 volume,
     397             :                                 DataContainer<data_t>(volumeDescriptor, backProjections[i % 2])));
     398             :                         }
     399             :                     }
     400             :                 }
     401             :             }
     402             :         }
     403             : 
     404             :         // Move the source far far back, so rays become nearly parallel
     405          48 :         stc = SourceToCenterOfRotation{volSize * 2000};
     406             : 
     407          56 :         GIVEN("Off center rays along axes")
     408             :         {
     409           8 :             std::array<real_t, numCases> offsetx = {-0.25, 0.0, -0.25, 0.0};
     410           8 :             std::array<real_t, numCases> offsety = {0.0, -0.25, 0.0, -0.25};
     411             : 
     412          16 :             backProjections[0] << 0, 0.25, 0.75, 0, 0, 0, 0.25, 0.75, 0, 0, 0, 0.25, 0.75, 0, 0, 0,
     413          16 :                 0.25, 0.75, 0, 0, 0, 0.25, 0.75, 0, 0;
     414             : 
     415          16 :             backProjections[1] << 0, 0, 0, 0, 0, 0.25, 0.25, 0.25, 0.25, 0.25, 0.75, 0.75, 0.75,
     416          16 :                 0.75, 0.75, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
     417             : 
     418          40 :             for (std::size_t i = 0; i < numCases; i++) {
     419          36 :                 WHEN("An axis-aligned ray with a fixed angle, which does not pass through the "
     420             :                      "center of a pixel")
     421             :                 {
     422           8 :                     INFO("An axis-aligned ray with an angle of ", angles[i],
     423             :                          " radians not passing through the center of a pixel");
     424           8 :                     geom.emplace_back(stc, ctr, Degree{angles[i]}, std::move(volData),
     425           4 :                                       std::move(sinoData), PrincipalPointOffset{0},
     426             :                                       RotationOffset2D{offsetx[i], offsety[i]});
     427             : 
     428           8 :                     PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
     429           8 :                     DataContainer<data_t> sino(sinoDescriptor);
     430             : 
     431           8 :                     TestType fast(volumeDescriptor, sinoDescriptor);
     432           8 :                     TestType slow(volumeDescriptor, sinoDescriptor, false);
     433             : 
     434           8 :                     THEN("The result of projecting through a pixel is the interpolated value "
     435             :                          "between "
     436             :                          "the two pixels closest to the ray")
     437             :                     {
     438          24 :                         for (index_t j = 0; j < volSize; j++) {
     439          20 :                             volume = 0;
     440          20 :                             if (i % 2 == 0)
     441          20 :                                 volume(volSize / 2, j) = 1;
     442             :                             else
     443           0 :                                 volume(j, volSize / 2) = 1;
     444             : 
     445          20 :                             fast.apply(volume, sino);
     446          20 :                             REQUIRE_EQ(sino[0], Approx(0.75));
     447             : 
     448          20 :                             slow.apply(volume, sino);
     449          20 :                             REQUIRE_EQ(sino[0], Approx(0.75));
     450             :                         }
     451             : 
     452           8 :                         AND_THEN("The slow backprojection yields the exact adjoint, the fast "
     453             :                                  "backprojection "
     454             :                                  "also yields the exact adjoint for a very distant x-ray source")
     455             :                         {
     456           4 :                             sino[0] = 1;
     457           4 :                             slow.applyAdjoint(sino, volume);
     458           4 :                             REQUIRE_UNARY(isCwiseApprox(
     459             :                                 volume,
     460             :                                 DataContainer<data_t>(volumeDescriptor, backProjections[i % 2])));
     461             : 
     462           4 :                             fast.applyAdjoint(sino, volume);
     463             : 
     464             :                             // Using doubles significantly increases interpolation accuracy
     465             :                             //  For example: when using floats, points very near the border (less
     466             :                             //  than
     467             :                             // ~1/500th of a pixel away from the border) are rounded to actually lie
     468             :                             // on the border. This then yields more accurate results when using
     469             :                             // floats in some of the axis-aligned test cases, despite the lower
     470             :                             // interpolation accuracy.
     471             :                             //
     472             :                             //  => different requirements for floats and doubles, looser
     473             :                             //  requirements
     474             :                             // for doubles
     475             :                             if constexpr (std::is_same_v<data_t, float>)
     476           2 :                                 REQUIRE_UNARY(isCwiseApprox(
     477             :                                     volume, DataContainer<data_t>(volumeDescriptor,
     478             :                                                                   backProjections[i % 2])));
     479             :                             else
     480           2 :                                 REQUIRE_UNARY(isApprox(
     481             :                                     volume,
     482             :                                     DataContainer<data_t>(volumeDescriptor, backProjections[i % 2]),
     483             :                                     static_cast<real_t>(0.001)));
     484             :                         }
     485             :                     }
     486             :                 }
     487             :             }
     488             :         }
     489             : 
     490          52 :         GIVEN("Rays going running allong right bolume boundary")
     491             :         {
     492           8 :             backProjections[0] << 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0,
     493           8 :                 0, 0, 1;
     494             : 
     495           8 :             WHEN("A y-axis-aligned ray runs along the right volume boundary")
     496             :             {
     497           4 :                 geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData),
     498             :                                   PrincipalPointOffset{0}, RotationOffset2D{volSize * 0.5, 0});
     499             : 
     500           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
     501           8 :                 DataContainer<data_t> sino(sinoDescriptor);
     502             : 
     503           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
     504           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
     505             : 
     506           8 :                 THEN("The result of projecting through a pixel is exactly the pixel's value (we "
     507             :                      "mirror "
     508             :                      "values at the border for the purpose of interpolation)")
     509             :                 {
     510          24 :                     for (index_t j = 0; j < volSize; j++) {
     511          20 :                         volume = 0;
     512          20 :                         volume(volSize - 1, j) = 1;
     513             : 
     514          20 :                         fast.apply(volume, sino);
     515          20 :                         REQUIRE_EQ(sino[0], Approx(1));
     516             : 
     517          20 :                         slow.apply(volume, sino);
     518          20 :                         REQUIRE_EQ(sino[0], Approx(1));
     519             :                     }
     520             : 
     521           8 :                     AND_THEN(
     522             :                         "The slow backprojection yields the exact adjoint, the fast backprojection "
     523             :                         "also yields the exact adjoint for a very distant x-ray source")
     524             :                     {
     525           4 :                         sino[0] = 1;
     526           4 :                         slow.applyAdjoint(sino, volume);
     527           4 :                         REQUIRE_UNARY(isCwiseApprox(
     528             :                             volume, DataContainer<data_t>(volumeDescriptor, backProjections[0])));
     529             : 
     530           4 :                         fast.applyAdjoint(sino, volume);
     531             :                         // Using doubles significantly increases interpolation accuracy
     532             :                         //  For example: when using floats, points very near the border (less than
     533             :                         // ~1/500th of a pixel away from the border) are rounded to actually lie on
     534             :                         // the border. This then yields more accurate results when using floats in
     535             :                         // some of the axis-aligned test cases, despite the lower interpolation
     536             :                         // accuracy.
     537             :                         //
     538             :                         //  => different requirements for floats and doubles, looser requirements
     539             :                         // for doubles
     540             :                         if constexpr (std::is_same_v<data_t, float>)
     541           2 :                             REQUIRE_UNARY(isCwiseApprox(
     542             :                                 volume, DataContainer<data_t>(volumeDescriptor,
     543             :                                                               (backProjections[0] / 2).eval())));
     544             :                         else
     545           2 :                             REQUIRE_UNARY(
     546             :                                 isApprox(volume,
     547             :                                          DataContainer<data_t>(volumeDescriptor,
     548             :                                                                (backProjections[0] / 2).eval()),
     549             :                                          static_cast<real_t>(0.001)));
     550             :                     }
     551             :                 }
     552             :             }
     553             :         }
     554             : 
     555          52 :         GIVEN("Rays going running along left bolume boundary")
     556             :         {
     557           8 :             backProjections[0] << 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0,
     558           8 :                 0, 0, 0;
     559             : 
     560           8 :             WHEN("A y-axis-aligned ray runs along the left volume boundary")
     561             :             {
     562           4 :                 geom.emplace_back(stc, ctr, Radian{0}, std::move(volData), std::move(sinoData),
     563             :                                   PrincipalPointOffset{0}, RotationOffset2D{-volSize * 0.5, 0});
     564             : 
     565           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
     566           8 :                 DataContainer<data_t> sino(sinoDescriptor);
     567             : 
     568           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
     569           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
     570             : 
     571             :                 // TestType fast(volumeDescriptor, sinoDescriptor, geom);
     572             :                 // TestType slow(volumeDescriptor, sinoDescriptor, geom, false);
     573           8 :                 THEN("The result of projecting through a pixel is exactly the pixel's value (we "
     574             :                      "mirror "
     575             :                      "values at the border for the purpose of interpolation)")
     576             :                 {
     577          24 :                     for (index_t j = 0; j < volSize; j++) {
     578          20 :                         volume = 0;
     579          20 :                         volume(0, j) = 1;
     580             : 
     581          20 :                         fast.apply(volume, sino);
     582          20 :                         REQUIRE_EQ(sino[0], Approx(1));
     583             : 
     584          20 :                         slow.apply(volume, sino);
     585          20 :                         REQUIRE_EQ(sino[0], Approx(1));
     586             :                     }
     587             : 
     588           8 :                     AND_THEN(
     589             :                         "The slow backprojection yields the exact adjoint, the fast backprojection "
     590             :                         "also yields the exact adjoint for a very distant x-ray source")
     591             :                     {
     592           4 :                         sino[0] = 1;
     593           4 :                         slow.applyAdjoint(sino, volume);
     594           4 :                         REQUIRE_UNARY(isCwiseApprox(
     595             :                             volume, DataContainer<data_t>(volumeDescriptor, backProjections[0])));
     596             : 
     597           4 :                         fast.applyAdjoint(sino, volume);
     598             :                         /** Using doubles significantly increases interpolation accuracy
     599             :                          *  For example: when using floats, points very near the border (less than
     600             :                          * ~1/500th of a pixel away from the border) are rounded to actually lie on
     601             :                          * the border. This then yields more accurate results when using floats in
     602             :                          * some of the axis-aligned test cases, despite the lower interpolation
     603             :                          * accuracy.
     604             :                          *
     605             :                          *  => different requirements for floats and doubles, looser requirements
     606             :                          * for doubles
     607             :                          */
     608             :                         if constexpr (std::is_same_v<data_t, float>)
     609           2 :                             REQUIRE_UNARY(isCwiseApprox(
     610             :                                 volume, DataContainer<data_t>(volumeDescriptor,
     611             :                                                               (backProjections[0] / 2).eval())));
     612             :                         else
     613           2 :                             REQUIRE_UNARY(
     614             :                                 isApprox(volume,
     615             :                                          DataContainer<data_t>(volumeDescriptor,
     616             :                                                                (backProjections[0] / 2).eval()),
     617             :                                          static_cast<real_t>(0.001)));
     618             :                     }
     619             :                 }
     620             :             }
     621             :         }
     622             :     }
     623             : 
     624         112 :     GIVEN("A 4x4 volume and detector of size 1 with 1 angle at -30 degree")
     625             :     {
     626             :         // Domain setup
     627          32 :         const index_t volSize = 4;
     628             : 
     629          64 :         IndexVector_t volumeDims(2);
     630          32 :         volumeDims << volSize, volSize;
     631             : 
     632          64 :         VolumeDescriptor volumeDescriptor(volumeDims);
     633          64 :         DataContainer<data_t> volume(volumeDescriptor);
     634             : 
     635             :         // range setup
     636          32 :         const index_t detectorSize = 1;
     637          32 :         const index_t numImgs = 1;
     638             : 
     639          64 :         IndexVector_t sinoDims(2);
     640          32 :         sinoDims << detectorSize, numImgs;
     641             : 
     642             :         // Setup geometry
     643          32 :         const auto stc = SourceToCenterOfRotation{20 * volSize};
     644          32 :         const auto ctr = CenterOfRotationToDetector{volSize};
     645          96 :         auto volData = VolumeData2D{Size2D{volumeDims}};
     646          96 :         auto sinoData = SinogramData2D{Size2D{sinoDims}};
     647             : 
     648          64 :         std::vector<Geometry> geom;
     649             : 
     650          32 :         const real_t sqrt3r = std::sqrt(static_cast<real_t>(3));
     651          32 :         const data_t sqrt3d = std::sqrt(static_cast<data_t>(3));
     652          32 :         const data_t halfd = static_cast<data_t>(0.5);
     653          32 :         const data_t thirdd = static_cast<data_t>(1.0 / 3);
     654             : 
     655          48 :         GIVEN("An angle of -30 degrees")
     656             :         {
     657          16 :             const auto angle = Degree{-30};
     658             : 
     659          20 :             WHEN("Projections goes through center of volume")
     660             :             {
     661           4 :                 geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData));
     662             : 
     663           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
     664           8 :                 DataContainer<data_t> sino(sinoDescriptor);
     665           4 :                 sino = 0;
     666             : 
     667           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
     668           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
     669             : 
     670           4 :                 real_t weight = 2 / sqrt3r;
     671             : 
     672           8 :                 THEN("Ray intersects the correct pixels")
     673             :                 {
     674           4 :                     volume = 1;
     675           4 :                     volume(3, 0) = 0;
     676           4 :                     volume(1, 1) = 0;
     677           4 :                     volume(1, 2) = 0;
     678           4 :                     volume(1, 3) = 0;
     679             : 
     680           4 :                     volume(2, 0) = 0;
     681           4 :                     volume(2, 1) = 0;
     682           4 :                     volume(2, 2) = 0;
     683           4 :                     volume(0, 3) = 0;
     684             : 
     685           4 :                     fast.apply(volume, sino);
     686           8 :                     DataContainer<data_t> zero(sinoDescriptor);
     687           4 :                     zero = 0;
     688           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     689             : 
     690           4 :                     slow.apply(volume, sino);
     691           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     692             : 
     693           8 :                     AND_THEN("The correct weighting is applied")
     694             :                     {
     695           4 :                         volume(3, 0) = 1;
     696           4 :                         volume(1, 1) = 1;
     697           4 :                         volume(1, 2) = 1;
     698           4 :                         volume(1, 3) = 1;
     699             : 
     700           4 :                         fast.apply(volume, sino);
     701           4 :                         REQUIRE_EQ(sino[0], Approx(2 * weight));
     702             : 
     703           4 :                         slow.apply(volume, sino);
     704           4 :                         REQUIRE_EQ(sino[0], Approx(2 * weight));
     705             : 
     706           4 :                         sino[0] = 1;
     707             : 
     708           8 :                         Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
     709           8 :                         slowExpected << 0, 0, (3 - sqrt3d) / 2, (sqrt3d - 1) / 2, 0,
     710           4 :                             (sqrt3d - 1) / (2 * sqrt3d), (sqrt3d + 1) / (2 * sqrt3d), 0, 0,
     711           4 :                             (sqrt3d + 1) / (2 * sqrt3d), (sqrt3d - 1) / (2 * sqrt3d), 0,
     712           8 :                             (sqrt3d - 1) / 2, (3 - sqrt3d) / 2, 0, 0;
     713             : 
     714           4 :                         slowExpected *= weight;
     715           4 :                         slow.applyAdjoint(sino, volume);
     716           4 :                         REQUIRE_UNARY(isCwiseApprox(
     717             :                             volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
     718             : 
     719           4 :                         fast.applyAdjoint(sino, volume);
     720          20 :                         for (real_t i = 0.5; i < volSize; i += 1) {
     721          80 :                             for (real_t j = 0.5; j < volSize; j += 1) {
     722             :                                 const real_t angle =
     723          64 :                                     std::abs(std::atan((sqrt3r * volSize * 10
     724          64 :                                                         - static_cast<real_t>(volSize / 2.0) + j)
     725             :                                                        / (volSize * 10
     726          64 :                                                           + static_cast<real_t>(volSize / 2.0) - i))
     727             :                                              - pi_t / 3);
     728          64 :                                 const real_t len = volSize * 21 * std::tan(angle);
     729          64 :                                 if (len < 1) {
     730          32 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j),
     731             :                                                Approx(1 - len).epsilon(0.005));
     732             :                                 } else {
     733          32 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
     734             :                                 }
     735             :                             }
     736             :                         }
     737             :                     }
     738             :                 }
     739             :             }
     740             : 
     741          20 :             WHEN("Ray enters through the right border")
     742             :             {
     743             :                 // In this case the ray exits through a border along the main ray direction, but
     744             :                 // enters through a border not along the main direction First pixel should be
     745             :                 // weighted differently
     746           4 :                 geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
     747             :                                   PrincipalPointOffset{0}, RotationOffset2D{sqrt3r, 0});
     748             : 
     749           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
     750           8 :                 DataContainer<data_t> sino(sinoDescriptor);
     751           4 :                 sino = 0;
     752             : 
     753           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
     754           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
     755             : 
     756             :                 // TestType fast(volumeDescriptor, sinoDescriptor, geom);
     757             :                 // TestType slow(volumeDescriptor, sinoDescriptor, geom, false);
     758             : 
     759           8 :                 THEN("Ray intersects the correct pixels")
     760             :                 {
     761           4 :                     volume = 1;
     762           4 :                     volume(3, 1) = 0;
     763           4 :                     volume(3, 2) = 0;
     764           4 :                     volume(3, 3) = 0;
     765           4 :                     volume(2, 3) = 0;
     766           4 :                     volume(2, 2) = 0;
     767             : 
     768           4 :                     fast.apply(volume, sino);
     769           8 :                     DataContainer<data_t> zero(sinoDescriptor);
     770           4 :                     zero = 0;
     771           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     772             : 
     773           4 :                     slow.apply(volume, sino);
     774           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     775             : 
     776           8 :                     AND_THEN("The correct weighting is applied")
     777             :                     {
     778           4 :                         volume(3, 1) = 1;
     779           4 :                         volume(2, 2) = 1;
     780           4 :                         volume(2, 3) = 1;
     781             : 
     782           4 :                         fast.apply(volume, sino);
     783           4 :                         REQUIRE_EQ(sino[0], Approx((4 - 2 * sqrt3d) * (sqrt3d - 1)
     784             :                                                    + (2 / sqrt3d) * (3 - 8 * sqrt3d / 6))
     785             :                                                 .epsilon(0.005));
     786             : 
     787           4 :                         slow.apply(volume, sino);
     788           4 :                         REQUIRE_EQ(sino[0], Approx((4 - 2 * sqrt3d) * (sqrt3d - 1)
     789             :                                                    + (2 / sqrt3d) * (3 - 8 * sqrt3d / 6))
     790             :                                                 .epsilon(0.005));
     791             : 
     792           4 :                         sino[0] = 1;
     793             : 
     794           8 :                         Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
     795           8 :                         slowExpected << 0, 0, 0, 0, 0, 0, 0, (4 - 2 * sqrt3d) * (sqrt3d - 1), 0, 0,
     796           4 :                             (2 / sqrt3d) * (1 + halfd - 5 * sqrt3d / 6),
     797             :                             (4 - 2 * sqrt3d) * (2 - sqrt3d)
     798           4 :                                 + (2 / sqrt3d) * (5 * sqrt3d / 6 - halfd),
     799           4 :                             0, 0, (2 / sqrt3d) * (1 + halfd - sqrt3d / 2),
     800           4 :                             (2 / sqrt3d) * (sqrt3d / 2 - halfd);
     801             : 
     802           4 :                         slow.applyAdjoint(sino, volume);
     803           4 :                         REQUIRE_UNARY(isCwiseApprox(
     804             :                             volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
     805             : 
     806           4 :                         fast.applyAdjoint(sino, volume);
     807          20 :                         for (real_t i = 0.5; i < volSize; i += 1) {
     808          80 :                             for (real_t j = 0.5; j < volSize; j += 1) {
     809             :                                 const real_t angle =
     810          64 :                                     std::abs(std::atan((40 * sqrt3r - 2 + j) / (42 + sqrt3r - i))
     811             :                                              - pi_t / 3);
     812          64 :                                 const real_t len = 84 * std::tan(angle);
     813          64 :                                 if (len < 1) {
     814          20 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j),
     815             :                                                Approx(1 - len).epsilon(0.01));
     816             :                                 } else {
     817          44 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
     818             :                                 }
     819             :                             }
     820             :                         }
     821             :                     }
     822             :                 }
     823             :             }
     824             : 
     825          20 :             WHEN("Ray exits through the left border")
     826             :             {
     827             :                 // In this case the ray enters through a border along the main ray direction, but
     828             :                 // exits through a border not along the main direction Last pixel should be weighted
     829             :                 // differently
     830           4 :                 geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
     831             :                                   PrincipalPointOffset{0}, RotationOffset2D{-sqrt3r, 0});
     832             : 
     833           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
     834           8 :                 DataContainer<data_t> sino(sinoDescriptor);
     835           4 :                 sino = 0;
     836             : 
     837           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
     838           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
     839             : 
     840           8 :                 THEN("Ray intersects the correct pixels")
     841             :                 {
     842           4 :                     volume = 1;
     843           4 :                     volume(0, 0) = 0;
     844           4 :                     volume(1, 0) = 0;
     845           4 :                     volume(0, 1) = 0;
     846           4 :                     volume(1, 1) = 0;
     847           4 :                     volume(0, 2) = 0;
     848             : 
     849           4 :                     fast.apply(volume, sino);
     850           8 :                     DataContainer<data_t> zero(sinoDescriptor);
     851           4 :                     zero = 0;
     852           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     853             : 
     854           4 :                     slow.apply(volume, sino);
     855           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     856             : 
     857           8 :                     AND_THEN("The correct weighting is applied")
     858             :                     {
     859           4 :                         volume(1, 0) = 1;
     860           4 :                         volume(0, 1) = 1;
     861             : 
     862           4 :                         fast.apply(volume, sino);
     863           4 :                         REQUIRE_EQ(sino[0], Approx((sqrt3d - 1) + (5.0 / 3.0 - 1 / sqrt3d)
     864             :                                                    + (4 - 2 * sqrt3d) * (2 - sqrt3d))
     865             :                                                 .epsilon(0.005));
     866             : 
     867           4 :                         slow.apply(volume, sino);
     868           4 :                         REQUIRE_EQ(sino[0], Approx((sqrt3d - 1) + (5.0 / 3.0 - 1 / sqrt3d)
     869             :                                                    + (4 - 2 * sqrt3d) * (2 - sqrt3d))
     870             :                                                 .epsilon(0.005));
     871             : 
     872           4 :                         sino[0] = 1;
     873             : 
     874           8 :                         Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
     875           8 :                         slowExpected << 1 - 1 / sqrt3d, sqrt3d - 1, 0, 0,
     876           4 :                             (5 * thirdd - 1 / sqrt3d) + (4 - 2 * sqrt3d) * (2 - sqrt3d),
     877           4 :                             sqrt3d - 5 * thirdd, 0, 0, (sqrt3d - 1) * (4 - 2 * sqrt3d), 0, 0, 0, 0,
     878           8 :                             0, 0, 0;
     879             : 
     880           4 :                         slow.applyAdjoint(sino, volume);
     881           4 :                         REQUIRE_UNARY(isCwiseApprox(
     882             :                             volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
     883             : 
     884           4 :                         fast.applyAdjoint(sino, volume);
     885             : 
     886          20 :                         for (real_t i = 0.5; i < volSize; i += 1) {
     887          80 :                             for (real_t j = 0.5; j < volSize; j += 1) {
     888             :                                 const real_t angle =
     889          64 :                                     std::abs(std::atan((40 * sqrt3r - 2 + j) / (42 - sqrt3r - i))
     890             :                                              - pi_t / 3);
     891          64 :                                 const real_t len = 84 * std::tan(angle);
     892             : 
     893          64 :                                 if (len < 1) {
     894          24 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j),
     895             :                                                Approx(1 - len).epsilon(0.002));
     896             :                                 } else {
     897          40 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
     898             :                                 }
     899             :                             }
     900             :                         }
     901             :                     }
     902             :                 }
     903             :             }
     904             : 
     905          20 :             WHEN("Ray only intersects a single pixel")
     906             :             {
     907             :                 // This is a special case that is handled separately in both forward and
     908             :                 // backprojection
     909           4 :                 geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
     910             :                                   PrincipalPointOffset{0}, RotationOffset2D{-2 - sqrt3r / 2, 0});
     911             : 
     912           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
     913           8 :                 DataContainer<data_t> sino(sinoDescriptor);
     914           4 :                 sino = 0;
     915             : 
     916           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
     917           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
     918             : 
     919           8 :                 THEN("Ray intersects the correct pixels")
     920             :                 {
     921           4 :                     volume = 1;
     922           4 :                     volume(0, 0) = 0;
     923             : 
     924           4 :                     fast.apply(volume, sino);
     925           8 :                     DataContainer<data_t> zero(sinoDescriptor);
     926           4 :                     zero = 0;
     927           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     928             : 
     929           4 :                     slow.apply(volume, sino);
     930           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
     931             : 
     932           8 :                     AND_THEN("The correct weighting is applied")
     933             :                     {
     934           4 :                         volume(0, 0) = 1;
     935             : 
     936           4 :                         fast.apply(volume, sino);
     937           4 :                         REQUIRE_EQ(sino[0], Approx(1 / sqrt3d).epsilon(0.005));
     938             : 
     939           4 :                         slow.apply(volume, sino);
     940           4 :                         REQUIRE_EQ(sino[0], Approx(1 / sqrt3d).epsilon(0.005));
     941             : 
     942           4 :                         sino[0] = 1;
     943             : 
     944           8 :                         Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
     945           4 :                         slowExpected << 1 / sqrt3d, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
     946             : 
     947           4 :                         slow.applyAdjoint(sino, volume);
     948           4 :                         REQUIRE_UNARY(isCwiseApprox(
     949             :                             volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
     950             : 
     951           4 :                         fast.applyAdjoint(sino, volume);
     952             : 
     953          20 :                         for (real_t i = 0.5; i < volSize; i += 1) {
     954          80 :                             for (real_t j = 0.5; j < volSize; j += 1) {
     955          64 :                                 const real_t angle = std::abs(
     956          64 :                                     std::atan((40 * sqrt3r - 2 + j) / (40 - sqrt3r / 2 - i))
     957             :                                     - pi_t / 3);
     958          64 :                                 const real_t len = 84 * std::tan(angle);
     959             : 
     960          64 :                                 if (len < 1) {
     961           8 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j),
     962             :                                                Approx(1 - len).epsilon(0.002));
     963             :                                 } else {
     964          56 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
     965             :                                 }
     966             :                             }
     967             :                         }
     968             :                     }
     969             :                 }
     970             :             }
     971             :         }
     972             : 
     973          48 :         GIVEN("An angle of -120 degrees")
     974             :         {
     975          16 :             const auto angle = Degree{-120};
     976             : 
     977          20 :             WHEN("Ray goes through center of volume")
     978             :             {
     979             :                 // In this case the ray enters and exits the volume through the borders along the
     980             :                 // main direction Weighting for all interpolated values should be the same
     981           4 :                 geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData));
     982           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
     983           8 :                 DataContainer<data_t> sino(sinoDescriptor);
     984           4 :                 sino = 0;
     985             : 
     986           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
     987           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
     988             : 
     989           4 :                 real_t weight = 2 / sqrt3r;
     990           8 :                 THEN("Ray intersects the correct pixels")
     991             :                 {
     992           4 :                     sino[0] = 1;
     993           4 :                     slow.applyAdjoint(sino, volume);
     994             : 
     995           4 :                     volume = 1;
     996           4 :                     volume(0, 0) = 0;
     997           4 :                     volume(0, 1) = 0;
     998           4 :                     volume(1, 1) = 0;
     999           4 :                     volume(1, 2) = 0;
    1000             : 
    1001           4 :                     volume(2, 1) = 0;
    1002           4 :                     volume(2, 2) = 0;
    1003           4 :                     volume(3, 2) = 0;
    1004           4 :                     volume(3, 3) = 0;
    1005             : 
    1006           4 :                     fast.apply(volume, sino);
    1007           8 :                     DataContainer<data_t> zero(sinoDescriptor);
    1008           4 :                     zero = 0;
    1009           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
    1010             : 
    1011           4 :                     slow.apply(volume, sino);
    1012           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
    1013             : 
    1014           8 :                     AND_THEN("The correct weighting is applied")
    1015             :                     {
    1016           4 :                         volume(3, 3) = 1;
    1017           4 :                         volume(0, 1) = 1;
    1018           4 :                         volume(1, 1) = 1;
    1019           4 :                         volume(2, 1) = 1;
    1020             : 
    1021           4 :                         fast.apply(volume, sino);
    1022           4 :                         REQUIRE_EQ(sino[0], Approx(2 * weight));
    1023             : 
    1024           4 :                         slow.apply(volume, sino);
    1025           4 :                         REQUIRE_EQ(sino[0], Approx(2 * weight));
    1026             : 
    1027           4 :                         sino[0] = 1;
    1028             : 
    1029           8 :                         Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
    1030             : 
    1031           8 :                         slowExpected << (sqrt3d - 1) / 2, 0, 0, 0, (3 - sqrt3d) / 2,
    1032           4 :                             (sqrt3d + 1) / (2 * sqrt3d), (sqrt3d - 1) / (2 * sqrt3d), 0, 0,
    1033           4 :                             (sqrt3d - 1) / (2 * sqrt3d), (sqrt3d + 1) / (2 * sqrt3d),
    1034           8 :                             (3 - sqrt3d) / 2, 0, 0, 0, (sqrt3d - 1) / 2;
    1035             : 
    1036           4 :                         slowExpected *= weight;
    1037           4 :                         slow.applyAdjoint(sino, volume);
    1038           4 :                         REQUIRE_UNARY(isCwiseApprox(
    1039             :                             volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
    1040             : 
    1041           4 :                         fast.applyAdjoint(sino, volume);
    1042          20 :                         for (real_t i = 0.5; i < volSize; i += 1) {
    1043          80 :                             for (real_t j = 0.5; j < volSize; j += 1) {
    1044          64 :                                 const real_t angle = std::abs(
    1045          64 :                                     std::atan((sqrt3r * 40 + 2 - i) / (42 - j)) - pi_t / 3);
    1046          64 :                                 const real_t len = volSize * 21 * std::tan(angle);
    1047          64 :                                 if (len < 1) {
    1048          32 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j),
    1049             :                                                Approx(1 - len).epsilon(0.005));
    1050             :                                 } else {
    1051          32 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
    1052             :                                 }
    1053             :                             }
    1054             :                         }
    1055             :                     }
    1056             :                 }
    1057             :             }
    1058             : 
    1059          20 :             WHEN("Ray enters through the top border")
    1060             :             {
    1061             :                 // In this case the ray exits through a border along the main ray direction, but
    1062             :                 // enters through a border not along the main direction First pixel should be
    1063             :                 // weighted differently
    1064           4 :                 geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
    1065           4 :                                   PrincipalPointOffset{0}, RotationOffset2D{0, std::sqrt(3.f)});
    1066           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
    1067           8 :                 DataContainer<data_t> sino(sinoDescriptor);
    1068           4 :                 sino = 0;
    1069             : 
    1070           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
    1071           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
    1072             : 
    1073           8 :                 THEN("Ray intersects the correct pixels")
    1074             :                 {
    1075           4 :                     volume = 1;
    1076           4 :                     volume(0, 2) = 0;
    1077           4 :                     volume(0, 3) = 0;
    1078           4 :                     volume(1, 2) = 0;
    1079           4 :                     volume(1, 3) = 0;
    1080           4 :                     volume(2, 3) = 0;
    1081             : 
    1082           4 :                     fast.apply(volume, sino);
    1083           8 :                     DataContainer<data_t> zero(sinoDescriptor);
    1084           4 :                     zero = 0;
    1085           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
    1086             : 
    1087           4 :                     slow.apply(volume, sino);
    1088           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
    1089             : 
    1090           8 :                     AND_THEN("The correct weighting is applied")
    1091             :                     {
    1092           4 :                         volume(2, 3) = 1;
    1093           4 :                         volume(1, 2) = 1;
    1094           4 :                         volume(1, 3) = 1;
    1095             : 
    1096           4 :                         fast.apply(volume, sino);
    1097           4 :                         REQUIRE_EQ(sino[0], Approx((4 - 2 * sqrt3d) + (2 / sqrt3d)));
    1098             : 
    1099           4 :                         slow.apply(volume, sino);
    1100           4 :                         REQUIRE_EQ(sino[0], Approx((4 - 2 * sqrt3d) + (2 / sqrt3d)));
    1101             : 
    1102           4 :                         sino[0] = 1;
    1103             : 
    1104           8 :                         Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
    1105             : 
    1106           8 :                         slowExpected << 0, 0, 0, 0, 0, 0, 0, 0,
    1107           4 :                             (2 / sqrt3d) * (1 + halfd - sqrt3d / 2),
    1108           4 :                             (2 / sqrt3d) * (1 + halfd - 5 * sqrt3d / 6), 0, 0,
    1109           4 :                             (2 / sqrt3d) * (sqrt3d / 2 - halfd),
    1110             :                             (4 - 2 * sqrt3d) * (2 - sqrt3d)
    1111           4 :                                 + (2 / sqrt3d) * (5 * sqrt3d / 6 - halfd),
    1112           8 :                             (4 - 2 * sqrt3d) * (sqrt3d - 1), 0;
    1113             : 
    1114           4 :                         slow.applyAdjoint(sino, volume);
    1115           4 :                         REQUIRE_UNARY(isCwiseApprox(
    1116             :                             volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
    1117             : 
    1118           4 :                         fast.applyAdjoint(sino, volume);
    1119          20 :                         for (real_t i = 0.5; i < volSize; i += 1) {
    1120          80 :                             for (real_t j = 0.5; j < volSize; j += 1) {
    1121             :                                 const real_t angle =
    1122          64 :                                     std::abs(std::atan((sqrt3r * 40 + 2 - i) / (42 + sqrt3r - j))
    1123             :                                              - pi_t / 3);
    1124          64 :                                 const real_t len = volSize * 21 * std::tan(angle);
    1125          64 :                                 if (len < 1) {
    1126          20 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j),
    1127             :                                                Approx(1 - len).epsilon(0.01));
    1128             :                                 } else {
    1129          44 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
    1130             :                                 }
    1131             :                             }
    1132             :                         }
    1133             :                     }
    1134             :                 }
    1135             :             }
    1136             : 
    1137          20 :             WHEN("Ray exits through the bottom border")
    1138             :             {
    1139             :                 // In this case the ray enters through a border along the main ray direction, but
    1140             :                 // exits through a border not along the main direction Last pixel should be weighted
    1141             :                 // differently
    1142           4 :                 geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
    1143           4 :                                   PrincipalPointOffset{0}, RotationOffset2D{0, -std::sqrt(3.f)});
    1144             : 
    1145           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
    1146           8 :                 DataContainer<data_t> sino(sinoDescriptor);
    1147           4 :                 sino = 0;
    1148             : 
    1149           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
    1150           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
    1151             : 
    1152           8 :                 THEN("Ray intersects the correct pixels")
    1153             :                 {
    1154           4 :                     volume = 1;
    1155           4 :                     volume(1, 0) = 0;
    1156           4 :                     volume(2, 0) = 0;
    1157           4 :                     volume(3, 0) = 0;
    1158           4 :                     volume(2, 1) = 0;
    1159           4 :                     volume(3, 1) = 0;
    1160             : 
    1161           4 :                     fast.apply(volume, sino);
    1162           8 :                     DataContainer<data_t> zero(sinoDescriptor);
    1163           4 :                     zero = 0;
    1164           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
    1165             : 
    1166           4 :                     slow.apply(volume, sino);
    1167           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
    1168             : 
    1169           8 :                     AND_THEN("The correct weighting is applied")
    1170             :                     {
    1171           4 :                         volume(2, 0) = 1;
    1172           4 :                         volume(3, 1) = 1;
    1173             : 
    1174           4 :                         fast.apply(volume, sino);
    1175           4 :                         REQUIRE_EQ(sino[0], Approx((sqrt3d - 1) + (5.0 / 3.0 - 1 / sqrt3d)
    1176             :                                                    + (4 - 2 * sqrt3d) * (2 - sqrt3d))
    1177             :                                                 .epsilon(0.005));
    1178             : 
    1179           4 :                         slow.apply(volume, sino);
    1180           4 :                         REQUIRE_EQ(sino[0], Approx((sqrt3d - 1) + (5.0 / 3.0 - 1 / sqrt3d)
    1181             :                                                    + (4 - 2 * sqrt3d) * (2 - sqrt3d))
    1182             :                                                 .epsilon(0.005));
    1183             : 
    1184           4 :                         sino[0] = 1;
    1185             : 
    1186           8 :                         Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
    1187             : 
    1188           8 :                         slowExpected << 0, (sqrt3d - 1) * (4 - 2 * sqrt3d),
    1189           4 :                             (5 * thirdd - 1 / sqrt3d) + (4 - 2 * sqrt3d) * (2 - sqrt3d),
    1190           4 :                             1 - 1 / sqrt3d, 0, 0, sqrt3d - 5 * thirdd, sqrt3d - 1, 0, 0, 0, 0, 0, 0,
    1191           8 :                             0, 0;
    1192             : 
    1193           4 :                         slow.applyAdjoint(sino, volume);
    1194           4 :                         REQUIRE_UNARY(isCwiseApprox(
    1195             :                             volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
    1196             : 
    1197           4 :                         fast.applyAdjoint(sino, volume);
    1198             : 
    1199          20 :                         for (real_t i = 0.5; i < volSize; i += 1) {
    1200          80 :                             for (real_t j = 0.5; j < volSize; j += 1) {
    1201             :                                 const real_t angle =
    1202          64 :                                     std::abs(std::atan((sqrt3r * 40 + 2 - i) / (42 - sqrt3r - j))
    1203             :                                              - pi_t / 3);
    1204          64 :                                 const real_t len = 84 * std::tan(angle);
    1205             : 
    1206          64 :                                 if (len < 1) {
    1207          24 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j),
    1208             :                                                Approx(1 - len).epsilon(0.002));
    1209             :                                 } else {
    1210          40 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
    1211             :                                 }
    1212             :                             }
    1213             :                         }
    1214             :                     }
    1215             :                 }
    1216             :             }
    1217             : 
    1218          20 :             WHEN("Ray only intersects a single pixel")
    1219             :             {
    1220             :                 // This is a special case that is handled separately in both forward and
    1221             :                 // backprojection
    1222           4 :                 geom.emplace_back(stc, ctr, angle, std::move(volData), std::move(sinoData),
    1223             :                                   PrincipalPointOffset{0},
    1224           4 :                                   RotationOffset2D{0, -2 - std::sqrt(3.f) / 2});
    1225             : 
    1226           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
    1227           8 :                 DataContainer<data_t> sino(sinoDescriptor);
    1228           4 :                 sino = 0;
    1229             : 
    1230           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
    1231           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
    1232             : 
    1233           8 :                 THEN("Ray intersects the correct pixels")
    1234             :                 {
    1235           4 :                     volume = 1;
    1236           4 :                     volume(3, 0) = 0;
    1237             : 
    1238           4 :                     fast.apply(volume, sino);
    1239           8 :                     DataContainer<data_t> zero(sinoDescriptor);
    1240           4 :                     zero = 0;
    1241           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
    1242             : 
    1243           4 :                     slow.apply(volume, sino);
    1244           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
    1245             : 
    1246           8 :                     AND_THEN("The correct weighting is applied")
    1247             :                     {
    1248           4 :                         volume(3, 0) = 1;
    1249             : 
    1250           4 :                         fast.apply(volume, sino);
    1251           4 :                         REQUIRE_EQ(sino[0], Approx(1 / sqrt3d).epsilon(0.005));
    1252             : 
    1253           4 :                         slow.apply(volume, sino);
    1254           4 :                         REQUIRE_EQ(sino[0], Approx(1 / sqrt3d).epsilon(0.005));
    1255             : 
    1256           4 :                         sino[0] = 1;
    1257             : 
    1258           8 :                         Eigen::Matrix<data_t, Eigen::Dynamic, 1> slowExpected(volSize * volSize);
    1259           4 :                         slowExpected << 0, 0, 0, 1 / sqrt3d, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
    1260             : 
    1261           4 :                         slow.applyAdjoint(sino, volume);
    1262           4 :                         REQUIRE_UNARY(isCwiseApprox(
    1263             :                             volume, DataContainer<data_t>(volumeDescriptor, slowExpected)));
    1264             : 
    1265           4 :                         fast.applyAdjoint(sino, volume);
    1266             : 
    1267          20 :                         for (real_t i = 0.5; i < volSize; i += 1) {
    1268          80 :                             for (real_t j = 0.5; j < volSize; j += 1) {
    1269          64 :                                 const real_t angle = std::abs(
    1270          64 :                                     std::atan((sqrt3r * 40 + 2 - i) / (40 - sqrt3r / 2 - j))
    1271             :                                     - pi_t / 3);
    1272          64 :                                 const real_t len = 84 * std::tan(angle);
    1273             : 
    1274          64 :                                 if (len < 1) {
    1275           8 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j),
    1276             :                                                Approx(1 - len).epsilon(0.002));
    1277             :                                 } else {
    1278          56 :                                     REQUIRE_EQ(volume((index_t) i, (index_t) j), Approx(0));
    1279             :                                 }
    1280             :                             }
    1281             :                         }
    1282             :                     }
    1283             :                 }
    1284             :             }
    1285             :         }
    1286             :     }
    1287          80 : }
    1288             : 
    1289          34 : TEST_CASE_TEMPLATE("JosephsMethodCUDA: 2D setup with a multiple rays", TestType,
    1290             :                    JosephsMethodCUDA<float>, JosephsMethodCUDA<double>)
    1291             : {
    1292             :     // Turn logger of
    1293           4 :     Logger::setLevel(Logger::LogLevel::OFF);
    1294             : 
    1295             :     using data_t = decltype(return_data_t(std::declval<TestType>()));
    1296             : 
    1297           8 :     GIVEN("a 5x5 Volume and detector of size 1, with 4 angles")
    1298             :     {
    1299             :         // Domain setup
    1300           4 :         const index_t volSize = 5;
    1301             : 
    1302           8 :         IndexVector_t volumeDims(2);
    1303           4 :         volumeDims << volSize, volSize;
    1304             : 
    1305           8 :         VolumeDescriptor volumeDescriptor(volumeDims);
    1306           8 :         DataContainer<data_t> volume(volumeDescriptor);
    1307             : 
    1308             :         // range setup
    1309           4 :         const index_t detectorSize = 1;
    1310           4 :         const index_t numImgs = 4;
    1311             : 
    1312           8 :         IndexVector_t sinoDims(2);
    1313           4 :         sinoDims << detectorSize, numImgs;
    1314             : 
    1315             :         // Setup geometry
    1316           4 :         auto stc = SourceToCenterOfRotation{20 * volSize};
    1317           4 :         auto ctr = CenterOfRotationToDetector{volSize};
    1318          12 :         auto volData = VolumeData2D{Size2D{volumeDims}};
    1319          12 :         auto sinoData = SinogramData2D{Size2D{sinoDims}};
    1320             : 
    1321           8 :         std::vector<Geometry> geom;
    1322             : 
    1323           8 :         WHEN("Both x- and y-axis-aligned rays are present")
    1324             :         {
    1325           4 :             geom.emplace_back(stc, ctr, Degree{0}, VolumeData2D{Size2D{volumeDims}},
    1326             :                               SinogramData2D{Size2D{sinoDims}});
    1327           4 :             geom.emplace_back(stc, ctr, Degree{90}, VolumeData2D{Size2D{volumeDims}},
    1328             :                               SinogramData2D{Size2D{sinoDims}});
    1329           4 :             geom.emplace_back(stc, ctr, Degree{180}, VolumeData2D{Size2D{volumeDims}},
    1330             :                               SinogramData2D{Size2D{sinoDims}});
    1331           4 :             geom.emplace_back(stc, ctr, Degree{270}, VolumeData2D{Size2D{volumeDims}},
    1332             :                               SinogramData2D{Size2D{sinoDims}});
    1333             : 
    1334           8 :             PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
    1335           8 :             DataContainer<data_t> sino(sinoDescriptor);
    1336           4 :             sino = 0;
    1337             : 
    1338           8 :             TestType fast(volumeDescriptor, sinoDescriptor);
    1339           8 :             TestType slow(volumeDescriptor, sinoDescriptor, false);
    1340             : 
    1341           8 :             THEN("Values are accumulated correctly along each ray's path")
    1342             :             {
    1343           4 :                 volume = 0;
    1344             : 
    1345             :                 // set only values along the rays' path to one to make sure interpolation is done
    1346             :                 // correctly
    1347          24 :                 for (index_t i = 0; i < volSize; i++) {
    1348          20 :                     volume(i, volSize / 2) = 1;
    1349          20 :                     volume(volSize / 2, i) = 1;
    1350             :                 }
    1351             : 
    1352           4 :                 slow.apply(volume, sino);
    1353          20 :                 for (index_t i = 0; i < numImgs; i++)
    1354          16 :                     REQUIRE_EQ(sino[i], Approx(5.0));
    1355             : 
    1356           4 :                 fast.apply(volume, sino);
    1357          20 :                 for (index_t i = 0; i < numImgs; i++)
    1358          16 :                     REQUIRE_EQ(sino[i], Approx(5.0));
    1359             : 
    1360           8 :                 AND_THEN("Both fast and slow backprojections yield the exact adjoint")
    1361             :                 {
    1362           4 :                     Eigen::Matrix<data_t, Eigen::Dynamic, 1> cmp(volSize * volSize);
    1363             : 
    1364           8 :                     cmp << 0, 0, 10, 0, 0, 0, 0, 10, 0, 0, 10, 10, 20, 10, 10, 0, 0, 10, 0, 0, 0, 0,
    1365           8 :                         10, 0, 0;
    1366             : 
    1367           4 :                     slow.applyAdjoint(sino, volume);
    1368           4 :                     REQUIRE_UNARY(
    1369             :                         isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor, cmp)));
    1370             : 
    1371           4 :                     fast.applyAdjoint(sino, volume);
    1372           4 :                     REQUIRE_UNARY(
    1373             :                         isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor, cmp)));
    1374             :                 }
    1375             :             }
    1376             :         }
    1377             :     }
    1378           4 : }
    1379          74 : TEST_CASE_TEMPLATE("JosephsMethodCUDA: 3D setup with a single ray", TestType,
    1380             :                    JosephsMethodCUDA<float>, JosephsMethodCUDA<double>)
    1381             : {
    1382             :     // Turn logger of
    1383          44 :     Logger::setLevel(Logger::LogLevel::OFF);
    1384             : 
    1385             :     using data_t = decltype(return_data_t(std::declval<TestType>()));
    1386             : 
    1387          80 :     GIVEN("A 3x3x3 Volume")
    1388             :     {
    1389             :         // Domain setup
    1390          36 :         const index_t volSize = 3;
    1391             : 
    1392          72 :         IndexVector_t volumeDims(3);
    1393          36 :         volumeDims << volSize, volSize, volSize;
    1394             : 
    1395          72 :         VolumeDescriptor volumeDescriptor(volumeDims);
    1396          72 :         DataContainer<data_t> volume(volumeDescriptor);
    1397             : 
    1398             :         // range setup
    1399          36 :         const index_t detectorSize = 1;
    1400          36 :         const index_t numImgs = 1;
    1401             : 
    1402          72 :         IndexVector_t sinoDims(3);
    1403          36 :         sinoDims << detectorSize, detectorSize, numImgs;
    1404             : 
    1405             :         // Setup geometry
    1406          36 :         const auto stc = SourceToCenterOfRotation{20 * volSize};
    1407          36 :         const auto ctr = CenterOfRotationToDetector{volSize};
    1408         108 :         auto volData = VolumeData3D{Size3D{volumeDims}};
    1409         108 :         auto sinoData = SinogramData3D{Size3D{sinoDims}};
    1410             : 
    1411          72 :         std::vector<Geometry> geom;
    1412             : 
    1413          44 :         GIVEN("Single ray along z-axis")
    1414             :         {
    1415             : 
    1416           8 :             geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
    1417             :                               RotationAngles3D{Gamma{0}});
    1418             : 
    1419          16 :             PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
    1420          16 :             DataContainer<data_t> sino(sinoDescriptor);
    1421           8 :             sino = 0;
    1422             : 
    1423          16 :             TestType fast(volumeDescriptor, sinoDescriptor);
    1424          16 :             TestType slow(volumeDescriptor, sinoDescriptor, false);
    1425             : 
    1426          12 :             WHEN("Sinogram conatainer is not zero initialized and we project through an empty "
    1427             :                  "volume")
    1428             :             {
    1429           4 :                 volume = 0;
    1430           4 :                 sino = 1;
    1431             : 
    1432           8 :                 THEN("Result is zero")
    1433             :                 {
    1434           4 :                     fast.apply(volume, sino);
    1435           4 :                     DataContainer<data_t> zero(sinoDescriptor);
    1436           4 :                     zero = 0;
    1437           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
    1438             : 
    1439           4 :                     sino = 1;
    1440           4 :                     slow.apply(volume, sino);
    1441           4 :                     REQUIRE_UNARY(isCwiseApprox(sino, zero));
    1442             :                 }
    1443             :             }
    1444             : 
    1445          12 :             WHEN("Volume container is not zero initialized and we backproject from an empty "
    1446             :                  "sinogram")
    1447             :             {
    1448           4 :                 sino = 0;
    1449           4 :                 volume = 1;
    1450             : 
    1451           8 :                 THEN("Result is zero")
    1452             :                 {
    1453           4 :                     fast.applyAdjoint(sino, volume);
    1454           4 :                     DataContainer<data_t> zero(volumeDescriptor);
    1455           4 :                     zero = 0;
    1456           4 :                     REQUIRE_UNARY(isCwiseApprox(volume, zero));
    1457             : 
    1458           4 :                     volume = 1;
    1459           4 :                     slow.applyAdjoint(sino, volume);
    1460           4 :                     REQUIRE_UNARY(isCwiseApprox(volume, zero));
    1461             :                 }
    1462             :             }
    1463             :         }
    1464             : 
    1465          36 :         const index_t numCases = 6;
    1466          36 :         std::array<real_t, numCases> beta = {0.0, 0.0, 0.0, 0.0, pi_t / 2, 3 * pi_t / 2};
    1467          36 :         std::array<real_t, numCases> gamma = {0.0,          pi_t,     pi_t / 2,
    1468             :                                               3 * pi_t / 2, pi_t / 2, 3 * pi_t / 2};
    1469          72 :         std::array<std::string, numCases> al = {"z", "-z", "x", "-x", "y", "-y"};
    1470             : 
    1471         504 :         Eigen::Matrix<data_t, Eigen::Dynamic, 1> backProjections[numCases];
    1472         252 :         for (auto& backPr : backProjections)
    1473         216 :             backPr.resize(volSize * volSize * volSize);
    1474             : 
    1475             :         // clang-format off
    1476          72 :         backProjections[0] << 0, 0, 0, 0, 1, 0, 0, 0, 0,
    1477          36 :                               0, 0, 0, 0, 1, 0, 0, 0, 0,
    1478          72 :                               0, 0, 0, 0, 1, 0, 0, 0, 0;
    1479             : 
    1480          72 :         backProjections[1] << 0, 0, 0, 0, 0, 0, 0, 0, 0,
    1481          36 :                               0, 0, 0, 1, 1, 1, 0, 0, 0,
    1482          72 :                               0, 0, 0, 0, 0, 0, 0, 0, 0;
    1483             : 
    1484          72 :         backProjections[2] << 0, 0, 0, 0, 0, 0, 0, 0, 0,
    1485          36 :                               0, 1, 0, 0, 1, 0, 0, 1, 0,
    1486          72 :                               0, 0, 0, 0, 0, 0, 0, 0, 0;
    1487             :         // clang-format on
    1488             : 
    1489         252 :         for (std::size_t i = 0; i < numCases; i++) {
    1490         220 :             WHEN("An axis-aligned ray passes through the center of a pixel")
    1491             :             {
    1492           8 :                 INFO("A ", al[i], "-axis-aligned ray passes through the center of a pixel");
    1493           4 :                 geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
    1494             :                                   RotationAngles3D{Gamma{gamma[i]}, Beta{beta[i]}});
    1495             : 
    1496           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
    1497           8 :                 DataContainer<data_t> sino(sinoDescriptor);
    1498           4 :                 sino = 0;
    1499             : 
    1500           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
    1501           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
    1502             : 
    1503           8 :                 THEN("The result of projecting through a voxel is exactly the voxel value")
    1504             :                 {
    1505          16 :                     for (index_t j = 0; j < volSize; j++) {
    1506          12 :                         volume = 0;
    1507          12 :                         if (i / 2 == 0)
    1508          12 :                             volume(volSize / 2, volSize / 2, j) = 1;
    1509           0 :                         else if (i / 2 == 1)
    1510           0 :                             volume(j, volSize / 2, volSize / 2) = 1;
    1511           0 :                         else if (i / 2 == 2)
    1512           0 :                             volume(volSize / 2, j, volSize / 2) = 1;
    1513             : 
    1514          12 :                         fast.apply(volume, sino);
    1515             :                         // Using doubles significantly increases interpolation accuracy
    1516             :                         // For example: when using floats, points very near the border (less than
    1517             :                         // ~1/500th of a pixel away from the border) are rounded to actually lie on
    1518             :                         // the border. This then yields more accurate results when using floats in
    1519             :                         // some of the axis-aligned test cases, despite the lower interpolation
    1520             :                         // accuracy.
    1521             :                         //  => different requirements for floats and doubles, looser requirements
    1522             :                         // for doubles
    1523          12 :                         REQUIRE_EQ(sino[0], Approx(1.0));
    1524             : 
    1525          12 :                         slow.apply(volume, sino);
    1526          12 :                         REQUIRE_EQ(sino[0], Approx(1.0));
    1527             :                     }
    1528             : 
    1529           8 :                     AND_THEN("The backprojection sets the values of all hit voxels to the detector "
    1530             :                              "value")
    1531             :                     {
    1532           4 :                         fast.applyAdjoint(sino, volume);
    1533           4 :                         REQUIRE_UNARY(
    1534             :                             isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor,
    1535             :                                                                         backProjections[i / 2])));
    1536             : 
    1537           4 :                         slow.applyAdjoint(sino, volume);
    1538           4 :                         REQUIRE_UNARY(
    1539             :                             isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor,
    1540             :                                                                         backProjections[i / 2])));
    1541             :                     }
    1542             :                 }
    1543             :             }
    1544             :         }
    1545             : 
    1546          36 :         std::array<real_t, numCases> offsetx = {-0.25, -0.25, 0.0, 0.0, 0.0, 0.0};
    1547          36 :         std::array<real_t, numCases> offsety = {0.0, 0.0, -0.25, -0.25, 0.0, 0.0};
    1548          36 :         std::array<real_t, numCases> offsetz = {0.0, 0.0, 0.0, 0.0, -0.25, -0.25};
    1549             : 
    1550             :         // clang-format off
    1551          72 :         backProjections[0] << 0, 0, 0, 0.25, 0.75, 0, 0, 0, 0,
    1552          36 :                               0, 0, 0, 0.25, 0.75, 0, 0, 0, 0,
    1553          72 :                               0, 0, 0, 0.25, 0.75, 0, 0, 0, 0;
    1554             : 
    1555          72 :         backProjections[1] <<  0   , 0   , 0   , 0   , 0   , 0   , 0, 0, 0,
    1556          36 :                                0.25, 0.25, 0.25, 0.75, 0.75, 0.75, 0, 0, 0,
    1557          72 :                                0   , 0   , 0   , 0   , 0   , 0   , 0, 0, 0;
    1558             : 
    1559          72 :         backProjections[2] << 0, 0.25, 0, 0, 0.25, 0, 0, 0.25, 0,
    1560          36 :                               0, 0.75, 0, 0, 0.75, 0, 0, 0.75, 0,
    1561          72 :                               0, 0   , 0, 0, 0   , 0, 0, 0   , 0;
    1562             :         // clang-format on
    1563             : 
    1564         252 :         for (std::size_t i = 0; i < numCases; i++) {
    1565         220 :             WHEN("An axis-aligned ray not passing through the center of a pixel")
    1566             :             {
    1567           8 :                 INFO("A ", al[i], "-axis-aligned ray not passing through the center of a pixel");
    1568             :                 // x-ray source must be very far from the volume center to make testing of the fast
    1569             :                 // backprojection simpler
    1570           8 :                 geom.emplace_back(SourceToCenterOfRotation{volSize * 2000}, ctr, std::move(volData),
    1571           4 :                                   std::move(sinoData),
    1572             :                                   RotationAngles3D{Gamma{gamma[i]}, Beta{beta[i]}},
    1573             :                                   PrincipalPointOffset2D{0, 0},
    1574             :                                   RotationOffset3D{offsetx[i], offsety[i], offsetz[i]});
    1575             : 
    1576           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
    1577           8 :                 DataContainer<data_t> sino(sinoDescriptor);
    1578           4 :                 sino = 0;
    1579             : 
    1580           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
    1581           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
    1582             : 
    1583           8 :                 THEN("The result of projecting through a voxel is the interpolated value between "
    1584             :                      "the four voxels nearest to the ray")
    1585             :                 {
    1586          16 :                     for (index_t j = 0; j < volSize; j++) {
    1587          12 :                         volume = 0;
    1588          12 :                         if (i / 2 == 0)
    1589          12 :                             volume(volSize / 2, volSize / 2, j) = 1;
    1590           0 :                         else if (i / 2 == 1)
    1591           0 :                             volume(j, volSize / 2, volSize / 2) = 1;
    1592           0 :                         else if (i / 2 == 2)
    1593           0 :                             volume(volSize / 2, j, volSize / 2) = 1;
    1594             : 
    1595          12 :                         fast.apply(volume, sino);
    1596             : 
    1597             :                         // Using doubles significantly increases interpolation accuracy
    1598             :                         //  For example: when using floats, points very near the border (less than
    1599             :                         // ~1/500th of a pixel away from the border) are rounded to actually lie on
    1600             :                         // the border. This then yields more accurate results when using floats in
    1601             :                         // some of the axis-aligned test cases, despite the lower interpolation
    1602             :                         // accuracy.
    1603             :                         //
    1604             :                         //  => different requirements for floats and doubles, looser requirements
    1605             :                         // for doubles
    1606             : 
    1607          12 :                         REQUIRE_EQ(sino[0], Approx(0.75));
    1608             : 
    1609          12 :                         slow.apply(volume, sino);
    1610          12 :                         REQUIRE_EQ(sino[0], Approx(0.75));
    1611             :                     }
    1612             : 
    1613           8 :                     AND_THEN("The slow backprojection yields the exact adjoint, the fast "
    1614             :                              "backprojection is also exact for a very distant x-ray source")
    1615             :                     {
    1616           4 :                         sino[0] = 1;
    1617           4 :                         fast.applyAdjoint(sino, volume);
    1618             :                         // Using doubles significantly increases interpolation accuracy
    1619             :                         //  For example: when using floats, points very near the border (less than
    1620             :                         // ~1/500th of a pixel away from the border) are rounded to actually lie on
    1621             :                         // the border. This then yields more accurate results when using floats in
    1622             :                         // some of the axis-aligned test cases, despite the lower interpolation
    1623             :                         // accuracy.
    1624             :                         //
    1625             :                         //  => different requirements for floats and doubles, looser requirements
    1626             :                         // for doubles
    1627             :                         if constexpr (std::is_same_v<data_t, float>)
    1628           2 :                             REQUIRE_UNARY(isCwiseApprox(
    1629             :                                 volume,
    1630             :                                 DataContainer<data_t>(volumeDescriptor, backProjections[i / 2])));
    1631             :                         else
    1632           2 :                             REQUIRE_UNARY(isApprox(
    1633             :                                 volume,
    1634             :                                 DataContainer<data_t>(volumeDescriptor, backProjections[i / 2]),
    1635             :                                 static_cast<real_t>(0.005)));
    1636             : 
    1637           4 :                         slow.applyAdjoint(sino, volume);
    1638           4 :                         REQUIRE_UNARY(
    1639             :                             isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor,
    1640             :                                                                         backProjections[i / 2])));
    1641             :                     }
    1642             :                 }
    1643             :             }
    1644             :         }
    1645             : 
    1646          36 :         offsetx[0] = -volSize / 2.0;
    1647          36 :         offsetx[1] = volSize / 2.0;
    1648          36 :         offsetx[2] = 0.0;
    1649          36 :         offsetx[3] = 0.0;
    1650          36 :         offsetx[4] = volSize / 2.0;
    1651          36 :         offsetx[5] = -volSize / 2.0;
    1652             : 
    1653          36 :         offsety[0] = 0.0;
    1654          36 :         offsety[1] = 0.0;
    1655          36 :         offsety[2] = -volSize / 2.0;
    1656          36 :         offsety[3] = volSize / 2.0;
    1657          36 :         offsety[4] = volSize / 2.0;
    1658          36 :         offsety[5] = -volSize / 2.0;
    1659             : 
    1660             :         // clang-format off
    1661          72 :         backProjections[0] << 0, 0, 0, 1, 0, 0, 0, 0, 0,
    1662          36 :                               0, 0, 0, 1, 0, 0, 0, 0, 0,
    1663          72 :                               0, 0, 0, 1, 0, 0, 0, 0, 0;
    1664             : 
    1665          72 :         backProjections[1] << 0, 0, 0, 0, 0, 1, 0, 0, 0,
    1666          36 :                               0, 0, 0, 0, 0, 1, 0, 0, 0,
    1667          72 :                               0, 0, 0, 0, 0, 1, 0, 0, 0;
    1668             : 
    1669          72 :         backProjections[2] << 0, 1, 0, 0, 0, 0, 0, 0, 0,
    1670          36 :                               0, 1, 0, 0, 0, 0, 0, 0, 0,
    1671          72 :                               0, 1, 0, 0, 0, 0, 0, 0, 0;
    1672             : 
    1673          72 :         backProjections[3] << 0, 0, 0, 0, 0, 0, 0, 1, 0,
    1674          36 :                               0, 0, 0, 0, 0, 0, 0, 1, 0,
    1675          72 :                               0, 0, 0, 0, 0, 0, 0, 1, 0;
    1676             : 
    1677          72 :         backProjections[4] << 0, 0, 0, 0, 0, 0, 0, 0, 1,
    1678          36 :                               0, 0, 0, 0, 0, 0, 0, 0, 1,
    1679          72 :                               0, 0, 0, 0, 0, 0, 0, 0, 1;
    1680             : 
    1681          72 :         backProjections[5] << 1, 0, 0, 0, 0, 0, 0, 0, 0,
    1682          36 :                               1, 0, 0, 0, 0, 0, 0, 0, 0,
    1683          72 :                               1, 0, 0, 0, 0, 0, 0, 0, 0;
    1684             :         // clang-format on
    1685             : 
    1686          36 :         al[0] = "left border";
    1687          36 :         al[1] = "right border";
    1688          36 :         al[2] = "top border";
    1689          36 :         al[3] = "bottom border";
    1690          36 :         al[4] = "top right edge";
    1691          36 :         al[5] = "bottom left edge";
    1692             : 
    1693         252 :         for (std::size_t i = 0; i < numCases; i++) {
    1694         220 :             WHEN("A z-axis-aligned ray runs along the corners and edges of the volume")
    1695             :             {
    1696           8 :                 INFO("A z-axis-aligned ray runs along the ", al[i], " of the volume");
    1697             :                 // x-ray source must be very far from the volume center to make testing of the fast
    1698             :                 // backprojection simpler
    1699           8 :                 geom.emplace_back(SourceToCenterOfRotation{volSize * 2000}, ctr, std::move(volData),
    1700           4 :                                   std::move(sinoData), RotationAngles3D{Gamma{0}},
    1701             :                                   PrincipalPointOffset2D{0, 0},
    1702             :                                   RotationOffset3D{offsetx[i], offsety[i], 0});
    1703           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
    1704           8 :                 DataContainer<data_t> sino(sinoDescriptor);
    1705           4 :                 sino = 0;
    1706             : 
    1707           8 :                 TestType fast(volumeDescriptor, sinoDescriptor);
    1708           8 :                 TestType slow(volumeDescriptor, sinoDescriptor, false);
    1709             : 
    1710           8 :                 THEN("The result of projecting through a voxel is exactly the voxel's value (we "
    1711             :                      "mirror values at the border for the purpose of interpolation)")
    1712             :                 {
    1713          16 :                     for (index_t j = 0; j < volSize; j++) {
    1714          12 :                         volume = 0;
    1715          12 :                         switch (i) {
    1716          12 :                             case 0:
    1717          12 :                                 volume(0, volSize / 2, j) = 1;
    1718          12 :                                 break;
    1719           0 :                             case 1:
    1720           0 :                                 volume(volSize - 1, volSize / 2, j) = 1;
    1721           0 :                                 break;
    1722           0 :                             case 2:
    1723           0 :                                 volume(volSize / 2, 0, j) = 1;
    1724           0 :                                 break;
    1725           0 :                             case 3:
    1726           0 :                                 volume(volSize / 2, volSize - 1, j) = 1;
    1727           0 :                                 break;
    1728           0 :                             case 4:
    1729           0 :                                 volume(volSize - 1, volSize - 1, j) = 1;
    1730           0 :                                 break;
    1731           0 :                             case 5:
    1732           0 :                                 volume(0, 0, j) = 1;
    1733           0 :                                 break;
    1734           0 :                             default:
    1735           0 :                                 break;
    1736             :                         }
    1737             : 
    1738          12 :                         fast.apply(volume, sino);
    1739          12 :                         REQUIRE_EQ(sino[0], Approx(1));
    1740             : 
    1741          12 :                         slow.apply(volume, sino);
    1742          12 :                         REQUIRE_EQ(sino[0], Approx(1));
    1743             :                     }
    1744             : 
    1745           8 :                     AND_THEN("The slow backprojection yields the exact adjoint, the fast "
    1746             :                              "backprojection is also exact for a very distant x-ray source")
    1747             :                     {
    1748           4 :                         sino[0] = 1;
    1749           4 :                         fast.applyAdjoint(sino, volume);
    1750             : 
    1751             :                         // Using doubles significantly increases interpolation accuracy
    1752             :                         //  For example: when using floats, points very near the border (less than
    1753             :                         // ~1/500th of a pixel away from the border) are rounded to actually lie on
    1754             :                         // the border. This then yields more accurate results when using floats in
    1755             :                         // some of the axis-aligned test cases, despite the lower interpolation
    1756             :                         // accuracy.
    1757             :                         //
    1758             :                         //  => different requirements for floats and doubles, looser requirements
    1759             :                         // for doubles
    1760             :                         if constexpr (std::is_same_v<data_t, float>)
    1761           2 :                             REQUIRE_UNARY(isCwiseApprox(
    1762             :                                 volume, DataContainer<data_t>(
    1763             :                                             volumeDescriptor,
    1764             :                                             (backProjections[i] / (i > 3 ? 4 : 2)).eval())));
    1765             :                         else
    1766           2 :                             REQUIRE_UNARY(
    1767             :                                 isApprox(volume,
    1768             :                                          DataContainer<data_t>(
    1769             :                                              volumeDescriptor,
    1770             :                                              (backProjections[i] / (i > 3 ? 4 : 2)).eval()),
    1771             :                                          static_cast<real_t>(0.005)));
    1772             : 
    1773           4 :                         slow.applyAdjoint(sino, volume);
    1774           4 :                         REQUIRE_UNARY(isCwiseApprox(
    1775             :                             volume, DataContainer<data_t>(volumeDescriptor, backProjections[i])));
    1776             :                     }
    1777             :                 }
    1778             :             }
    1779             :         }
    1780             : 
    1781          36 :         const data_t sqrt3d = std::sqrt(static_cast<data_t>(3));
    1782          36 :         const data_t thirdd = static_cast<data_t>(1.0 / 3);
    1783             : 
    1784          72 :         Eigen::Matrix<data_t, Eigen::Dynamic, 1> backProj(volSize * volSize * volSize);
    1785             : 
    1786          52 :         GIVEN("A 30 degree rotate around z axis")
    1787             :         {
    1788          16 :             auto gamma = Gamma{Degree{30}};
    1789             : 
    1790          20 :             WHEN("A ray goes through the center of the volume")
    1791             :             {
    1792             :                 // In this case the ray enters and exits the volume along the main direction
    1793           4 :                 geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
    1794             :                                   RotationAngles3D{gamma});
    1795             : 
    1796           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
    1797           8 :                 DataContainer<data_t> sino(sinoDescriptor);
    1798             : 
    1799           8 :                 TestType op(volumeDescriptor, sinoDescriptor, false);
    1800             : 
    1801           8 :                 THEN("The ray intersects the correct voxels")
    1802             :                 {
    1803           4 :                     volume = 1;
    1804           4 :                     volume(1, 1, 1) = 0;
    1805           4 :                     volume(2, 1, 0) = 0;
    1806           4 :                     volume(1, 1, 0) = 0;
    1807           4 :                     volume(0, 1, 2) = 0;
    1808           4 :                     volume(1, 1, 2) = 0;
    1809             : 
    1810           4 :                     op.apply(volume, sino);
    1811           4 :                     REQUIRE_EQ(sino[0], Approx(0).epsilon(1e-5));
    1812             : 
    1813           8 :                     AND_THEN("The correct weighting is applied")
    1814             :                     {
    1815           4 :                         volume(1, 1, 1) = 1;
    1816           4 :                         volume(2, 1, 0) = 3;
    1817           4 :                         volume(1, 1, 2) = 2;
    1818             : 
    1819           4 :                         op.apply(volume, sino);
    1820           4 :                         REQUIRE_EQ(sino[0], Approx(6 / sqrt3d + 2.0 / 3).epsilon(0.001));
    1821             : 
    1822           4 :                         sino[0] = 1;
    1823             :                         // clang-format off
    1824           8 :                             backProj << 0, 0, 0, 0         , 2 / sqrt3d - 2 * thirdd, 2 * thirdd, 0, 0, 0,
    1825           4 :                                         0, 0, 0, 0         , 2 / sqrt3d             , 0         , 0, 0, 0,
    1826           8 :                                         0, 0, 0, 2 * thirdd, 2 / sqrt3d - 2 * thirdd, 0         , 0, 0, 0;
    1827             :                         // clang-format on
    1828             : 
    1829           4 :                         op.applyAdjoint(sino, volume);
    1830           4 :                         REQUIRE_UNARY(isCwiseApprox(
    1831             :                             volume, DataContainer<data_t>(volumeDescriptor, backProj)));
    1832             :                     }
    1833             :                 }
    1834             :             }
    1835             : 
    1836          20 :             WHEN("A ray enters through the right border")
    1837             :             {
    1838             :                 // In this case the ray enters through a border orthogonal to a non-main direction
    1839           4 :                 geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
    1840             :                                   RotationAngles3D{gamma}, PrincipalPointOffset2D{0, 0},
    1841             :                                   RotationOffset3D{1, 0, 0});
    1842             : 
    1843           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
    1844           8 :                 DataContainer<data_t> sino(sinoDescriptor);
    1845             : 
    1846           8 :                 TestType op(volumeDescriptor, sinoDescriptor, false);
    1847             : 
    1848           8 :                 THEN("The ray intersects the correct voxels")
    1849             :                 {
    1850           4 :                     volume = 1;
    1851           4 :                     volume(2, 1, 1) = 0;
    1852           4 :                     volume(2, 1, 0) = 0;
    1853           4 :                     volume(2, 1, 2) = 0;
    1854           4 :                     volume(1, 1, 2) = 0;
    1855             : 
    1856           4 :                     op.apply(volume, sino);
    1857           4 :                     REQUIRE_EQ(sino[0], Approx(0).epsilon(1e-5));
    1858             : 
    1859           8 :                     AND_THEN("The correct weighting is applied")
    1860             :                     {
    1861           4 :                         volume(2, 1, 0) = 4;
    1862           4 :                         volume(1, 1, 2) = 3;
    1863           4 :                         volume(2, 1, 1) = 1;
    1864             : 
    1865           4 :                         op.apply(volume, sino);
    1866           4 :                         REQUIRE_EQ(sino[0], Approx((sqrt3d + 1) * (1 - 1 / sqrt3d) + 3 - sqrt3d / 2
    1867             :                                                    + 2 / sqrt3d)
    1868             :                                                 .epsilon(0.001));
    1869             : 
    1870           4 :                         sino[0] = 1;
    1871             :                         // clang-format off
    1872           8 :                         backProj << 0, 0, 0, 0, 0, ((sqrt3d + 1) / 4) * (1 - 1 / sqrt3d), 0, 0, 0,
    1873           4 :                             0, 0, 0, 0, 0, 2 / sqrt3d + 1 - sqrt3d / 2, 0, 0, 0, 0, 0, 0, 0,
    1874           8 :                             2 * thirdd, 2 / sqrt3d - 2 * thirdd, 0, 0, 0;
    1875             :                         // clang-format on
    1876             : 
    1877           4 :                         op.applyAdjoint(sino, volume);
    1878           4 :                         REQUIRE_UNARY(isCwiseApprox(
    1879             :                             volume, DataContainer<data_t>(volumeDescriptor, backProj)));
    1880             :                     }
    1881             :                 }
    1882             :             }
    1883             : 
    1884          20 :             WHEN("A ray exits through the left border")
    1885             :             {
    1886             :                 // In this case the ray exit through a border orthogonal to a non-main direction
    1887           4 :                 geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
    1888             :                                   RotationAngles3D{gamma}, PrincipalPointOffset2D{0, 0},
    1889             :                                   RotationOffset3D{-1, 0, 0});
    1890             : 
    1891           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
    1892           8 :                 DataContainer<data_t> sino(sinoDescriptor);
    1893             : 
    1894           8 :                 TestType op(volumeDescriptor, sinoDescriptor, false);
    1895             : 
    1896           8 :                 THEN("The ray intersects the correct voxels")
    1897             :                 {
    1898           4 :                     volume = 1;
    1899           4 :                     volume(0, 1, 0) = 0;
    1900           4 :                     volume(1, 1, 0) = 0;
    1901           4 :                     volume(0, 1, 1) = 0;
    1902           4 :                     volume(0, 1, 2) = 0;
    1903             : 
    1904           4 :                     op.apply(volume, sino);
    1905           4 :                     REQUIRE_EQ(sino[0], Approx(0).epsilon(1e-5));
    1906             : 
    1907           8 :                     AND_THEN("The correct weighting is applied")
    1908             :                     {
    1909           4 :                         volume(0, 1, 2) = 4;
    1910           4 :                         volume(1, 1, 0) = 3;
    1911           4 :                         volume(0, 1, 1) = 1;
    1912             : 
    1913           4 :                         op.apply(volume, sino);
    1914           4 :                         REQUIRE_EQ(sino[0], Approx((sqrt3d + 1) * (1 - 1 / sqrt3d) + 3 - sqrt3d / 2
    1915             :                                                    + 2 / sqrt3d)
    1916             :                                                 .epsilon(0.001));
    1917             : 
    1918           4 :                         sino[0] = 1;
    1919             :                         // clang-format off
    1920           8 :                         backProj << 0, 0, 0, 2 / sqrt3d - 2 * thirdd, 2 * thirdd, 0, 0, 0, 0, 0, 0,
    1921           4 :                             0, 2 / sqrt3d + 1 - sqrt3d / 2, 0, 0, 0, 0, 0, 0, 0, 0,
    1922           8 :                             ((sqrt3d + 1) / 4) * (1 - 1 / sqrt3d), 0, 0, 0, 0, 0;
    1923             :                         // clang-format on
    1924             : 
    1925           4 :                         op.applyAdjoint(sino, volume);
    1926           4 :                         REQUIRE_UNARY(isCwiseApprox(
    1927             :                             volume, DataContainer<data_t>(volumeDescriptor, backProj)));
    1928             :                     }
    1929             :                 }
    1930             :             }
    1931             : 
    1932          20 :             WHEN("A ray with an angle of 30 degrees only intersects a single voxel")
    1933             :             {
    1934             :                 // special case - no interior voxels, entry and exit voxels are the same
    1935           4 :                 geom.emplace_back(stc, ctr, std::move(volData), std::move(sinoData),
    1936             :                                   RotationAngles3D{gamma}, PrincipalPointOffset2D{0, 0},
    1937             :                                   RotationOffset3D{-2, 0, 0});
    1938             : 
    1939           8 :                 PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
    1940           8 :                 DataContainer<data_t> sino(sinoDescriptor);
    1941             : 
    1942           8 :                 TestType op(volumeDescriptor, sinoDescriptor, false);
    1943             : 
    1944           8 :                 THEN("The ray intersects the correct voxels")
    1945             :                 {
    1946           4 :                     volume = 1;
    1947           4 :                     volume(0, 1, 0) = 0;
    1948             : 
    1949           4 :                     op.apply(volume, sino);
    1950           4 :                     REQUIRE_EQ(sino[0], Approx(0).epsilon(1e-5));
    1951             : 
    1952           8 :                     AND_THEN("The correct weighting is applied")
    1953             :                     {
    1954           4 :                         volume(0, 1, 0) = 1;
    1955             : 
    1956           4 :                         op.apply(volume, sino);
    1957           4 :                         REQUIRE_EQ(sino[0], Approx(sqrt3d - 1));
    1958             : 
    1959           4 :                         sino[0] = 1;
    1960             :                         // clang-format off
    1961           8 :                         backProj << 0, 0, 0, sqrt3d - 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
    1962           8 :                             0, 0, 0, 0, 0, 0, 0, 0, 0;
    1963             :                         // clang-format off
    1964             : 
    1965           4 :                         op.applyAdjoint(sino, volume);
    1966           4 :                         REQUIRE_UNARY(
    1967             :                             isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor, backProj)));
    1968             :                     }
    1969             :                 }
    1970             :             }
    1971             :         }
    1972             :     }
    1973             : 
    1974          52 :     GIVEN("A 5x5x5 Volume")
    1975             :     {
    1976             :         // Domain setup
    1977           8 :         const index_t volSize = 5;
    1978             : 
    1979          16 :         IndexVector_t volumeDims(3);
    1980           8 :         volumeDims << volSize, volSize, volSize;
    1981             : 
    1982          16 :         VolumeDescriptor volumeDescriptor(volumeDims);
    1983          16 :         DataContainer<data_t> volume(volumeDescriptor);
    1984             : 
    1985             :         // range setup
    1986           8 :         const index_t detectorSize = 1;
    1987           8 :         const index_t numImgs = 1;
    1988             : 
    1989          16 :         IndexVector_t sinoDims(3);
    1990           8 :         sinoDims << detectorSize, detectorSize, numImgs;
    1991             : 
    1992             :         // Setup geometry
    1993           8 :         const auto stc = SourceToCenterOfRotation{20 * volSize};
    1994           8 :         const auto ctr = CenterOfRotationToDetector{volSize};
    1995          24 :         auto volData = VolumeData3D{Size3D{volumeDims}};
    1996          24 :         auto sinoData = SinogramData3D{Size3D{sinoDims}};
    1997             : 
    1998          16 :         std::vector<Geometry> geom;
    1999             : 
    2000          16 :         GIVEN("Rays not intersecting the volume")
    2001             :         {
    2002           8 :             constexpr index_t numCases = 9;
    2003           8 :             std::array<real_t, numCases> alpha = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    2004           8 :             std::array<real_t, numCases> beta = {0.0, 0.0,      0.0,      0.0,     0.0,
    2005             :                                                  0.0, pi_t / 2, pi_t / 2, pi_t / 2};
    2006           8 :             std::array<real_t, numCases> gamma = {0.0,      0.0,      0.0,      pi_t / 2, pi_t / 2,
    2007             :                                                   pi_t / 2, pi_t / 2, pi_t / 2, pi_t / 2};
    2008           8 :             std::array<real_t, numCases> offsetx = {-volSize, 0.0,      -volSize, 0.0,     0.0,
    2009             :                                                     0.0,      -volSize, 0.0,      -volSize};
    2010           8 :             std::array<real_t, numCases> offsety = {0.0,      -volSize, -volSize, -volSize, 0.0,
    2011             :                                                     -volSize, 0.0,      0.0,      0.0};
    2012           8 :             std::array<real_t, numCases> offsetz = {0.0,      0.0, 0.0,      0.0,     -volSize,
    2013             :                                                     -volSize, 0.0, -volSize, -volSize};
    2014          16 :             std::array<std::string, numCases> neg = {"x",       "y", "x and y", "y",      "z",
    2015             :                                                      "y and z", "x", "z",       "x and z"};
    2016          16 :             std::array<std::string, numCases> ali = {"z", "z", "z", "x", "x", "x", "y", "y", "y"};
    2017             : 
    2018          80 :             for (std::size_t i = 0; i < numCases; i++) {
    2019          76 :                 WHEN("Tracing along a fixed axis-aligned ray with negative coordinate of origin")
    2020             :                 {
    2021           8 :                     INFO("Tracing along a ", ali[i], "-axis-aligned ray with negative ", neg[i],
    2022             :                          "-coodinate of origin");
    2023           8 :                     geom.emplace_back(
    2024           4 :                         stc, ctr, std::move(volData), std::move(sinoData),
    2025             :                         RotationAngles3D{Gamma{gamma[i]}, Beta{beta[i]}, Alpha{alpha[i]}},
    2026             :                         PrincipalPointOffset2D{0, 0},
    2027           4 :                         RotationOffset3D{-offsetx[i], -offsety[i], -offsetz[i]});
    2028             : 
    2029           8 :                     PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
    2030           8 :                     DataContainer<data_t> sino(sinoDescriptor);
    2031           4 :                     sino = 0;
    2032             : 
    2033           8 :                     TestType fast(volumeDescriptor, sinoDescriptor);
    2034           8 :                     TestType slow(volumeDescriptor, sinoDescriptor, false);
    2035             : 
    2036           8 :                     THEN("Result of forward projection is zero")
    2037             :                     {
    2038           4 :                         fast.apply(volume, sino);
    2039           8 :                         DataContainer<data_t> zero(sinoDescriptor);
    2040           4 :                         zero = 0;
    2041           4 :                         REQUIRE_UNARY(isCwiseApprox(sino, zero));
    2042             : 
    2043           4 :                         slow.apply(volume, sino);
    2044           4 :                         REQUIRE_UNARY(isCwiseApprox(sino, zero));
    2045             : 
    2046           8 :                         AND_THEN("Result of backprojection is zero")
    2047             :                         {
    2048           4 :                             fast.applyAdjoint(sino, volume);
    2049           4 :                             DataContainer<data_t> zero(volumeDescriptor);
    2050           4 :                             zero = 0;
    2051           4 :                             REQUIRE_UNARY(isCwiseApprox(volume, zero));
    2052             : 
    2053           4 :                             slow.applyAdjoint(sino, volume);
    2054           4 :                             REQUIRE_UNARY(isCwiseApprox(volume, zero));
    2055             :                         }
    2056             :                     }
    2057             :                 }
    2058             :             }
    2059             :         }
    2060             :     }
    2061          44 : }
    2062          34 : TEST_CASE_TEMPLATE("JosephsMethodCUDA: 3D setup with a multiple rays",TestType, JosephsMethodCUDA<float>,
    2063             :                    JosephsMethodCUDA<double>)
    2064             : {
    2065             :     // Turn logger of
    2066           4 :     Logger::setLevel(Logger::LogLevel::OFF);
    2067             : 
    2068             :     using data_t = decltype(return_data_t(std::declval<TestType>()));
    2069             : 
    2070           8 :     GIVEN("A 3x3x3 volume with 6 projection angles")
    2071             :     {
    2072             :         // Domain setup
    2073           4 :         const index_t volSize = 3;
    2074             : 
    2075           8 :         IndexVector_t volumeDims(3);
    2076           4 :         volumeDims << volSize, volSize, volSize;
    2077             : 
    2078           8 :         VolumeDescriptor volumeDescriptor(volumeDims);
    2079           8 :         DataContainer<data_t> volume(volumeDescriptor);
    2080             : 
    2081             :         // range setup
    2082           4 :         const index_t detectorSize = 1;
    2083           4 :         const index_t numImgs = 6;
    2084             : 
    2085           8 :         IndexVector_t sinoDims(3);
    2086           4 :         sinoDims << detectorSize, detectorSize, numImgs;
    2087             : 
    2088             :         // Setup geometry
    2089           4 :         const auto stc = SourceToCenterOfRotation{20 * volSize};
    2090           4 :         const auto ctr = CenterOfRotationToDetector{volSize};
    2091          12 :         auto volData = VolumeData3D{Size3D{volumeDims}};
    2092          12 :         auto sinoData = SinogramData3D{Size3D{sinoDims}};
    2093             : 
    2094           8 :         std::vector<Geometry> geom;
    2095             : 
    2096           8 :         WHEN("x-, y and z-axis-aligned rays are present")
    2097             :         {
    2098           4 :             real_t beta[numImgs] = {0.0, 0.0, 0.0, 0.0, pi_t / 2, 3 * pi_t / 2};
    2099           4 :             real_t gamma[numImgs] = {0.0, pi_t, pi_t / 2, 3 * pi_t / 2, pi_t / 2, 3 * pi_t / 2};
    2100             : 
    2101          28 :             for (index_t i = 0; i < numImgs; i++)
    2102          24 :                 geom.emplace_back(stc, ctr, VolumeData3D{Size3D{volumeDims}},
    2103             :                                   SinogramData3D{Size3D{sinoDims}},
    2104             :                                   RotationAngles3D{Gamma{gamma[i]}, Beta{beta[i]}});
    2105             : 
    2106           8 :             PlanarDetectorDescriptor sinoDescriptor(sinoDims, geom);
    2107           8 :             DataContainer<data_t> sino(sinoDescriptor);
    2108             : 
    2109           8 :             TestType slow(volumeDescriptor, sinoDescriptor, false);
    2110           8 :             TestType fast(volumeDescriptor, sinoDescriptor);
    2111             : 
    2112           8 :             THEN("Values are accumulated correctly along each ray's path")
    2113             :             {
    2114           4 :                 volume = 0;
    2115             : 
    2116             :                 // set only values along the rays' path to one to make sure interpolation is done
    2117             :                 // correctly
    2118          16 :                 for (index_t i = 0; i < volSize; i++) {
    2119          12 :                     volume(i, volSize / 2, volSize / 2) = 1;
    2120          12 :                     volume(volSize / 2, i, volSize / 2) = 1;
    2121          12 :                     volume(volSize / 2, volSize / 2, i) = 1;
    2122             :                 }
    2123             : 
    2124           4 :                 slow.apply(volume, sino);
    2125          28 :                 for (index_t i = 0; i < numImgs; i++)
    2126          24 :                     REQUIRE_EQ(sino[i], Approx(3.0));
    2127             : 
    2128           4 :                 fast.apply(volume, sino);
    2129          28 :                 for (index_t i = 0; i < numImgs; i++)
    2130          24 :                     REQUIRE_EQ(sino[i], Approx(3.0));
    2131             : 
    2132           8 :                 AND_THEN("Both fast and slow backprojections yield the exact adjoint")
    2133             :                 {
    2134           4 :                     Eigen::Matrix<data_t, Eigen::Dynamic, 1> cmp(volSize * volSize * volSize);
    2135             : 
    2136             :                     // clang-format off
    2137           8 :                     cmp << 0, 0, 0, 0,  6, 0, 0, 0, 0,
    2138           4 :                            0, 6, 0, 6, 18, 6, 0, 6, 0,
    2139           8 :                            0, 0, 0, 0,  6, 0, 0, 0, 0;
    2140             :                     // clang-format on
    2141             : 
    2142           4 :                     slow.applyAdjoint(sino, volume);
    2143           4 :                     REQUIRE_UNARY(
    2144             :                         isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor, cmp)));
    2145             : 
    2146           4 :                     fast.applyAdjoint(sino, volume);
    2147           4 :                     REQUIRE_UNARY(
    2148             :                         isCwiseApprox(volume, DataContainer<data_t>(volumeDescriptor, cmp)));
    2149             :                 }
    2150             :             }
    2151             :         }
    2152             :     }
    2153           4 : }

Generated by: LCOV version 1.15