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