Line data    Source code

       1             : /**
       2             :  * @file test_PlanarDetectorDescriptor.cpp
       3             :  *
       4             :  * @brief Test for PlanarDetectorDescriptor
       5             :  *
       6             :  * @author David Frank - initial code
       7             :  */
       8             : 
       9             : #include "doctest/doctest.h"
      10             : 
      11             : #include "PlanarDetectorDescriptor.h"
      12             : #include "VolumeDescriptor.h"
      13             : 
      14             : using namespace elsa;
      15             : using namespace elsa::geometry;
      16             : using namespace doctest;
      17             : 
      18             : TEST_SUITE_BEGIN("core");
      19             : 
      20             : TEST_CASE("PlanarDetectorDescriptor: Testing 2D PlanarDetectorDescriptor")
      21          10 : {
      22          10 :     GIVEN("Given a 5x5 Volume and a single 5 wide detector pose")
      23          10 :     {
      24           4 :         IndexVector_t volSize(2);
      25           4 :         volSize << 5, 5;
      26           4 :         VolumeDescriptor ddVol(volSize);
      27             : 
      28           4 :         IndexVector_t sinoSize(2);
      29           4 :         sinoSize << 5, 1;
      30             : 
      31           4 :         real_t s2c = 10;
      32           4 :         real_t c2d = 4;
      33             : 
      34           4 :         Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Radian{0},
      35           4 :                    VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
      36             : 
      37           4 :         PlanarDetectorDescriptor desc(sinoSize, {g});
      38             : 
      39           4 :         WHEN("Retrieving the single geometry pose")
      40           4 :         {
      41           1 :             auto geom = desc.getGeometryAt(0);
      42           1 :             auto geomList = desc.getGeometry();
      43             : 
      44           1 :             CHECK_EQ(desc.getNumberOfGeometryPoses(), 1);
      45           1 :             CHECK_EQ(geomList.size(), 1);
      46             : 
      47           1 :             THEN("Geometry is equal")
      48           1 :             {
      49           1 :                 CHECK_EQ((geom), g);
      50           1 :                 CHECK_EQ(geomList[0], g);
      51           1 :             }
      52           1 :         }
      53             : 
      54           4 :         WHEN("Generating rays for detector pixels 0, 2 and 4")
      55           4 :         {
      56           3 :             for (real_t detPixel : std::initializer_list<real_t>{0, 2, 4}) {
      57           3 :                 RealVector_t pixel(1);
      58           3 :                 pixel << detPixel + 0.5f;
      59             : 
      60             :                 // Check that ray for IndexVector_t is equal to previous one
      61           3 :                 auto ray = desc.computeRayFromDetectorCoord(pixel, 0);
      62             : 
      63             :                 // Create variables, which make typing quicker
      64           3 :                 auto ro = ray.origin();
      65           3 :                 auto rd = ray.direction();
      66             : 
      67             :                 // Check that ray origin is camera center
      68           3 :                 auto c = g.getCameraCenter();
      69           3 :                 CHECK_EQ((ro - c).sum(), Approx(0));
      70             : 
      71             :                 // compute intersection manually
      72           3 :                 real_t t = Approx(rd[0]) == 0 ? (s2c + c2d) : ((pixel[0] - ro[0]) / rd[0]);
      73             : 
      74           3 :                 auto detCoordY = ro[1] + t * rd[1];
      75             : 
      76           3 :                 CHECK_EQ(detCoordY, Approx(ddVol.getLocationOfOrigin()[1] + c2d));
      77           3 :             }
      78           1 :         }
      79             : 
      80           4 :         WHEN("Center axis voxels are projected to middle detector pixel with correct scaling")
      81           4 :         {
      82           5 :             for (real_t slice : std::initializer_list<real_t>{0, 1, 2, 3, 4}) {
      83           5 :                 RealVector_t voxelCoord{{2.f, slice}};
      84             : 
      85             :                 // move to voxel center
      86           5 :                 voxelCoord = voxelCoord.array() + 0.5f;
      87             : 
      88             :                 // Check that detector Pixel is the center one
      89           5 :                 auto [pixel, scaling] = desc.projectAndScaleVoxelOnDetector(voxelCoord, 0);
      90           5 :                 real_t pixelIndex = pixel[0] - 0.5f;
      91             : 
      92           5 :                 CHECK_EQ(pixelIndex, Approx(2));
      93             : 
      94             :                 // verify scaling
      95           5 :                 auto correctScaling = g.getSourceDetectorDistance() / (s2c - 2.f + slice);
      96           5 :                 CHECK_EQ(scaling, Approx(correctScaling));
      97           5 :             }
      98           1 :         }
      99             : 
     100           4 :         WHEN("All voxels are projected to correct detector pixel with correct scaling")
     101           4 :         {
     102           5 :             for (real_t x : std::initializer_list<real_t>{0, 1, 2, 3, 4}) {
     103          25 :                 for (real_t y : std::initializer_list<real_t>{0, 1, 2, 3, 4}) {
     104          25 :                     RealVector_t voxelCoord{{x, y}};
     105             : 
     106          25 :                     auto centerAxisOffset = x - 2.f;
     107             : 
     108             :                     // move to voxel center
     109          25 :                     voxelCoord = voxelCoord.array() + 0.5f;
     110             : 
     111          25 :                     auto [pixel, scaling] = desc.projectAndScaleVoxelOnDetector(voxelCoord, 0);
     112          25 :                     real_t pixelIndex = pixel[0] - 0.5f;
     113             : 
     114          25 :                     auto zAxisDistance = s2c - 2.f + y;
     115          25 :                     auto s2v = std::sqrt(centerAxisOffset * centerAxisOffset
     116          25 :                                          + zAxisDistance * zAxisDistance);
     117             : 
     118          25 :                     auto correctScaling = g.getSourceDetectorDistance() / s2v;
     119             :                     // verify scaling
     120          25 :                     CHECK_EQ(scaling, Approx(correctScaling));
     121             : 
     122             :                     // verify detector pixel
     123          25 :                     auto scaled_center_offset =
     124          25 :                         centerAxisOffset * g.getSourceDetectorDistance() / zAxisDistance;
     125          25 :                     CHECK_EQ(pixelIndex, Approx(2 + scaled_center_offset));
     126          25 :                 }
     127           5 :             }
     128           1 :         }
     129           4 :     }
     130             : 
     131          10 :     GIVEN("Given a 5x5 Volume and a multiple 5 wide detector pose")
     132          10 :     {
     133           6 :         IndexVector_t volSize(2);
     134           6 :         volSize << 5, 5;
     135           6 :         VolumeDescriptor ddVol(volSize);
     136             : 
     137           6 :         IndexVector_t sinoSize(2);
     138           6 :         sinoSize << 5, 4;
     139             : 
     140           6 :         real_t s2c = 10;
     141           6 :         real_t c2d = 4;
     142             : 
     143           6 :         Geometry g1(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{0},
     144           6 :                     VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
     145           6 :         Geometry g2(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{90},
     146           6 :                     VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
     147           6 :         Geometry g3(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{180},
     148           6 :                     VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
     149           6 :         Geometry g4(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d}, Degree{270},
     150           6 :                     VolumeData2D{Size2D{volSize}}, SinogramData2D{Size2D{sinoSize}});
     151             : 
     152           6 :         PlanarDetectorDescriptor desc(sinoSize, {g1, g2, g3, g4});
     153             : 
     154           6 :         WHEN("Retreiving geometry poses")
     155           6 :         {
     156           4 :             auto geom = desc.getGeometryAt(0);
     157           4 :             auto geomList = desc.getGeometry();
     158             : 
     159           4 :             CHECK_EQ(desc.getNumberOfGeometryPoses(), 4);
     160           4 :             CHECK_EQ(geomList.size(), 4);
     161             : 
     162           4 :             THEN("Geometry is equal")
     163           4 :             {
     164           1 :                 CHECK_EQ((desc.getGeometryAt(0)), g1);
     165           1 :                 CHECK_EQ(geomList[0], g1);
     166           1 :             }
     167           4 :             THEN("Geometry is equal")
     168           4 :             {
     169           1 :                 CHECK_EQ((desc.getGeometryAt(1)), g2);
     170           1 :                 CHECK_EQ(geomList[1], g2);
     171           1 :             }
     172           4 :             THEN("Geometry is equal")
     173           4 :             {
     174           1 :                 CHECK_EQ((desc.getGeometryAt(2)), g3);
     175           1 :                 CHECK_EQ(geomList[2], g3);
     176           1 :             }
     177           4 :             THEN("Geometry is equal")
     178           4 :             {
     179           1 :                 CHECK_EQ((desc.getGeometryAt(3)), g4);
     180           1 :                 CHECK_EQ(geomList[3], g4);
     181           1 :             }
     182           4 :         }
     183             : 
     184           6 :         WHEN("Check for multiple poses, that all the overloads compute the same rays")
     185           6 :         {
     186           4 :             for (index_t pose : {0, 1, 2, 3}) {
     187             : 
     188          12 :                 for (index_t detPixel : {0, 2, 4}) {
     189          12 :                     IndexVector_t pixel(2);
     190          12 :                     pixel << detPixel, pose;
     191             : 
     192          12 :                     RealVector_t pixelReal(1);
     193          12 :                     pixelReal << static_cast<real_t>(detPixel) + 0.5f;
     194             : 
     195          12 :                     auto ray1 =
     196          12 :                         desc.computeRayFromDetectorCoord(desc.getIndexFromCoordinate(pixel));
     197          12 :                     auto ray2 = desc.computeRayFromDetectorCoord(pixel);
     198          12 :                     auto ray3 = desc.computeRayFromDetectorCoord(pixelReal, pose);
     199             : 
     200          12 :                     auto ro1 = ray1.origin();
     201          12 :                     auto rd1 = ray1.direction();
     202             : 
     203          12 :                     auto ro2 = ray2.origin();
     204          12 :                     auto rd2 = ray2.direction();
     205             : 
     206          12 :                     auto ro3 = ray3.origin();
     207          12 :                     auto rd3 = ray3.direction();
     208             : 
     209          12 :                     CHECK_EQ(ro1, ro2);
     210          12 :                     CHECK_EQ(ro1, ro3);
     211             : 
     212          12 :                     CHECK_EQ(rd1, rd2);
     213          12 :                     CHECK_EQ(rd1, rd3);
     214             : 
     215             :                     // Shouldn't be necessary, but whatever
     216          12 :                     CHECK_EQ(ro2, ro3);
     217          12 :                     CHECK_EQ(rd2, rd3);
     218          12 :                 }
     219           4 :             }
     220           1 :         }
     221             : 
     222           6 :         WHEN("Center voxels are projected to middle detector pixel with correct scaling")
     223           6 :         {
     224           4 :             for (index_t pose : {0, 1, 2, 3}) {
     225           4 :                 RealVector_t voxelCoord{{2.f, 2.f}};
     226             : 
     227             :                 // move to voxel center
     228           4 :                 voxelCoord = voxelCoord.array() + 0.5f;
     229             : 
     230             :                 // Check that detector Pixel is the center one
     231           4 :                 auto [pixel, scaling] = desc.projectAndScaleVoxelOnDetector(voxelCoord, pose);
     232           4 :                 real_t pixelIndex = pixel[0] - 0.5f;
     233             : 
     234           4 :                 CHECK_EQ(pixelIndex, Approx(2));
     235             : 
     236             :                 // verify scaling
     237           4 :                 auto correctScaling = (s2c + c2d) / (s2c);
     238           4 :                 CHECK_EQ(scaling, Approx(correctScaling));
     239           4 :             }
     240           1 :         }
     241           6 :     }
     242          10 : }
     243             : 
     244             : TEST_CASE("PlanarDetectorDescriptor: Testing 3D PlanarDetectorDescriptor")
     245           4 : {
     246           4 :     GIVEN("Given a 5x5x5 Volume and a single 5x5 wide detector pose")
     247           4 :     {
     248           4 :         IndexVector_t volSize(3);
     249           4 :         volSize << 5, 5, 5;
     250           4 :         VolumeDescriptor ddVol(volSize);
     251             : 
     252           4 :         IndexVector_t sinoSize(3);
     253           4 :         sinoSize << 5, 5, 1;
     254             : 
     255           4 :         real_t s2c = 10;
     256           4 :         real_t c2d = 4;
     257             : 
     258           4 :         Geometry g(SourceToCenterOfRotation{s2c}, CenterOfRotationToDetector{c2d},
     259           4 :                    VolumeData3D{Size3D{volSize}}, SinogramData3D{Size3D{sinoSize}},
     260           4 :                    RotationAngles3D{Gamma{0}});
     261             : 
     262           4 :         PlanarDetectorDescriptor desc(sinoSize, {g});
     263             : 
     264           4 :         WHEN("Retrieving the single geometry pose")
     265           4 :         {
     266           1 :             auto geom = desc.getGeometryAt(0);
     267           1 :             auto geomList = desc.getGeometry();
     268             : 
     269           1 :             CHECK_EQ(desc.getNumberOfGeometryPoses(), 1);
     270           1 :             CHECK_EQ(geomList.size(), 1);
     271             : 
     272           1 :             THEN("Geometry is equal")
     273           1 :             {
     274           1 :                 CHECK_EQ((geom), g);
     275           1 :                 CHECK_EQ(geomList[0], g);
     276           1 :             }
     277           1 :         }
     278             : 
     279           4 :         WHEN("Generating rays for detector pixels 0, 2 and 4 for each dim")
     280           4 :         {
     281           3 :             for (index_t detPixel1 : {0, 2, 4}) {
     282           9 :                 for (index_t detPixel2 : {0, 2, 4}) {
     283           9 :                     RealVector_t pixel(2);
     284           9 :                     pixel << static_cast<real_t>(detPixel1) + 0.5f,
     285           9 :                         static_cast<real_t>(detPixel2) + 0.5f;
     286             : 
     287             :                     // Check that ray for IndexVector_t is equal to previous one
     288           9 :                     auto ray = desc.computeRayFromDetectorCoord(pixel, 0);
     289             : 
     290             :                     // Create variables, which make typing quicker
     291           9 :                     auto ro = ray.origin();
     292           9 :                     auto rd = ray.direction();
     293             : 
     294             :                     // Check that ray origin is camera center
     295           9 :                     auto c = g.getCameraCenter();
     296           9 :                     CHECK_EQ((ro - c).sum(), Approx(0));
     297             : 
     298           9 :                     auto o = ddVol.getLocationOfOrigin();
     299           9 :                     RealVector_t detCoordWorld(3);
     300           9 :                     detCoordWorld << pixel[0] - o[0], pixel[1] - o[1], c2d;
     301           9 :                     RealVector_t rotD = g.getRotationMatrix().transpose() * detCoordWorld + o;
     302             : 
     303           9 :                     real_t factor = 0;
     304           9 :                     if (std::abs(rd[0]) > 0)
     305           6 :                         factor = (rotD[0] - ro[0]) / rd[0];
     306           3 :                     else if (std::abs(rd[1]) > 0)
     307           2 :                         factor = (rotD[1] - ro[1]) / rd[1];
     308           1 :                     else if (std::abs(rd[2]) > 0)
     309           1 :                         factor = (rotD[2] - ro[2] / rd[2]);
     310             : 
     311           9 :                     CHECK_EQ((ro[0] + factor * rd[0]), Approx(rotD[0]));
     312           9 :                     CHECK_EQ((ro[1] + factor * rd[1]), Approx(rotD[1]));
     313           9 :                     CHECK_EQ((ro[2] + factor * rd[2]), Approx(rotD[2]));
     314           9 :                 }
     315           3 :             }
     316           1 :         }
     317             : 
     318           4 :         WHEN("Center voxels are projected to middle detector pixel with correct scaling")
     319           4 :         {
     320           5 :             for (real_t slice : std::initializer_list<real_t>{0, 1, 2, 3, 4}) {
     321           5 :                 RealVector_t voxelCoord{{2.f, 2.f, slice}};
     322             : 
     323             :                 // move to voxel center
     324           5 :                 voxelCoord = voxelCoord.array() + 0.5f;
     325             : 
     326             :                 // Check that detector Pixel is the center one
     327           5 :                 auto [pixel, scaling] = desc.projectAndScaleVoxelOnDetector(voxelCoord, 0);
     328           5 :                 real_t pixelIndex = pixel[0] - 0.5f;
     329             : 
     330           5 :                 CHECK_EQ(pixelIndex, Approx(2));
     331             : 
     332             :                 // verify scaling
     333           5 :                 auto correctScaling = (s2c + c2d) / (s2c - 2.f + slice);
     334           5 :                 CHECK_EQ(scaling, Approx(correctScaling));
     335           5 :             }
     336           1 :         }
     337             : 
     338           4 :         WHEN("All voxels are projected to correct detector pixel with correct scaling")
     339           4 :         {
     340           5 :             for (real_t x : std::initializer_list<real_t>{0, 1, 2, 3, 4}) {
     341          25 :                 for (real_t y : std::initializer_list<real_t>{0, 1, 2, 3, 4}) {
     342         125 :                     for (real_t z : std::initializer_list<real_t>{0, 1, 2, 3, 4}) {
     343         125 :                         RealVector_t voxelCoord{{x, y, z}};
     344             : 
     345         125 :                         RealVector_t centeredVoxelCoord = voxelCoord.array() - 2;
     346         125 :                         auto centerAxisOffset = centeredVoxelCoord.head(2).norm();
     347             : 
     348             :                         // move to voxel center
     349         125 :                         voxelCoord = voxelCoord.array() + 0.5f;
     350             : 
     351         125 :                         auto [pixel, scaling] = desc.projectAndScaleVoxelOnDetector(voxelCoord, 0);
     352         125 :                         RealVector_t pixelIndex = pixel.head(2).array() - 0.5f;
     353             : 
     354         125 :                         auto zAxisDistance = s2c - 2.f + z;
     355         125 :                         auto s2v = std::sqrt(centerAxisOffset * centerAxisOffset
     356         125 :                                              + zAxisDistance * zAxisDistance);
     357             : 
     358         125 :                         auto correctScaling = g.getSourceDetectorDistance() / s2v;
     359             :                         // verify scaling
     360         125 :                         CHECK_EQ(scaling, Approx(correctScaling));
     361             : 
     362             :                         // verify detector pixel
     363         125 :                         auto scaled_x_offset =
     364         125 :                             centeredVoxelCoord[0] * g.getSourceDetectorDistance() / zAxisDistance;
     365         125 :                         auto scaled_y_offset =
     366         125 :                             centeredVoxelCoord[1] * g.getSourceDetectorDistance() / zAxisDistance;
     367         125 :                         CHECK_EQ(pixelIndex[0], Approx(2 + scaled_x_offset));
     368         125 :                         CHECK_EQ(pixelIndex[1], Approx(2 + scaled_y_offset));
     369         125 :                     }
     370          25 :                 }
     371           5 :             }
     372           1 :         }
     373           4 :     }
     374           4 : }
     375             : 
     376             : TEST_SUITE_END();