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           4 : TEST_CASE("CircleTrajectoryGenerator: Create a Circular Trajectory")
      22             : {
      23             :     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           6 :     GIVEN("A 2D descriptor and 256 angles")
      31             :     {
      32           2 :         index_t numberOfAngles = 256;
      33           4 :         IndexVector_t volSize(2);
      34           2 :         volSize << s, s;
      35           4 :         VolumeDescriptor desc{volSize};
      36             : 
      37           3 :         WHEN("We create a half circular trajectory for this scenario")
      38             :         {
      39           1 :             index_t halfCircular = 180;
      40           1 :             real_t diffCenterSource{s * 100};
      41           1 :             real_t diffCenterDetector{s};
      42             : 
      43             :             auto sdesc = CircleTrajectoryGenerator::createTrajectory(
      44           2 :                 numberOfAngles, desc, halfCircular, diffCenterSource, diffCenterDetector);
      45             : 
      46             :             // Check that the detector size is correct
      47           1 :             REQUIRE_EQ(sdesc->getNumberOfCoefficientsPerDimension()[0], expectedDetectorSize);
      48             : 
      49           2 :             THEN("Every geomList in our list has the same camera center and the same projection "
      50             :                  "matrix")
      51             :             {
      52           1 :                 const real_t sourceToCenter = diffCenterSource;
      53           1 :                 const real_t centerToDetector = diffCenterDetector;
      54             : 
      55           2 :                 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             :                     Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
      60             :                                      CenterOfRotationToDetector{centerToDetector},
      61         512 :                                      Radian{currAngle}, VolumeData2D{volSize},
      62         768 :                                      SinogramData2D{sdesc->getSpacingPerDimension(),
      63        1280 :                                                     sdesc->getLocationOfOrigin()});
      64             : 
      65         256 :                     auto geom = sdesc->getGeometryAt(i);
      66         256 :                     CHECK(geom);
      67             : 
      68             :                     const auto centerNorm =
      69         256 :                         (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
      70             :                     const auto projMatNorm =
      71         256 :                         (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
      72             :                     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             :                 }
      79             :             }
      80             :         }
      81             : 
      82           3 :         WHEN("We create a full circular trajectory for this scenario")
      83             :         {
      84           1 :             index_t fullyCircular = 359;
      85           1 :             real_t diffCenterSource{s * 100};
      86           1 :             real_t diffCenterDetector{s};
      87             : 
      88             :             auto sdesc = CircleTrajectoryGenerator::createTrajectory(
      89           2 :                 numberOfAngles, desc, fullyCircular, diffCenterSource, diffCenterDetector);
      90             : 
      91             :             // Check that the detector size is correct
      92           1 :             REQUIRE_EQ(sdesc->getNumberOfCoefficientsPerDimension()[0], expectedDetectorSize);
      93             : 
      94           2 :             THEN("Every geomList in our list has the same camera center and the same projection "
      95             :                  "matrix")
      96             :             {
      97           1 :                 const real_t sourceToCenter = diffCenterSource;
      98           1 :                 const real_t centerToDetector = diffCenterDetector;
      99             : 
     100           2 :                 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             :                     Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
     106             :                                      CenterOfRotationToDetector{centerToDetector},
     107         512 :                                      Radian{currAngle}, VolumeData2D{volSize},
     108         768 :                                      SinogramData2D{sdesc->getSpacingPerDimension(),
     109        1280 :                                                     sdesc->getLocationOfOrigin()});
     110             : 
     111         256 :                     auto geom = sdesc->getGeometryAt(i);
     112         256 :                     CHECK(geom);
     113             : 
     114             :                     const auto centerNorm =
     115         256 :                         (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
     116             :                     const auto projMatNorm =
     117         256 :                         (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
     118             :                     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             :                 }
     125             :             }
     126             :         }
     127             :     }
     128             : 
     129           6 :     GIVEN("A 3D descriptor and 256 angles")
     130             :     {
     131           2 :         index_t numberOfAngles = 256;
     132           4 :         IndexVector_t volSize(3);
     133           2 :         volSize << s, s, s;
     134           4 :         VolumeDescriptor desc{volSize};
     135             : 
     136           3 :         WHEN("We create a half circular trajectory for this scenario")
     137             :         {
     138           1 :             index_t halfCircular = 180;
     139           1 :             real_t diffCenterSource{s * 100};
     140           1 :             real_t diffCenterDetector{s};
     141             : 
     142             :             auto sdesc = CircleTrajectoryGenerator::createTrajectory(
     143           2 :                 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           2 :             THEN("Every geomList in our list has the same camera center and the same projection "
     150             :                  "matrix")
     151             :             {
     152           1 :                 const real_t sourceToCenter = diffCenterSource;
     153           1 :                 const real_t centerToDetector = diffCenterDetector;
     154             : 
     155           2 :                 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             :                     Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
     161             :                                      CenterOfRotationToDetector{centerToDetector},
     162         512 :                                      VolumeData3D{volSize},
     163         768 :                                      SinogramData3D{sdesc->getSpacingPerDimension(),
     164         512 :                                                     sdesc->getLocationOfOrigin()},
     165        1024 :                                      RotationAngles3D{Gamma{angle}});
     166             : 
     167         256 :                     auto geom = sdesc->getGeometryAt(i);
     168         256 :                     CHECK(geom);
     169             : 
     170             :                     const auto centerNorm =
     171         256 :                         (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
     172             :                     const auto projMatNorm =
     173         256 :                         (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
     174             :                     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             :                 }
     181             :             }
     182             :         }
     183           3 :         WHEN("We create a full circular trajectory for this scenario")
     184             :         {
     185           1 :             const index_t fullyCircular = 359;
     186           1 :             real_t diffCenterSource{s * 100};
     187           1 :             real_t diffCenterDetector{s};
     188             : 
     189             :             auto sdesc = CircleTrajectoryGenerator::createTrajectory(
     190           2 :                 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           2 :             THEN("Every geomList in our list has the same camera center and the same projection "
     197             :                  "matrix")
     198             :             {
     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             :                     Geometry tmpGeom(SourceToCenterOfRotation{sourceToCenter},
     208             :                                      CenterOfRotationToDetector{centerToDetector},
     209         512 :                                      VolumeData3D{volSize},
     210         768 :                                      SinogramData3D{sdesc->getSpacingPerDimension(),
     211         512 :                                                     sdesc->getLocationOfOrigin()},
     212        1024 :                                      RotationAngles3D{Gamma{angle}});
     213             : 
     214         256 :                     auto geom = sdesc->getGeometryAt(i);
     215         256 :                     CHECK(geom);
     216             : 
     217             :                     const auto centerNorm =
     218         256 :                         (tmpGeom.getCameraCenter() - geom->getCameraCenter()).norm();
     219             :                     const auto projMatNorm =
     220         256 :                         (tmpGeom.getProjectionMatrix() - geom->getProjectionMatrix()).norm();
     221             :                     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             :                 }
     228             :             }
     229             :         }
     230             :     }
     231           4 : }