mxlib
c++ tools for analyzing astronomical data and other tasks by Jared R. Males. [git repo]
Loading...
Searching...
No Matches
imageXCorrDiscrete_test.cpp
Go to the documentation of this file.
1/** \file imageXCorrDiscrete_test.cpp
2 * \brief Tests discrete image cross-correlation.
3 */
4#include "../../catch2/catch.hpp"
5
6#include <cmath>
7#include <vector>
8#include <Eigen/Dense>
9
10#define MX_NO_ERROR_REPORTS
11
12#include "../../../include/math/func/gaussian.hpp"
13#include "../../../include/improc/imageXCorrDiscrete.hpp"
14#include "../../../include/improc/eigenCube.hpp"
15
16/** \cond
17 * Explicit instantiation compile-checks discrete correlation with the canonical double-precision Eigen image and
18 * emits its header-defined methods for coverage accounting in this test translation unit.
19 */
20template class mx::improc::imageXCorrDiscrete<mx::improc::eigenImage<double>>;
21/** \endcond */
22
23/** Recovering Gaussian displacement and centroid
24 *
25 * Verify imageXCorrDiscrete interpolation and imageCenterOfLight calculation.
26 *
27 */
28/**
29 * \ingroup imageXCorrDiscrete_unit_tests
30 */
31TEST_CASE( "Recover Gaussian displacement and centroid", "[improc::imageXCorrDiscrete]" )
32{
33 GIVEN( "two Gaussians" )
34 {
35 WHEN( "the target is shifted from a centered reference" )
36 {
37 mx::improc::eigenImage<double> im0, im2;
38 im0.resize( 64, 64 );
39 im2.resize( 64, 64 );
40
41 mx::math::func::gaussian2D<double>( im0.data(), im0.rows(), im0.cols(), 0., 1.0, 31.5, 31.5, 2 );
42 mx::math::func::gaussian2D<double>( im2.data(), im2.rows(), im2.cols(), 0., 1.0, 31.5 + 4, 31.5 + 4, 2 );
43
44 double x, y;
45 mx::improc::imageXCorrDiscrete<mx::improc::eigenImage<double>> xcf( 5 );
46
47 mx::improc::eigenImage<double> refIm = im0.block( 10, 10, im0.rows() - 20, im0.cols() - 20 );
48 REQUIRE( xcf.refIm( refIm ) == 0 );
49 xcf.m_peakMethod = mx::improc::xcorrPeakMethod::centroid;
50
51 REQUIRE( xcf( x, y, im2 ) == 0 );
52
53 int xPeak = -1;
54 int yPeak = -1;
55 xcf.ccIm().maxCoeff( &xPeak, &yPeak );
56 REQUIRE( xPeak == 9 );
57 REQUIRE( yPeak == 9 );
58 REQUIRE_THAT( x, Catch::Matchers::WithinAbs( 4, 1e-8 ) );
59 REQUIRE_THAT( y, Catch::Matchers::WithinAbs( 4, 1e-8 ) );
60 REQUIRE( std::isfinite( x ) );
61 REQUIRE( std::isfinite( y ) );
62 }
63
64 WHEN( "the target has an asymmetric signed shift" )
65 {
66 mx::improc::eigenImage<double> centered( 64, 64 );
67 mx::improc::eigenImage<double> shifted( 64, 64 );
68 mx::math::func::gaussian2D<double>(
69 centered.data(), centered.rows(), centered.cols(), 0., 1.0, 31.5, 31.5, 2 );
70 mx::math::func::gaussian2D<double>(
71 shifted.data(), shifted.rows(), shifted.cols(), 0., 1.0, 35.5, 28.5, 2 );
72
73 mx::improc::imageXCorrDiscrete<mx::improc::eigenImage<double>> xcf( 5 );
74 mx::improc::eigenImage<double> reference = centered.block( 10, 10, 44, 44 );
75 REQUIRE( xcf.refIm( reference ) == 0 );
76 xcf.m_peakMethod = mx::improc::xcorrPeakMethod::centroid;
77
78 double x = 0;
79 double y = 0;
80 REQUIRE( xcf( x, y, shifted ) == 0 );
81
82 int xPeak = -1;
83 int yPeak = -1;
84 xcf.ccIm().maxCoeff( &xPeak, &yPeak );
85 REQUIRE( xPeak == 9 );
86 REQUIRE( yPeak == 2 );
87 REQUIRE_THAT( x, Catch::Matchers::WithinAbs( 4, 1e-8 ) );
88 REQUIRE_THAT( y, Catch::Matchers::WithinAbs( -3, 1e-8 ) );
89 REQUIRE( std::isfinite( x ) );
90 REQUIRE( std::isfinite( y ) );
91 }
92
93 WHEN( "the target has zero shift" )
94 {
95 mx::improc::eigenImage<double> target( 64, 64 );
96 mx::math::func::gaussian2D<double>(
97 target.data(), target.rows(), target.cols(), 0., 1.0, 31.5, 31.5, 2 );
98
99 mx::improc::imageXCorrDiscrete<mx::improc::eigenImage<double>> xcf( 5 );
100 mx::improc::eigenImage<double> reference = target.block( 10, 10, 44, 44 );
101 REQUIRE( xcf.refIm( reference ) == 0 );
102 xcf.m_peakMethod = mx::improc::xcorrPeakMethod::centroid;
103
104 double x = 0;
105 double y = 0;
106 REQUIRE( xcf( x, y, target ) == 0 );
107
108 int xPeak = -1;
109 int yPeak = -1;
110 xcf.ccIm().maxCoeff( &xPeak, &yPeak );
111 REQUIRE( xPeak == 5 );
112 REQUIRE( yPeak == 5 );
113 REQUIRE_THAT( x, Catch::Matchers::WithinAbs( 0, 1e-8 ) );
114 REQUIRE_THAT( y, Catch::Matchers::WithinAbs( 0, 1e-8 ) );
115 REQUIRE( std::isfinite( x ) );
116 REQUIRE( std::isfinite( y ) );
117 }
118 WHEN( "geometric quarter" )
119 {
120 mx::improc::eigenImage<double> im;
121 im.resize( 64, 64 );
122
123 mx::math::func::gaussian2D<double>( im.data(), im.rows(), im.cols(), 0., 1.0, 15.5, 15.5, 2 );
124
125 double x, y;
126 mx::improc::imageCenterOfLight( x, y, im );
127
128 REQUIRE_THAT( x, Catch::Matchers::WithinAbs( 15.5, 1e-8 ) );
129 REQUIRE_THAT( y, Catch::Matchers::WithinAbs( 15.5, 1e-8 ) );
130 }
131
132 WHEN( "automatic lags span only valid target blocks" )
133 {
134 mx::improc::eigenImage<double> target( 9, 10 );
135 mx::math::func::gaussian2D<double>(
136 target.data(), target.rows(), target.cols(), 0., 1.0, 4.0, 4.5, 1.5 );
137
138 mx::improc::eigenImage<double> reference = target.block( 2, 2, 4, 5 );
139 mx::improc::imageXCorrDiscrete<mx::improc::eigenImage<double>> xcf;
140 REQUIRE( xcf.refIm( reference ) == 0 );
141
142 double x, y;
143 REQUIRE( xcf( x, y, target ) == 0 );
144 REQUIRE( xcf.ccIm().rows() == 5 );
145 REQUIRE( xcf.ccIm().cols() == 5 );
146 int xPeak = -1;
147 int yPeak = -1;
148 xcf.ccIm().maxCoeff( &xPeak, &yPeak );
149 REQUIRE( xPeak == 2 );
150 REQUIRE( yPeak == 2 );
151 REQUIRE( std::isfinite( x ) );
152 REQUIRE( std::isfinite( y ) );
153 }
154
155 WHEN( "interpolation uses the lag range allowed by the image dimensions" )
156 {
157 mx::improc::eigenImage<double> centered( 64, 64 );
158 mx::improc::eigenImage<double> shifted( 64, 64 );
159 mx::math::func::gaussian2D<double>(
160 centered.data(), centered.rows(), centered.cols(), 0., 1.0, 31.5, 31.5, 2 );
161 mx::math::func::gaussian2D<double>(
162 shifted.data(), shifted.rows(), shifted.cols(), 0., 1.0, 32.5, 31.5, 2 );
163
164 mx::improc::imageXCorrDiscrete<mx::improc::eigenImage<double>> xcf( 5 );
165 mx::improc::eigenImage<double> reference = centered.block( 2, 2, 60, 60 );
166 REQUIRE( xcf.refIm( reference ) == 0 );
167 xcf.m_peakMethod = mx::improc::xcorrPeakMethod::interp;
168
169 double x = 0;
170 double y = 0;
171 REQUIRE( xcf( x, y, shifted ) == 0 );
172 REQUIRE( xcf.ccIm().rows() == 5 );
173 REQUIRE( xcf.ccIm().cols() == 5 );
174 REQUIRE( xcf.magIm().rows() == 41 );
175 REQUIRE( xcf.magIm().cols() == 41 );
176 REQUIRE_THAT( x, Catch::Matchers::WithinAbs( 1, 0.11 ) );
177 REQUIRE_THAT( y, Catch::Matchers::WithinAbs( 0, 0.11 ) );
178 REQUIRE( std::isfinite( x ) );
179 REQUIRE( std::isfinite( y ) );
180 }
181
182 WHEN( "interpolation falls back safely for a three-sample lag range" )
183 {
184 mx::improc::eigenImage<double> centered( 64, 64 );
185 mx::improc::eigenImage<double> shifted( 64, 64 );
186 mx::math::func::gaussian2D<double>(
187 centered.data(), centered.rows(), centered.cols(), 0., 1.0, 31.5, 31.5, 2 );
188 mx::math::func::gaussian2D<double>(
189 shifted.data(), shifted.rows(), shifted.cols(), 0., 1.0, 32.5, 31.5, 2 );
190
191 mx::improc::imageXCorrDiscrete<mx::improc::eigenImage<double>> xcf( 5 );
192 mx::improc::eigenImage<double> reference = centered.block( 1, 1, 62, 62 );
193 REQUIRE( xcf.refIm( reference ) == 0 );
194 xcf.m_peakMethod = mx::improc::xcorrPeakMethod::interp;
195
196 double x = 0;
197 double y = 0;
198 REQUIRE( xcf( x, y, centered ) == 0 );
199 REQUIRE( xcf.ccIm().rows() == 3 );
200 REQUIRE( xcf.ccIm().cols() == 3 );
201 REQUIRE( xcf.magIm().size() == 0 );
202 REQUIRE_THAT( x, Catch::Matchers::WithinAbs( 0, 1e-8 ) );
203 REQUIRE_THAT( y, Catch::Matchers::WithinAbs( 0, 1e-8 ) );
204
205 REQUIRE( xcf( x, y, shifted ) == 0 );
206 REQUIRE( xcf.magIm().size() == 0 );
207 REQUIRE_THAT( x, Catch::Matchers::WithinAbs( 1, 1e-8 ) );
208 REQUIRE_THAT( y, Catch::Matchers::WithinAbs( 0, 1e-8 ) );
209 REQUIRE( std::isfinite( x ) );
210 REQUIRE( std::isfinite( y ) );
211 }
212
213 WHEN( "interpolation falls back safely for a boundary correlation peak" )
214 {
215 mx::improc::eigenImage<double> centered( 64, 64 );
216 mx::improc::eigenImage<double> shiftedNegative( 64, 64 );
217 mx::improc::eigenImage<double> shiftedPositive( 64, 64 );
218 mx::math::func::gaussian2D<double>(
219 centered.data(), centered.rows(), centered.cols(), 0., 1.0, 31.5, 31.5, 2 );
220 mx::math::func::gaussian2D<double>( shiftedNegative.data(),
221 shiftedNegative.rows(),
222 shiftedNegative.cols(),
223 0.,
224 1.0,
225 29.5,
226 31.5,
227 2 );
228 mx::math::func::gaussian2D<double>( shiftedPositive.data(),
229 shiftedPositive.rows(),
230 shiftedPositive.cols(),
231 0.,
232 1.0,
233 33.5,
234 31.5,
235 2 );
236
237 mx::improc::imageXCorrDiscrete<mx::improc::eigenImage<double>> xcf( 5 );
238 mx::improc::eigenImage<double> reference = centered.block( 2, 2, 60, 60 );
239 REQUIRE( xcf.refIm( reference ) == 0 );
240 xcf.m_peakMethod = mx::improc::xcorrPeakMethod::interp;
241
242 double x = 0;
243 double y = 0;
244 REQUIRE( xcf( x, y, shiftedNegative ) == 0 );
245 REQUIRE( xcf.ccIm().rows() == 5 );
246 REQUIRE( xcf.ccIm().cols() == 5 );
247 REQUIRE( xcf.magIm().size() == 0 );
248 REQUIRE_THAT( x, Catch::Matchers::WithinAbs( -2, 1e-8 ) );
249 REQUIRE_THAT( y, Catch::Matchers::WithinAbs( 0, 1e-8 ) );
250
251 REQUIRE( xcf( x, y, shiftedPositive ) == 0 );
252 REQUIRE( xcf.magIm().size() == 0 );
253 REQUIRE_THAT( x, Catch::Matchers::WithinAbs( 2, 1e-8 ) );
254 REQUIRE_THAT( y, Catch::Matchers::WithinAbs( 0, 1e-8 ) );
255 REQUIRE( std::isfinite( x ) );
256 REQUIRE( std::isfinite( y ) );
257 }
258
259 WHEN( "interpolation falls back safely for a negative correlation surface" )
260 {
261 mx::improc::eigenImage<double> target( 7, 7 );
262 mx::improc::eigenImage<double> reference( 3, 3 );
263 for( int row = 0; row < target.rows(); ++row )
264 {
265 for( int col = 0; col < target.cols(); ++col )
266 {
267 target( row, col ) = row + 0.2 * col + 1.2 * std::sin( row + 0.3 * col );
268 }
269 }
270 for( int row = 0; row < reference.rows(); ++row )
271 {
272 for( int col = 0; col < reference.cols(); ++col )
273 {
274 reference( row, col ) = -( row + 0.2 * col );
275 }
276 }
277
278 mx::improc::imageXCorrDiscrete<mx::improc::eigenImage<double>> xcf( 5 );
279 REQUIRE( xcf.refIm( reference ) == 0 );
280 xcf.m_peakMethod = mx::improc::xcorrPeakMethod::interp;
281
282 double x = 0;
283 double y = 0;
284 REQUIRE( xcf( x, y, target ) == 0 );
285 REQUIRE( xcf.ccIm().rows() == 5 );
286 REQUIRE( xcf.ccIm().cols() == 5 );
287
288 int xPeak = -1;
289 int yPeak = -1;
290 REQUIRE( xcf.ccIm().maxCoeff( &xPeak, &yPeak ) < 0 );
291 REQUIRE( xPeak == 1 );
292 REQUIRE( yPeak == 3 );
293 REQUIRE( xcf.magIm().size() == 0 );
294 REQUIRE_THAT( x, Catch::Matchers::WithinAbs( -1, 1e-8 ) );
295 REQUIRE_THAT( y, Catch::Matchers::WithinAbs( 1, 1e-8 ) );
296 REQUIRE( std::isfinite( x ) );
297 REQUIRE( std::isfinite( y ) );
298 }
299
300 WHEN( "a mask differs from the reference in exactly one dimension" )
301 {
302 mx::improc::eigenImage<double> reference( 5, 6 );
303 reference.setOnes();
304
305 mx::improc::imageXCorrDiscrete<mx::improc::eigenImage<double>> xcf;
306
307 mx::improc::eigenImage<double> wrongRows( 4, 6 );
308 wrongRows.setOnes();
309 REQUIRE( xcf.maskIm( wrongRows ) == 0 );
310 REQUIRE( xcf.refIm( reference ) == -1 );
311
312 mx::improc::eigenImage<double> wrongCols( 5, 7 );
313 wrongCols.setOnes();
314 REQUIRE( xcf.maskIm( wrongCols ) == 0 );
315 REQUIRE( xcf.refIm( reference ) == -1 );
316 }
317 }
318}
mx::improc::imageXCorrDiscrete
Find the optimum shift to align two images using the discrete cross correlation.
Definition imageXCorrDiscrete.hpp:73
mx::improc::imageXCorrDiscrete::maskIm
int maskIm(const ccImT &mask)
Set the mask image.
Definition imageXCorrDiscrete.hpp:256
mx::improc::imageXCorrDiscrete::refIm
int refIm(const ccImT &im0)
Set the reference image.
Definition imageXCorrDiscrete.hpp:271
mx::improc::imageXCorrDiscrete::magIm
const ccImT & magIm()
Get a reference to the magnified image.
Definition imageXCorrDiscrete.hpp:315
mx::improc::imageXCorrDiscrete::ccIm
const ccImT & ccIm()
Get a reference to the cross correlation image.
Definition imageXCorrDiscrete.hpp:309
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::improc::imageCenterOfLight
int imageCenterOfLight(typename imageT::Scalar &x, typename imageT::Scalar &y, const imageT &im)
Calculate the center of light of an image.
Definition imageUtils.hpp:374
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("Recover Gaussian displacement and centroid", "[improc::imageXCorrDiscrete]")
Definition imageXCorrDiscrete_test.cpp:31
imageXCorrDiscrete.hpp
A class to register images using the discrete cross correlation.