Line data Source code
1 : /** \file imageXCorrDiscrete.hpp
2 : * \brief A class to register images using the discrete cross correlation.
3 : * \ingroup image_processing_files
4 : * \author Jared R. Males (jaredmales@gmail.com)
5 : *
6 : */
7 :
8 : //***********************************************************************//
9 : // Copyright 2020 Jared R. Males (jaredmales@gmail.com)
10 : //
11 : // This file is part of mxlib.
12 : //
13 : // mxlib is free software: you can redistribute it and/or modify
14 : // it under the terms of the GNU General Public License as published by
15 : // the Free Software Foundation, either version 3 of the License, or
16 : // (at your option) any later version.
17 : //
18 : // mxlib is distributed in the hope that it will be useful,
19 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
20 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 : // GNU General Public License for more details.
22 : //
23 : // You should have received a copy of the GNU General Public License
24 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
25 : //***********************************************************************//
26 :
27 : #ifndef imageXCorrDiscrete_hpp
28 : #define imageXCorrDiscrete_hpp
29 :
30 : #include <algorithm>
31 : #include <cmath>
32 : #include <limits>
33 :
34 : #include "../mxlib.hpp"
35 : #include "../math/fit/fitGaussian.hpp"
36 : #include "imageUtils.hpp"
37 :
38 : namespace mx
39 : {
40 : namespace improc
41 : {
42 :
43 : enum class xcorrPeakMethod
44 : {
45 : centroid,
46 : gaussfit,
47 : interp
48 : };
49 :
50 : /// Find the optimum shift to align two images using the discrete cross correlation.
51 : /** The reference image must be smaller than the target image. Both the reference image and the section of
52 : * target image being analyzed are mean subtracted and variance normalized. An optional mask can be supplied,
53 : * which limits the pixels used for the mean and variance calculation.
54 : *
55 : * Typical usage will be to set the mask, then the reference, then repeatedly call operator() to
56 : * determine the shifts for a sequence of imaages. No new heap allocations take place on these calls
57 : * to operator(), and the reference image is not re-normalized on each call.
58 : *
59 : * The shift is reported in pixels such that if the mxlib imageShift function is used
60 : * to shift the input image by the negative of the shifts, it will align with the
61 : * reference at the center of the array.
62 : *
63 : * Three peak finding methods are provided. xcorrPeakMethod::centroid uses center of light,
64 : * xcorrPeakMethod::gaussfit uses Gaussian centroiding, and xcorrPeakMethod::interp uses interpolation
65 : * to find the peak to a given tolerance.
66 : *
67 : * \tparam _ccImT is the Eigen-like array type used for image processing. See typedefs.
68 : *
69 : * \ingroup image_reg
70 : */
71 : template <class _ccImT>
72 : class imageXCorrDiscrete
73 : {
74 : public:
75 : typedef _ccImT ccImT; ///< the Eigen-like array type used for image processing
76 : typedef typename _ccImT::Scalar Scalar; ///< the scalar type of the image type
77 :
78 : protected:
79 : public:
80 : /** \name Working Memory
81 : * @{
82 : */
83 :
84 : ccImT m_refIm; ///< The normalized reference image.
85 :
86 : ccImT m_maskIm; ///< Mask image to use, may be needed for proper normalization even if refIm has 0 mask applied.
87 :
88 : bool m_haveMask{ false };
89 :
90 : ccImT m_normIm; ///< The normalized image.
91 :
92 : ccImT m_ccIm; ///< The cross-correlation image
93 :
94 : ccImT m_magIm; ///< The magnified image, or empty when interpolation falls back to the discrete peak.
95 :
96 : ///@}
97 :
98 : math::fit::fitGaussian2D<mx::math::fit::gaussian2D_gen_fitter<Scalar>> m_fitter;
99 :
100 : int m_maxLag{ 0 }; ///< The maximum lag to consider in the initial cross-correlation.
101 :
102 : Scalar m_tol{ 0.1 }; ///< The tolerance of the interpolated-magnified image, in pixels.
103 :
104 : Scalar m_magSize{ 0 }; ///< Magnified size of the ccIm when using interp. Set as function of m_tol and m_maxLag.
105 :
106 : public:
107 : xcorrPeakMethod m_peakMethod{ xcorrPeakMethod::centroid };
108 :
109 : public:
110 : /// Default c'tor
111 : imageXCorrDiscrete();
112 :
113 : /// Construct seeting maxLag.
114 : explicit imageXCorrDiscrete( int maxLag );
115 :
116 : /// Get the current maximum lag
117 : /**
118 : * \returns the current value of m_maxLag
119 : */
120 : int maxLag();
121 :
122 : /// Set the maximum lag
123 : void maxLag( int ml /**< [in] the new maximum lag */ );
124 :
125 : /// Get the tolerance of the interpolated-magnified image, in pixels.
126 : /**
127 : * \returns the current value of m_tol.
128 : */
129 : Scalar tol();
130 :
131 : /// Set the tolerance of the interpolated-magnified image, in pixels.
132 : void tol( Scalar nt /**< [in] The new value of the interpolation tolerance. */ );
133 :
134 : /// Set the size of the cross-correlation images.
135 : /** This resizes all working memory.
136 : *
137 : * \returns 0 on success
138 : * \returns -1 on error
139 : */
140 : int resize( int nrows, ///< [in] the number of rows in the images to register
141 : int ncols ///< [in] the number of columns in the images to register
142 : );
143 :
144 : /// Set the mask image
145 : /**
146 : * \returns 0 on success
147 : * \returns -1 on error
148 : */
149 : int maskIm( const ccImT &mask /**< [in] the new mask image */ );
150 :
151 : /// Get a reference to the mask image.
152 : /**
153 : * \returns a const referance to the mask image.
154 : */
155 : const ccImT &maskIm();
156 :
157 : /// Set the reference image
158 : /** Normalizes the reference image by mean subtraction and variance division. Applies
159 : * the mask first if supplied.
160 : *
161 : * \returns 0 on success
162 : * \returns -1 on error
163 : */
164 : int refIm( const ccImT &im0 );
165 :
166 : /// Get a reference to the reference image.
167 : /**
168 : * \returns a const referent to m_refIm.
169 : */
170 : const ccImT &refIm();
171 :
172 : /// Get a reference to the normalized image.
173 : /**
174 : * \returns a const referent to m_normIm.
175 : */
176 : const ccImT &normIm();
177 :
178 : /// Get a reference to the cross correlation image.
179 : /**
180 : * \returns a const referent to m_ccIm.
181 : */
182 : const ccImT &ccIm();
183 :
184 : /// Get a reference to the magnified image.
185 : /**
186 : * \returns a const referent to m_magIm.
187 : */
188 : const ccImT &magIm();
189 :
190 : /// Conduct the cross correlation to a specified tolerance
191 : /** Cubic interpolation requires at least five correlation samples in each
192 : * dimension and cannot interpolate a maximum on the outer boundary. It
193 : * also pads unsupported magnified borders with zero, which would dominate
194 : * a non-positive correlation surface. In these cases, this falls back to
195 : * the discrete correlation maximum and leaves `magIm()` empty.
196 : *
197 : * \returns 0 on success
198 : * \returns -1 on error
199 : */
200 : template <class imT>
201 : int operator()( Scalar &xShift, ///< [out] the x shift of im w.r.t. im0, in pixels
202 : Scalar &yShift, ///< [out] the y shift of im w.r.t. im0, in pixels
203 : const imT &im ///< [in] the image to cross-correlate with the reference
204 : );
205 : };
206 :
207 : template <class ccImT>
208 2 : imageXCorrDiscrete<ccImT>::imageXCorrDiscrete()
209 : {
210 2 : }
211 :
212 : template <class ccImT>
213 7 : imageXCorrDiscrete<ccImT>::imageXCorrDiscrete( int mL )
214 : {
215 7 : maxLag( mL );
216 7 : }
217 :
218 : template <class ccImT>
219 0 : int imageXCorrDiscrete<ccImT>::maxLag()
220 : {
221 0 : return m_maxLag;
222 : }
223 :
224 : template <class ccImT>
225 7 : void imageXCorrDiscrete<ccImT>::maxLag( int ml )
226 : {
227 7 : m_maxLag = ml;
228 7 : tol( m_tol );
229 7 : }
230 :
231 : template <class ccImT>
232 0 : typename ccImT::Scalar imageXCorrDiscrete<ccImT>::tol()
233 : {
234 0 : return m_tol;
235 : }
236 :
237 : template <class ccImT>
238 7 : void imageXCorrDiscrete<ccImT>::tol( Scalar nt )
239 : {
240 7 : m_tol = nt;
241 7 : m_magSize = 0;
242 :
243 7 : if( !math::isFinite( nt ) || nt <= 0 || m_maxLag <= 0 )
244 : {
245 0 : return;
246 : }
247 :
248 7 : m_magSize = ceil( ( ( 2. * m_maxLag + 1 ) - 1.0 ) / nt ) + 1;
249 :
250 7 : Scalar mag = ( m_magSize - 1.0 ) / ( ( 2. * m_maxLag + 1 ) - 1.0 );
251 :
252 7 : m_tol = 1.0 / mag;
253 : }
254 :
255 : template <class ccImT>
256 2 : int imageXCorrDiscrete<ccImT>::maskIm( const ccImT &mask )
257 : {
258 2 : m_maskIm = mask;
259 2 : m_haveMask = true;
260 :
261 2 : return 0;
262 : }
263 :
264 : template <class ccImT>
265 0 : const ccImT &imageXCorrDiscrete<ccImT>::maskIm()
266 : {
267 0 : return m_maskIm;
268 : }
269 :
270 : template <class ccImT>
271 10 : int imageXCorrDiscrete<ccImT>::refIm( const ccImT &im )
272 : {
273 10 : ccImT im0;
274 10 : if( m_haveMask )
275 : {
276 2 : if( im.rows() != m_maskIm.rows() || im.cols() != m_maskIm.cols() )
277 : {
278 4 : internal::mxlib_error_report( error_t::sizeerr, "reference and mask are not the same size" );
279 2 : return -1;
280 : }
281 :
282 0 : im0 = im * m_maskIm;
283 : }
284 : else
285 : {
286 8 : im0 = im;
287 : }
288 :
289 8 : Scalar m = imageMean<Scalar>( im0 );
290 8 : Scalar v = imageVariance( im0, m );
291 8 : m_refIm = ( im0 - m ) / sqrt( v );
292 :
293 8 : return 0;
294 10 : }
295 :
296 : template <class ccImT>
297 0 : const ccImT &imageXCorrDiscrete<ccImT>::refIm()
298 : {
299 0 : return m_refIm;
300 : }
301 :
302 : template <class ccImT>
303 0 : const ccImT &imageXCorrDiscrete<ccImT>::normIm()
304 : {
305 0 : return m_normIm;
306 : }
307 :
308 : template <class ccImT>
309 15 : const ccImT &imageXCorrDiscrete<ccImT>::ccIm()
310 : {
311 15 : return m_ccIm;
312 : }
313 :
314 : template <class ccImT>
315 7 : const ccImT &imageXCorrDiscrete<ccImT>::magIm()
316 : {
317 7 : return m_magIm;
318 : }
319 :
320 : template <class ccImT>
321 : template <class imT>
322 10 : int imageXCorrDiscrete<ccImT>::operator()( Scalar &xShift, Scalar &yShift, const imT &im )
323 : {
324 10 : if( im.rows() <= m_refIm.rows() )
325 : {
326 0 : internal::mxlib_error_report(error_t::sizeerr,"reference must be smaller than target image (rows)" );
327 0 : return -1;
328 : }
329 :
330 10 : if( im.cols() <= m_refIm.cols() )
331 : {
332 0 : internal::mxlib_error_report(error_t::sizeerr, "reference must be smaller than target image (cols)" );
333 0 : return -1;
334 : }
335 : // 16/4 15/4 16/3 15/3
336 10 : int maxLag_r = ( im.rows() - m_refIm.rows() ) / 2; // 6 5 6 6
337 10 : int maxLag_c = ( im.cols() - m_refIm.cols() ) / 2;
338 :
339 10 : if( maxLag_r > m_maxLag && m_maxLag != 0 )
340 3 : maxLag_r = m_maxLag;
341 10 : if( maxLag_c > m_maxLag && m_maxLag != 0 )
342 3 : maxLag_c = m_maxLag;
343 :
344 10 : m_ccIm.resize( 2 * maxLag_r + 1, 2 * maxLag_c + 1 );
345 :
346 74 : for( int rL = -maxLag_r; rL <= maxLag_r; ++rL )
347 : {
348 570 : for( int cL = -maxLag_c; cL <= maxLag_c; ++cL )
349 : {
350 : // 16/4 15/4 16/3 15/3
351 506 : int r0 = 0.5 * im.rows() + rL - 0.5 * m_refIm.rows(); // 0-15 0-13 0-14 0-14
352 506 : int c0 = 0.5 * im.cols() + cL - 0.5 * m_refIm.cols();
353 :
354 506 : if( m_haveMask )
355 : {
356 0 : m_normIm = im.block( r0, c0, m_refIm.rows(), m_refIm.cols() ) * m_maskIm;
357 : }
358 : else
359 : {
360 506 : m_normIm = im.block( r0, c0, m_refIm.rows(), m_refIm.cols() );
361 : }
362 :
363 506 : Scalar m = imageMean<Scalar>( m_normIm );
364 506 : Scalar sv = sqrt( imageVariance( m_normIm, m ) );
365 506 : if( sv <= 0 )
366 : {
367 0 : m_ccIm( maxLag_r + rL, maxLag_c + cL ) = 0;
368 : }
369 : else
370 : {
371 506 : m_normIm = ( m_normIm - m ) / sv;
372 506 : m_ccIm( maxLag_r + rL, maxLag_c + cL ) = ( m_normIm * m_refIm ).sum();
373 : // ds9Interface ds9( m_normIm, 1);
374 : }
375 : }
376 : }
377 :
378 10 : if( m_peakMethod == xcorrPeakMethod::gaussfit )
379 : {
380 : int xLag0, yLag0;
381 0 : Scalar pk = m_ccIm.maxCoeff( &xLag0, &yLag0 );
382 0 : Scalar mn = m_ccIm.minCoeff();
383 0 : m_fitter.setArray( m_ccIm.data(), m_ccIm.rows(), m_ccIm.cols() );
384 0 : m_fitter.setGuess( mn, pk, xLag0, yLag0, 0.2 * m_ccIm.rows(), 0.2 * m_ccIm.cols(), 0 );
385 0 : m_fitter.fit();
386 :
387 0 : xShift = m_fitter.x0() - maxLag_r;
388 0 : yShift = m_fitter.y0() - maxLag_c;
389 : }
390 10 : else if( m_peakMethod == xcorrPeakMethod::interp )
391 : {
392 6 : if( !math::isFinite( m_tol ) || m_tol <= 0 )
393 : {
394 0 : internal::mxlib_error_report( error_t::invalidarg, "interpolation tolerance must be finite and positive" );
395 5 : return -1;
396 : }
397 :
398 : int xPeak, yPeak;
399 6 : const Scalar peakValue = m_ccIm.maxCoeff( &xPeak, &yPeak );
400 10 : if( m_ccIm.rows() < 5 || m_ccIm.cols() < 5 || !( peakValue > 0 ) || !math::isFinite( peakValue ) ||
401 10 : xPeak == 0 || yPeak == 0 || xPeak == m_ccIm.rows() - 1 || yPeak == m_ccIm.cols() - 1 )
402 : {
403 5 : m_magIm.resize( 0, 0 );
404 5 : m_magSize = 0;
405 5 : xShift = xPeak - maxLag_r;
406 5 : yShift = yPeak - maxLag_c;
407 5 : return 0;
408 : }
409 :
410 1 : const Scalar magRowsScalar = ceil( ( m_ccIm.rows() - 1.0 ) / m_tol ) + 1.0;
411 1 : const Scalar magColsScalar = ceil( ( m_ccIm.cols() - 1.0 ) / m_tol ) + 1.0;
412 2 : if( !math::isFinite( magRowsScalar ) || !math::isFinite( magColsScalar ) || magRowsScalar <= 1 ||
413 3 : magColsScalar <= 1 || magRowsScalar > std::numeric_limits<int>::max() ||
414 1 : magColsScalar > std::numeric_limits<int>::max() )
415 : {
416 0 : internal::mxlib_error_report( error_t::sizeerr, "invalid interpolated cross-correlation dimensions" );
417 0 : return -1;
418 : }
419 :
420 1 : const int magRows = static_cast<int>( magRowsScalar );
421 1 : const int magCols = static_cast<int>( magColsScalar );
422 1 : m_magSize = std::max( magRows, magCols );
423 1 : m_magIm.resize( magRows, magCols );
424 1 : imageMagnify( m_magIm, m_ccIm, cubicConvolTransform<Scalar>() );
425 : int x, y;
426 1 : m_magIm.maxCoeff( &x, &y );
427 1 : const Scalar xScale = ( m_ccIm.rows() - 1.0 ) / ( m_magIm.rows() - 1.0 );
428 1 : const Scalar yScale = ( m_ccIm.cols() - 1.0 ) / ( m_magIm.cols() - 1.0 );
429 1 : xShift = x * xScale - maxLag_r;
430 1 : yShift = y * yScale - maxLag_c;
431 : }
432 : else
433 : {
434 : int xPeak, yPeak;
435 4 : m_ccIm.maxCoeff( &xPeak, &yPeak );
436 :
437 4 : int centroidHalfWidth = std::min( xPeak, static_cast<int>( m_ccIm.rows() ) - 1 - xPeak );
438 4 : centroidHalfWidth = std::min( centroidHalfWidth, yPeak );
439 4 : centroidHalfWidth =
440 4 : std::min( centroidHalfWidth, static_cast<int>( m_ccIm.cols() ) - 1 - yPeak );
441 :
442 4 : if( centroidHalfWidth == 0 )
443 : {
444 0 : xShift = xPeak - maxLag_r;
445 0 : yShift = yPeak - maxLag_c;
446 0 : return 0;
447 : }
448 :
449 4 : const int centroidWidth = 2 * centroidHalfWidth + 1;
450 4 : const int x0 = xPeak - centroidHalfWidth;
451 4 : const int y0 = yPeak - centroidHalfWidth;
452 4 : m_magIm = m_ccIm.block( x0, y0, centroidWidth, centroidWidth );
453 4 : m_magIm -= m_magIm.minCoeff();
454 :
455 4 : const Scalar centroidWeight = m_magIm.sum();
456 4 : if( !( centroidWeight > 0 ) || !math::isFinite( centroidWeight ) )
457 : {
458 0 : xShift = xPeak - maxLag_r;
459 0 : yShift = yPeak - maxLag_c;
460 0 : return 0;
461 : }
462 :
463 : Scalar x, y;
464 4 : imageCenterOfLight( x, y, m_magIm );
465 4 : xShift = x + x0 - maxLag_r;
466 4 : yShift = y + y0 - maxLag_c;
467 : }
468 :
469 5 : return 0;
470 : }
471 :
472 : } // namespace improc
473 : } // namespace mx
474 :
475 : #endif // imageXCorrDiscrete_hpp
|
