Line data Source code
1 : /**
2 : * @file test_FourierTransform.cpp
3 : *
4 : * @brief Tests for the fourier transform operator
5 : *
6 : * @author Jonas Jelten
7 : */
8 :
9 : #include "ExpressionPredicates.h"
10 : #include "doctest/doctest.h"
11 : #include "elsaDefines.h"
12 : #include "testHelpers.h"
13 : #include "FourierTransform.h"
14 : #include "VolumeDescriptor.h"
15 : #include "DataHandlerCPU.h"
16 : #include <type_traits>
17 :
18 : using namespace elsa;
19 :
20 : TEST_SUITE_BEGIN("core");
21 :
22 : TEST_CASE_TEMPLATE("FourierTransform: Testing construction", data_t, float, double)
23 4 : {
24 4 : GIVEN("a descriptor")
25 4 : {
26 4 : IndexVector_t numCoeff(3);
27 4 : numCoeff << 13, 45, 28;
28 4 : VolumeDescriptor dd(numCoeff);
29 :
30 4 : WHEN("instantiating a FourierTransform operator")
31 4 : {
32 2 : FourierTransform fftOp(dd);
33 :
34 2 : THEN("the DataDescriptors are as expected")
35 2 : {
36 2 : REQUIRE(fftOp.getDomainDescriptor() == dd);
37 2 : REQUIRE(fftOp.getRangeDescriptor() == dd);
38 2 : }
39 2 : }
40 :
41 4 : WHEN("cloning a FourierTransform operator")
42 4 : {
43 2 : FourierTransform fftOp(dd);
44 2 : auto fftOpClone = fftOp.clone();
45 :
46 2 : THEN("everything matches")
47 2 : {
48 2 : REQUIRE(fftOpClone.get() != &fftOp);
49 2 : REQUIRE(*fftOpClone == fftOp);
50 2 : }
51 2 : }
52 4 : }
53 4 : }
54 :
55 : TEST_CASE_TEMPLATE("FourierTransform: Testing 1d", data_t, float, double)
56 6 : {
57 6 : GIVEN("a volume descriptor")
58 6 : {
59 6 : IndexVector_t numCoeff{1};
60 6 : numCoeff << 4;
61 6 : VolumeDescriptor dd{numCoeff};
62 :
63 6 : WHEN("performing 1d fft transform")
64 6 : {
65 6 : using namespace std::complex_literals;
66 6 : FourierTransform<complex<data_t>> fftOp{dd};
67 :
68 6 : Vector_t<complex<data_t>> input{4};
69 :
70 : // TODO: somehow magically do complex conversions with eigens , operator
71 6 : if constexpr (std::is_same_v<data_t, float>) {
72 3 : input << 0.f + 0.if, 42.f + 0.if, 42.f + 0.if, 0.f + 0.if;
73 3 : } else {
74 3 : input << 0. + 0.i, 42. + 0.i, 42. + 0.i, 0. + 0.i;
75 3 : }
76 6 : Vector_t<complex<data_t>> expected{4};
77 :
78 : // TODO: also check other datahandlers
79 :
80 6 : THEN("the forward transformed values are correct")
81 6 : {
82 2 : DataContainer<complex<data_t>> inputdc{dd, input, DataHandlerType::CPU};
83 2 : auto output = fftOp.apply(inputdc);
84 :
85 2 : if constexpr (std::is_same_v<data_t, float>) {
86 1 : expected << 84.f + 0.if, -42.f - 42.if, 0.f + 0.if, -42.f + 42.if;
87 1 : } else {
88 1 : expected << 84. + 0.i, -42. - 42.i, 0. + 0.i, -42. + 42.i;
89 1 : }
90 :
91 10 : for (index_t i = 0; i < output.getSize(); ++i) {
92 8 : REQUIRE_UNARY(checkApproxEq(output[i], expected(i)));
93 8 : }
94 2 : }
95 :
96 6 : THEN("the inverse transformed values are correct")
97 6 : {
98 2 : DataContainer<complex<data_t>> inputdc{dd, input, DataHandlerType::CPU};
99 2 : auto output = fftOp.applyAdjoint(inputdc);
100 :
101 2 : if constexpr (std::is_same_v<data_t, float>) {
102 1 : expected << 21.f + 0.if, -10.5f + 10.5if, 0.f + 0.if, -10.5f - 10.5if;
103 1 : } else {
104 1 : expected << 21. + 0.i, -10.5 + 10.5i, 0. + 0.i, -10.5 - 10.5i;
105 1 : }
106 :
107 10 : for (index_t i = 0; i < output.getSize(); ++i) {
108 8 : REQUIRE_UNARY(checkApproxEq(output[i], expected(i)));
109 8 : }
110 2 : }
111 :
112 6 : THEN("the forward-inverse transformed values are correct")
113 6 : {
114 2 : DataContainer<complex<data_t>> inputdc{dd, input, DataHandlerType::CPU};
115 :
116 2 : auto output = fftOp.applyAdjoint(fftOp.apply(inputdc));
117 :
118 2 : if constexpr (std::is_same_v<data_t, float>) {
119 1 : expected << 0.f + 0.if, 42.f + 0.if, 42.f + 0.if, 0.f + 0.if;
120 1 : } else {
121 1 : expected << 0. + 0.i, 42. + 0.i, 42. + 0.i, 0. + 0.i;
122 1 : }
123 :
124 10 : for (index_t i = 0; i < output.getSize(); ++i) {
125 8 : REQUIRE_UNARY(checkApproxEq(output[i], expected(i)));
126 8 : }
127 2 : }
128 6 : }
129 6 : }
130 6 : }
131 :
132 : TEST_CASE_TEMPLATE("FourierTransform: 2d test", data_t, float, double)
133 10 : {
134 10 : using namespace std::complex_literals;
135 :
136 10 : GIVEN("some data")
137 10 : {
138 :
139 10 : WHEN("applying the fft")
140 10 : {
141 2 : IndexVector_t size{2};
142 2 : size << 4, 4;
143 2 : VolumeDescriptor dd{size};
144 :
145 2 : FourierTransform<complex<data_t>> fftOp{dd};
146 :
147 2 : DataContainer<data_t> testdata{dd};
148 2 : testdata = 0.0;
149 2 : testdata(2, 2) = 42.0f;
150 2 : testdata(0, 2) = 42.0f;
151 :
152 2 : auto input = testdata.asComplex();
153 :
154 2 : auto output = fftOp.apply(input);
155 :
156 2 : Vector_t<complex<data_t>> expected{4 * 4};
157 : // transposed because vector!
158 2 : expected << 84, 0, 84, 0, -84, 0, -84, 0, 84, 0, 84, 0, -84, 0, -84, 0;
159 :
160 2 : THEN("the result is as expected")
161 2 : {
162 2 : REQUIRE(output.getSize() == input.getSize());
163 34 : for (index_t i = 0; i < output.getSize(); ++i) {
164 32 : REQUIRE_UNARY(checkApproxEq(output[i], expected(i)));
165 32 : }
166 2 : }
167 2 : }
168 :
169 10 : WHEN("applying the adjoint of fft")
170 10 : {
171 2 : IndexVector_t size{2};
172 2 : size << 4, 4;
173 2 : VolumeDescriptor dd{size};
174 :
175 2 : FourierTransform<complex<data_t>> fftOp{dd};
176 :
177 2 : DataContainer<data_t> testdata{dd};
178 2 : testdata = 0.0;
179 2 : testdata(2, 2) = 42.0f;
180 2 : testdata(0, 2) = 42.0f;
181 :
182 2 : auto input = testdata.asComplex();
183 :
184 2 : auto output = fftOp.apply(input);
185 :
186 2 : Vector_t<complex<data_t>> expected{4 * 4};
187 : // transposed because vector!
188 2 : expected << 84, 0, 84, 0, -84, 0, -84, 0, 84, 0, 84, 0, -84, 0, -84, 0;
189 :
190 2 : THEN("the result is as expected")
191 2 : {
192 2 : REQUIRE(output.getSize() == input.getSize());
193 34 : for (index_t i = 0; i < output.getSize(); ++i) {
194 32 : REQUIRE_UNARY(checkApproxEq(output[i], expected(i)));
195 32 : }
196 2 : }
197 2 : }
198 :
199 10 : WHEN("applying the fft and inverse fft")
200 10 : {
201 :
202 2 : IndexVector_t size{4};
203 2 : size << 5, 10, 15, 20;
204 2 : VolumeDescriptor dd{size};
205 2 : FourierTransform<complex<data_t>> fftOp{dd};
206 :
207 2 : auto [input, randVec] =
208 2 : generateRandomContainer<complex<data_t>>(dd, DataHandlerType::CPU);
209 :
210 2 : auto mid = fftOp.apply(input);
211 2 : auto output = fftOp.applyAdjoint(mid);
212 :
213 2 : THEN("the back and forward fourier transform result in the original values")
214 2 : {
215 2 : REQUIRE(output.getSize() == input.getSize());
216 30002 : for (index_t i = 0; i < output.getSize(); ++i) {
217 30000 : REQUIRE_UNARY(checkApproxEq(output[i], randVec(i)));
218 30000 : }
219 2 : }
220 2 : }
221 :
222 10 : WHEN("ensuring the linearity of a fft")
223 10 : {
224 :
225 4 : IndexVector_t size{2};
226 4 : size << 20, 10;
227 4 : VolumeDescriptor dd{size};
228 :
229 4 : auto [inputA, randVecA] =
230 4 : generateRandomContainer<complex<data_t>>(dd, DataHandlerType::CPU);
231 4 : auto [inputB, randVecB] =
232 4 : generateRandomContainer<complex<data_t>>(dd, DataHandlerType::CPU);
233 :
234 4 : FourierTransform<complex<data_t>> fftOp{dd};
235 :
236 4 : auto fftA = fftOp.apply(inputA);
237 4 : auto fftB = fftOp.apply(inputB);
238 :
239 4 : THEN("fft(A + B) = fft(A) + fft(B)")
240 4 : {
241 2 : DataContainer<complex<data_t>> fftA_plus_fftB = fftA + fftB;
242 2 : DataContainer<complex<data_t>> A_plus_B = inputA + inputB;
243 :
244 2 : auto fftAB = fftOp.apply(A_plus_B);
245 :
246 2 : REQUIRE(fftAB.getSize() == fftA_plus_fftB.getSize());
247 402 : for (index_t i = 0; i < fftAB.getSize(); ++i) {
248 400 : REQUIRE_UNARY(checkApproxEq(fftA_plus_fftB[i], fftAB[i]));
249 400 : }
250 2 : }
251 :
252 4 : THEN("f(x * A) = x * f(A)")
253 4 : {
254 2 : complex<data_t> x = 42;
255 2 : DataContainer<complex<data_t>> fftA_scaled = (fftA * x);
256 :
257 2 : auto fftAx = fftOp.apply(inputA * x);
258 :
259 2 : REQUIRE(fftA_scaled.getSize() == fftAx.getSize());
260 402 : for (index_t i = 0; i < fftA_scaled.getSize(); ++i) {
261 400 : REQUIRE_UNARY(checkApproxEq(fftA_scaled[i], fftAx[i]));
262 400 : }
263 2 : }
264 4 : }
265 :
266 : // TODO impulse function response
267 10 : }
268 10 : }
269 : TEST_SUITE_END();
|
