LCOV - code coverage report
Current view: top level - solvers/tests - test_SQS.cpp (source / functions) Hit Total Coverage
Test: test_coverage.info.cleaned Lines: 150 150 100.0 %
Date: 2022-02-28 03:37:41 Functions: 14 14 100.0 %

          Line data    Source code
       1             : /**
       2             :  * @file test_SQS.cpp
       3             :  *
       4             :  * @brief Tests for the SQS class
       5             :  *
       6             :  * @author Michael Loipführer - initial code
       7             :  */
       8             : 
       9             : #include "doctest/doctest.h"
      10             : 
      11             : #include <iostream>
      12             : #include "JacobiPreconditioner.h"
      13             : #include "SQS.h"
      14             : #include "WLSProblem.h"
      15             : #include "WLSSubsetProblem.h"
      16             : #include "SubsetSampler.h"
      17             : #include "PlanarDetectorDescriptor.h"
      18             : #include "Identity.h"
      19             : #include "Logger.h"
      20             : #include "VolumeDescriptor.h"
      21             : #include "CircleTrajectoryGenerator.h"
      22             : #include "SiddonsMethod.h"
      23             : #include "PhantomGenerator.h"
      24             : #include "TypeCasts.hpp"
      25             : #include "testHelpers.h"
      26             : 
      27             : using namespace elsa;
      28             : 
      29             : template <template <typename> typename T, typename data_t>
      30             : constexpr data_t return_data_t(const T<data_t>&);
      31             : 
      32             : TEST_SUITE_BEGIN("solvers");
      33             : 
      34          12 : TYPE_TO_STRING(SQS<float>);
      35          12 : TYPE_TO_STRING(SQS<double>);
      36             : 
      37          22 : TEST_CASE_TEMPLATE("SQS: Solving a simple linear problem", TestType, SQS<float>, SQS<double>)
      38             : {
      39             :     // Set seed for Eigen Matrices!
      40           4 :     srand((unsigned int) 666);
      41             : 
      42             :     using data_t = decltype(return_data_t(std::declval<TestType>()));
      43             :     // eliminate the timing info from console for the tests
      44           4 :     Logger::setLevel(Logger::LogLevel::OFF);
      45             : 
      46           8 :     GIVEN("a linear problem")
      47             :     {
      48           8 :         IndexVector_t numCoeff(2);
      49           4 :         numCoeff << 13, 24;
      50           8 :         VolumeDescriptor dd{numCoeff};
      51             : 
      52           8 :         Eigen::Matrix<data_t, -1, 1> bVec(dd.getNumberOfCoefficients());
      53           4 :         bVec.setRandom();
      54           8 :         DataContainer<data_t> dcB{dd, bVec};
      55             : 
      56           4 :         bVec.setRandom();
      57           4 :         bVec = bVec.cwiseAbs();
      58           8 :         Scaling<data_t> scalingOp{dd, DataContainer<data_t>{dd, bVec}};
      59             : 
      60             :         // using WLS problem here for ease of use
      61           8 :         WLSProblem prob{scalingOp, dcB};
      62             : 
      63           4 :         data_t epsilon = std::numeric_limits<data_t>::epsilon();
      64             : 
      65           6 :         WHEN("setting up a SQS solver")
      66             :         {
      67           4 :             TestType solver{prob, true, epsilon};
      68             : 
      69           4 :             THEN("the clone works correctly")
      70             :             {
      71           4 :                 auto sqsClone = solver.clone();
      72             : 
      73           2 :                 REQUIRE_NE(sqsClone.get(), &solver);
      74           2 :                 REQUIRE_EQ(*sqsClone, solver);
      75             : 
      76           4 :                 AND_THEN("it works as expected")
      77             :                 {
      78           4 :                     auto solution = solver.solve(1000);
      79             : 
      80           2 :                     DataContainer<data_t> resultsDifference = scalingOp.apply(solution) - dcB;
      81             : 
      82             :                     // should have converged for the given number of iterations
      83           2 :                     REQUIRE_UNARY(checkApproxEq(resultsDifference.squaredL2Norm(),
      84             :                                                 epsilon * epsilon * dcB.squaredL2Norm()));
      85             :                 }
      86             :             }
      87             :         }
      88             : 
      89           6 :         WHEN("setting up a preconditioned SQS solver")
      90             :         {
      91           4 :             auto preconditioner = JacobiPreconditioner<data_t>(scalingOp, false);
      92           4 :             TestType solver{prob, preconditioner, true, epsilon};
      93             : 
      94           4 :             THEN("the clone works correctly")
      95             :             {
      96           4 :                 auto sqsClone = solver.clone();
      97             : 
      98           2 :                 REQUIRE_NE(sqsClone.get(), &solver);
      99           2 :                 REQUIRE_EQ(*sqsClone, solver);
     100             : 
     101           4 :                 AND_THEN("it works as expected")
     102             :                 {
     103             :                     // with a good preconditioner we should need fewer iterations than without
     104           4 :                     auto solution = solver.solve(500);
     105             : 
     106           2 :                     DataContainer<data_t> resultsDifference = scalingOp.apply(solution) - dcB;
     107             : 
     108             :                     // should have converged for the given number of iterations
     109           2 :                     REQUIRE_UNARY(checkApproxEq(resultsDifference.squaredL2Norm(),
     110             :                                                 epsilon * epsilon * dcB.squaredL2Norm()));
     111             :                 }
     112             :             }
     113             :         }
     114             :     }
     115           4 : }
     116             : 
     117          22 : TEST_CASE_TEMPLATE("SQS: Solving a Tikhonov problem", TestType, SQS<float>, SQS<double>)
     118             : {
     119             :     // Set seed for Eigen Matrices!
     120           4 :     srand((unsigned int) 666);
     121             : 
     122             :     using data_t = decltype(return_data_t(std::declval<TestType>()));
     123             :     // eliminate the timing info from console for the tests
     124           4 :     Logger::setLevel(Logger::LogLevel::OFF);
     125             : 
     126           8 :     GIVEN("a Tikhonov problem")
     127             :     {
     128           8 :         IndexVector_t numCoeff(2);
     129           4 :         numCoeff << 13, 24;
     130           8 :         VolumeDescriptor dd(numCoeff);
     131             : 
     132           8 :         Eigen::Matrix<data_t, -1, 1> bVec(dd.getNumberOfCoefficients());
     133           4 :         bVec.setRandom();
     134           8 :         DataContainer dcB(dd, bVec);
     135             : 
     136           4 :         bVec.setRandom();
     137           4 :         bVec = bVec.cwiseProduct(bVec);
     138           8 :         Scaling<data_t> scalingOp{dd, DataContainer<data_t>{dd, bVec}};
     139             : 
     140           4 :         auto lambda = static_cast<data_t>(0.1);
     141           8 :         Scaling<data_t> lambdaOp{dd, lambda};
     142             : 
     143             :         // using WLS problem here for ease of use
     144           8 :         WLSProblem<data_t> prob{scalingOp + lambdaOp, dcB};
     145             : 
     146           4 :         data_t epsilon = std::numeric_limits<data_t>::epsilon();
     147             : 
     148           6 :         WHEN("setting up a SQS solver")
     149             :         {
     150           4 :             TestType solver{prob, true, epsilon};
     151             : 
     152           4 :             THEN("the clone works correctly")
     153             :             {
     154           4 :                 auto sqsClone = solver.clone();
     155             : 
     156           2 :                 REQUIRE_NE(sqsClone.get(), &solver);
     157           2 :                 REQUIRE_EQ(*sqsClone, solver);
     158             : 
     159           4 :                 AND_THEN("it works as expected")
     160             :                 {
     161           4 :                     auto solution = solver.solve(dd.getNumberOfCoefficients());
     162             : 
     163           2 :                     DataContainer<data_t> resultsDifference =
     164             :                         (scalingOp + lambdaOp).apply(solution) - dcB;
     165             : 
     166             :                     // should have converged for the given number of iterations
     167             :                     // does not converge to the optimal solution because of the regularization term
     168           2 :                     REQUIRE_UNARY(checkApproxEq(resultsDifference.squaredL2Norm(),
     169             :                                                 epsilon * epsilon * dcB.squaredL2Norm()));
     170             :                 }
     171             :             }
     172             :         }
     173             : 
     174           6 :         WHEN("setting up a preconditioned SQS solver")
     175             :         {
     176           4 :             auto preconditioner = JacobiPreconditioner<data_t>(scalingOp + lambdaOp, false);
     177           4 :             TestType solver{prob, preconditioner, true, epsilon};
     178             : 
     179           4 :             THEN("the clone works correctly")
     180             :             {
     181           4 :                 auto sqsClone = solver.clone();
     182             : 
     183           2 :                 REQUIRE_NE(sqsClone.get(), &solver);
     184           2 :                 REQUIRE_EQ(*sqsClone, solver);
     185             : 
     186           4 :                 AND_THEN("it works as expected")
     187             :                 {
     188             :                     // a perfect preconditioner should allow for convergence in a single step
     189           4 :                     auto solution = solver.solve(dd.getNumberOfCoefficients());
     190             : 
     191           2 :                     DataContainer<data_t> resultsDifference =
     192             :                         (scalingOp + lambdaOp).apply(solution) - dcB;
     193             : 
     194             :                     // should have converged for the given number of iterations
     195           2 :                     REQUIRE_UNARY(checkApproxEq(resultsDifference.squaredL2Norm(),
     196             :                                                 epsilon * epsilon * dcB.squaredL2Norm()));
     197             :                 }
     198             :             }
     199             :         }
     200             :     }
     201           4 : }
     202             : 
     203           1 : TEST_CASE("SQS: Solving a simple phantom reconstruction")
     204             : {
     205             :     // Set seed for Eigen Matrices!
     206           1 :     srand((unsigned int) 666);
     207             : 
     208             :     // eliminate the timing info from console for the tests
     209           1 :     Logger::setLevel(Logger::LogLevel::OFF);
     210             : 
     211           2 :     GIVEN("a Phantom reconstruction problem")
     212             :     {
     213           2 :         IndexVector_t size(2);
     214           1 :         size << 16, 16; // TODO: determine optimal phantom size for efficient testing
     215           2 :         auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
     216           1 :         auto& volumeDescriptor = phantom.getDataDescriptor();
     217             : 
     218           1 :         index_t numAngles{90}, arc{360};
     219             :         auto sinoDescriptor = CircleTrajectoryGenerator::createTrajectory(
     220           2 :             numAngles, phantom.getDataDescriptor(), arc, static_cast<real_t>(size(0)) * 100.0f,
     221           2 :             static_cast<real_t>(size(0)));
     222             : 
     223           2 :         SiddonsMethod projector(downcast<VolumeDescriptor>(volumeDescriptor), *sinoDescriptor);
     224             : 
     225           2 :         auto sinogram = projector.apply(phantom);
     226             : 
     227           2 :         WLSProblem problem(projector, sinogram);
     228           1 :         real_t epsilon = std::numeric_limits<real_t>::epsilon();
     229             : 
     230           2 :         WHEN("setting up a SQS solver")
     231             :         {
     232           2 :             SQS solver{problem, true, epsilon};
     233             : 
     234           2 :             THEN("the clone works correctly")
     235             :             {
     236           2 :                 auto sqsClone = solver.clone();
     237             : 
     238           1 :                 REQUIRE_NE(sqsClone.get(), &solver);
     239           1 :                 REQUIRE_EQ(*sqsClone, solver);
     240             : 
     241           2 :                 AND_THEN("it works as expected")
     242             :                 {
     243           2 :                     auto reconstruction = solver.solve(40);
     244             : 
     245           1 :                     DataContainer resultsDifference = reconstruction - phantom;
     246             : 
     247             :                     // should have converged for the given number of iterations
     248             :                     // does not converge to the optimal solution because of the regularization term
     249           1 :                     REQUIRE_UNARY(checkApproxEq(resultsDifference.squaredL2Norm(), 0.034f));
     250             :                 }
     251             :             }
     252             :         }
     253             :     }
     254           1 : }
     255             : 
     256           2 : TEST_CASE("SQS: Solving a simple phantom problem using ordered subsets")
     257             : {
     258             :     // Set seed for Eigen Matrices!
     259           2 :     srand((unsigned int) 666);
     260             : 
     261             :     // eliminate the timing info from console for the tests
     262           2 :     Logger::setLevel(Logger::LogLevel::OFF);
     263             : 
     264           4 :     GIVEN("a Phantom reconstruction problem")
     265             :     {
     266             : 
     267           4 :         IndexVector_t size(2);
     268           2 :         size << 16, 16; // TODO: determine optimal phantom size for efficient testing
     269           4 :         auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
     270           2 :         auto& volumeDescriptor = phantom.getDataDescriptor();
     271             : 
     272           2 :         index_t numAngles{20}, arc{180};
     273             :         auto sinoDescriptor = CircleTrajectoryGenerator::createTrajectory(
     274           4 :             numAngles, phantom.getDataDescriptor(), arc, static_cast<real_t>(size(0)) * 100.0f,
     275           4 :             static_cast<real_t>(size(0)));
     276             : 
     277             :         SiddonsMethod projector(static_cast<const VolumeDescriptor&>(volumeDescriptor),
     278           4 :                                 *sinoDescriptor);
     279             : 
     280           4 :         auto sinogram = projector.apply(phantom);
     281             : 
     282           2 :         real_t epsilon = std::numeric_limits<real_t>::epsilon();
     283             : 
     284           3 :         WHEN("setting up a SQS solver with ROUND_ROBIN subsampling")
     285             :         {
     286           1 :             index_t nSubsets{4};
     287             :             SubsetSampler<PlanarDetectorDescriptor, real_t> subsetSampler(
     288             :                 static_cast<const VolumeDescriptor&>(volumeDescriptor),
     289           2 :                 static_cast<const PlanarDetectorDescriptor&>(*sinoDescriptor), nSubsets);
     290             : 
     291             :             WLSSubsetProblem<real_t> problem(
     292           2 :                 *subsetSampler.getProjector<SiddonsMethod<real_t>>(),
     293           2 :                 subsetSampler.getPartitionedData(sinogram),
     294           3 :                 subsetSampler.getSubsetProjectors<SiddonsMethod<real_t>>());
     295           2 :             SQS solver{problem, true, epsilon};
     296             : 
     297           2 :             THEN("the clone works correctly")
     298             :             {
     299           2 :                 auto sqsClone = solver.clone();
     300             : 
     301           1 :                 REQUIRE(sqsClone.get() != &solver);
     302           1 :                 REQUIRE(*sqsClone == solver);
     303             : 
     304           2 :                 AND_THEN("it works as expected")
     305             :                 {
     306           2 :                     auto reconstruction = solver.solve(10);
     307             : 
     308           1 :                     DataContainer resultsDifference = reconstruction - phantom;
     309             : 
     310             :                     // should have converged for the given number of iterations
     311             :                     // does not converge to the optimal solution because of the regularization term
     312           1 :                     REQUIRE_UNARY(checkApproxEq(resultsDifference.squaredL2Norm(),
     313             :                                                 epsilon * epsilon * phantom.squaredL2Norm(), 0.1));
     314             :                 }
     315             :             }
     316             :         }
     317           3 :         WHEN("setting up a SQS solver with ROTATIONAL_CLUSTERING subsampling")
     318             :         {
     319           1 :             index_t nSubsets{4};
     320             :             SubsetSampler<PlanarDetectorDescriptor, real_t> subsetSampler(
     321             :                 static_cast<const VolumeDescriptor&>(volumeDescriptor),
     322           1 :                 static_cast<const PlanarDetectorDescriptor&>(*sinoDescriptor), nSubsets,
     323             :                 SubsetSampler<PlanarDetectorDescriptor,
     324           2 :                               real_t>::SamplingStrategy::ROTATIONAL_CLUSTERING);
     325             : 
     326             :             WLSSubsetProblem<real_t> problem(
     327           2 :                 *subsetSampler.getProjector<SiddonsMethod<real_t>>(),
     328           2 :                 subsetSampler.getPartitionedData(sinogram),
     329           3 :                 subsetSampler.getSubsetProjectors<SiddonsMethod<real_t>>());
     330           2 :             SQS solver{problem, true, epsilon};
     331             : 
     332           2 :             THEN("the clone works correctly")
     333             :             {
     334           2 :                 auto sqsClone = solver.clone();
     335             : 
     336           1 :                 REQUIRE_NE(sqsClone.get(), &solver);
     337           1 :                 REQUIRE_EQ(*sqsClone, solver);
     338             : 
     339           2 :                 AND_THEN("it works as expected")
     340             :                 {
     341           2 :                     auto reconstruction = solver.solve(10);
     342             : 
     343           1 :                     DataContainer resultsDifference = reconstruction - phantom;
     344             : 
     345             :                     // should have converged for the given number of iterations
     346           1 :                     REQUIRE_UNARY(checkApproxEq(resultsDifference.squaredL2Norm(),
     347             :                                                 epsilon * epsilon * phantom.squaredL2Norm(), 0.1));
     348             :                 }
     349             :             }
     350             :         }
     351             :     }
     352           2 : }
     353             : 
     354             : TEST_SUITE_END();

Generated by: LCOV version 1.15