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 : TYPE_TO_STRING(SQS<float>);
35 : TYPE_TO_STRING(SQS<double>);
36 :
37 : TEST_CASE_TEMPLATE("SQS: Solving a simple linear problem", TestType, SQS<float>, SQS<double>)
38 4 : {
39 : // Set seed for Eigen Matrices!
40 4 : srand((unsigned int) 666);
41 :
42 4 : 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 4 : GIVEN("a linear problem")
47 4 : {
48 4 : IndexVector_t numCoeff(2);
49 4 : numCoeff << 8, 14;
50 4 : VolumeDescriptor dd{numCoeff};
51 :
52 4 : Eigen::Matrix<data_t, -1, 1> bVec(dd.getNumberOfCoefficients());
53 4 : bVec.setRandom();
54 4 : DataContainer<data_t> dcB{dd, bVec};
55 :
56 4 : bVec.setRandom();
57 4 : bVec = bVec.cwiseAbs();
58 4 : Scaling<data_t> scalingOp{dd, DataContainer<data_t>{dd, bVec}};
59 :
60 : // using WLS problem here for ease of use
61 4 : WLSProblem prob{scalingOp, dcB};
62 :
63 4 : data_t epsilon = std::numeric_limits<data_t>::epsilon();
64 :
65 4 : WHEN("setting up a SQS solver")
66 4 : {
67 2 : TestType solver{prob, true, epsilon};
68 :
69 2 : THEN("the clone works correctly")
70 2 : {
71 2 : auto sqsClone = solver.clone();
72 :
73 2 : REQUIRE_NE(sqsClone.get(), &solver);
74 2 : REQUIRE_EQ(*sqsClone, solver);
75 :
76 2 : AND_THEN("it works as expected")
77 2 : {
78 2 : auto solution = solver.solve(300);
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 2 : epsilon * epsilon * dcB.squaredL2Norm()));
85 2 : }
86 2 : }
87 2 : }
88 :
89 4 : WHEN("setting up a preconditioned SQS solver")
90 4 : {
91 2 : auto preconditioner = JacobiPreconditioner<data_t>(scalingOp, false);
92 2 : TestType solver{prob, preconditioner, true, epsilon};
93 :
94 2 : THEN("the clone works correctly")
95 2 : {
96 2 : auto sqsClone = solver.clone();
97 :
98 2 : REQUIRE_NE(sqsClone.get(), &solver);
99 2 : REQUIRE_EQ(*sqsClone, solver);
100 :
101 2 : AND_THEN("it works as expected")
102 2 : {
103 : // with a good preconditioner we should need fewer iterations than without
104 2 : auto solution = solver.solve(200);
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 2 : epsilon * epsilon * dcB.squaredL2Norm()));
111 2 : }
112 2 : }
113 2 : }
114 4 : }
115 4 : }
116 :
117 : TEST_CASE_TEMPLATE("SQS: Solving a Tikhonov problem", TestType, SQS<float>, SQS<double>)
118 4 : {
119 : // Set seed for Eigen Matrices!
120 4 : srand((unsigned int) 666);
121 :
122 4 : 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 4 : GIVEN("a Tikhonov problem")
127 4 : {
128 4 : IndexVector_t numCoeff(2);
129 4 : numCoeff << 13, 24;
130 4 : VolumeDescriptor dd(numCoeff);
131 :
132 4 : Eigen::Matrix<data_t, -1, 1> bVec(dd.getNumberOfCoefficients());
133 4 : bVec.setRandom();
134 4 : DataContainer dcB(dd, bVec);
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 SQS solver")
149 4 : {
150 2 : TestType solver{prob, true, epsilon};
151 :
152 2 : THEN("the clone works correctly")
153 2 : {
154 2 : auto sqsClone = solver.clone();
155 :
156 2 : REQUIRE_NE(sqsClone.get(), &solver);
157 2 : REQUIRE_EQ(*sqsClone, 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()));
170 2 : }
171 2 : }
172 2 : }
173 :
174 4 : WHEN("setting up a preconditioned SQS solver")
175 4 : {
176 2 : auto preconditioner = JacobiPreconditioner<data_t>(scalingOp + lambdaOp, false);
177 2 : TestType solver{prob, preconditioner, true, epsilon};
178 :
179 2 : THEN("the clone works correctly")
180 2 : {
181 2 : auto sqsClone = solver.clone();
182 :
183 2 : REQUIRE_NE(sqsClone.get(), &solver);
184 2 : REQUIRE_EQ(*sqsClone, 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()));
197 2 : }
198 2 : }
199 2 : }
200 4 : }
201 4 : }
202 :
203 : TEST_CASE("SQS: 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 1 : IndexVector_t size(2);
214 1 : size << 16, 16; // TODO: determine optimal phantom size for efficient testing
215 1 : auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
216 1 : auto& volumeDescriptor = phantom.getDataDescriptor();
217 :
218 1 : index_t numAngles{90}, arc{360};
219 1 : auto sinoDescriptor = CircleTrajectoryGenerator::createTrajectory(
220 1 : numAngles, phantom.getDataDescriptor(), arc, static_cast<real_t>(size(0)) * 100.0f,
221 1 : static_cast<real_t>(size(0)));
222 :
223 1 : SiddonsMethod projector(downcast<VolumeDescriptor>(volumeDescriptor), *sinoDescriptor);
224 :
225 1 : auto sinogram = projector.apply(phantom);
226 :
227 1 : WLSProblem problem(projector, sinogram);
228 1 : real_t epsilon = std::numeric_limits<real_t>::epsilon();
229 :
230 1 : WHEN("setting up a SQS solver")
231 1 : {
232 1 : SQS solver{problem, true, epsilon};
233 :
234 1 : THEN("the clone works correctly")
235 1 : {
236 1 : auto sqsClone = solver.clone();
237 :
238 1 : REQUIRE_NE(sqsClone.get(), &solver);
239 1 : REQUIRE_EQ(*sqsClone, solver);
240 :
241 1 : AND_THEN("it works as expected")
242 1 : {
243 1 : 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 1 : }
251 1 : }
252 1 : }
253 1 : }
254 1 : }
255 :
256 : TEST_CASE("SQS: Solving a simple phantom problem using ordered subsets")
257 2 : {
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 2 : GIVEN("a Phantom reconstruction problem")
265 2 : {
266 :
267 2 : IndexVector_t size(2);
268 2 : size << 16, 16; // TODO: determine optimal phantom size for efficient testing
269 2 : auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
270 2 : auto& volumeDescriptor = phantom.getDataDescriptor();
271 :
272 2 : index_t numAngles{20}, arc{180};
273 2 : auto sinoDescriptor = CircleTrajectoryGenerator::createTrajectory(
274 2 : numAngles, phantom.getDataDescriptor(), arc, static_cast<real_t>(size(0)) * 100.0f,
275 2 : static_cast<real_t>(size(0)));
276 :
277 2 : SiddonsMethod projector(static_cast<const VolumeDescriptor&>(volumeDescriptor),
278 2 : *sinoDescriptor);
279 :
280 2 : auto sinogram = projector.apply(phantom);
281 :
282 2 : real_t epsilon = std::numeric_limits<real_t>::epsilon();
283 :
284 2 : WHEN("setting up a SQS solver with ROUND_ROBIN subsampling")
285 2 : {
286 1 : index_t nSubsets{4};
287 1 : SubsetSampler<PlanarDetectorDescriptor, real_t> subsetSampler(
288 1 : static_cast<const VolumeDescriptor&>(volumeDescriptor),
289 1 : static_cast<const PlanarDetectorDescriptor&>(*sinoDescriptor), nSubsets);
290 :
291 1 : WLSSubsetProblem<real_t> problem(
292 1 : *subsetSampler.getProjector<SiddonsMethod<real_t>>(),
293 1 : subsetSampler.getPartitionedData(sinogram),
294 1 : subsetSampler.getSubsetProjectors<SiddonsMethod<real_t>>());
295 1 : SQS solver{problem, true, epsilon};
296 :
297 1 : THEN("the clone works correctly")
298 1 : {
299 1 : auto sqsClone = solver.clone();
300 :
301 1 : REQUIRE(sqsClone.get() != &solver);
302 1 : REQUIRE(*sqsClone == solver);
303 :
304 1 : AND_THEN("it works as expected")
305 1 : {
306 1 : 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 1 : epsilon * epsilon * phantom.squaredL2Norm(), 0.1));
314 1 : }
315 1 : }
316 1 : }
317 2 : WHEN("setting up a SQS solver with ROTATIONAL_CLUSTERING subsampling")
318 2 : {
319 1 : index_t nSubsets{4};
320 1 : SubsetSampler<PlanarDetectorDescriptor, real_t> subsetSampler(
321 1 : static_cast<const VolumeDescriptor&>(volumeDescriptor),
322 1 : static_cast<const PlanarDetectorDescriptor&>(*sinoDescriptor), nSubsets,
323 1 : SubsetSampler<PlanarDetectorDescriptor,
324 1 : real_t>::SamplingStrategy::ROTATIONAL_CLUSTERING);
325 :
326 1 : WLSSubsetProblem<real_t> problem(
327 1 : *subsetSampler.getProjector<SiddonsMethod<real_t>>(),
328 1 : subsetSampler.getPartitionedData(sinogram),
329 1 : subsetSampler.getSubsetProjectors<SiddonsMethod<real_t>>());
330 1 : SQS solver{problem, true, epsilon};
331 :
332 1 : THEN("the clone works correctly")
333 1 : {
334 1 : auto sqsClone = solver.clone();
335 :
336 1 : REQUIRE_NE(sqsClone.get(), &solver);
337 1 : REQUIRE_EQ(*sqsClone, solver);
338 :
339 1 : AND_THEN("it works as expected")
340 1 : {
341 1 : 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 1 : epsilon * epsilon * phantom.squaredL2Norm(), 0.1));
348 1 : }
349 1 : }
350 1 : }
351 2 : }
352 2 : }
353 :
354 : TEST_SUITE_END();
|
