Line data    Source code

       1             : /**
       2             :  * @file test_CircleTrajectoryGenerator.cpp
       3             :  *
       4             :  * @brief Test for CircleTrajectoryGenerator class
       5             :  *
       6             :  * @author David Frank - initial code
       7             :  * @author Nikola Dinev - fixes
       8             :  * @author Tobias Lasser - modernization, fixes
       9             :  */
      10             : 
      11             : #include "doctest/doctest.h"
      12             : 
      13             : #include "CircleTrajectoryGenerator.h"
      14             : #include "Logger.h"
      15             : #include "VolumeDescriptor.h"
      16             : #include "testHelpers.h"
      17             : 
      18             : using namespace elsa;
      19             : using namespace doctest;
      20             : 
      21             : TEST_CASE("CircleTrajectoryGenerator: Create a Circular Trajectory")
      22           4 : {
      23           4 :     using namespace geometry;
      24             : 
      25           4 :     const index_t s = 64;
      26             : 
      27             :     // Detector size is the volume size scalled by the square root of 2
      28           4 :     const auto expectedDetectorSize = static_cast<index_t>(s * std::sqrt(2));
      29             : 
      30           4 :     GIVEN("A 2D descriptor and 256 angles")
      31           4 :     {
      32           2 :         index_t numberOfAngles = 256;
      33           2 :         IndexVector_t volSize(2);
      34           2 :         volSize << s, s;
      35           2 :         VolumeDescriptor desc{volSize};
      36             : 
      37           2 :         WHEN("We create a half circular trajectory for this scenario")
      38           2 :         {
      39           1 :             index_t halfCircular = 180;
      40           1 :             real_t diffCenterSource{s * 100};
      41           1 :             real_t diffCenterDetector{s};
      42             : 
      43           1 :             auto sdesc = CircleTrajectoryGenerator::createTrajectory(
      44           1 :                 numberOfAngles, desc, halfCircular, diffCenterSource, diffCenterDetector);
      45             : 
      46             :             // Check that the detector size is correct
      47           1 :             REQUIRE_EQ(sdesc->getNumberOfCoefficientsPerDimension()[0], expectedDetectorSize);
      48             : 
      49           1 :             THEN("Every geomList in our list has the same camera center and the same projection "
      50           1 :                  "matrix")
      51           1 :             {
      52           1 :                 const real_t sourceToCenter = diffCenterSource;
      53           1 :                 const real_t centerToDetector = diffCenterDetector;
      54             : 
      55           1 :                 real_t angle = static_cast<real_t>(1.0) * static_cast<real_t>(halfCircular)
      56           1 :                                / real_t(numberOfAngles - 1);
      57         257 :                 for (index_t i = 0; i < numberOfAngles; ++i) {
      58         256 :                     real_t currAngle = static_cast<real_t>(i) * angle * pi_t / 180.0f;
      59         256 :                     Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
      60         256 :                                      CenterOfRotationToDetector{centerToDetector},
      61         256 :                                      Radian{currAngle}, VolumeData2D{volSize},
      62         256 :                                      SinogramData2D{sdesc->getSpacingPerDimension(),
      63         256 :                                                     sdesc->getLocationOfOrigin()});
      64             : 
      65         256 :                     auto geom = sdesc->getGeometryAt(i);
      66         256 :                     CHECK(geom);
      67             : 
      68         256 :                     const auto centerNorm =
      69         256 :                         (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
      70         256 :                     const auto projMatNorm =
      71         256 :                         (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
      72         256 :                     const auto invProjMatNorm =
      73         256 :                         (tmpGeom.getInverseProjectionMatrix() - geom->getInverseProjectionMatrix())
      74         256 :                             .norm();
      75         256 :                     REQUIRE_UNARY(checkApproxEq(centerNorm, 0));
      76         256 :                     REQUIRE_UNARY(checkApproxEq(projMatNorm, 0, 0.0000001));
      77         256 :                     REQUIRE_UNARY(checkApproxEq(invProjMatNorm, 0, 0.0000001));
      78         256 :                 }
      79           1 :             }
      80           1 :         }
      81             : 
      82           2 :         WHEN("We create a full circular trajectory for this scenario")
      83           2 :         {
      84           1 :             index_t fullyCircular = 359;
      85           1 :             real_t diffCenterSource{s * 100};
      86           1 :             real_t diffCenterDetector{s};
      87             : 
      88           1 :             auto sdesc = CircleTrajectoryGenerator::createTrajectory(
      89           1 :                 numberOfAngles, desc, fullyCircular, diffCenterSource, diffCenterDetector);
      90             : 
      91             :             // Check that the detector size is correct
      92           1 :             REQUIRE_EQ(sdesc->getNumberOfCoefficientsPerDimension()[0], expectedDetectorSize);
      93             : 
      94           1 :             THEN("Every geomList in our list has the same camera center and the same projection "
      95           1 :                  "matrix")
      96           1 :             {
      97           1 :                 const real_t sourceToCenter = diffCenterSource;
      98           1 :                 const real_t centerToDetector = diffCenterDetector;
      99             : 
     100           1 :                 real_t angle = static_cast<real_t>(1.0) * static_cast<real_t>(fullyCircular)
     101           1 :                                / static_cast<real_t>(numberOfAngles - 1);
     102         257 :                 for (index_t i = 0; i < numberOfAngles; ++i) {
     103         256 :                     real_t currAngle = static_cast<real_t>(i) * angle * pi_t / 180.0f;
     104             : 
     105         256 :                     Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
     106         256 :                                      CenterOfRotationToDetector{centerToDetector},
     107         256 :                                      Radian{currAngle}, VolumeData2D{volSize},
     108         256 :                                      SinogramData2D{sdesc->getSpacingPerDimension(),
     109         256 :                                                     sdesc->getLocationOfOrigin()});
     110             : 
     111         256 :                     auto geom = sdesc->getGeometryAt(i);
     112         256 :                     CHECK(geom);
     113             : 
     114         256 :                     const auto centerNorm =
     115         256 :                         (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
     116         256 :                     const auto projMatNorm =
     117         256 :                         (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
     118         256 :                     const auto invProjMatNorm =
     119         256 :                         (tmpGeom.getInverseProjectionMatrix() - geom->getInverseProjectionMatrix())
     120         256 :                             .norm();
     121         256 :                     REQUIRE_UNARY(checkApproxEq(centerNorm, 0));
     122         256 :                     REQUIRE_UNARY(checkApproxEq(projMatNorm, 0, 0.0000001));
     123         256 :                     REQUIRE_UNARY(checkApproxEq(invProjMatNorm, 0, 0.0000001));
     124         256 :                 }
     125           1 :             }
     126           1 :         }
     127           2 :     }
     128             : 
     129           4 :     GIVEN("A 3D descriptor and 256 angles")
     130           4 :     {
     131           2 :         index_t numberOfAngles = 256;
     132           2 :         IndexVector_t volSize(3);
     133           2 :         volSize << s, s, s;
     134           2 :         VolumeDescriptor desc{volSize};
     135             : 
     136           2 :         WHEN("We create a half circular trajectory for this scenario")
     137           2 :         {
     138           1 :             index_t halfCircular = 180;
     139           1 :             real_t diffCenterSource{s * 100};
     140           1 :             real_t diffCenterDetector{s};
     141             : 
     142           1 :             auto sdesc = CircleTrajectoryGenerator::createTrajectory(
     143           1 :                 numberOfAngles, desc, halfCircular, diffCenterSource, diffCenterDetector);
     144             : 
     145             :             // Check that the detector size is correct
     146           1 :             REQUIRE_EQ(sdesc->getNumberOfCoefficientsPerDimension()[0], expectedDetectorSize);
     147           1 :             REQUIRE_EQ(sdesc->getNumberOfCoefficientsPerDimension()[1], expectedDetectorSize);
     148             : 
     149           1 :             THEN("Every geomList in our list has the same camera center and the same projection "
     150           1 :                  "matrix")
     151           1 :             {
     152           1 :                 const real_t sourceToCenter = diffCenterSource;
     153           1 :                 const real_t centerToDetector = diffCenterDetector;
     154             : 
     155           1 :                 real_t angleInc = 1.0f * static_cast<real_t>(halfCircular)
     156           1 :                                   / static_cast<real_t>(numberOfAngles - 1);
     157         257 :                 for (index_t i = 0; i < numberOfAngles; ++i) {
     158         256 :                     real_t angle = static_cast<real_t>(i) * angleInc * pi_t / 180.0f;
     159             : 
     160         256 :                     Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
     161         256 :                                      CenterOfRotationToDetector{centerToDetector},
     162         256 :                                      VolumeData3D{volSize},
     163         256 :                                      SinogramData3D{sdesc->getSpacingPerDimension(),
     164         256 :                                                     sdesc->getLocationOfOrigin()},
     165         256 :                                      RotationAngles3D{Gamma{angle}});
     166             : 
     167         256 :                     auto geom = sdesc->getGeometryAt(i);
     168         256 :                     CHECK(geom);
     169             : 
     170         256 :                     const auto centerNorm =
     171         256 :                         (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
     172         256 :                     const auto projMatNorm =
     173         256 :                         (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
     174         256 :                     const auto invProjMatNorm =
     175         256 :                         (tmpGeom.getInverseProjectionMatrix() - geom->getInverseProjectionMatrix())
     176         256 :                             .norm();
     177         256 :                     REQUIRE(checkApproxEq(centerNorm, 0));
     178         256 :                     REQUIRE(checkApproxEq(projMatNorm, 0, 0.0000001));
     179         256 :                     REQUIRE(checkApproxEq(invProjMatNorm, 0, 0.0000001));
     180         256 :                 }
     181           1 :             }
     182           1 :         }
     183           2 :         WHEN("We create a full circular trajectory for this scenario")
     184           2 :         {
     185           1 :             const index_t fullyCircular = 359;
     186           1 :             real_t diffCenterSource{s * 100};
     187           1 :             real_t diffCenterDetector{s};
     188             : 
     189           1 :             auto sdesc = CircleTrajectoryGenerator::createTrajectory(
     190           1 :                 numberOfAngles, desc, fullyCircular, diffCenterSource, diffCenterDetector);
     191             : 
     192             :             // Check that the detector size is correct
     193           1 :             REQUIRE_EQ(sdesc->getNumberOfCoefficientsPerDimension()[0], expectedDetectorSize);
     194           1 :             REQUIRE_EQ(sdesc->getNumberOfCoefficientsPerDimension()[1], expectedDetectorSize);
     195             : 
     196           1 :             THEN("Every geomList in our list has the same camera center and the same projection "
     197           1 :                  "matrix")
     198           1 :             {
     199           1 :                 const auto sourceToCenter = diffCenterSource;
     200           1 :                 const auto centerToDetector = diffCenterDetector;
     201             : 
     202           1 :                 real_t angleInc = 1.0f * static_cast<real_t>(fullyCircular)
     203           1 :                                   / static_cast<real_t>(numberOfAngles - 1);
     204         257 :                 for (index_t i = 0; i < numberOfAngles; ++i) {
     205         256 :                     real_t angle = static_cast<real_t>(i) * angleInc * pi_t / 180.0f;
     206             : 
     207         256 :                     Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
     208         256 :                                      CenterOfRotationToDetector{centerToDetector},
     209         256 :                                      VolumeData3D{volSize},
     210         256 :                                      SinogramData3D{sdesc->getSpacingPerDimension(),
     211         256 :                                                     sdesc->getLocationOfOrigin()},
     212         256 :                                      RotationAngles3D{Gamma{angle}});
     213             : 
     214         256 :                     auto geom = sdesc->getGeometryAt(i);
     215         256 :                     CHECK(geom);
     216             : 
     217         256 :                     const auto centerNorm =
     218         256 :                         (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
     219         256 :                     const auto projMatNorm =
     220         256 :                         (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
     221         256 :                     const auto invProjMatNorm =
     222         256 :                         (tmpGeom.getInverseProjectionMatrix() - geom->getInverseProjectionMatrix())
     223         256 :                             .norm();
     224         256 :                     REQUIRE_UNARY(checkApproxEq(centerNorm, 0));
     225         256 :                     REQUIRE_UNARY(checkApproxEq(projMatNorm, 0, 0.0000001));
     226         256 :                     REQUIRE_UNARY(checkApproxEq(invProjMatNorm, 0, 0.0000001));
     227         256 :                 }
     228           1 :             }
     229           1 :         }
     230           2 :     }
     231           4 : }