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 10 : TYPE_TO_STRING(FGM<float>);
33 10 : TYPE_TO_STRING(FGM<double>);
34 :
35 19 : TEST_CASE_TEMPLATE("FGM: Solving a simple linear problem", TestType, FGM<float>, FGM<double>)
36 : {
37 : // Set seed for Eigen Matrices!
38 4 : srand((unsigned int) 666);
39 :
40 : 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 8 : GIVEN("a linear problem")
45 : {
46 8 : IndexVector_t numCoeff(2);
47 4 : numCoeff << 13, 24;
48 8 : VolumeDescriptor dd{numCoeff};
49 :
50 8 : Eigen::Matrix<data_t, -1, 1> bVec(dd.getNumberOfCoefficients());
51 4 : bVec.setRandom();
52 8 : DataContainer<data_t> dcB{dd, bVec};
53 :
54 4 : bVec.setRandom();
55 4 : bVec = bVec.cwiseAbs();
56 8 : Scaling<data_t> scalingOp{dd, DataContainer<data_t>{dd, bVec}};
57 :
58 : // using WLS problem here for ease of use
59 8 : WLSProblem prob{scalingOp, dcB};
60 :
61 4 : data_t epsilon = std::numeric_limits<data_t>::epsilon();
62 :
63 6 : WHEN("setting up a FGM solver")
64 : {
65 4 : TestType solver{prob, epsilon};
66 :
67 4 : THEN("the clone works correctly")
68 : {
69 4 : auto fgmClone = solver.clone();
70 :
71 2 : REQUIRE_NE(fgmClone.get(), &solver);
72 2 : REQUIRE_EQ(*fgmClone, solver);
73 :
74 4 : AND_THEN("it works as expected")
75 : {
76 4 : 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 : epsilon * epsilon * dcB.squaredL2Norm(), 0.5f));
83 : }
84 : }
85 : }
86 :
87 6 : WHEN("setting up a preconditioned FGM solver")
88 : {
89 4 : auto preconditionerInverse = JacobiPreconditioner<data_t>(scalingOp, true);
90 4 : TestType solver{prob, preconditionerInverse, epsilon};
91 :
92 4 : THEN("the clone works correctly")
93 : {
94 4 : auto fgmClone = solver.clone();
95 :
96 2 : REQUIRE_NE(fgmClone.get(), &solver);
97 2 : REQUIRE_EQ(*fgmClone, solver);
98 :
99 4 : AND_THEN("it works as expected")
100 : {
101 : // with a good preconditioner we should need fewer iterations than without
102 4 : 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 : epsilon * epsilon * dcB.squaredL2Norm(), 0.1f));
109 : }
110 : }
111 : }
112 : }
113 4 : }
114 :
115 19 : TEST_CASE_TEMPLATE("FGM: Solving a Tikhonov problem", TestType, FGM<float>, FGM<double>)
116 : {
117 : // Set seed for Eigen Matrices!
118 4 : srand((unsigned int) 666);
119 :
120 : 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 8 : GIVEN("a Tikhonov problem")
125 : {
126 8 : IndexVector_t numCoeff(2);
127 4 : numCoeff << 13, 24;
128 8 : VolumeDescriptor dd(numCoeff);
129 :
130 8 : Eigen::Matrix<data_t, -1, 1> bVec(dd.getNumberOfCoefficients());
131 4 : bVec.setRandom();
132 8 : DataContainer dcB(dd, bVec);
133 :
134 : // the regularization term
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 FGM solver")
149 : {
150 4 : TestType solver{prob, epsilon};
151 :
152 4 : THEN("the clone works correctly")
153 : {
154 4 : auto fgmClone = solver.clone();
155 :
156 2 : REQUIRE_NE(fgmClone.get(), &solver);
157 2 : REQUIRE_EQ(*fgmClone, 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(), 0.1f));
170 : }
171 : }
172 : }
173 :
174 6 : WHEN("setting up a preconditioned FGM solver")
175 : {
176 4 : auto preconditionerInverse = JacobiPreconditioner<data_t>(scalingOp + lambdaOp, true);
177 4 : TestType solver{prob, preconditionerInverse, epsilon};
178 :
179 4 : THEN("the clone works correctly")
180 : {
181 4 : auto fgmClone = solver.clone();
182 :
183 2 : REQUIRE_NE(fgmClone.get(), &solver);
184 2 : REQUIRE_EQ(*fgmClone, 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(), 0.1f));
197 : }
198 : }
199 : }
200 : }
201 4 : }
202 :
203 1 : TEST_CASE("FGM: 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 :
214 2 : IndexVector_t size(2);
215 1 : size << 16, 16; // TODO: determine optimal phantom size for efficient testing
216 2 : auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
217 1 : auto& volumeDescriptor = phantom.getDataDescriptor();
218 :
219 1 : index_t numAngles{30}, arc{180};
220 : auto sinoDescriptor = CircleTrajectoryGenerator::createTrajectory(
221 2 : numAngles, phantom.getDataDescriptor(), arc, static_cast<real_t>(size(0)) * 100.0f,
222 2 : static_cast<real_t>(size(0)));
223 :
224 2 : SiddonsMethod projector(downcast<VolumeDescriptor>(volumeDescriptor), *sinoDescriptor);
225 :
226 2 : auto sinogram = projector.apply(phantom);
227 :
228 2 : WLSProblem problem(projector, sinogram);
229 1 : real_t epsilon = std::numeric_limits<real_t>::epsilon();
230 :
231 2 : WHEN("setting up a FGM solver")
232 : {
233 2 : FGM solver{problem, epsilon};
234 :
235 2 : THEN("the clone works correctly")
236 : {
237 2 : auto fgmClone = solver.clone();
238 :
239 1 : REQUIRE_NE(fgmClone.get(), &solver);
240 1 : REQUIRE_EQ(*fgmClone, solver);
241 :
242 2 : AND_THEN("it works as expected")
243 : {
244 2 : 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 : epsilon * epsilon * phantom.squaredL2Norm(), 0.1));
252 : }
253 : }
254 : }
255 : }
256 1 : }
257 :
258 : TEST_SUITE_END();
|
