mxlib
c++ tools for analyzing astronomical data and other tasks by Jared R. Males. [git repo]
Loading...
Searching...
No Matches
imageXCorrFFT_test.cpp
Go to the documentation of this file.
1/** \file imageXCorrFFT_test.cpp
2 * \brief Tests FFT-based image cross-correlation.
3 */
4#include "../../catch2/catch.hpp"
5
6#include <vector>
7#include <Eigen/Dense>
8
9#define MX_NO_ERROR_REPORTS
10
11#include "../../../include/math/func/gaussian.hpp"
12#include "../../../include/improc/imageXCorrFFT.hpp"
13#include "../../../include/improc/eigenCube.hpp"
14
15/** \cond
16 * Explicit instantiation compile-checks FFT correlation with the canonical double-precision Eigen image and emits
17 * its header-defined methods for coverage accounting in this test translation unit.
18 */
19template class mx::improc::imageXCorrFFT<mx::improc::eigenImage<double>>;
20/** \endcond */
21
22/** centroiding Gaussians with center of light
23 *
24 * Verify center of light calculation
25 *
26 */
27/**
28 * \ingroup imageXCorrFFT_unit_tests
29 */
30TEST_CASE( "Image cross-correlation with FFT using center of light", "[improc::imageXCorrFFT]" )
31{
32 GIVEN( "two Gaussians" )
33 {
34 WHEN( "ref at geometric center, equal even sizes, shift=(+4,+4)" )
35 {
36 double xshift = 4;
37 double yshift = 4;
38 int xsz = 64;
39 int ysz = 64;
40
41 mx::improc::eigenImage<double> im0, im2;
42 im0.resize( xsz, ysz );
43 im2.resize( xsz, ysz );
44
45 double xcen = 0.5 * ( im0.rows() - 1 );
46 double ycen = 0.5 * ( im0.cols() - 1 );
47
48 mx::math::func::gaussian2D<double>( im0.data(), im0.rows(), im0.cols(), 0., 1.0, xcen, ycen, 2 );
49 mx::math::func::gaussian2D<double>( im2.data(),
50 im2.rows(),
51 im2.cols(),
52 0.,
53 1.0,
54 xcen + xshift,
55 ycen + yshift,
56 2 );
57
58 double x, y, peak;
59 mx::improc::imageXCorrFFT<mx::improc::eigenImage<double>> xcf;
60
61 mx::improc::eigenImage<double> refIm = im0; //.block(10,10, im0.rows()-11, im0.cols()-11);
62 xcf.peakMethod( mx::improc::xcorrPeakMethod::centerOfLight );
63 xcf.maxLag( 20 );
64 xcf.refIm( refIm );
65 xcf( x, y, peak, im2 );
66
67 REQUIRE( x == Approx( xshift ) );
68 REQUIRE( y == Approx( yshift ) );
69 REQUIRE( peak > 0 );
70 }
71 WHEN( "ref at geometric center, equal odd sizes, shift=(+4,+4)" )
72 {
73 double xshift = 4;
74 double yshift = 4;
75 int xsz = 65;
76 int ysz = 65;
77
78 mx::improc::eigenImage<double> im0, im2;
79 im0.resize( xsz, ysz );
80 im2.resize( xsz, ysz );
81
82 double xcen = 0.5 * ( im0.rows() - 1 );
83 double ycen = 0.5 * ( im0.cols() - 1 );
84
85 mx::math::func::gaussian2D<double>( im0.data(), im0.rows(), im0.cols(), 0., 1.0, xcen, ycen, 2 );
86
87 mx::math::func::gaussian2D<double>( im2.data(),
88 im2.rows(),
89 im2.cols(),
90 0.,
91 1.0,
92 xcen + xshift,
93 ycen + yshift,
94 2 );
95
96 double x, y, peak;
97 mx::improc::imageXCorrFFT<mx::improc::eigenImage<double>> xcf;
98 xcf.peakMethod( mx::improc::xcorrPeakMethod::centerOfLight );
99 xcf.maxLag( 20 );
100
101 mx::improc::eigenImage<double> refIm = im0; //.block(10,10, im0.rows()-11, im0.cols()-11);
102 xcf.refIm( refIm );
103
104 xcf( x, y, peak, im2 );
105
106 REQUIRE( x == Approx( xshift ) );
107 REQUIRE( y == Approx( yshift ) );
108 REQUIRE( peak > 0 );
109 }
110 WHEN( "ref at geometric center, equal even sizes, shift=(-3.6,+2.25)" )
111 {
112 double xshift = -3.6;
113 double yshift = 2.25;
114 int xsz = 64;
115 int ysz = 64;
116
117 mx::improc::eigenImage<double> im0, im2;
118 im0.resize( xsz, ysz );
119 im2.resize( xsz, ysz );
120
121 double xcen = 0.5 * ( im0.rows() - 1 );
122 double ycen = 0.5 * ( im0.cols() - 1 );
123
124 mx::math::func::gaussian2D<double>( im0.data(), im0.rows(), im0.cols(), 0., 1.0, xcen, ycen, 2 );
125 mx::math::func::gaussian2D<double>( im2.data(),
126 im2.rows(),
127 im2.cols(),
128 0.,
129 1.0,
130 xcen + xshift,
131 ycen + yshift,
132 2 );
133
134 double x, y, peak;
135 mx::improc::imageXCorrFFT<mx::improc::eigenImage<double>> xcf;
136 xcf.peakMethod( mx::improc::xcorrPeakMethod::centerOfLight );
137 xcf.maxLag( 20 );
138
139 mx::improc::eigenImage<double> refIm = im0; //.block(10,10, im0.rows()-11, im0.cols()-11);
140 xcf.refIm( refIm );
141
142 xcf( x, y, peak, im2 );
143
144 REQUIRE( x == Approx( xshift ) );
145 REQUIRE( y == Approx( yshift ) );
146 REQUIRE( peak > 0 );
147 }
148 WHEN( "ref at geometric center, equal odd sizes, shift=(+1.3,-0.6)" )
149 {
150 double xshift = 1.3;
151 double yshift = -0.6;
152 int xsz = 65;
153 int ysz = 65;
154
155 mx::improc::eigenImage<double> im0, im2;
156 im0.resize( xsz, ysz );
157 im2.resize( xsz, ysz );
158
159 double xcen = 0.5 * ( im0.rows() - 1 );
160 double ycen = 0.5 * ( im0.cols() - 1 );
161
162 mx::math::func::gaussian2D<double>( im0.data(), im0.rows(), im0.cols(), 0., 1.0, xcen, ycen, 2 );
163 mx::math::func::gaussian2D<double>( im2.data(),
164 im2.rows(),
165 im2.cols(),
166 0.,
167 1.0,
168 xcen + xshift,
169 ycen + yshift,
170 2 );
171
172 double x, y, peak;
173 mx::improc::imageXCorrFFT<mx::improc::eigenImage<double>> xcf;
174 xcf.peakMethod( mx::improc::xcorrPeakMethod::centerOfLight );
175 xcf.maxLag( 20 );
176
177 mx::improc::eigenImage<double> refIm = im0; //.block(10,10, im0.rows()-11, im0.cols()-11);
178 xcf.refIm( refIm );
179
180 xcf( x, y, peak, im2 );
181
182 REQUIRE( x == Approx( xshift ) );
183 REQUIRE( y == Approx( yshift ) );
184 REQUIRE( peak > 0 );
185 }
186 }
187}
188
189/** centroiding by magnification
190 *
191 * Verify magnification peak finding
192 *
193 */
194/**
195 * \ingroup imageXCorrFFT_unit_tests
196 */
197TEST_CASE( "Image cross-correlation with FFT using magnification peak finding", "[improc::imageXCorrFFT]" )
198{
199 GIVEN( "two Gaussians" )
200 {
201 WHEN( "ref at geometric center, equal even sizes, shift=(+4,+4)" )
202 {
203 double xshift = 4;
204 double yshift = 4;
205 int xsz = 64;
206 int ysz = 64;
207
208 mx::improc::eigenImage<double> im0, im2;
209 im0.resize( xsz, ysz );
210 im2.resize( xsz, ysz );
211
212 double xcen = 0.5 * ( im0.rows() - 1 );
213 double ycen = 0.5 * ( im0.cols() - 1 );
214
215 mx::math::func::gaussian2D<double>( im0.data(), im0.rows(), im0.cols(), 0., 1.0, xcen, ycen, 2 );
216 mx::math::func::gaussian2D<double>( im2.data(),
217 im2.rows(),
218 im2.cols(),
219 0.,
220 1.0,
221 xcen + xshift,
222 ycen + yshift,
223 2 );
224
225 double x, y, peak;
226 mx::improc::imageXCorrFFT<mx::improc::eigenImage<double>> xcf;
227 xcf.peakMethod( mx::improc::xcorrPeakMethod::interpPeak );
228 xcf.maxLag( 5 );
229
230 mx::improc::eigenImage<double> refIm = im0; //.block(10,10, im0.rows()-11, im0.cols()-11);
231
232 xcf.refIm( refIm );
233 xcf( x, y, peak, im2 );
234
235 REQUIRE( x == Approx( xshift ) );
236 REQUIRE( y == Approx( yshift ) );
237 REQUIRE( peak > 0 );
238 }
239 WHEN( "ref at geometric center, equal odd sizes, shift=(+4,+4)" )
240 {
241 double xshift = 4;
242 double yshift = 4;
243 int xsz = 65;
244 int ysz = 65;
245
246 mx::improc::eigenImage<double> im0, im2;
247 im0.resize( xsz, ysz );
248 im2.resize( xsz, ysz );
249
250 double xcen = 0.5 * ( im0.rows() - 1 );
251 double ycen = 0.5 * ( im0.cols() - 1 );
252
253 mx::math::func::gaussian2D<double>( im0.data(), im0.rows(), im0.cols(), 0., 1.0, xcen, ycen, 2 );
254 mx::math::func::gaussian2D<double>( im2.data(),
255 im2.rows(),
256 im2.cols(),
257 0.,
258 1.0,
259 xcen + xshift,
260 ycen + yshift,
261 2 );
262
263 double x, y, peak;
264 mx::improc::imageXCorrFFT<mx::improc::eigenImage<double>> xcf;
265
266 mx::improc::eigenImage<double> refIm = im0; //.block(10,10, im0.rows()-11, im0.cols()-11);
267
268 xcf.peakMethod( mx::improc::xcorrPeakMethod::interpPeak );
269 xcf.maxLag( 5 );
270
271 xcf.refIm( refIm );
272 xcf( x, y, peak, im2 );
273
274 REQUIRE( x == Approx( xshift ) );
275 REQUIRE( y == Approx( yshift ) );
276 REQUIRE( peak > 0 );
277 }
278 WHEN( "ref at geometric center, equal even sizes, shift=(-3.6,+2.3)" )
279 {
280 double xshift = -3.6;
281 double yshift = 2.3;
282 int xsz = 64;
283 int ysz = 64;
284
285 mx::improc::eigenImage<double> im0, im2;
286 im0.resize( xsz, ysz );
287 im2.resize( xsz, ysz );
288
289 double xcen = 0.5 * ( im0.rows() - 1 );
290 double ycen = 0.5 * ( im0.cols() - 1 );
291
292 mx::math::func::gaussian2D<double>( im0.data(), im0.rows(), im0.cols(), 0., 1.0, xcen, ycen, 2 );
293 mx::math::func::gaussian2D<double>( im2.data(),
294 im2.rows(),
295 im2.cols(),
296 0.,
297 1.0,
298 xcen + xshift,
299 ycen + yshift,
300 2 );
301
302 double x, y, peak;
303 mx::improc::imageXCorrFFT<mx::improc::eigenImage<double>> xcf;
304 xcf.peakMethod( mx::improc::xcorrPeakMethod::interpPeak );
305 xcf.tol( 0.1 );
306
307 mx::improc::eigenImage<double> refIm = im0;
308 xcf.refIm( refIm );
309
310 xcf( x, y, peak, im2 );
311
312 REQUIRE( x == Approx( xshift ) );
313 REQUIRE( y == Approx( yshift ) );
314 REQUIRE( peak > 0 );
315 }
316 WHEN( "ref at geometric center, equal odd sizes, shift=(+1.3,-0.6)" )
317 {
318 double xshift = 1.3;
319 double yshift = -0.6;
320 int xsz = 65;
321 int ysz = 65;
322
323 mx::improc::eigenImage<double> im0, im2;
324 im0.resize( xsz, ysz );
325 im2.resize( xsz, ysz );
326
327 double xcen = 0.5 * ( im0.rows() - 1 );
328 double ycen = 0.5 * ( im0.cols() - 1 );
329
330 mx::math::func::gaussian2D<double>( im0.data(), im0.rows(), im0.cols(), 0., 1.0, xcen, ycen, 2 );
331 mx::math::func::gaussian2D<double>( im2.data(),
332 im2.rows(),
333 im2.cols(),
334 0.,
335 1.0,
336 xcen + xshift,
337 ycen + yshift,
338 2 );
339
340 double x, y, peak;
341 mx::improc::imageXCorrFFT<mx::improc::eigenImage<double>> xcf;
342 xcf.peakMethod( mx::improc::xcorrPeakMethod::interpPeak );
343 xcf.maxLag( 5 );
344
345 mx::improc::eigenImage<double> refIm = im0; //.block(10,10, im0.rows()-11, im0.cols()-11);
346 xcf.refIm( refIm );
347
348 xcf( x, y, peak, im2 );
349
350 REQUIRE( x == Approx( xshift ) );
351 REQUIRE( y == Approx( yshift ) );
352 REQUIRE( peak > 0 );
353 }
354 }
355}
356
357/// Verify imageXCorrFFT displacement recovery with Gaussian peak fitting.
358/**
359 * \ingroup imageXCorrFFT_unit_tests
360 */
361TEST_CASE( "Image cross-correlation with FFT using Gaussian peak fit", "[improc::imageXCorrFFT]" )
362{
363 GIVEN( "two Gaussians" )
364 {
365 WHEN( "ref at geometric center, equal even sizes, shift=(+4,+4)" )
366 {
367 double xshift = 4;
368 double yshift = 4;
369 int xsz = 64;
370 int ysz = 64;
371
372 mx::improc::eigenImage<double> im0, im2;
373 im0.resize( xsz, ysz );
374 im2.resize( xsz, ysz );
375
376 double xcen = 0.5 * ( im0.rows() - 1 );
377 double ycen = 0.5 * ( im0.cols() - 1 );
378
379 mx::math::func::gaussian2D<double>( im0.data(), im0.rows(), im0.cols(), 0., 1.0, xcen, ycen, 2 );
380 mx::math::func::gaussian2D<double>( im2.data(),
381 im2.rows(),
382 im2.cols(),
383 0.,
384 1.0,
385 xcen + xshift,
386 ycen + yshift,
387 2 );
388
389 double x, y, peak;
390 mx::improc::imageXCorrFFT<mx::improc::eigenImage<double>> xcf;
391
392 mx::improc::eigenImage<double> refIm = im0;
393 xcf.peakMethod( mx::improc::xcorrPeakMethod::gaussFit );
394 xcf.maxLag( 32 );
395 xcf.refIm( refIm );
396 xcf( x, y, peak, im2 );
397
398 REQUIRE( x == Approx( xshift ) );
399 REQUIRE( y == Approx( yshift ) );
400 REQUIRE( peak > 0 );
401 }
402 WHEN( "ref at geometric center, equal odd sizes, shift=(+4,+4)" )
403 {
404 double xshift = 4;
405 double yshift = 4;
406 int xsz = 65;
407 int ysz = 65;
408
409 mx::improc::eigenImage<double> im0, im2;
410 im0.resize( xsz, ysz );
411 im2.resize( xsz, ysz );
412
413 double xcen = 0.5 * ( im0.rows() - 1 );
414 double ycen = 0.5 * ( im0.cols() - 1 );
415
416 mx::math::func::gaussian2D<double>( im0.data(), im0.rows(), im0.cols(), 0., 1.0, xcen, ycen, 2 );
417 mx::math::func::gaussian2D<double>( im2.data(),
418 im2.rows(),
419 im2.cols(),
420 0.,
421 1.0,
422 xcen + xshift,
423 ycen + yshift,
424 2 );
425
426 double x, y, peak;
427 mx::improc::imageXCorrFFT<mx::improc::eigenImage<double>> xcf;
428 xcf.peakMethod( mx::improc::xcorrPeakMethod::gaussFit );
429 xcf.maxLag( 32 );
430
431 mx::improc::eigenImage<double> refIm = im0; //.block(10,10, im0.rows()-11, im0.cols()-11);
432 xcf.refIm( refIm );
433
434 xcf( x, y, peak, im2 );
435
436 REQUIRE( x == Approx( xshift ) );
437 REQUIRE( y == Approx( yshift ) );
438 REQUIRE( peak > 0 );
439 }
440 WHEN( "ref at geometric center, equal even sizes, shift=(-3.6,+2.25)" )
441 {
442 double xshift = -3.6;
443 double yshift = 2.25;
444 int xsz = 64;
445 int ysz = 64;
446
447 mx::improc::eigenImage<double> im0, im2;
448 im0.resize( xsz, ysz );
449 im2.resize( xsz, ysz );
450
451 double xcen = 0.5 * ( im0.rows() - 1 );
452 double ycen = 0.5 * ( im0.cols() - 1 );
453
454 mx::math::func::gaussian2D<double>( im0.data(), im0.rows(), im0.cols(), 0., 1.0, xcen, ycen, 2 );
455 mx::math::func::gaussian2D<double>( im2.data(),
456 im2.rows(),
457 im2.cols(),
458 0.,
459 1.0,
460 xcen + xshift,
461 ycen + yshift,
462 2 );
463
464 double x, y, peak;
465 mx::improc::imageXCorrFFT<mx::improc::eigenImage<double>> xcf;
466
467 mx::improc::eigenImage<double> refIm = im0; //.block(10,10, im0.rows()-11, im0.cols()-11);
468 xcf.peakMethod( mx::improc::xcorrPeakMethod::gaussFit );
469 xcf.maxLag( 32 );
470 xcf.refIm( refIm );
471
472 xcf( x, y, peak, im2 );
473
474 REQUIRE( x == Approx( xshift ) );
475 REQUIRE( y == Approx( yshift ) );
476 REQUIRE( peak > 0 );
477 }
478 WHEN( "ref at geometric center, equal odd sizes, shift=(+1.3,-0.6)" )
479 {
480 double xshift = 1.3;
481 double yshift = -0.6;
482 int xsz = 65;
483 int ysz = 65;
484
485 mx::improc::eigenImage<double> im0, im2;
486 im0.resize( xsz, ysz );
487 im2.resize( xsz, ysz );
488
489 double xcen = 0.5 * ( im0.rows() - 1 );
490 double ycen = 0.5 * ( im0.cols() - 1 );
491
492 mx::math::func::gaussian2D<double>( im0.data(), im0.rows(), im0.cols(), 0., 1.0, xcen, ycen, 2 );
493 mx::math::func::gaussian2D<double>( im2.data(),
494 im2.rows(),
495 im2.cols(),
496 0.,
497 1.0,
498 xcen + xshift,
499 ycen + yshift,
500 2 );
501
502 double x, y, peak;
503 mx::improc::imageXCorrFFT<mx::improc::eigenImage<double>> xcf;
504 xcf.peakMethod( mx::improc::xcorrPeakMethod::gaussFit );
505 xcf.maxLag( 32 );
506
507 mx::improc::eigenImage<double> refIm = im0; //.block(10,10, im0.rows()-11, im0.cols()-11);
508 xcf.refIm( refIm );
509
510 xcf( x, y, peak, im2 );
511
512 REQUIRE( x == Approx( xshift ) );
513 REQUIRE( y == Approx( yshift ) );
514 REQUIRE( peak > 0 );
515 }
516 }
517}
mx::improc::imageXCorrFFT
Find the optimum shift to align two images using the FFT cross correlation.
Definition imageXCorrFFT.hpp:81
mx::improc::imageXCorrFFT::maxLag
void maxLag(int ml)
Set the maximum lag.
Definition imageXCorrFFT.hpp:926
mx::improc::imageXCorrFFT::tol
void tol(realT nt)
Set the tolerance of the interpolated-magnified image, in pixels.
Definition imageXCorrFFT.hpp:945
mx::improc::imageXCorrFFT::refIm
int refIm(const realImageT &im0, realT padFactor)
Set the reference image.
Definition imageXCorrFFT.hpp:809
eigenCube.hpp
An image cube with an Eigen API.
gaussian.hpp
Declarations for utilities related to the Gaussian function.
mx::improc::eigenImage
Eigen::Array< scalarT, -1, -1 > eigenImage
Definition of the eigenImage type, which is an alias for Eigen::Array.
Definition eigenImage.hpp:44
mx::math::func::gaussian2D
realT gaussian2D(const realT x, const realT y, const realT G0, const realT G, const realT x0, const realT y0, const realT sigma)
Find value at position (x,y) of the 2D arbitrarily-centered symmetric Gaussian.
Definition gaussian.hpp:126
TEST_CASE
TEST_CASE("Image cross-correlation with FFT using center of light", "[improc::imageXCorrFFT]")
Definition imageXCorrFFT_test.cpp:30
imageXCorrFFT.hpp
A class to register images using the Fourier cross correlation with a peak fit.