1 /************************************************************************/
2 /* */
3 /* Copyright 2004 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35
36 #include "vigra/unittest.hxx"
37 #include <stdlib.h>
38 #include <algorithm>
39 #include <functional>
40 #include <vigra/stdimage.hxx>
41 #include <vigra/stdimagefunctions.hxx>
42 #include <vigra/functorexpression.hxx>
43 #include <vigra/fftw3.hxx>
44 #include <vigra/impex.hxx>
45 #include <vigra/inspectimage.hxx>
46 #include <vigra/gaborfilter.hxx>
47 #include <vigra/multi_fft.hxx>
48 #include <vigra/multi_pointoperators.hxx>
49 #include <vigra/convolution.hxx>
50 #include "test.hxx"
51
52 using namespace std;
53 using namespace vigra;
54 using namespace vigra::functor;
55
56 struct Compare1
57 {
58 double s;
59
60 typedef FFTWComplexImage::value_type value_type;
61 typedef FFTWComplexImage::value_type result_type;
62
Compare1Compare163 Compare1(double is)
64 : s(is)
65 {}
66
operator ()Compare167 result_type operator()(value_type v1, value_type v2) const
68 { return v1 - v2/s; }
69 };
70
71 struct FFTWComplexTest
72 {
73 FFTWComplex<> clx0, clx1, clx2, clx3;
74
75 template <class VECTOR>
printVectorFFTWComplexTest76 void printVector(VECTOR const & v)
77 {
78 std::cerr << "(";
79 for(int i=0; i<v.size(); ++i)
80 std::cerr << v[i] << ", ";
81 std::cerr << ")\n";
82 }
83
FFTWComplexTestFFTWComplexTest84 FFTWComplexTest()
85 : clx0(), clx1(2.0), clx2(0.0, -2.0), clx3(2.0, -2.0)
86 {}
87
testConstructionFFTWComplexTest88 void testConstruction()
89 {
90 shouldEqual(clx0.re() , 0.0);
91 shouldEqual(clx0.im() , 0.0);
92 shouldEqual(clx1.re() , 2.0);
93 shouldEqual(clx1.im() , 0.0);
94 shouldEqual(clx2.re() , 0.0);
95 shouldEqual(clx2.im() , -2.0);
96 shouldEqual(clx3.re() , 2.0);
97 shouldEqual(real(clx3) , 2.0);
98 shouldEqual(clx3.im() , -2.0);
99 shouldEqual(imag(clx3) , -2.0);
100 shouldEqual(clx0[0] , 0.0);
101 shouldEqual(clx0[1] , 0.0);
102 shouldEqual(clx1[0] , 2.0);
103 shouldEqual(clx1[1] , 0.0);
104 shouldEqual(clx2[0] , 0.0);
105 shouldEqual(clx2[1] , -2.0);
106 shouldEqual(clx3[0] , 2.0);
107 shouldEqual(clx3[1] , -2.0);
108 }
109
testComparisonFFTWComplexTest110 void testComparison()
111 {
112 should(clx0 == FFTWComplex<>(0.0, 0.0));
113 should(clx1 == FFTWComplex<>(2.0, 0.0));
114 should(clx2 == FFTWComplex<>(0.0, -2.0));
115 should(clx3 == FFTWComplex<>(2.0, -2.0));
116 should(clx0 != clx1);
117 should(clx0 != clx2);
118 should(clx1 != clx2);
119 should(clx1 != clx3);
120 }
121
testArithmeticFFTWComplexTest122 void testArithmetic()
123 {
124 shouldEqual(abs(clx0) , 0.0);
125 shouldEqual(abs(clx1) , 2.0);
126 shouldEqual(abs(clx2) , 2.0);
127 shouldEqual(abs(clx3) , sqrt(8.0));
128 should(conj(clx3) == FFTWComplex<>(2.0, 2.0));
129
130 shouldEqual(clx1.phase() , 0.0);
131 shouldEqual(arg(clx1) , 0.0);
132 shouldEqualTolerance(arg(clx3), -0.25*M_PI, 1e-16);
133 shouldEqual(sin(clx2.phase()) , -1.0);
134 shouldEqualTolerance(sin(clx3.phase()), -sqrt(2.0)/2.0, 1.0e-7);
135
136 should(FFTWComplex<>(2.0, -2.0) == clx1 + clx2);
137
138 should(clx0 == clx1 - clx1);
139 should(clx0 == clx2 - clx2);
140 should(FFTWComplex<>(2.0, 2.0) == clx1 - clx2);
141
142 should(2.0*clx1 == FFTWComplex<>(4.0, 0.0));
143 should(2.0*clx2 == FFTWComplex<>(0.0, -4.0));
144 should(clx1*clx2 == FFTWComplex<>(0.0, -4.0));
145 should(clx3*conj(clx3) == clx3.squaredMagnitude());
146
147 should(clx1/2.0 == FFTWComplex<>(1.0, 0.0));
148 should(clx2/2.0 == FFTWComplex<>(0.0, -1.0));
149 should(clx2/clx1 == FFTWComplex<>(0.0, -1.0));
150 should(clx3*conj(clx3) == clx3.squaredMagnitude());
151
152 std::complex<FFTWComplex<>::value_type> c(clx3.re(), clx3.im()),
153 c1(clx1.re(), clx1.im());
154
155 shouldEqual(cos(clx3), FFTWComplex<>(cos(c)));
156 shouldEqual(cosh(clx3), FFTWComplex<>(cosh(c)));
157 shouldEqual(exp(clx3), FFTWComplex<>(exp(c)));
158 shouldEqual(log(clx3), FFTWComplex<>(log(c)));
159 shouldEqual(log10(clx3), FFTWComplex<>(log10(c)));
160 shouldEqual(sin(clx3), FFTWComplex<>(sin(c)));
161 shouldEqual(sinh(clx3), FFTWComplex<>(sinh(c)));
162 shouldEqual(sqrt(clx3), FFTWComplex<>(sqrt(c)));
163 shouldEqual(tan(clx3), FFTWComplex<>(tan(c)));
164 shouldEqual(tanh(clx3), FFTWComplex<>(tanh(c)));
165 shouldEqual(pow(clx3, 2), FFTWComplex<>(pow(c, 2)));
166 shouldEqual(pow(clx3, 2.0), FFTWComplex<>(pow(c, 2.0)));
167 shouldEqual(pow(clx3, clx1), FFTWComplex<>(pow(c, c1)));
168 shouldEqual(pow(2.0, clx3), FFTWComplex<>(pow(2.0, c)));
169 }
170
testAccessorFFTWComplexTest171 void testAccessor()
172 {
173 FFTWComplexImage img(2,2);
174 img = clx2;
175 img(0,0) = clx3;
176 FFTWComplexImage::Iterator i = img.upperLeft();
177
178 FFTWComplexImage::Accessor get = img.accessor();
179 FFTWRealAccessor<> real;
180 FFTWImaginaryAccessor<> imag;
181 FFTWSquaredMagnitudeAccessor<> smag;
182 FFTWMagnitudeAccessor<> mag;
183 FFTWPhaseAccessor<> phase;
184 FFTWWriteRealAccessor<> writeReal;
185
186 should(get(i) == clx3);
187 should(get(i, Diff2D(1,1)) == clx2);
188 should(real(i) == 2.0);
189 should(real(i, Diff2D(1,1)) == 0.0);
190 should(imag(i) == -2.0);
191 should(imag(i, Diff2D(1,1)) == -2.0);
192 shouldEqual(smag(i) , 8.0);
193 shouldEqual(smag(i, Diff2D(1,1)) , 4.0);
194 shouldEqual(mag(i) , sqrt(8.0));
195 shouldEqual(mag(i, Diff2D(1,1)) , 2.0);
196 shouldEqualTolerance(sin(phase(i)), -sqrt(2.0)/2.0, 1.0e-7);
197 shouldEqual(sin(phase(i, Diff2D(1,1))), -1.0);
198
199 writeReal.set(2.0, i);
200 writeReal.set(2.0, i, Diff2D(1,1));
201 should(get(i) == clx1);
202 should(get(i, Diff2D(1,1)) == clx1);
203 }
204
testForwardBackwardTransFFTWComplexTest205 void testForwardBackwardTrans()
206 {
207 const int w=256, h=256;
208
209 FFTWComplexImage in(w, h);
210 for (int y=0; y<in.height(); y++)
211 for (int x=0; x<in.width(); x++)
212 {
213 in(x,y)= rand()/(double)RAND_MAX;
214 }
215
216 FFTWComplexImage out(w, h);
217
218 fourierTransform(srcImageRange(in), destImage(out));
219 fourierTransformInverse(srcImageRange(out), destImage(out));
220
221 FindAverage<FFTWComplex<>::value_type> average;
222 inspectImage(srcImageRange(out, FFTWImaginaryAccessor<>()), average);
223
224 shouldEqualTolerance(average(), 0.0, 1e-14);
225
226 combineTwoImages(srcImageRange(in), srcImage(out), destImage(out),
227 Compare1((double)w*h));
228
229 average = FindAverage<FFTWMagnitudeAccessor<>::value_type>();
230 inspectImage(srcImageRange(out, FFTWMagnitudeAccessor<>()), average);
231
232 shouldEqualTolerance(average(), 0.0, 1e-14);
233
234 for (int y=0; y<in.height(); y++)
235 for (int x=0; x<in.width(); x++)
236 {
237 in(x,y)[1]= rand()/(double)RAND_MAX;
238 }
239
240 fourierTransform(srcImageRange(in), destImage(out));
241 fourierTransformInverse(srcImageRange(out), destImage(out));
242
243 combineTwoImages(srcImageRange(in), srcImage(out), destImage(out),
244 Compare1((double)w*h));
245
246 FindAverage<FFTWComplex<> > caverage;
247 inspectImage(srcImageRange(out), caverage);
248
249 shouldEqualTolerance(caverage().magnitude(), 0.0, 1e-14);
250 }
251
testRearrangeQuadrantsFFTWComplexTest252 void testRearrangeQuadrants()
253 {
254 double t4[] = { 0, 1, 2, 2,
255 1, 1, 2, 2,
256 3, 3, 4, 4,
257 3, 3, 4, 4};
258 double t4res[] = { 4, 4, 3, 3,
259 4, 4, 3, 3,
260 2, 2, 0, 1,
261 2, 2, 1, 1};
262
263 FFTWComplexImage in4(4,4), out4(4,4);
264 copy(t4, t4+16, in4.begin());
265 moveDCToCenter(srcImageRange(in4), destImage(out4));
266 moveDCToUpperLeft(srcImageRange(out4), destImage(in4));
267 for(int i=0; i<16; ++i)
268 {
269 should(out4.begin()[i] == t4res[i]);
270 should(in4.begin()[i] == t4[i]);
271 }
272
273 double t5[] = { 0, 1, 1, 2, 2,
274 1, 1, 1, 2, 2,
275 1, 1, 1, 2, 2,
276 3, 3, 3, 4, 4,
277 3, 3, 3, 4, 4};
278 double t5res[] = { 4, 4, 3, 3, 3,
279 4, 4, 3, 3, 3,
280 2, 2, 0, 1, 1,
281 2, 2, 1, 1, 1,
282 2, 2, 1, 1, 1};
283
284 FFTWComplexImage in5(5,5), out5(5,5);
285 copy(t5, t5+25, in5.begin());
286 moveDCToCenter(srcImageRange(in5), destImage(out5));
287 moveDCToUpperLeft(srcImageRange(out5), destImage(in5));
288 for(int i=0; i<25; ++i)
289 {
290 should(out5.begin()[i] == t5res[i]);
291 should(in5.begin()[i] == t5[i]);
292 }
293 }
294 };
295
296 struct MultiFFTTest
297 {
298 typedef double R;
299 typedef FFTWComplex<R> C;
300 typedef MultiArray<2, R, FFTWAllocator<R> > DArray2;
301 typedef MultiArray<2, C, FFTWAllocator<C> > CArray2;
302 typedef MultiArrayShape<2>::type Shape2;
303 typedef MultiArray<3, R, FFTWAllocator<R> > DArray3;
304 typedef MultiArray<3, C, FFTWAllocator<C> > CArray3;
305 typedef MultiArrayShape<3>::type Shape3;
306
testFFTShiftMultiFFTTest307 void testFFTShift()
308 {
309 Shape3 s(5,5,5);
310 MultiArray<3, unsigned int> a(s), ref(s), iref(s);
311
312 for(int z=0; z<5; ++z)
313 {
314 for(int y=0; y<5; ++y)
315 {
316 for(int x=0; x<5; ++x)
317 {
318 unsigned int v = ((z > 2) ? 4 : 0) + ((y > 2) ? 2 : 0) + ((x > 2) ? 1 : 0);
319 a(x,y,z) = iref(x,y,z) = v;
320 v = ((z < 2) ? 4 : 0) + ((y < 2) ? 2 : 0) + ((x < 2) ? 1 : 0);
321 ref(x,y,z) = v;
322 }
323 }
324 }
325
326 moveDCToCenter(a);
327 should(a == ref);
328 moveDCToUpperLeft(a);
329 should(a == iref);
330
331 MultiArrayView<3, unsigned int> b = a.subarray(Shape3(1,1,1), s);
332 moveDCToCenter(b);
333 should(b == ref.subarray(Shape3(0,0,0), Shape3(4,4,4)));
334 moveDCToUpperLeft(b);
335 should(b == iref.subarray(Shape3(1,1,1), s));
336 moveDCToCenter(b);
337 moveDCToCenter(b);
338 should(b == iref.subarray(Shape3(1,1,1), s));
339
340 MultiArrayShape<1>::type s1_e(10), s1_5(5), s1_6(6);
341 MultiArray<1, double> e1(s1_e), k1_5(s1_5), k1_6(s1_6);
342
343 for(int k=0; k<k1_5.size(); ++k)
344 k1_5(k) = k+1;
345 detail::fftEmbedKernel(k1_5, e1);
346
347 double k1_5ref[] = {3, 4, 5, 0, 0, 0, 0, 0, 1, 2};
348 shouldEqualSequence(e1.data(), e1.data()+e1.size(), k1_5ref);
349
350 for(int k=0; k<k1_6.size(); ++k)
351 k1_6(k) = k+1;
352 detail::fftEmbedKernel(k1_6, e1);
353
354 double k1_6ref[] = {4, 5, 6, 0, 0, 0, 0, 1, 2, 3};
355 shouldEqualSequence(e1.data(), e1.data()+e1.size(), k1_6ref);
356
357 detail::fftEmbedArray(k1_5, e1);
358 double a1_5ref[] = {3, 2, 1, 2, 3, 4, 5, 4, 3, 2 };
359 shouldEqualSequence(e1.data(), e1.data()+e1.size(), a1_5ref);
360
361 detail::fftEmbedArray(k1_6, e1);
362 double a1_6ref[] = {3, 2, 1, 2, 3, 4, 5, 6, 5, 4 };
363 shouldEqualSequence(e1.data(), e1.data()+e1.size(), a1_6ref);
364
365 MultiArrayShape<2>::type s2_e(8, 8), s2_4(4, 4), s2_5(5, 5);
366 MultiArray<2, double> e2(s2_e), k2_4(s2_4), k2_5(s2_5);
367
368 for(int k=0; k<k2_4.size(); ++k)
369 k2_4(k) = k+1;
370 detail::fftEmbedKernel(k2_4, e2);
371
372 double k2_4ref[] = {11,12, 0, 0, 0, 0, 9,10,
373 15,16, 0, 0, 0, 0,13,14,
374 0, 0, 0, 0, 0, 0, 0, 0,
375 0, 0, 0, 0, 0, 0, 0, 0,
376 0, 0, 0, 0, 0, 0, 0, 0,
377 0, 0, 0, 0, 0, 0, 0, 0,
378 3, 4, 0, 0, 0, 0, 1, 2,
379 7, 8, 0, 0, 0, 0, 5, 6};
380 shouldEqualSequence(e2.data(), e2.data()+e2.size(), k2_4ref);
381
382 for(int k=0; k<k2_5.size(); ++k)
383 k2_5(k) = k+1;
384 detail::fftEmbedKernel(k2_5, e2);
385
386 double k2_5ref[] = {13,14,15, 0, 0, 0,11,12,
387 18,19,20, 0, 0, 0,16,17,
388 23,24,25, 0, 0, 0,21,22,
389 0, 0, 0, 0, 0, 0, 0, 0,
390 0, 0, 0, 0, 0, 0, 0, 0,
391 0, 0, 0, 0, 0, 0, 0, 0,
392 3, 4, 5, 0, 0, 0, 1, 2,
393 8, 9,10, 0, 0, 0, 6, 7};
394 shouldEqualSequence(e2.data(), e2.data()+e2.size(), k2_5ref);
395
396 detail::fftEmbedArray(k2_4, e2);
397 double a2_4ref[] = {11,10, 9,10,11,12,11,10,
398 7, 6, 5, 6, 7, 8, 7, 6,
399 3, 2, 1, 2, 3, 4, 3, 2,
400 7, 6, 5, 6, 7, 8, 7, 6,
401 11,10, 9,10,11,12,11,10,
402 15,14,13,14,15,16,15,14,
403 11,10, 9,10,11,12,11,10,
404 7, 6, 5, 6, 7, 8, 7, 6};
405 shouldEqualSequence(e2.data(), e2.data()+e2.size(), a2_4ref);
406
407 detail::fftEmbedArray(k2_5, e2);
408 double a2_5ref[] = { 7, 6, 7, 8, 9,10, 9, 8,
409 2, 1, 2, 3, 4, 5, 4, 3,
410 7, 6, 7, 8, 9,10, 9, 8,
411 12,11,12,13,14,15,14,13,
412 17,16,17,18,19,20,19,18,
413 22,21,22,23,24,25,24,23,
414 17,16,17,18,19,20,19,18,
415 12,11,12,13,14,15,14,13};
416 shouldEqualSequence(e2.data(), e2.data()+e2.size(), a2_5ref);
417 }
418
testFFT2DMultiFFTTest419 void testFFT2D()
420 {
421 Shape2 s(256, 256);
422
423 FFTWComplexImage in(s[0], s[1]);
424 for (int y=0; y<in.height(); y++)
425 for (int x=0; x<in.width(); x++)
426 {
427 in(x,y)= rand()/(double)RAND_MAX;
428 }
429
430 FFTWComplexImage out(in.size());
431 CArray2 aout(s);
432
433 fourierTransform(srcImageRange(in), destImage(out));
434 fourierTransform(MultiArrayView<2, C>(s, const_cast<C*>(in.data())),
435 aout);
436
437 shouldEqualSequence(aout.data(), aout.data()+aout.size(), out.data());
438
439 DArray2 rin(s);
440 copyImage(srcImageRange(in, FFTWRealAccessor<>()), destImage(rin));
441
442 fourierTransform(rin, aout);
443 shouldEqualSequence(aout.data(), aout.data()+aout.size(), out.data());
444 }
445
testFFT3DMultiFFTTest446 void testFFT3D()
447 {
448 Shape3 s(32, 24, 16);
449 CArray3 r(s), ir(s), irr(s);
450
451 fourierTransform(MultiArrayView<3, double>(s, f3data), r);
452
453 DArray3 re(s);
454 copyMultiArray(srcMultiArrayRange(r, FFTWRealAccessor<>()), destMultiArray(re));
455 shouldEqualSequenceTolerance(re.data(), re.data()+re.size(), f3ref, 1e-10);
456
457 FindMinMax<double> minmax;
458 inspectMultiArray(srcMultiArrayRange(r, FFTWImaginaryAccessor<>()), minmax);
459 shouldEqualTolerance(minmax.min, 0.0, 1e-10);
460 shouldEqualTolerance(minmax.max, 0.0, 1e-10);
461
462 fourierTransformInverse(r, ir);
463
464 copyMultiArray(srcMultiArrayRange(ir, FFTWRealAccessor<>()), destMultiArray(re));
465 shouldEqualSequenceTolerance(re.data(), re.data()+re.size(), f3data, 1e-10);
466
467 inspectMultiArray(srcMultiArrayRange(ir, FFTWImaginaryAccessor<>()), minmax);
468 shouldEqualTolerance(minmax.min, 0.0, 1e-10);
469 shouldEqualTolerance(minmax.max, 0.0, 1e-10);
470 }
471
testPaddingMultiFFTTest472 void testPadding()
473 {
474 shouldEqual(0, detail::FFTWPaddingSize<0>::find(0));
475 shouldEqual(1, detail::FFTWPaddingSize<0>::find(1));
476 shouldEqual(3, detail::FFTWPaddingSize<0>::find(3));
477 shouldEqual(256, detail::FFTWPaddingSize<0>::find(255));
478 shouldEqual(256, detail::FFTWPaddingSize<0>::find(256));
479 shouldEqual(260, detail::FFTWPaddingSize<0>::find(257));
480 shouldEqual(0, detail::FFTWPaddingSize<0>::findEven(0));
481 shouldEqual(2, detail::FFTWPaddingSize<0>::findEven(1));
482 shouldEqual(4, detail::FFTWPaddingSize<0>::findEven(3));
483 shouldEqual(256, detail::FFTWPaddingSize<0>::findEven(255));
484 shouldEqual(256, detail::FFTWPaddingSize<0>::findEven(256));
485 shouldEqual(260, detail::FFTWPaddingSize<0>::findEven(257));
486
487 Shape3 s(113, 256, 257);
488 shouldEqual(Shape3(117, 256, 260), fftwBestPaddedShape(s));
489 shouldEqual(Shape3(120, 256, 260), fftwBestPaddedShapeR2C(s));
490 }
491
testConvolveFFTMultiFFTTest492 void testConvolveFFT()
493 {
494 typedef MultiArrayView<2, double> MV;
495 ImageImportInfo info("ghouse.gif");
496 Shape2 s(info.width(), info.height());
497 DArray2 in(s), out(s), ref(s), out2(s), ref2(s);
498 importImage(info, destImage(in));
499
500 double scale = 2.0;
501 gaussianSmoothing(srcImageRange(in), destImage(ref), scale);
502 gaussianSmoothing(srcImageRange(in), destImage(ref2), 2.0*scale);
503
504 Kernel2D<double> gauss;
505 gauss.initGaussian(scale);
506 MV kernel(Shape2(gauss.width(), gauss.height()), &gauss[gauss.upperLeft()]);
507 convolveFFT(in, kernel, out);
508
509 shouldEqualSequenceTolerance(out.data(), out.data()+out.size(),
510 ref.data(), 1e-14);
511
512 Kernel2D<double> gauss2;
513 gauss2.initGaussian(2.0*scale);
514 MV kernel2(Shape2(gauss2.width(), gauss2.height()), &gauss2[gauss2.upperLeft()]);
515
516 MV kernels[] = { kernel, kernel2 };
517 MV outs[] = { out, out2 };
518 convolveFFTMany(in, kernels, kernels+2, outs);
519
520 shouldEqualSequenceTolerance(out.data(), out.data()+out.size(),
521 ref.data(), 1e-14);
522 shouldEqualSequenceTolerance(out2.data(), out2.data()+out2.size(),
523 ref2.data(), 1e-14);
524 }
525
testConvolveFFTComplexMultiFFTTest526 void testConvolveFFTComplex()
527 {
528 typedef MultiArrayView<2, double> MV;
529 typedef MultiArrayView<2, FFTWComplex<double> > MVC;
530
531 ImageImportInfo info("ghouse.gif");
532 Shape2 s(info.width(), info.height());
533 DArray2 in(s), out(s), ref(s), ref2(s);
534 importImage(info, destImage(in));
535
536 double scale = 2.0;
537 Kernel2D<double> gauss, gauss2;
538 gauss.initGaussian(scale);
539 gauss2.initGaussian(2.0*scale);
540
541 convolveImage(srcImageRange(in), destImage(ref), kernel2d(gauss));
542 convolveImage(srcImageRange(in), destImage(ref2), kernel2d(gauss2));
543
544 MV kernel(Shape2(gauss.width(), gauss.height()), &gauss[gauss.upperLeft()]);
545 MV kernel2(Shape2(gauss2.width(), gauss2.height()), &gauss2[gauss2.upperLeft()]);
546
547 //test complex double 2D convolution with spatial domain kernels
548
549 MultiArray<2, FFTWComplex<double> > inc(in);
550 MultiArray<2, FFTWComplex<double> > outc(inc.shape()), outc2(inc.shape());
551 MultiArray<2, FFTWComplex<double> > kernelc(kernel), kernelc2(kernel2);
552 convolveFFTComplex(inc, kernelc, outc, false);
553
554 copyMultiArray(srcMultiArrayRange(outc, FFTWRealAccessor<double>()), destMultiArray(out));
555
556 shouldEqualSequenceTolerance(out.data(), out.data()+out.size(),
557 ref.data(), 1e-14);
558
559 outc.init(0.0);
560
561 MVC kernels[] = { kernelc, kernelc2 };
562 MVC outs[] = { outc, outc2 };
563 convolveFFTComplexMany(inc, kernels, kernels+2, outs, false);
564
565 copyMultiArray(srcMultiArrayRange(outc, FFTWRealAccessor<double>()), destMultiArray(out));
566 shouldEqualSequenceTolerance(out.data(), out.data()+out.size(),
567 ref.data(), 1e-14);
568 copyMultiArray(srcMultiArrayRange(outc2, FFTWRealAccessor<double>()), destMultiArray(out));
569 shouldEqualSequenceTolerance(out.data(), out.data()+out.size(),
570 ref2.data(), 1e-14);
571
572 //test complex double 2D convolution with Fourier domain kernels
573
574 gauss.setBorderTreatment(BORDER_TREATMENT_WRAP);
575 gauss2.setBorderTreatment(BORDER_TREATMENT_WRAP);
576
577 convolveImage(srcImageRange(in), destImage(ref), kernel2d(gauss));
578 convolveImage(srcImageRange(in), destImage(ref2), kernel2d(gauss2));
579
580 MultiArray<2, FFTWComplex<double> > kernelf(in.shape()), kernelf2(in.shape());
581
582 detail::fftEmbedKernel(kernelc, kernelf);
583 fourierTransform(kernelf, kernelf);
584 moveDCToCenter(kernelf);
585
586 detail::fftEmbedKernel(kernelc2, kernelf2);
587 fourierTransform(kernelf2, kernelf2);
588 moveDCToCenter(kernelf2);
589
590 outc.init(0.0);
591 convolveFFTComplex(inc, kernelf, outc, true);
592
593 copyMultiArray(srcMultiArrayRange(outc, FFTWRealAccessor<double>()), destMultiArray(out));
594
595 shouldEqualSequenceTolerance(out.data(), out.data()+out.size(),
596 ref.data(), 1e-14);
597
598 outc.init(0.0);
599 outc2.init(0.0);
600
601 MVC kernelsf[] = { kernelf, kernelf2 };
602 convolveFFTComplexMany(inc, kernelsf, kernelsf+2, outs, true);
603
604 copyMultiArray(srcMultiArrayRange(outc, FFTWRealAccessor<double>()), destMultiArray(out));
605 shouldEqualSequenceTolerance(out.data(), out.data()+out.size(),
606 ref.data(), 1e-14);
607 copyMultiArray(srcMultiArrayRange(outc2, FFTWRealAccessor<double>()), destMultiArray(out));
608 shouldEqualSequenceTolerance(out.data(), out.data()+out.size(),
609 ref2.data(), 1e-14);
610 #if 0
611 // test complex float 3D convolution
612 // compile-only test, currently disabled because it requires libfftwf
613
614 Shape3 sc(100,100,100);
615 vigra::MultiArray<3, FFTWComplex<float> > incf(sc);
616 vigra::MultiArray<3, FFTWComplex<float> > outcf(sc);
617 vigra::MultiArray<3, FFTWComplex<float> > kernelcf(Shape3(10,10,10));
618 convolveFFTComplex(incf, kernelcf, outcf, false);
619 #endif
620 }
621
testConvolveFourierKernelMultiFFTTest622 void testConvolveFourierKernel()
623 {
624 ImageImportInfo info("ghouse.gif");
625 Shape2 s(info.width(), info.height());
626 DArray2 in(s), out(s), out2(s), out3(s), out4(s), ref(s), tmp(s);
627 importImage(info, destImage(in));
628
629 double scale = 2.0;
630 gaussianSmoothing(srcImageRange(in), destImage(ref), scale);
631
632 Shape2 paddedShape = fftwBestPaddedShapeR2C(s + Shape2(16)),
633 kernelShape(paddedShape);
634 kernelShape[0] = kernelShape[0] / 2 + 1;
635 Shape2 center = div(kernelShape, Shape2::value_type(2));
636 center[0] = 0;
637
638 CArray2 kernel(kernelShape);
639
640 for(int y=0; y<kernelShape[1]; ++y)
641 {
642 for(int x=0; x<kernelShape[0]; ++x)
643 {
644 double xx = 2.0 * M_PI * (x - center[0]) / paddedShape[0];
645 double yy = 2.0 * M_PI * (y - center[1]) / paddedShape[1];
646 double r2 = sq(xx) + sq(yy);
647 kernel(x,y) = std::exp(-0.5 * sq(scale) * r2);
648 }
649 }
650 convolveFFT(in, kernel, out);
651
652 shouldEqualSequenceTolerance(out.data(), out.data()+out.size(),
653 ref.data(), 1e-2);
654
655 Kernel1D<double> gauss, grad;
656 gauss.initGaussian(scale);
657 grad.initGaussianDerivative(scale, 1);
658
659 separableConvolveX(srcImageRange(in), destImage(tmp), kernel1d(grad));
660 separableConvolveY(srcImageRange(tmp), destImage(ref), kernel1d(gauss));
661
662 CArray2 kernel2(kernelShape);
663 for(int y=0; y<kernelShape[1]; ++y)
664 {
665 for(int x=0; x<kernelShape[0]; ++x)
666 {
667 double xx = 2.0 * M_PI * (x - center[0]) / paddedShape[0];
668 double yy = 2.0 * M_PI * (y - center[1]) / paddedShape[1];
669 double r2 = sq(xx) + sq(yy);
670 kernel2(x,y) = C(0, xx*std::exp(-0.5 * sq(scale) * r2));
671 }
672 }
673 convolveFFT(in, kernel2, out2);
674
675 ref -= out2;
676
677 FindMinMax<double> minmax;
678 FindAverage<double> average;
679 inspectImage(srcImageRange(ref), minmax);
680 inspectImage(srcImageRange(ref), average);
681
682 should(std::max(minmax.max, -minmax.min) < 0.2);
683 should(average.average() < 0.001);
684
685 MultiArrayView<2, C> kernels[] = { kernel, kernel2 };
686 MultiArrayView<2, R> outs[] = { out3, out4 };
687 convolveFFTMany(in, kernels, kernels+2, outs);
688
689 shouldEqualSequenceTolerance(out.data(), out.data()+out.size(),
690 out3.data(), 1e-15);
691 shouldEqualSequenceTolerance(out2.data(), out2.data()+out2.size(),
692 out4.data(), 1e-15);
693 }
694 };
695
696 struct FFTWTestSuite
697 : public test_suite
698 {
699
FFTWTestSuiteFFTWTestSuite700 FFTWTestSuite()
701 : test_suite("FFTWTest")
702 {
703
704 add( testCase(&FFTWComplexTest::testConstruction));
705 add( testCase(&FFTWComplexTest::testComparison));
706 add( testCase(&FFTWComplexTest::testArithmetic));
707 add( testCase(&FFTWComplexTest::testAccessor));
708 add( testCase(&FFTWComplexTest::testForwardBackwardTrans));
709 add( testCase(&FFTWComplexTest::testRearrangeQuadrants));
710
711 add( testCase(&MultiFFTTest::testFFTShift));
712 add( testCase(&MultiFFTTest::testFFT2D));
713 add( testCase(&MultiFFTTest::testFFT3D));
714 add( testCase(&MultiFFTTest::testPadding));
715 add( testCase(&MultiFFTTest::testConvolveFFT));
716 add( testCase(&MultiFFTTest::testConvolveFFTComplex));
717 add( testCase(&MultiFFTTest::testConvolveFourierKernel));
718 }
719 };
720
721 int initializing;
722
723 struct CompareFunctor
724 {
725 double sumDifference_;
726
CompareFunctorCompareFunctor727 CompareFunctor(): sumDifference_(0) {}
728
operator ()CompareFunctor729 void operator()(const fftw_real &a, const fftw_real &b)
730 { sumDifference_+= abs(a-b); }
731
operator ()CompareFunctor732 double operator()()
733 { return sumDifference_; }
734 };
735
736 struct GaborTests
737 {
738 typedef MultiArrayView<2, fftw_real> RView;
739
740 ImageImportInfo info;
741 int w, h;
742 BasicImage<fftw_real> image;
743
GaborTestsGaborTests744 GaborTests()
745 : info("ghouse.gif"),
746 w(info.width()), h(info.height()),
747 image(w, h)
748 {
749 importImage(info, destImage(image));
750 }
751
752 template<class Iterator, class Accessor>
checkImageGaborTests753 void checkImage(triple<Iterator, Iterator, Accessor> src, const char *filename)
754 {
755 if (initializing)
756 {
757 exportImage(src, ImageExportInfo(filename));
758 cout << "wrote " << filename << endl;
759 }
760 else
761 {
762 ImageImportInfo info(filename);
763 FImage image(info.width(), info.height());
764 importImage(info, destImage(image));
765
766 CompareFunctor cmp;
767
768 inspectTwoImages(src, srcImage(image), cmp);
769 cout << "difference to " << filename << ": " << cmp() << endl;
770 shouldEqualTolerance(cmp(), 0.0, 1e-4);
771 }
772 }
773
testImagesGaborTests774 void testImages()
775 {
776 BasicImage<fftw_real> filter(w, h);
777 int directionCount= 8;
778 int dir= 1, scale= 1;
779 double angle = dir * M_PI / directionCount;
780 double centerFrequency = 3.0/8.0 / VIGRA_CSTD::pow(2.0,scale);
781
782 createGaborFilter(RView(filter),
783 angle, centerFrequency,
784 angularGaborSigma(directionCount, centerFrequency),
785 radialGaborSigma(centerFrequency));
786
787 checkImage(srcImageRange(filter), "filter.xv");
788
789 cout << "Applying filter...\n";
790 FFTWComplexImage result(w, h);
791 applyFourierFilter(srcImageRange(image), srcImage(filter), destImage(result));
792 checkImage(srcImageRange(result, FFTWMagnitudeAccessor<>()), "gaborresult.xv");
793
794 BasicImage<fftw_real> realPart(w, h);
795 applyFourierFilter(srcImageRange(image), srcImage(filter), destImage(realPart));
796
797 CompareFunctor cmp;
798 inspectTwoImages(srcImageRange(result, FFTWRealAccessor<>()),
799 srcImage(realPart), cmp);
800 cout << "difference between real parts: " << cmp() << endl;
801 shouldEqualTolerance(cmp(), 0.0, 1e-4);
802
803 cout << "testing vector results..\n";
804 FVector2Image vectorResult(w,h);
805 applyFourierFilter(srcImageRange(image), srcImage(filter),
806 destImage(vectorResult));
807 CompareFunctor iCmp;
808 inspectTwoImages(srcImageRange(result, FFTWImaginaryAccessor<>()),
809 srcImage(vectorResult,
810 VectorComponentAccessor<FVector2Image::PixelType>(1)),
811 iCmp);
812 cout << "difference between imaginary parts: " << iCmp() << endl;
813 shouldEqualTolerance(cmp(), 0.0, 1e-4);
814
815 cout << "applying on ROI...\n";
816 BasicImage<fftw_real> bigImage(w+20, h+20);
817 BasicImage<fftw_real>::Iterator bigUL= bigImage.upperLeft() + Diff2D(5, 5);
818 copyImage(srcImageRange(image), destIter(bigUL));
819 applyFourierFilter(srcIterRange(bigUL, bigUL + Diff2D(w, h)),
820 srcImage(filter), destImage(result));
821 checkImage(srcImageRange(result, FFTWMagnitudeAccessor<>()), "gaborresult.xv");
822 #if 0
823 cout << "Creating plans with measurement...\n";
824 fftwnd_plan forwardPlan=
825 fftw2d_create_plan(h, w, FFTW_FORWARD, FFTW_MEASURE );
826 fftwnd_plan backwardPlan=
827 fftw2d_create_plan(h, w, FFTW_BACKWARD, FFTW_MEASURE | FFTW_IN_PLACE);
828
829 cout << "Applying again...\n";
830 applyFourierFilter(srcImageRange(image), srcImage(filter), destImage(result),
831 forwardPlan, backwardPlan);
832
833 checkImage(srcImageRange(result, FFTWMagnitudeAccessor()), "gaborresult.xv");
834 #endif
835 }
836
testFamilyGaborTests837 void testFamily()
838 {
839 cout << "testing 8x2 GaborFilterFamily...\n";
840 GaborFilterFamily<FImage> filters(w, h, 8, 2);
841 ImageArray<FFTWComplexImage> results((unsigned int)filters.size(), filters.imageSize());
842 applyFourierFilterFamily(srcImageRange(image), filters, results);
843 checkImage(srcImageRange(results[filters.filterIndex(1,1)], FFTWMagnitudeAccessor<>()),
844 "gaborresult.xv");
845
846 ImageArray<FImage> realParts((unsigned int)filters.size(), filters.imageSize());
847 applyFourierFilterFamily(srcImageRange(image), filters, realParts);
848
849 CompareFunctor cmp;
850 inspectTwoImages(srcImageRange(results[3], FFTWRealAccessor<>()),
851 srcImage(realParts[3]), cmp);
852 cout << "difference between real parts: " << cmp() << endl;
853 shouldEqualTolerance(cmp(), 0.0, 1e-4);
854 }
855 };
856
857 struct GaborTestSuite
858 : public test_suite
859 {
GaborTestSuiteGaborTestSuite860 GaborTestSuite()
861 : test_suite("FFTWWithGaborTest")
862 {
863 add(testCase(&GaborTests::testImages));
864 add(testCase(&GaborTests::testFamily));
865 }
866 };
867
main(int argc,char ** argv)868 int main(int argc, char **argv)
869 {
870 initializing= argc>1 && std::string(argv[1]) == "init"; // global variable,
871 // initializing means: don't compare with image files but create them
872
873 FFTWTestSuite fftwTest;
874
875 int failed = fftwTest.run(testsToBeExecuted(argc, argv));
876
877 std::cout << fftwTest.report() << std::endl;
878
879 GaborTestSuite gaborTest;
880
881 failed = failed || gaborTest.run(testsToBeExecuted(argc, argv));
882
883 std::cout << gaborTest.report() << std::endl;
884
885 return (failed != 0);
886 }
887