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