Line data    Source code

       1             : /**
       2             :  * @file test_SubsetSampler.cpp
       3             :  *
       4             :  * @brief Tests for the SubsetSampler class
       5             :  *
       6             :  * @author Michael Loipführer - initial code
       7             :  */
       8             : #include <SphereTrajectoryGenerator.h>
       9             : #include <numeric>
      10             : #include "doctest/doctest.h"
      11             : 
      12             : #include "SubsetSampler.h"
      13             : #include "Logger.h"
      14             : #include "CircleTrajectoryGenerator.h"
      15             : #include "SiddonsMethod.h"
      16             : #include "PlanarDetectorDescriptor.h"
      17             : #include "PhantomGenerator.h"
      18             : #include "SiddonsMethod.h"
      19             : #include "JosephsMethod.h"
      20             : 
      21             : using namespace elsa;
      22             : using namespace doctest;
      23             : 
      24           6 : TEST_CASE("SubsetSampler: Testing subset sampling strategies")
      25             : {
      26           6 :     Logger::setLevel(Logger::LogLevel::WARN);
      27             : 
      28           9 :     GIVEN("A circular trajectory with 32 angles")
      29             :     {
      30           6 :         IndexVector_t size(2);
      31           3 :         size << 32, 32;
      32           6 :         VolumeDescriptor volumeDescriptor{size};
      33           3 :         index_t numAngles{32}, arc{360};
      34             :         auto sinoDescriptor = CircleTrajectoryGenerator::createTrajectory(
      35           6 :             numAngles, volumeDescriptor, arc, static_cast<real_t>(size(0)) * 100.0f,
      36           6 :             static_cast<real_t>(size(0)));
      37           6 :         const auto numCoeffsPerDim = sinoDescriptor->getNumberOfCoefficientsPerDimension();
      38           3 :         const index_t nDimensions = sinoDescriptor->getNumberOfDimensions();
      39           3 :         const auto numElements = numCoeffsPerDim[nDimensions - 1];
      40             : 
      41           4 :         WHEN("performing round robin sampling into 4 subsets")
      42             :         {
      43           1 :             const auto nSubsets = 4;
      44           2 :             std::vector<index_t> indices(static_cast<std::size_t>(numElements));
      45           1 :             std::iota(indices.begin(), indices.end(), 0);
      46             :             const std::vector<std::vector<index_t>> mapping =
      47           2 :                 SubsetSampler<PlanarDetectorDescriptor, real_t>::splitRoundRobin(indices, nSubsets);
      48           2 :             THEN("The mapping is correct with every subset having 8 elements")
      49             :             {
      50           5 :                 for (std::size_t i = 0; i < static_cast<std::size_t>(nSubsets); ++i) {
      51           4 :                     REQUIRE_EQ(mapping[i].size(), 8);
      52          36 :                     for (std::size_t j = 0; j < mapping[i].size(); j++) {
      53          32 :                         REQUIRE_EQ(mapping[i][j], j * nSubsets + i);
      54             :                     }
      55             :                 }
      56             :             }
      57             :         }
      58           4 :         WHEN("performing round robin sampling into 5 subsets")
      59             :         {
      60           1 :             const auto nSubsets = 5;
      61           2 :             std::vector<index_t> indices(static_cast<std::size_t>(numElements));
      62           1 :             std::iota(indices.begin(), indices.end(), 0);
      63             :             const auto mapping =
      64           2 :                 SubsetSampler<PlanarDetectorDescriptor, real_t>::splitRoundRobin(indices, nSubsets);
      65           2 :             THEN("The mapping is correct with every subset having 6 elements apart from the first "
      66             :                  "two")
      67             :             {
      68           6 :                 for (std::size_t i = 0; i < static_cast<std::size_t>(nSubsets); ++i) {
      69           5 :                     if (i <= 1) {
      70           2 :                         REQUIRE_EQ(mapping[i].size(), 7);
      71             :                     } else {
      72           3 :                         REQUIRE_EQ(mapping[i].size(), 6);
      73             :                     }
      74          37 :                     for (std::size_t j = 0; j < mapping[i].size(); j++) {
      75          32 :                         REQUIRE_EQ(mapping[i][j], j * nSubsets + i);
      76             :                     }
      77             :                 }
      78             :             }
      79             :         }
      80           4 :         WHEN("performing equi rotation sampling into 4 subsets")
      81             :         {
      82           1 :             const auto nSubsets = 4;
      83             :             const std::vector<std::vector<index_t>> mapping =
      84             :                 SubsetSampler<PlanarDetectorDescriptor, real_t>::splitRotationalClustering(
      85           2 :                     static_cast<const PlanarDetectorDescriptor&>(*sinoDescriptor), nSubsets);
      86           2 :             THEN("The mapping is correct with every subset having 8 elements")
      87             :             {
      88           1 :                 REQUIRE_EQ(mapping[0], std::vector<index_t>{15, 11, 7, 3, 0, 27, 23, 19});
      89           1 :                 REQUIRE_EQ(mapping[1], std::vector<index_t>{14, 10, 6, 2, 30, 26, 22, 18});
      90           1 :                 REQUIRE_EQ(mapping[2], std::vector<index_t>{13, 9, 5, 1, 29, 25, 21, 17});
      91           1 :                 REQUIRE_EQ(mapping[3], std::vector<index_t>{12, 8, 4, 31, 28, 24, 20, 16});
      92             :             }
      93             :         }
      94             :     }
      95             : 
      96           9 :     GIVEN("A spherical trajectory with 32 angles and 4 circles")
      97             :     {
      98           6 :         IndexVector_t size(3);
      99           3 :         size << 32, 32, 32;
     100           6 :         VolumeDescriptor volumeDescriptor{size};
     101           3 :         index_t numPoses{32}, numCircles{4};
     102             :         auto sinoDescriptor = SphereTrajectoryGenerator::createTrajectory(
     103             :             numPoses, volumeDescriptor, numCircles,
     104           3 :             geometry::SourceToCenterOfRotation(static_cast<real_t>(size(0)) * 100.0f),
     105           9 :             geometry::CenterOfRotationToDetector(static_cast<real_t>(size(0))));
     106           6 :         const auto numCoeffsPerDim = sinoDescriptor->getNumberOfCoefficientsPerDimension();
     107           3 :         const index_t nDimensions = sinoDescriptor->getNumberOfDimensions();
     108           3 :         const auto numElements = numCoeffsPerDim[nDimensions - 1];
     109             : 
     110           4 :         WHEN("performing round robin sampling into 4 subsets")
     111             :         {
     112           1 :             const auto nSubsets = 4;
     113           2 :             std::vector<index_t> indices(static_cast<std::size_t>(numElements));
     114           1 :             std::iota(indices.begin(), indices.end(), 0);
     115             :             const auto mapping =
     116           2 :                 SubsetSampler<PlanarDetectorDescriptor, real_t>::splitRoundRobin(indices, nSubsets);
     117           2 :             THEN("The mapping is correct with every subset having 8 elements")
     118             :             {
     119           5 :                 for (std::size_t i = 0; i < static_cast<std::size_t>(nSubsets); ++i) {
     120           4 :                     REQUIRE_EQ(mapping[i].size(), 8);
     121          36 :                     for (std::size_t j = 0; j < mapping[i].size(); j++) {
     122          32 :                         REQUIRE_EQ(mapping[i][j], j * nSubsets + i);
     123             :                     }
     124             :                 }
     125             :             }
     126             :         }
     127           4 :         WHEN("performing round robin sampling into 5 subsets")
     128             :         {
     129           1 :             const auto nSubsets = 5;
     130           2 :             std::vector<index_t> indices(static_cast<std::size_t>(numElements));
     131           1 :             std::iota(indices.begin(), indices.end(), 0);
     132             :             const auto mapping =
     133           2 :                 SubsetSampler<PlanarDetectorDescriptor, real_t>::splitRoundRobin(indices, nSubsets);
     134           2 :             THEN("The mapping is correct with every subset having 6 elements apart from the first "
     135             :                  "two")
     136             :             {
     137           6 :                 for (std::size_t i = 0; i < static_cast<std::size_t>(nSubsets); ++i) {
     138           5 :                     if (i <= 1) {
     139           2 :                         REQUIRE_EQ(mapping[i].size(), 7);
     140             :                     } else {
     141           3 :                         REQUIRE_EQ(mapping[i].size(), 6);
     142             :                     }
     143          37 :                     for (std::size_t j = 0; j < mapping[i].size(); j++) {
     144          32 :                         REQUIRE_EQ(mapping[i][j], j * nSubsets + i);
     145             :                     }
     146             :                 }
     147             :             }
     148             :         }
     149           4 :         WHEN("performing equi rotation sampling into 4 subsets")
     150             :         {
     151           1 :             const auto nSubsets = 4;
     152             :             const auto mapping =
     153             :                 SubsetSampler<PlanarDetectorDescriptor, real_t>::splitRotationalClustering(
     154           2 :                     static_cast<const PlanarDetectorDescriptor&>(*sinoDescriptor), nSubsets);
     155           2 :             THEN("The mapping is correct with every subset having 8 elements")
     156             :             {
     157           1 :                 REQUIRE_EQ(mapping[0], std::vector<index_t>{0, 12, 1, 22, 30, 18, 25, 26});
     158           1 :                 REQUIRE_EQ(mapping[1], std::vector<index_t>{4, 13, 21, 31, 23, 15, 28, 27});
     159           1 :                 REQUIRE_EQ(mapping[2], std::vector<index_t>{5, 20, 10, 14, 24, 7, 17, 8});
     160           1 :                 REQUIRE_EQ(mapping[3], std::vector<index_t>{11, 3, 6, 19, 29, 9, 16, 2});
     161             :             }
     162             :         }
     163             :     }
     164           6 : }
     165             : 
     166           6 : TEST_CASE(
     167             :     "SubsetSampler: 2D Integration test with PlanarDetectorDescriptor and circular trajectory")
     168             : {
     169           6 :     Logger::setLevel(Logger::LogLevel::WARN);
     170             : 
     171          12 :     IndexVector_t size(2);
     172           6 :     size << 16, 16;
     173          12 :     auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
     174           6 :     auto& volumeDescriptor = phantom.getDataDescriptor();
     175             : 
     176           6 :     index_t numAngles{20}, arc{360};
     177             :     auto sinoDescriptor = CircleTrajectoryGenerator::createTrajectory(
     178          12 :         numAngles, phantom.getDataDescriptor(), arc, static_cast<real_t>(size(0)) * 100.0f,
     179          12 :         static_cast<real_t>(size(0)));
     180             : 
     181             :     SiddonsMethod projector(static_cast<const VolumeDescriptor&>(volumeDescriptor),
     182          12 :                             *sinoDescriptor);
     183             : 
     184          12 :     auto sinogram = projector.apply(phantom);
     185             :     const auto coeffsPerDimSinogram =
     186          12 :         sinogram.getDataDescriptor().getNumberOfCoefficientsPerDimension();
     187             : 
     188          12 :     GIVEN("A small phantom problem")
     189             :     {
     190           9 :         WHEN("Setting up a subset sampler with 4 subsets")
     191             :         {
     192           3 :             index_t nSubsets{4};
     193             :             SubsetSampler<PlanarDetectorDescriptor, real_t> subsetSampler(
     194             :                 static_cast<const VolumeDescriptor&>(volumeDescriptor),
     195           6 :                 static_cast<const PlanarDetectorDescriptor&>(*sinoDescriptor), nSubsets);
     196           4 :             THEN("the clone works as expected")
     197             :             {
     198             : 
     199           2 :                 auto subsetSamplerClone = subsetSampler.clone();
     200             : 
     201           1 :                 REQUIRE_NE(subsetSamplerClone.get(), &subsetSampler);
     202           1 :                 REQUIRE_EQ(*subsetSamplerClone, subsetSampler);
     203             :             }
     204           4 :             AND_THEN("The full data has the correct dimensions")
     205             :             {
     206           1 :                 const auto data = subsetSampler.getPartitionedData(sinogram);
     207           1 :                 REQUIRE_EQ(data.getDataDescriptor().getNumberOfDimensions(),
     208             :                            sinogram.getDataDescriptor().getNumberOfDimensions());
     209           1 :                 REQUIRE_EQ(data.getDataDescriptor().getNumberOfCoefficientsPerDimension(),
     210             :                            coeffsPerDimSinogram);
     211             :             }
     212             : 
     213           4 :             AND_THEN("It returns the correct data blocks and projectors for each subsets")
     214             :             {
     215           1 :                 const auto data = subsetSampler.getPartitionedData(sinogram);
     216           5 :                 for (index_t i = 0; i < nSubsets; i++) {
     217           8 :                     const auto subset = data.getBlock(i);
     218             :                     const auto coeffsPerDimension =
     219           8 :                         subset.getDataDescriptor().getNumberOfCoefficientsPerDimension();
     220             : 
     221           4 :                     REQUIRE_EQ(coeffsPerDimension[0], coeffsPerDimSinogram[0]);
     222             :                 }
     223             : 
     224           1 :                 subsetSampler.getProjector<SiddonsMethod<real_t>>();
     225           1 :                 REQUIRE_EQ(subsetSampler.getSubsetProjectors<SiddonsMethod<real_t>>().size(),
     226             :                            nSubsets);
     227             :             }
     228             :         }
     229           9 :         WHEN("Setting up a subset sampler with 4 subsets and ROTATIONAL_CLUSTERING sampling")
     230             :         {
     231           3 :             index_t nSubsets{4};
     232             :             SubsetSampler<PlanarDetectorDescriptor, real_t> subsetSampler(
     233             :                 static_cast<const VolumeDescriptor&>(volumeDescriptor),
     234           3 :                 static_cast<const PlanarDetectorDescriptor&>(*sinoDescriptor), nSubsets,
     235             :                 SubsetSampler<PlanarDetectorDescriptor,
     236           6 :                               real_t>::SamplingStrategy::ROTATIONAL_CLUSTERING);
     237           4 :             THEN("the clone works as expected")
     238             :             {
     239             : 
     240           2 :                 auto subsetSamplerClone = subsetSampler.clone();
     241             : 
     242           1 :                 REQUIRE_NE(subsetSamplerClone.get(), &subsetSampler);
     243           1 :                 REQUIRE_EQ(*subsetSamplerClone, subsetSampler);
     244             :             }
     245           4 :             AND_THEN("The full data has the correct dimensions")
     246             :             {
     247           1 :                 const auto data = subsetSampler.getPartitionedData(sinogram);
     248           1 :                 REQUIRE_EQ(data.getDataDescriptor().getNumberOfDimensions(),
     249             :                            sinogram.getDataDescriptor().getNumberOfDimensions());
     250           1 :                 REQUIRE_EQ(data.getDataDescriptor().getNumberOfCoefficientsPerDimension(),
     251             :                            coeffsPerDimSinogram);
     252             :             }
     253             : 
     254           4 :             AND_THEN("It returns the correct data blocks and projectors for each subsets")
     255             :             {
     256           1 :                 const auto data = subsetSampler.getPartitionedData(sinogram);
     257           5 :                 for (index_t i = 0; i < nSubsets; i++) {
     258           8 :                     const auto subset = data.getBlock(i);
     259             :                     const auto coeffsPerDimension =
     260           8 :                         subset.getDataDescriptor().getNumberOfCoefficientsPerDimension();
     261             : 
     262           4 :                     REQUIRE_EQ(coeffsPerDimension[0], coeffsPerDimSinogram[0]);
     263             :                 }
     264             : 
     265           1 :                 subsetSampler.getProjector<SiddonsMethod<real_t>>();
     266           1 :                 REQUIRE_EQ(subsetSampler.getSubsetProjectors<SiddonsMethod<real_t>>().size(),
     267             :                            nSubsets);
     268             :             }
     269             :         }
     270             :     }
     271           6 : }
     272             : 
     273           6 : TEST_CASE(
     274             :     "SubsetSampler: 3D Integration test with PlanarDetectorDescriptor and spherical trajectory")
     275             : {
     276           6 :     Logger::setLevel(Logger::LogLevel::WARN);
     277             : 
     278          12 :     IndexVector_t size(3);
     279           6 :     size << 8, 8, 8;
     280          12 :     auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
     281           6 :     auto& volumeDescriptor = phantom.getDataDescriptor();
     282             : 
     283           6 :     index_t numPoses{16}, numCircles{5};
     284             :     auto sinoDescriptor = SphereTrajectoryGenerator::createTrajectory(
     285             :         numPoses, phantom.getDataDescriptor(), numCircles,
     286           6 :         geometry::SourceToCenterOfRotation(static_cast<real_t>(size(0)) * 100.0f),
     287          18 :         geometry::CenterOfRotationToDetector(static_cast<real_t>(size(0))));
     288             : 
     289             :     SiddonsMethod projector(static_cast<const VolumeDescriptor&>(volumeDescriptor),
     290          12 :                             *sinoDescriptor);
     291             : 
     292          12 :     auto sinogram = projector.apply(phantom);
     293             :     const auto coeffsPerDimSinogram =
     294          12 :         sinogram.getDataDescriptor().getNumberOfCoefficientsPerDimension();
     295             : 
     296          12 :     GIVEN("A small phantom problem")
     297             :     {
     298             : 
     299           9 :         WHEN("Setting up a subset sampler with 4 subsets")
     300             :         {
     301           3 :             index_t nSubsets{4};
     302             :             SubsetSampler<PlanarDetectorDescriptor, real_t> subsetSampler(
     303             :                 static_cast<const VolumeDescriptor&>(volumeDescriptor),
     304           6 :                 static_cast<const PlanarDetectorDescriptor&>(*sinoDescriptor), nSubsets);
     305           4 :             THEN("the clone works as expected")
     306             :             {
     307             : 
     308           2 :                 auto subsetSamplerClone = subsetSampler.clone();
     309             : 
     310           1 :                 REQUIRE_NE(subsetSamplerClone.get(), &subsetSampler);
     311           1 :                 REQUIRE_EQ(*subsetSamplerClone, subsetSampler);
     312             :             }
     313           4 :             AND_THEN("The full data has the correct dimensions")
     314             :             {
     315           1 :                 const auto data = subsetSampler.getPartitionedData(sinogram);
     316           1 :                 REQUIRE_EQ(data.getDataDescriptor().getNumberOfDimensions(),
     317             :                            sinogram.getDataDescriptor().getNumberOfDimensions());
     318           1 :                 REQUIRE_EQ(data.getDataDescriptor().getNumberOfCoefficientsPerDimension(),
     319             :                            coeffsPerDimSinogram);
     320             :             }
     321             : 
     322           4 :             AND_THEN("It returns the correct data blocks and projectors for each subsets")
     323             :             {
     324           1 :                 const auto data = subsetSampler.getPartitionedData(sinogram);
     325           5 :                 for (index_t i = 0; i < nSubsets; i++) {
     326           8 :                     const auto subset = data.getBlock(i);
     327             :                     const auto coeffsPerDimension =
     328           8 :                         subset.getDataDescriptor().getNumberOfCoefficientsPerDimension();
     329             : 
     330           4 :                     REQUIRE_EQ(coeffsPerDimension[0], coeffsPerDimSinogram[0]);
     331             :                 }
     332             : 
     333           1 :                 subsetSampler.getProjector<SiddonsMethod<real_t>>();
     334           1 :                 REQUIRE_EQ(subsetSampler.getSubsetProjectors<SiddonsMethod<real_t>>().size(),
     335             :                            nSubsets);
     336             :             }
     337             :         }
     338           9 :         WHEN("Setting up a subset sampler with 4 subsets and ROTATIONAL_CLUSTERING sampling")
     339             :         {
     340           3 :             index_t nSubsets{4};
     341             :             SubsetSampler<PlanarDetectorDescriptor, real_t> subsetSampler(
     342             :                 static_cast<const VolumeDescriptor&>(volumeDescriptor),
     343           3 :                 static_cast<const PlanarDetectorDescriptor&>(*sinoDescriptor), nSubsets,
     344             :                 SubsetSampler<PlanarDetectorDescriptor,
     345           6 :                               real_t>::SamplingStrategy::ROTATIONAL_CLUSTERING);
     346           4 :             THEN("the clone works as expected")
     347             :             {
     348             : 
     349           2 :                 auto subsetSamplerClone = subsetSampler.clone();
     350             : 
     351           1 :                 REQUIRE_NE(subsetSamplerClone.get(), &subsetSampler);
     352           1 :                 REQUIRE_EQ(*subsetSamplerClone, subsetSampler);
     353             :             }
     354           4 :             AND_THEN("The full data has the correct dimensions")
     355             :             {
     356           1 :                 const auto data = subsetSampler.getPartitionedData(sinogram);
     357           1 :                 REQUIRE_EQ(data.getDataDescriptor().getNumberOfDimensions(),
     358             :                            sinogram.getDataDescriptor().getNumberOfDimensions());
     359           1 :                 REQUIRE_EQ(data.getDataDescriptor().getNumberOfCoefficientsPerDimension(),
     360             :                            coeffsPerDimSinogram);
     361             :             }
     362             : 
     363           4 :             AND_THEN("It returns the correct data blocks and projectors for each subsets")
     364             :             {
     365           1 :                 const auto data = subsetSampler.getPartitionedData(sinogram);
     366           5 :                 for (index_t i = 0; i < nSubsets; i++) {
     367           8 :                     const auto subset = data.getBlock(i);
     368             :                     const auto coeffsPerDimension =
     369           8 :                         subset.getDataDescriptor().getNumberOfCoefficientsPerDimension();
     370             : 
     371           4 :                     REQUIRE_EQ(coeffsPerDimension[0], coeffsPerDimSinogram[0]);
     372             :                 }
     373             : 
     374           1 :                 subsetSampler.getProjector<SiddonsMethod<real_t>>();
     375           1 :                 REQUIRE_EQ(subsetSampler.getSubsetProjectors<SiddonsMethod<real_t>>().size(),
     376             :                            nSubsets);
     377             :             }
     378             :         }
     379             :     }
     380           6 : }