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