Line data Source code
1 : /** \file imageXCorrFFT.hpp
2 : * \brief A class to register images using the Fourier cross correlation with a peak fit.
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 imageXCorrFFT_hpp
28 : #define imageXCorrFFT_hpp
29 :
30 : #include "../mxlib.hpp"
31 : #include "../math/ft/fftT.hpp"
32 : #include "../math/ft/mftT.hpp"
33 : #include "../math/ft/ftUtils.hpp"
34 : #include "../math/ft/fftwEnvironment.hpp"
35 :
36 : #include "../math/fit/fitGaussian.hpp"
37 :
38 : #include "eigenImage.hpp"
39 : #include "imageUtils.hpp"
40 : #include "imageTransforms.hpp"
41 : #include "imageXCorr.hpp"
42 :
43 : namespace mx
44 : {
45 : namespace improc
46 : {
47 :
48 : /// Find the optimum shift to align two images using the FFT cross correlation.
49 : /** The reference image must be the same size as the target image. Both the reference image and the
50 : * target image being analyzed can be masked, normalized by mean subtraction and variance division,
51 : * and windowed. Seperate masks and windows for the reference and target image are used.
52 : *
53 : * The reference mask, reference window, and normalization flag must be set before setting the reference
54 : * (unless using the defaults). If the reference mask, reference window, or normalization flag are changed
55 : * then the reference must be set again. After configuring the reference, then operator() can be called
56 : * repeatedly to determine the shifts for a sequence of imaages. The mask and window can be changed between
57 : * calls without requiring resetting the reference. No new heap allocations take place on these calls
58 : * to operator(), and the reference image is not re-normalized on each call.
59 : *
60 : * The shift is reported in pixels such that if the mxlib imageShift function is used
61 : * to shift the input image by the negative of the shifts, it will align with the
62 : * reference at the center of the array.
63 : *
64 : * Several peak finding methods are provided:
65 : * - xcorrPeakMethod::centerOfLight uses center of light
66 : * - xcorrPeakMethod::gaussFit uses Gaussian centroiding
67 : * - xcorrPeakMethod::interpPeak uses interpolation to find the peak to a given tolerance.
68 : * - xcorrPeakMethod::none will just return the coordinate of the max pixel
69 : *
70 : * This class uses the real-to-complex transforms of FFTW by default. This should offer
71 : * improved speed and memory use. For testing the complex-to-complex transforms can be
72 : * used by defining `XCFFT_C2C` before including this file.
73 : *
74 : * \tparam _realImageT is the Eigen-like array type used for image processing. Must have
75 : * a real floating point Scalar type (float, double, etc.). See typedefs.
76 : *
77 : * \ingroup image_reg
78 : */
79 : template <class _realImageT>
80 : class imageXCorrFFT
81 : {
82 : public:
83 : typedef _realImageT realImageT; ///< the Eigen-like array type used for image processing
84 :
85 : typedef typename _realImageT::Scalar realT; ///< the scalar type of the image type
86 :
87 : typedef std::complex<realT> complexT; ///< Complex floating point type.
88 :
89 : // defaut is the real-to-complex. Only use C2C if defined:
90 : #ifdef XCFFT_C2C
91 : typedef complexT fftOutT;
92 : typedef complexT fftInT;
93 : #else
94 : typedef realT fftOutT;
95 : typedef complexT fftInT;
96 : #endif
97 :
98 : typedef eigenImage<complexT> complexArrayT; ///< Complex eigen array type with Scalar==complexT
99 :
100 : protected:
101 : /** \name Configuration Member Data
102 : * @{
103 : */
104 : int m_rows{ 0 }; ///< The number of rows in the images
105 :
106 : int m_cols{ 0 }; ///< The number of columns in the images
107 :
108 : realT m_padFactor{ 1 }; ///< The padding factor for the CC
109 :
110 : realT m_padFactorR{ 1 }; /**< The actual padding factor for rows of the CC
111 : set based on the resultant array sizes.*/
112 :
113 : realT m_padFactorC{ 1 }; /**< The actual padding factor for cols of the CC
114 : set based on the resultant array sizes.*/
115 :
116 : int m_rowsPadded{ 0 }; ///< The number of rows in the padded CC
117 :
118 : int m_colsPadded{ 0 }; ///< The number of columns in the padded CC
119 :
120 : realImageT m_refMaskIm; /**< Mask image to use for the reference, may be needed for proper normalization
121 : even if refIm has 0 mask applied.*/
122 :
123 : bool m_haveRefMask{ false }; ///< Flag indicating that a referece mask has been provided.
124 :
125 : realImageT m_maskIm; /**< Mask image to use, may be needed for proper normalization even if
126 : refIm has 0 mask applied.*/
127 :
128 : bool m_haveMask{ false }; ///< Flag indicating that a mask has been provided.
129 :
130 : realImageT m_refWinIm; ///< Window image to use for the reference
131 :
132 : bool m_haveRefWindow{ false }; ///< Flag indicating that a window has been provided for the reference
133 :
134 : realImageT m_winIm; ///< Window image for the target image.
135 :
136 : bool m_haveWindow{ false }; ///< Flag indicating that a window has been provided.
137 :
138 : bool m_normalize{ true }; /**< Flag specifying whether images should be normalized prior to performing
139 : the cross-correlation*/
140 :
141 : realImageT m_refIm; ///< The normalized reference image.
142 :
143 : bool m_refValid{ false }; /**< Flag indicating whether the reference is valid. It can be invalidated
144 : if a mask, window, or normalize is set after the reference is set.*/
145 :
146 : int m_maxLag{ 5 }; ///< The maximum lag to consider in the initial cross-correlation. Default is 5.
147 :
148 : realT m_tol{ 0.1 }; /**< For xcorrPeakMethod::interpPeak, the tolerance of the
149 : interpolated-magnified image, in pixels.*/
150 :
151 : realT m_magSize{ 0 }; /**< Magnified size of the ccIm when using interp.
152 : Set as function of m_tol and m_maxLag.*/
153 :
154 : realT m_mag{ 1 }; /**< The magnification, in principle is 1/m_tol but will vary
155 : slightly based on integer array sizes*/
156 :
157 : /// The peak finding method to use
158 : xcorrPeakMethod m_peakMethod{ xcorrPeakMethod::interpPeak };
159 :
160 : ///@}
161 :
162 : /** \name Working Memory Member Data
163 : * @{
164 : */
165 :
166 : realImageT m_refACIm; ///< The auto-correlation image of the reference.
167 :
168 : realT m_refX0{ 0 }; /**< The x-shift of the reference image to itself using the selected
169 : algorithm, used as coordinate origin*/
170 :
171 : realT m_refY0{ 0 }; /**< The x-shift of the reference image to itself using the selected
172 : algorithm, used as coordinate origin*/
173 :
174 : realT m_refPeak{ 0 }; /**< the x-corr peak of the reference shift to itself */
175 :
176 : realImageT m_normIm; ///< The normalized image.
177 :
178 : realImageT m_ccIm; ///< The cross-correlation image
179 :
180 : realImageT m_magIm; ///< The magnified image, used if m_peakMethod == xcorrPeakMethod::interpPeak
181 :
182 : complexArrayT m_ftIm0; ///< Working memory for the FT of the reference image.
183 :
184 : public:
185 : complexArrayT m_ftWork; ///< Working memory for the FFT.
186 :
187 : #ifdef XCFFT_C2C
188 : complexArrayT m_ftWorkIn; ///< Working memory for the FFT input.
189 : #endif
190 :
191 : complexArrayT m_ftWorkPadded; ///< Working memory for the padded inverse FFT input.
192 :
193 : complexArrayT m_ftWorkPaddedOut; ///< Working memory for the FFT and MFT output.
194 :
195 : complexArrayT m_mtWorkPadded; ///< Working memory for the MFT for xcorrPeakMethod::mftOversamp
196 :
197 : math::ft::fftT<fftOutT, fftInT, 2> m_fft_fwd; ///< FFT object for the forward transform.
198 :
199 : math::ft::fftT<fftInT, fftOutT, 2> m_fft_bwd; ///< FFT object for the backward transform.
200 :
201 : math::ft::mftT<complexT, complexT, 2> m_mft_bwd; /**< MFT object used to oversample the
202 : peak for xcorrPeakMethod::mftOversamp*/
203 :
204 : ///@}
205 :
206 : /** \name Peak Finding Tools
207 : * @{
208 : */
209 : math::fit::fitGaussian2D<mx::math::fit::gaussian2D_gen_fitter<realT>> m_fitter;
210 :
211 : ///@}
212 :
213 : public:
214 : /** \name Construction and Destruction
215 : * @{
216 : */
217 :
218 : /// Default c'tor
219 : imageXCorrFFT();
220 :
221 : ///@}
222 :
223 : /** \name Configuration Member Access
224 : * @{
225 : */
226 :
227 : /// Get the number of rows in the input images
228 : /** This can only be set by a call to resize() or refIm().
229 : *
230 : * \returns the current value of m_rows
231 : */
232 : int rows();
233 :
234 : /// Get the number of columns in the input images
235 : /** This can only be set by a call to resize() or refIm().
236 : *
237 : * \returns the current value of m_cols
238 : */
239 : int cols();
240 :
241 : /// Set the padding factor
242 : /** This can also be set by a call to resize() or refIm().
243 : * The padding factor must be greater than or equal to 1.
244 : *
245 : * If m_rows and m_cols are both greater than 0 this calls resize(int, int, realT).
246 : *
247 : */
248 : int padFactor( realT os /**< [in] the new padding factor*/ );
249 :
250 : /// Get the padding factor
251 : /**
252 : * \returns the current value of m_padFactor
253 : */
254 : realT padFactor();
255 :
256 : /// Get the number of rows in the padded cross-correlation
257 : /** This can only be set by a call to resize() or refIm().
258 : *
259 : * \returns the current value of m_rowsPadded
260 : */
261 : int rowsPadded();
262 :
263 : /// Get the number of rows in the padded cross-correlation
264 : /** This can only be set by a call to resize() or refIm().
265 : *
266 : * \returns the current value of m_colsPadded
267 : */
268 : int colsPadded();
269 :
270 : /// Set the size of the cross-correlation problem based on input image size.
271 : /** This resizes all working memory and conducts fftw planning.
272 : * If there are no changes to size or padding, this returns immediately
273 : * making no changes.
274 : *
275 : * If there are changes to size or padding the reference is invalidated.
276 : *
277 : * If nrows or ncols is 0, this resizes all working memory to 0.
278 : *
279 : * \returns 0 on success
280 : * \returns -1 on error
281 : */
282 : int resize( int nrows, /**< [in] the number of rows in the images to register */
283 : int ncols, /**< [in] the number of columns in the images to register */
284 : realT padFactor /**< [in] the padding factor. Must be `>=` 1. */ );
285 :
286 : /// Set the size of the cross-correlation problem based on input image size.
287 : /** Calls resize(int, int, realT) with the current value of m_padFactor.
288 : *
289 : * \overload
290 : *
291 : * \returns 0 on success
292 : * \returns -1 on error
293 : */
294 : int resize( int nrows, /**< [in] the number of rows in the images to register */
295 : int ncols /**< [in] the number of columns in the images to register */ );
296 :
297 : /// Set the reference mask image
298 : /** This invalidates the reference
299 : * \returns 0 on success
300 : * \returns -1 on error
301 : */
302 : int refMaskIm( const realImageT &mask /**< [in] the new reference mask image */ );
303 :
304 : /// Get a reference to the reference mask image.
305 : /**
306 : * \returns a const reference to the reference mask image.
307 : */
308 : const realImageT &refMaskIm();
309 :
310 : /// Get the value of the haveRefMask flag
311 : /**
312 : * \returns the current value of m_haveRefMask
313 : */
314 : bool haveRefMask();
315 :
316 : /// Set the mask image
317 : /** This does not invalidate the reference
318 : *
319 : * \returns 0 on success
320 : * \returns -1 on error
321 : */
322 : int maskIm( const realImageT &mask /**< [in] the new mask image */ );
323 :
324 : /// Get a reference to the mask image.
325 : /**
326 : * \returns a const referance to the mask image.
327 : */
328 : const realImageT &maskIm();
329 :
330 : /// Get the value of the haveMask flag
331 : /**
332 : * \returns the current value of m_haveMask
333 : */
334 : bool haveMask();
335 :
336 : /// Set the reference window image
337 : /** This invalidates the reference.
338 : *
339 : * \returns 0 on success
340 : * \returns -1 on error
341 : */
342 : int refWinIm( const realImageT &win /**< [in] the new reference window image */ );
343 :
344 : /// Get a reference to the reference window image.
345 : /**
346 : * \returns a const referance to the reference window image.
347 : */
348 : const realImageT &refWinIm();
349 :
350 : /// Get the value of the haveRefWindow flag
351 : /**
352 : * \returns the current value of m_haveRefWindow
353 : */
354 : bool haveRefWindow();
355 :
356 : /// Set the window image
357 : /** This does not invalidate the reference.
358 : *
359 : * \returns 0 on success
360 : * \returns -1 on error
361 : */
362 : int winIm( const realImageT &win /**< [in] the new window image */ );
363 :
364 : /// Get a reference to the window image.
365 : /**
366 : * \returns a const referance to the window image.
367 : */
368 : const realImageT &winIm();
369 :
370 : /// Get the value of the haveWindow flag
371 : /**
372 : * \returns the current value of m_haveWindow
373 : */
374 : bool haveWindow();
375 :
376 : /// Set the normalize flag
377 : /** Changing this invalidates the reference.
378 : */
379 : void normalize( bool no /**< [in] the new normalize flag value */ );
380 :
381 : /// Get the current value of the normalize flag
382 : /**
383 : * \returns the current value of m_normalize
384 : */
385 : bool normalize();
386 :
387 : /// Set the reference image
388 : /** Performs the following steps
389 : * - Applies the mask if supplied.
390 : * - If desired, normalizes the reference image by mean subtraction and variance division.
391 : * - Applies the window if supplied.
392 : * - Stores the result as m_refIm
393 : * - Resizes working memory with a call to resize()
394 : * - Shifts to FFT order and Fourier transforms and conjugates the result.
395 : * - Calls operator() to measure the reference shifts
396 : * - Stores the resultant auto-correlation in m_refACIm
397 : *
398 : * \returns 0 on success
399 : * \returns -1 on error
400 : */
401 : int refIm( const realImageT &im0, /**< [in] The new reference image */
402 : realT padFactor /**< [in] [optional] the padding factor to use */ );
403 :
404 : /// Set the reference image
405 : /** Calls refIm(const realImageT &, realT) with the current value of m_padFactor.
406 : *
407 : * \overload
408 : *
409 : * \returns 0 on success
410 : * \returns -1 on error
411 : */
412 : int refIm( const realImageT &im0 /**< [in] The new reference image */ );
413 :
414 : /// Get a reference to the reference image.
415 : /**
416 : * \returns a const referent to m_refIm.
417 : */
418 : const realImageT &refIm();
419 :
420 : /// Get the value of the reference valid flag
421 : /**
422 : * \returns the value of m_refValid
423 : */
424 : bool refValid();
425 :
426 : /// Set the maximum lag
427 : void maxLag( int ml /**< [in] the new maximum lag */ );
428 :
429 : /// Get the current maximum lag
430 : /**
431 : * \returns the current value of m_maxLag
432 : */
433 : int maxLag();
434 :
435 : /// Set the tolerance of the interpolated-magnified image, in pixels.
436 : void tol( realT nt /**< [in] The new value of the interpolation tolerance. */ );
437 :
438 : /// Get the tolerance of the interpolated-magnified image, in pixels.
439 : /**
440 : * \returns the current value of m_tol.
441 : */
442 : realT tol();
443 :
444 : void peakMethod( xcorrPeakMethod xpm );
445 :
446 : xcorrPeakMethod peakMethod();
447 :
448 : ///@}
449 :
450 : /** \name Working Memory Member Access
451 : * @{
452 : */
453 :
454 : /// Get a reference to the reference auto-correlation image.
455 : /**
456 : * \returns a const referent to m_refACIm.
457 : */
458 : const realImageT &refACIm();
459 :
460 : /// Get the x-shift of the reference image to itself.
461 : /**
462 : * \returns m_refX0
463 : */
464 : realT refX0();
465 :
466 : /// Get the y-shift of the reference image to itself.
467 : /**
468 : * \returns m_refY0
469 : */
470 : realT refY0();
471 :
472 : /// Get a reference to the normalized image.
473 : /**
474 : * \returns a const referent to m_normIm.
475 : */
476 : const realImageT &normIm();
477 :
478 : /// Get a reference to the cross correlation image.
479 : /**
480 : * \returns a const referent to m_ccIm.
481 : */
482 : const realImageT &ccIm();
483 :
484 : /// Get a reference to the magnified image.
485 : /**
486 : * \returns a const referent to m_magIm.
487 : */
488 : const realImageT &magIm();
489 :
490 : ///@}
491 :
492 : /** \name Peak Finding
493 : * @{
494 : */
495 : protected:
496 : void findPeak( realT &xShift, /**< [out] the x shift of im w.r.t. im0, in pixels */
497 : realT &yShift, /**< [out] the y shift of im w.r.t. im0, in pixels */
498 : realT &xcPeak /**< [out] the x-corr peak */
499 : );
500 :
501 : ///@}
502 :
503 : /** \name Cross Correlation Operator
504 : * @{
505 : */
506 :
507 : public:
508 : /// Conduct the cross correlation to a specified tolerance
509 : /**
510 : * \returns 0 on success
511 : * \returns -1 on error
512 : */
513 : template <class imT>
514 : int operator()( realT &xShift, ///< [out] the x shift of im w.r.t. im0, in pixels
515 : realT &yShift, ///< [out] the y shift of im w.r.t. im0, in pixels
516 : realT &xcPeak, ///< [out] the xc peak
517 : const imT &im /**< [in] the image to cross-correlate with the reference*/ );
518 :
519 : /// Conduct the cross correlation to a specified tolerance
520 : /**
521 : * \overload
522 : *
523 : * \returns 0 on success
524 : * \returns -1 on error
525 : */
526 : template <class im0T, class imT>
527 : int operator()( realT &xShift, ///< [out] the x shift of im w.r.t. im0, in pixels
528 : realT &yShift, ///< [out] the y shift of im w.r.t. im0, in pixels
529 : realT &xcPeak, ///< [out] the xc peak
530 : im0T &im0, ///< [in] a new reference image
531 : imT &im ///< [in] the image to cross-correlate with the reference
532 : );
533 :
534 : ///@}
535 : };
536 :
537 : template <class realImageT>
538 12 : imageXCorrFFT<realImageT>::imageXCorrFFT()
539 : {
540 12 : maxLag( m_maxLag ); // this sets up m_magSize.
541 12 : }
542 :
543 : template <class realImageT>
544 : int imageXCorrFFT<realImageT>::rows()
545 : {
546 : return m_rows;
547 : }
548 :
549 : template <class realImageT>
550 : int imageXCorrFFT<realImageT>::cols()
551 : {
552 : return m_cols;
553 : }
554 :
555 : template <class realImageT>
556 : int imageXCorrFFT<realImageT>::padFactor( realT os )
557 : {
558 : if( os < 1 )
559 : {
560 : internal::mxlib_error_report( error_t::invalidarg, "padding factor can't be less than 1" );
561 : return -1;
562 : }
563 :
564 : if( m_rows > 0 && m_cols > 0 )
565 : {
566 : return resize( m_rows, m_cols, os );
567 : }
568 : else
569 : {
570 : m_padFactor = os;
571 : return 0;
572 : }
573 : }
574 :
575 : template <class realImageT>
576 : typename imageXCorrFFT<realImageT>::realT imageXCorrFFT<realImageT>::padFactor()
577 : {
578 : return m_padFactor;
579 : }
580 :
581 : template <class realImageT>
582 : int imageXCorrFFT<realImageT>::rowsPadded()
583 : {
584 : return m_rowsPadded;
585 : }
586 :
587 : template <class realImageT>
588 : int imageXCorrFFT<realImageT>::colsPadded()
589 : {
590 : return m_colsPadded;
591 : }
592 :
593 : template <class realImageT>
594 12 : int imageXCorrFFT<realImageT>::resize( int nrows, int ncols, realT padFactor )
595 : {
596 : // No-op if no change
597 12 : if( m_rows == nrows && m_cols == ncols && m_padFactor == padFactor )
598 : {
599 0 : return 0;
600 : }
601 :
602 12 : if( padFactor < 1 )
603 : {
604 0 : internal::mxlib_error_report( error_t::invalidarg, "padding factor can't be less than 1" );
605 0 : return -1;
606 : }
607 :
608 12 : m_rows = nrows;
609 :
610 12 : m_cols = ncols;
611 :
612 12 : m_padFactor = padFactor;
613 :
614 : /* Any change here invalidates the reference.
615 : * This is trivially true for rows or cols. Also must be
616 : * handled for padding which will change the reference shifts.
617 : */
618 12 : m_refValid = false;
619 :
620 12 : if( m_rows == 0 || m_cols == 0 )
621 : {
622 0 : m_ftIm0.resize( 0, 0 );
623 0 : m_ftWork.resize( 0, 0 );
624 : #ifdef XCFFT_C2C
625 : m_ftWorkIn.resize( 0, 0 );
626 : #endif
627 0 : m_rowsPadded = 0;
628 0 : m_colsPadded = 0;
629 0 : m_ccIm.resize( 0, 0 );
630 0 : m_ftWorkPadded.resize( 0, 0 );
631 0 : m_mtWorkPadded.resize( 0, 0 );
632 0 : m_ftWorkPaddedOut.resize( 0, 0 );
633 :
634 0 : return 0;
635 : }
636 :
637 : #ifdef XCFFT_C2C
638 :
639 : m_ftIm0.resize( m_rows, m_cols );
640 :
641 : m_ftWork.resize( m_rows, m_cols );
642 :
643 : m_ftWorkIn.resize( m_rows, m_cols );
644 :
645 : #else
646 :
647 12 : m_ftIm0.resize( (int)( 0.5 * m_rows ) + 1, m_cols );
648 :
649 12 : m_ftWork.resize( (int)( 0.5 * m_rows ) + 1, m_cols );
650 :
651 : #endif
652 :
653 : // fftw is row-major, eigen defaults to column-major
654 12 : m_fft_fwd.plan( m_cols, m_rows, math::ft::dir::forward, false );
655 :
656 12 : m_rowsPadded = ceil( m_rows * m_padFactor );
657 12 : m_colsPadded = ceil( m_cols * m_padFactor );
658 :
659 12 : m_padFactorR = static_cast<realT>( m_rowsPadded ) / static_cast<realT>( m_rows );
660 12 : m_padFactorC = static_cast<realT>( m_colsPadded ) / static_cast<realT>( m_cols );
661 :
662 12 : m_ccIm.resize( m_rowsPadded, m_colsPadded );
663 :
664 : #ifdef XCFFT_C2C
665 :
666 : m_ftWorkPadded.resize( m_rowsPadded, m_colsPadded );
667 : m_ftWorkPaddedOut.resize( m_rowsPadded, m_colsPadded );
668 :
669 : #else
670 :
671 12 : m_ftWorkPadded.resize( (int)( 0.5 * m_rowsPadded ) + 1, m_colsPadded );
672 :
673 12 : if( m_peakMethod == xcorrPeakMethod::mftOversamp )
674 : {
675 0 : m_ftWorkPaddedOut.resize( m_rowsPadded, m_colsPadded );
676 : }
677 : else
678 : {
679 12 : m_ftWorkPaddedOut.resize( 0, 0 );
680 : }
681 :
682 : #endif
683 :
684 12 : if( m_peakMethod == xcorrPeakMethod::mftOversamp )
685 : {
686 0 : m_mtWorkPadded.resize( m_rowsPadded, m_colsPadded );
687 : }
688 : else
689 : {
690 12 : m_mtWorkPadded.resize( 0, 0 );
691 : }
692 :
693 12 : m_fft_bwd.plan( m_colsPadded, m_rowsPadded, math::ft::dir::backward, false );
694 :
695 12 : return 0;
696 : } // resize(int, int, realT)
697 :
698 : template <class realImageT>
699 : int imageXCorrFFT<realImageT>::resize( int nrows, int ncols )
700 : {
701 : return resize( nrows, ncols, m_padFactor );
702 : }
703 :
704 : template <class realImageT>
705 : int imageXCorrFFT<realImageT>::refMaskIm( const realImageT &mask )
706 : {
707 : m_refMaskIm = mask;
708 : m_haveRefMask = true;
709 :
710 : m_refValid = false;
711 :
712 : return 0;
713 : }
714 :
715 : template <class realImageT>
716 : const realImageT &imageXCorrFFT<realImageT>::refMaskIm()
717 : {
718 : return m_refMaskIm;
719 : }
720 :
721 : template <class realImageT>
722 : bool imageXCorrFFT<realImageT>::haveRefMask()
723 : {
724 : return m_haveRefMask;
725 : }
726 :
727 : template <class realImageT>
728 : int imageXCorrFFT<realImageT>::maskIm( const realImageT &mask )
729 : {
730 : m_maskIm = mask;
731 : m_haveMask = true;
732 :
733 : return 0;
734 : }
735 :
736 : template <class realImageT>
737 : const realImageT &imageXCorrFFT<realImageT>::maskIm()
738 : {
739 : return m_maskIm;
740 : }
741 :
742 : template <class realImageT>
743 : bool imageXCorrFFT<realImageT>::haveMask()
744 : {
745 : return m_haveMask;
746 : }
747 :
748 : template <class realImageT>
749 : int imageXCorrFFT<realImageT>::refWinIm( const realImageT &win )
750 : {
751 : m_refWinIm = win;
752 : m_haveRefWindow = true;
753 :
754 : m_refValid = false;
755 :
756 : return 0;
757 : }
758 :
759 : template <class realImageT>
760 : const realImageT &imageXCorrFFT<realImageT>::refWinIm()
761 : {
762 : return m_refWinIm;
763 : }
764 :
765 : template <class realImageT>
766 : bool imageXCorrFFT<realImageT>::haveRefWindow()
767 : {
768 : return m_haveRefWindow;
769 : }
770 :
771 : template <class realImageT>
772 : int imageXCorrFFT<realImageT>::winIm( const realImageT &win )
773 : {
774 : m_winIm = win;
775 : m_haveWindow = true;
776 :
777 : return 0;
778 : }
779 :
780 : template <class realImageT>
781 : const realImageT &imageXCorrFFT<realImageT>::winIm()
782 : {
783 : return m_winIm;
784 : }
785 :
786 : template <class realImageT>
787 : bool imageXCorrFFT<realImageT>::haveWindow()
788 : {
789 : return m_haveWindow;
790 : }
791 :
792 : template <class realImageT>
793 : void imageXCorrFFT<realImageT>::normalize( bool no )
794 : {
795 : if( no != m_normalize )
796 : {
797 : m_refValid = false;
798 : m_normalize = no;
799 : }
800 : }
801 :
802 : template <class realImageT>
803 : bool imageXCorrFFT<realImageT>::normalize()
804 : {
805 : return m_normalize;
806 : }
807 :
808 : template <class realImageT>
809 12 : int imageXCorrFFT<realImageT>::refIm( const realImageT &im, realT padFactor )
810 : {
811 12 : realImageT im0;
812 :
813 : // Mask if needed
814 12 : if( m_haveRefMask )
815 : {
816 0 : if( im.rows() != m_refMaskIm.rows() && im.cols() != m_refMaskIm.cols() )
817 : {
818 0 : internal::mxlib_error_report( error_t::sizeerr, "reference and reference mask are not the same size" );
819 0 : return -1;
820 : }
821 :
822 : // Now normalize
823 0 : if( m_normalize )
824 : {
825 0 : realT m = imageMean<realT>( im, m_refMaskIm );
826 0 : realT v = imageVariance<realT>( im, m, m_refMaskIm );
827 0 : im0 = ( ( im - m ) * m_refMaskIm ) / sqrt( v );
828 : }
829 : else
830 : {
831 0 : im0 = im * m_refMaskIm;
832 : }
833 : }
834 : else
835 : {
836 : // Normalize
837 12 : if( m_normalize )
838 : {
839 12 : realT m = imageMean<realT>( im );
840 12 : realT v = imageVariance<realT>( im, m );
841 12 : im0 = ( im - m ) / sqrt( v );
842 : }
843 : else
844 : {
845 0 : im0 = im;
846 : }
847 : }
848 :
849 : // Window if needed
850 12 : if( m_haveRefWindow )
851 : {
852 0 : if( im0.rows() != m_refWinIm.rows() && im.cols() != m_refWinIm.cols() )
853 : {
854 0 : internal::mxlib_error_report( error_t::sizeerr, "reference and reference window are not the same size" );
855 0 : return -1;
856 : }
857 :
858 0 : im0 *= m_refWinIm;
859 : }
860 :
861 : // We save refIm as the un-shifted version
862 12 : m_refIm = im0;
863 :
864 : /* Setup the FFTW space
865 : * Note: we don't do this earlier in case we ever implement something (e.g. a filter)
866 : * that changes size before this point
867 : */
868 12 : resize( im0.rows(), im0.cols(), padFactor );
869 :
870 : // Now shift so center pixel is 0,0
871 12 : imageShiftWP( im0, m_refIm, 0.5 * m_rows + 1, 0.5 * m_cols + 1 );
872 :
873 : // Then FT
874 : #ifdef XCFFT_C2C
875 :
876 : m_ftWork = im0.template cast<complexT>();
877 :
878 : m_fft_fwd( m_ftIm0.data(), m_ftWork.data() );
879 : #else
880 12 : m_fft_fwd( m_ftIm0.data(), im0.data() );
881 : #endif
882 :
883 : // Conjugate and normalize for FFTW scaling.
884 786 : for( int c = 0; c < m_ftIm0.cols(); ++c )
885 : {
886 26316 : for( int r = 0; r < m_ftIm0.rows(); ++r )
887 : {
888 25542 : complexT val = std::conj( m_ftIm0( r, c ) );
889 25542 : val /= ( m_rows * m_cols );
890 25542 : m_ftIm0( r, c ) = val;
891 : }
892 : }
893 :
894 : // And finally find the reference shift
895 12 : m_refX0 = 0;
896 12 : m_refY0 = 0;
897 12 : m_refPeak = 0;
898 12 : m_refValid = true;
899 :
900 12 : operator()( m_refX0, m_refY0, m_refPeak, m_refIm );
901 :
902 12 : m_refACIm = m_ccIm;
903 :
904 12 : return 0;
905 12 : }
906 :
907 : template <class realImageT>
908 12 : int imageXCorrFFT<realImageT>::refIm( const realImageT &im )
909 : {
910 12 : return refIm( im, m_padFactor );
911 : }
912 :
913 : template <class realImageT>
914 : const realImageT &imageXCorrFFT<realImageT>::refIm()
915 : {
916 : return m_refIm;
917 : }
918 :
919 : template <class realImageT>
920 : bool imageXCorrFFT<realImageT>::refValid()
921 : {
922 : return m_refValid;
923 : }
924 :
925 : template <class realImageT>
926 23 : void imageXCorrFFT<realImageT>::maxLag( int ml )
927 : {
928 23 : m_maxLag = ml;
929 23 : tol( m_tol );
930 23 : }
931 :
932 : template <class realImageT>
933 : int imageXCorrFFT<realImageT>::maxLag()
934 : {
935 : return m_maxLag;
936 : }
937 :
938 : template <class realImageT>
939 : typename imageXCorrFFT<realImageT>::realT imageXCorrFFT<realImageT>::tol()
940 : {
941 : return m_tol;
942 : }
943 :
944 : template <class realImageT>
945 24 : void imageXCorrFFT<realImageT>::tol( realT nt )
946 : {
947 24 : m_tol = nt;
948 :
949 24 : m_magSize = ceil( ( ( 2. * m_maxLag + 1 ) - 1.0 ) / nt ) + 1;
950 :
951 24 : m_mag = ( m_magSize - 1.0 ) / ( ( 2. * m_maxLag + 1 ) - 1.0 );
952 :
953 24 : m_refValid = false;
954 :
955 24 : return;
956 : }
957 :
958 : template <class realImageT>
959 12 : void imageXCorrFFT<realImageT>::peakMethod( xcorrPeakMethod xpm )
960 : {
961 12 : if( xpm != m_peakMethod )
962 : {
963 8 : m_peakMethod = xpm;
964 :
965 8 : m_refValid = false;
966 : }
967 :
968 12 : return;
969 : }
970 :
971 : template <class realImageT>
972 : xcorrPeakMethod imageXCorrFFT<realImageT>::peakMethod()
973 : {
974 : return m_peakMethod;
975 : }
976 :
977 : template <class realImageT>
978 : const realImageT &imageXCorrFFT<realImageT>::refACIm()
979 : {
980 : return m_refACIm;
981 : }
982 :
983 : template <class realImageT>
984 : typename imageXCorrFFT<realImageT>::realT imageXCorrFFT<realImageT>::refX0()
985 : {
986 : return m_refX0;
987 : }
988 :
989 : template <class realImageT>
990 : typename imageXCorrFFT<realImageT>::realT imageXCorrFFT<realImageT>::refY0()
991 : {
992 : return m_refY0;
993 : }
994 :
995 : template <class realImageT>
996 : const realImageT &imageXCorrFFT<realImageT>::normIm()
997 : {
998 : return m_normIm;
999 : }
1000 :
1001 : template <class realImageT>
1002 : const realImageT &imageXCorrFFT<realImageT>::ccIm()
1003 : {
1004 : return m_ccIm;
1005 : }
1006 :
1007 : template <class realImageT>
1008 : const realImageT &imageXCorrFFT<realImageT>::magIm()
1009 : {
1010 : return m_magIm;
1011 : }
1012 :
1013 : template <class realImageT>
1014 24 : void imageXCorrFFT<realImageT>::findPeak( realT &xShift, realT &yShift, realT &xcPeak )
1015 : {
1016 : int xLag0, yLag0;
1017 :
1018 : /*realImageT sm = m_ccIm;
1019 : filterImage(sm, m_ccIm, gaussKernel<eigenImage<realT>, 2>(3));
1020 : m_ccIm = sm;*/
1021 :
1022 : // m_ccIm = m_ccIm.pow(3);
1023 24 : realT pk = m_ccIm.maxCoeff( &xLag0, &yLag0 );
1024 24 : realT mn = m_ccIm.minCoeff();
1025 :
1026 24 : if( xLag0 - m_maxLag < 0 )
1027 : {
1028 2 : m_maxLag = xLag0;
1029 : }
1030 24 : if( xLag0 + 2 * m_maxLag + 1 >= m_ccIm.rows() )
1031 : {
1032 14 : m_maxLag = ( m_ccIm.rows() - 1 - xLag0 ) / 2;
1033 : }
1034 24 : if( yLag0 - m_maxLag < 0 )
1035 : {
1036 0 : m_maxLag = yLag0;
1037 : }
1038 24 : if( yLag0 + 2 * m_maxLag + 1 >= m_ccIm.cols() )
1039 : {
1040 14 : m_maxLag = ( m_ccIm.cols() - 1 - yLag0 ) / 2;
1041 : }
1042 :
1043 24 : realT x0 = xLag0 - m_maxLag;
1044 24 : realT y0 = yLag0 - m_maxLag;
1045 :
1046 24 : if( m_peakMethod == xcorrPeakMethod::gaussFit )
1047 : {
1048 8 : m_magIm = m_ccIm.block( x0, y0, 2 * m_maxLag + 1, 2 * m_maxLag + 1 );
1049 :
1050 8 : m_fitter.setArray( m_magIm.data(), m_magIm.rows(), m_magIm.cols() );
1051 8 : m_magIm.maxCoeff( &xLag0, &yLag0 );
1052 :
1053 8 : m_fitter.setGuess( mn, pk - mn, m_maxLag, m_maxLag, 3, 3, 0 );
1054 8 : m_fitter.fit();
1055 :
1056 8 : x0 = ( m_fitter.x0() + x0 );
1057 8 : y0 = ( m_fitter.y0() + y0 );
1058 8 : xcPeak = m_fitter.G();
1059 : }
1060 16 : else if( m_peakMethod == xcorrPeakMethod::interpPeak )
1061 : {
1062 : int x, y;
1063 8 : m_magIm.resize( m_magSize, m_magSize );
1064 :
1065 8 : imageMagnify( m_magIm,
1066 8 : m_ccIm.block( x0, y0, 2 * m_maxLag + 1, 2 * m_maxLag + 1 ),
1067 8 : cubicConvolTransform<realT>() );
1068 :
1069 8 : xcPeak = m_magIm.maxCoeff( &x, &y );
1070 8 : x0 = ( x * ( 1.0 / m_mag ) + x0 );
1071 8 : y0 = ( y * ( 1.0 / m_mag ) + y0 );
1072 : }
1073 8 : else if( m_peakMethod == xcorrPeakMethod::centerOfLight )
1074 : {
1075 : realT x, y;
1076 :
1077 8 : m_magIm = m_ccIm.block( x0, y0, 2 * m_maxLag + 1, 2 * m_maxLag + 1 );
1078 8 : m_magIm -= m_magIm.minCoeff(); // Must sum to > 0.
1079 8 : xcPeak = m_magIm.maxCoeff();
1080 8 : imageCenterOfLight( x, y, m_magIm );
1081 :
1082 8 : x0 = ( x + x0 );
1083 8 : y0 = ( y + y0 );
1084 : }
1085 0 : else if( m_peakMethod == xcorrPeakMethod::mftOversamp )
1086 : {
1087 : int x, y;
1088 0 : m_ccIm.maxCoeff( &x, &y );
1089 :
1090 0 : if( m_refX0 == 0 && m_refY0 == 0 ) // reference finding
1091 : {
1092 0 : x0 = x;
1093 0 : y0 = y;
1094 : }
1095 : else
1096 : {
1097 0 : x = x / m_padFactorR - m_refX0;
1098 0 : y = y / m_padFactorC - m_refY0;
1099 :
1100 0 : realT mag = 1.0 / m_tol; // m_mag is the effective mag for interpPeak...
1101 :
1102 0 : realT tmagR = mag * m_padFactorR; // the total mag
1103 0 : realT tmagC = mag * m_padFactorC; // the total mag
1104 :
1105 0 : realT shiftx0 = -1 * ( m_refX0 * m_padFactorR ) * ( mag - 1 ) - x * tmagR;
1106 0 : realT shifty0 = -1 * ( m_refY0 * m_padFactorC ) * ( mag - 1 ) - y * tmagC;
1107 :
1108 0 : m_mft_bwd.plan( m_rowsPadded, m_colsPadded, math::ft::dir::backward, shiftx0, shifty0, mag );
1109 :
1110 : #ifdef XCFFT_C2C
1111 :
1112 : m_mft_bwd( m_ftWorkPaddedOut, m_ftWorkPadded );
1113 :
1114 : #else
1115 :
1116 0 : m_mtWorkPadded.resize( m_rowsPadded, m_colsPadded );
1117 0 : m_mtWorkPadded.setZero();
1118 0 : math::ft::augmentR2CFFTOutput( m_mtWorkPadded, m_ftWork );
1119 :
1120 0 : m_mft_bwd( m_ftWorkPaddedOut, m_mtWorkPadded );
1121 :
1122 : #endif
1123 :
1124 0 : m_magIm = m_ftWorkPaddedOut.real();
1125 :
1126 : int xf, yf;
1127 0 : xcPeak = m_magIm.maxCoeff( &xf, &yf );
1128 :
1129 0 : xShift = x + ( xf - m_refX0 * m_padFactorR ) / tmagR;
1130 0 : yShift = y + ( yf - m_refY0 * m_padFactorC ) / tmagC;
1131 :
1132 0 : return;
1133 : }
1134 : }
1135 0 : else if( m_peakMethod == xcorrPeakMethod::none )
1136 : {
1137 : int x, y;
1138 0 : xcPeak = m_ccIm.maxCoeff( &x, &y );
1139 0 : x0 = x;
1140 0 : y0 = y;
1141 : }
1142 : else
1143 : {
1144 0 : throw( mx::exception( error_t::invalidconfig, "unknown peak finding method" ) );
1145 : }
1146 :
1147 : //--> unpad here, scaling the shifts
1148 24 : xShift = x0 / m_padFactorR - m_refX0;
1149 24 : yShift = y0 / m_padFactorC - m_refY0;
1150 : }
1151 :
1152 : template <class realImageT>
1153 : template <class imT>
1154 24 : int imageXCorrFFT<realImageT>::operator()( realT &xShift, realT &yShift, realT &xcPeak, const imT &im )
1155 : {
1156 24 : if( !m_refValid )
1157 : {
1158 0 : throw( mx::exception( error_t::invalidconfig, "reference image is not valid" ) );
1159 : }
1160 :
1161 24 : if( im.rows() != m_rows )
1162 : {
1163 0 : throw( mx::exception( error_t::sizeerr, "image must be same size as reference (rows)" ) );
1164 : }
1165 :
1166 24 : if( im.cols() != m_cols )
1167 : {
1168 0 : throw( mx::exception( error_t::sizeerr, "image must be same size as reference (rows)" ) );
1169 : }
1170 :
1171 : // Mask and normalize as needed
1172 24 : if( m_haveMask )
1173 : {
1174 0 : if( m_normalize )
1175 : {
1176 0 : realT m = imageMean<realT>( im, m_maskIm );
1177 0 : realT v = imageVariance<realT>( im, m, m_maskIm );
1178 :
1179 0 : m_normIm = ( im - m ) * m_maskIm / sqrt( v );
1180 : }
1181 : else
1182 : {
1183 0 : m_normIm = im * m_maskIm;
1184 : }
1185 : }
1186 : else
1187 : {
1188 :
1189 24 : if( m_normalize )
1190 : {
1191 24 : realT m = imageMean<realT>( im );
1192 24 : realT v = imageVariance<realT>( im, m );
1193 :
1194 24 : m_normIm = ( im - m ) / sqrt( v );
1195 : }
1196 : else
1197 : {
1198 0 : m_normIm = im;
1199 : }
1200 : }
1201 :
1202 24 : if( m_haveWindow )
1203 : {
1204 0 : m_normIm *= m_winIm;
1205 : }
1206 :
1207 : #ifdef XCFFT_C2C
1208 : m_ftWorkIn = m_normIm.template cast<complexT>();
1209 :
1210 : m_fft_fwd( m_ftWork.data(), m_ftWorkIn.data() );
1211 : #else
1212 24 : m_fft_fwd( m_ftWork.data(), m_normIm.data() );
1213 : #endif
1214 :
1215 : // m_ftIm0 is the conjugated, fftw-normalized, reference image in the Fourier domain
1216 : // So this is the FT of the cross-correlation:
1217 24 : m_ftWork *= m_ftIm0;
1218 :
1219 24 : m_ftWorkPadded.setZero();
1220 :
1221 : #ifdef XCFFT_C2C
1222 :
1223 : math::ft::padC2CFFTOutput( m_ftWorkPadded, m_ftWork );
1224 :
1225 : m_fft_bwd( m_ftWorkPaddedOut.data(), m_ftWorkPadded.data() );
1226 :
1227 : m_ccIm = m_ftWorkPaddedOut.real();
1228 :
1229 : #else
1230 :
1231 24 : math::ft::padR2CFFTOutput( m_ftWorkPadded, m_ftWork );
1232 :
1233 24 : m_fft_bwd( m_ccIm, m_ftWorkPadded );
1234 :
1235 : #endif
1236 :
1237 24 : findPeak( xShift, yShift, xcPeak );
1238 :
1239 24 : return 0;
1240 : }
1241 :
1242 : template <class realImageT>
1243 : template <class im0T, class imT>
1244 : int imageXCorrFFT<realImageT>::operator()( realT &xShift, realT &yShift, realT &xcPeak, im0T &im0, imT &im )
1245 : {
1246 : setReference( im0 );
1247 : return operator()( xShift, yShift, xcPeak, im );
1248 : }
1249 :
1250 : } // namespace improc
1251 : } // namespace mx
1252 :
1253 : #endif // imageXCorrFFT_hpp
|
