Line data    Source code

       1             : /**
       2             :  * @file test_FGM.cpp
       3             :  *
       4             :  * @brief Tests for the Fast Gradient Method class
       5             :  *
       6             :  * @author Michael Loipführer - initial code
       7             :  */
       8             : 
       9             : #include "doctest/doctest.h"
      10             : 
      11             : #include <iostream>
      12             : #include "FGM.h"
      13             : #include "WLSProblem.h"
      14             : #include "Identity.h"
      15             : #include "Logger.h"
      16             : #include "VolumeDescriptor.h"
      17             : #include "SiddonsMethod.h"
      18             : #include "CircleTrajectoryGenerator.h"
      19             : #include "PhantomGenerator.h"
      20             : #include "JacobiPreconditioner.h"
      21             : #include "TypeCasts.hpp"
      22             : #include "testHelpers.h"
      23             : 
      24             : using namespace elsa;
      25             : using namespace doctest;
      26             : 
      27             : TEST_SUITE_BEGIN("solvers");
      28             : 
      29             : template <template <typename> typename T, typename data_t>
      30             : constexpr data_t return_data_t(const T<data_t>&);
      31             : 
      32             : TYPE_TO_STRING(FGM<float>);
      33             : TYPE_TO_STRING(FGM<double>);
      34             : 
      35             : TEST_CASE_TEMPLATE("FGM: Solving a simple linear problem", TestType, FGM<float>, FGM<double>)
      36           4 : {
      37             :     // Set seed for Eigen Matrices!
      38           4 :     srand((unsigned int) 666);
      39             : 
      40           4 :     using data_t = decltype(return_data_t(std::declval<TestType>()));
      41             :     // eliminate the timing info from console for the tests
      42           4 :     Logger::setLevel(Logger::LogLevel::OFF);
      43             : 
      44           4 :     GIVEN("a linear problem")
      45           4 :     {
      46           4 :         IndexVector_t numCoeff(2);
      47           4 :         numCoeff << 13, 24;
      48           4 :         VolumeDescriptor dd{numCoeff};
      49             : 
      50           4 :         Eigen::Matrix<data_t, -1, 1> bVec(dd.getNumberOfCoefficients());
      51           4 :         bVec.setRandom();
      52           4 :         DataContainer<data_t> dcB{dd, bVec};
      53             : 
      54           4 :         bVec.setRandom();
      55           4 :         bVec = bVec.cwiseAbs();
      56           4 :         Scaling<data_t> scalingOp{dd, DataContainer<data_t>{dd, bVec}};
      57             : 
      58             :         // using WLS problem here for ease of use
      59           4 :         WLSProblem prob{scalingOp, dcB};
      60             : 
      61           4 :         data_t epsilon = std::numeric_limits<data_t>::epsilon();
      62             : 
      63           4 :         WHEN("setting up a FGM solver")
      64           4 :         {
      65           2 :             TestType solver{prob, epsilon};
      66             : 
      67           2 :             THEN("the clone works correctly")
      68           2 :             {
      69           2 :                 auto fgmClone = solver.clone();
      70             : 
      71           2 :                 REQUIRE_NE(fgmClone.get(), &solver);
      72           2 :                 REQUIRE_EQ(*fgmClone, solver);
      73             : 
      74           2 :                 AND_THEN("it works as expected")
      75           2 :                 {
      76           2 :                     auto solution = solver.solve(1000);
      77             : 
      78           2 :                     DataContainer<data_t> resultsDifference = scalingOp.apply(solution) - dcB;
      79             : 
      80             :                     // should have converged for the given number of iterations
      81           2 :                     REQUIRE_UNARY(checkApproxEq(resultsDifference.squaredL2Norm(),
      82           2 :                                                 epsilon * epsilon * dcB.squaredL2Norm(), 0.5f));
      83           2 :                 }
      84           2 :             }
      85           2 :         }
      86             : 
      87           4 :         WHEN("setting up a preconditioned FGM solver")
      88           4 :         {
      89           2 :             auto preconditionerInverse = JacobiPreconditioner<data_t>(scalingOp, true);
      90           2 :             TestType solver{prob, preconditionerInverse, epsilon};
      91             : 
      92           2 :             THEN("the clone works correctly")
      93           2 :             {
      94           2 :                 auto fgmClone = solver.clone();
      95             : 
      96           2 :                 REQUIRE_NE(fgmClone.get(), &solver);
      97           2 :                 REQUIRE_EQ(*fgmClone, solver);
      98             : 
      99           2 :                 AND_THEN("it works as expected")
     100           2 :                 {
     101             :                     // with a good preconditioner we should need fewer iterations than without
     102           2 :                     auto solution = solver.solve(1000);
     103             : 
     104           2 :                     DataContainer<data_t> resultsDifference = scalingOp.apply(solution) - dcB;
     105             : 
     106             :                     // should have converged for the given number of iterations
     107           2 :                     REQUIRE_UNARY(checkApproxEq(resultsDifference.squaredL2Norm(),
     108           2 :                                                 epsilon * epsilon * dcB.squaredL2Norm(), 0.1f));
     109           2 :                 }
     110           2 :             }
     111           2 :         }
     112           4 :     }
     113           4 : }
     114             : 
     115             : TEST_CASE_TEMPLATE("FGM: Solving a Tikhonov problem", TestType, FGM<float>, FGM<double>)
     116           4 : {
     117             :     // Set seed for Eigen Matrices!
     118           4 :     srand((unsigned int) 666);
     119             : 
     120           4 :     using data_t = decltype(return_data_t(std::declval<TestType>()));
     121             :     // eliminate the timing info from console for the tests
     122           4 :     Logger::setLevel(Logger::LogLevel::OFF);
     123             : 
     124           4 :     GIVEN("a Tikhonov problem")
     125           4 :     {
     126           4 :         IndexVector_t numCoeff(2);
     127           4 :         numCoeff << 13, 24;
     128           4 :         VolumeDescriptor dd(numCoeff);
     129             : 
     130           4 :         Eigen::Matrix<data_t, -1, 1> bVec(dd.getNumberOfCoefficients());
     131           4 :         bVec.setRandom();
     132           4 :         DataContainer dcB(dd, bVec);
     133             : 
     134             :         // the regularization term
     135             : 
     136           4 :         bVec.setRandom();
     137           4 :         bVec = bVec.cwiseProduct(bVec);
     138           4 :         Scaling<data_t> scalingOp{dd, DataContainer<data_t>{dd, bVec}};
     139             : 
     140           4 :         auto lambda = static_cast<data_t>(0.1);
     141           4 :         Scaling<data_t> lambdaOp{dd, lambda};
     142             : 
     143             :         // using WLS problem here for ease of use
     144           4 :         WLSProblem<data_t> prob{scalingOp + lambdaOp, dcB};
     145             : 
     146           4 :         data_t epsilon = std::numeric_limits<data_t>::epsilon();
     147             : 
     148           4 :         WHEN("setting up a FGM solver")
     149           4 :         {
     150           2 :             TestType solver{prob, epsilon};
     151             : 
     152           2 :             THEN("the clone works correctly")
     153           2 :             {
     154           2 :                 auto fgmClone = solver.clone();
     155             : 
     156           2 :                 REQUIRE_NE(fgmClone.get(), &solver);
     157           2 :                 REQUIRE_EQ(*fgmClone, solver);
     158             : 
     159           2 :                 AND_THEN("it works as expected")
     160           2 :                 {
     161           2 :                     auto solution = solver.solve(dd.getNumberOfCoefficients());
     162             : 
     163           2 :                     DataContainer<data_t> resultsDifference =
     164           2 :                         (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           2 :                                                 epsilon * epsilon * dcB.squaredL2Norm(), 0.1f));
     170           2 :                 }
     171           2 :             }
     172           2 :         }
     173             : 
     174           4 :         WHEN("setting up a preconditioned FGM solver")
     175           4 :         {
     176           2 :             auto preconditionerInverse = JacobiPreconditioner<data_t>(scalingOp + lambdaOp, true);
     177           2 :             TestType solver{prob, preconditionerInverse, epsilon};
     178             : 
     179           2 :             THEN("the clone works correctly")
     180           2 :             {
     181           2 :                 auto fgmClone = solver.clone();
     182             : 
     183           2 :                 REQUIRE_NE(fgmClone.get(), &solver);
     184           2 :                 REQUIRE_EQ(*fgmClone, solver);
     185             : 
     186           2 :                 AND_THEN("it works as expected")
     187           2 :                 {
     188             :                     // a perfect preconditioner should allow for convergence in a single step
     189           2 :                     auto solution = solver.solve(dd.getNumberOfCoefficients());
     190             : 
     191           2 :                     DataContainer<data_t> resultsDifference =
     192           2 :                         (scalingOp + lambdaOp).apply(solution) - dcB;
     193             : 
     194             :                     // should have converged for the given number of iterations
     195           2 :                     REQUIRE_UNARY(checkApproxEq(resultsDifference.squaredL2Norm(),
     196           2 :                                                 epsilon * epsilon * dcB.squaredL2Norm(), 0.1f));
     197           2 :                 }
     198           2 :             }
     199           2 :         }
     200           4 :     }
     201           4 : }
     202             : 
     203             : TEST_CASE("FGM: Solving a simple phantom reconstruction")
     204           1 : {
     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           1 :     GIVEN("a Phantom reconstruction problem")
     212           1 :     {
     213             : 
     214           1 :         IndexVector_t size(2);
     215           1 :         size << 16, 16; // TODO: determine optimal phantom size for efficient testing
     216           1 :         auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
     217           1 :         auto& volumeDescriptor = phantom.getDataDescriptor();
     218             : 
     219           1 :         index_t numAngles{30}, arc{180};
     220           1 :         auto sinoDescriptor = CircleTrajectoryGenerator::createTrajectory(
     221           1 :             numAngles, phantom.getDataDescriptor(), arc, static_cast<real_t>(size(0)) * 100.0f,
     222           1 :             static_cast<real_t>(size(0)));
     223             : 
     224           1 :         SiddonsMethod projector(downcast<VolumeDescriptor>(volumeDescriptor), *sinoDescriptor);
     225             : 
     226           1 :         auto sinogram = projector.apply(phantom);
     227             : 
     228           1 :         WLSProblem problem(projector, sinogram);
     229           1 :         real_t epsilon = std::numeric_limits<real_t>::epsilon();
     230             : 
     231           1 :         WHEN("setting up a FGM solver")
     232           1 :         {
     233           1 :             FGM solver{problem, epsilon};
     234             : 
     235           1 :             THEN("the clone works correctly")
     236           1 :             {
     237           1 :                 auto fgmClone = solver.clone();
     238             : 
     239           1 :                 REQUIRE_NE(fgmClone.get(), &solver);
     240           1 :                 REQUIRE_EQ(*fgmClone, solver);
     241             : 
     242           1 :                 AND_THEN("it works as expected")
     243           1 :                 {
     244           1 :                     auto reconstruction = solver.solve(15);
     245             : 
     246           1 :                     DataContainer resultsDifference = reconstruction - phantom;
     247             : 
     248             :                     // should have converged for the given number of iterations
     249             :                     // does not converge to the optimal solution because of the regularization term
     250           1 :                     REQUIRE(checkApproxEq(resultsDifference.squaredL2Norm(),
     251           1 :                                           epsilon * epsilon * phantom.squaredL2Norm(), 0.1));
     252           1 :                 }
     253           1 :             }
     254           1 :         }
     255           1 :     }
     256           1 : }
     257             : 
     258             : TEST_SUITE_END();