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 22 : TEST_CASE_TEMPLATE("FourierTransform: Testing construction", data_t, float, double)
23 : {
24 8 : GIVEN("a descriptor")
25 : {
26 8 : IndexVector_t numCoeff(3);
27 4 : numCoeff << 13, 45, 28;
28 8 : VolumeDescriptor dd(numCoeff);
29 :
30 6 : WHEN("instantiating a FourierTransform operator")
31 : {
32 4 : FourierTransform fftOp(dd);
33 :
34 4 : THEN("the DataDescriptors are as expected")
35 : {
36 2 : REQUIRE(fftOp.getDomainDescriptor() == dd);
37 2 : REQUIRE(fftOp.getRangeDescriptor() == dd);
38 : }
39 : }
40 :
41 6 : WHEN("cloning a FourierTransform operator")
42 : {
43 4 : FourierTransform fftOp(dd);
44 4 : auto fftOpClone = fftOp.clone();
45 :
46 4 : THEN("everything matches")
47 : {
48 2 : REQUIRE(fftOpClone.get() != &fftOp);
49 2 : REQUIRE(*fftOpClone == fftOp);
50 : }
51 : }
52 : }
53 4 : }
54 :
55 24 : TEST_CASE_TEMPLATE("FourierTransform: Testing 1d", data_t, float, double)
56 : {
57 12 : GIVEN("a volume descriptor")
58 : {
59 12 : IndexVector_t numCoeff{1};
60 6 : numCoeff << 4;
61 12 : VolumeDescriptor dd{numCoeff};
62 :
63 12 : WHEN("performing 1d fft transform")
64 : {
65 : using namespace std::complex_literals;
66 12 : FourierTransform<complex<data_t>> fftOp{dd};
67 :
68 12 : Vector_t<complex<data_t>> input{4};
69 :
70 : // TODO: somehow magically do complex conversions with eigens , operator
71 : if constexpr (std::is_same_v<data_t, float>) {
72 27 : input << 0.f + 0.if, 42.f + 0.if, 42.f + 0.if, 0.f + 0.if;
73 : } else {
74 27 : input << 0. + 0.i, 42. + 0.i, 42. + 0.i, 0. + 0.i;
75 : }
76 12 : Vector_t<complex<data_t>> expected{4};
77 :
78 : // TODO: also check other datahandlers
79 :
80 8 : THEN("the forward transformed values are correct")
81 : {
82 4 : DataContainer<complex<data_t>> inputdc{dd, input, DataHandlerType::CPU};
83 4 : auto output = fftOp.apply(inputdc);
84 :
85 : 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 : } else {
88 1 : expected << 84. + 0.i, -42. - 42.i, 0. + 0.i, -42. + 42.i;
89 : }
90 :
91 10 : for (index_t i = 0; i < output.getSize(); ++i) {
92 8 : REQUIRE_UNARY(checkApproxEq(output[i], expected(i)));
93 : }
94 : }
95 :
96 8 : THEN("the inverse transformed values are correct")
97 : {
98 4 : DataContainer<complex<data_t>> inputdc{dd, input, DataHandlerType::CPU};
99 4 : auto output = fftOp.applyAdjoint(inputdc);
100 :
101 : 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 : } else {
104 1 : expected << 21. + 0.i, -10.5 + 10.5i, 0. + 0.i, -10.5 - 10.5i;
105 : }
106 :
107 10 : for (index_t i = 0; i < output.getSize(); ++i) {
108 8 : REQUIRE_UNARY(checkApproxEq(output[i], expected(i)));
109 : }
110 : }
111 :
112 8 : THEN("the forward-inverse transformed values are correct")
113 : {
114 4 : DataContainer<complex<data_t>> inputdc{dd, input, DataHandlerType::CPU};
115 :
116 4 : auto output = fftOp.applyAdjoint(fftOp.apply(inputdc));
117 :
118 : 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 : } else {
121 1 : expected << 0. + 0.i, 42. + 0.i, 42. + 0.i, 0. + 0.i;
122 : }
123 :
124 10 : for (index_t i = 0; i < output.getSize(); ++i) {
125 8 : REQUIRE_UNARY(checkApproxEq(output[i], expected(i)));
126 : }
127 : }
128 : }
129 : }
130 6 : }
131 :
132 28 : TEST_CASE_TEMPLATE("FourierTransform: 2d test", data_t, float, double)
133 : {
134 : using namespace std::complex_literals;
135 :
136 20 : GIVEN("some data")
137 : {
138 :
139 12 : WHEN("applying the fft")
140 : {
141 4 : IndexVector_t size{2};
142 2 : size << 4, 4;
143 4 : VolumeDescriptor dd{size};
144 :
145 4 : FourierTransform<complex<data_t>> fftOp{dd};
146 :
147 4 : 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 4 : auto input = testdata.asComplex();
153 :
154 4 : auto output = fftOp.apply(input);
155 :
156 4 : 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 4 : THEN("the result is as expected")
161 : {
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 : }
166 : }
167 : }
168 :
169 12 : WHEN("applying the adjoint of fft")
170 : {
171 4 : IndexVector_t size{2};
172 2 : size << 4, 4;
173 4 : VolumeDescriptor dd{size};
174 :
175 4 : FourierTransform<complex<data_t>> fftOp{dd};
176 :
177 4 : 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 4 : auto input = testdata.asComplex();
183 :
184 4 : auto output = fftOp.apply(input);
185 :
186 4 : 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 4 : THEN("the result is as expected")
191 : {
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 : }
196 : }
197 : }
198 :
199 12 : WHEN("applying the fft and inverse fft")
200 : {
201 :
202 4 : IndexVector_t size{4};
203 2 : size << 5, 10, 15, 20;
204 4 : VolumeDescriptor dd{size};
205 4 : FourierTransform<complex<data_t>> fftOp{dd};
206 :
207 4 : auto [input, randVec] =
208 : generateRandomContainer<complex<data_t>>(dd, DataHandlerType::CPU);
209 :
210 4 : auto mid = fftOp.apply(input);
211 4 : auto output = fftOp.applyAdjoint(mid);
212 :
213 4 : THEN("the back and forward fourier transform result in the original values")
214 : {
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 : }
219 : }
220 : }
221 :
222 14 : WHEN("ensuring the linearity of a fft")
223 : {
224 :
225 8 : IndexVector_t size{2};
226 4 : size << 20, 10;
227 8 : VolumeDescriptor dd{size};
228 :
229 8 : auto [inputA, randVecA] =
230 : generateRandomContainer<complex<data_t>>(dd, DataHandlerType::CPU);
231 8 : auto [inputB, randVecB] =
232 : generateRandomContainer<complex<data_t>>(dd, DataHandlerType::CPU);
233 :
234 8 : FourierTransform<complex<data_t>> fftOp{dd};
235 :
236 8 : auto fftA = fftOp.apply(inputA);
237 8 : auto fftB = fftOp.apply(inputB);
238 :
239 6 : THEN("fft(A + B) = fft(A) + fft(B)")
240 : {
241 4 : DataContainer<complex<data_t>> fftA_plus_fftB = fftA + fftB;
242 4 : DataContainer<complex<data_t>> A_plus_B = inputA + inputB;
243 :
244 4 : 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 : }
250 : }
251 :
252 6 : THEN("f(x * A) = x * f(A)")
253 : {
254 2 : complex<data_t> x = 42;
255 4 : DataContainer<complex<data_t>> fftA_scaled = (fftA * x);
256 :
257 4 : 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 : }
263 : }
264 : }
265 :
266 : // TODO impulse function response
267 : }
268 10 : }
269 : TEST_SUITE_END();
|
