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