Cross-correlation between images

Question

I'd like to compute the cross correlation using de Fast Fourier Transform, for cloud motion tracking following the steps of the image below.

def roi_image(image):
    image = cv.imread(image, 0)
    roi = image[700:900, 1900:2100]
    return roi

def FouTransf(image):
    img_f32 = np.float32(image)
    d_ft = cv.dft(img_f32, flags = cv.DFT_COMPLEX_OUTPUT)
    d_ft_shift = np.fft.fftshift(d_ft)

    rows, cols = image.shape
    opt_rows = cv.getOptimalDFTSize(rows)
    opt_cols = cv.getOptimalDFTSize(cols)
    opt_img = np.zeros((opt_rows, opt_cols))
    opt_img[:rows, :cols] = image 
    crow, ccol = opt_rows / 2 , opt_cols / 2
    mask = np.zeros((opt_rows, opt_cols, 2), np.uint8)
    mask[int(crow-50):int(crow+50), int(ccol-50):int(ccol+50)] = 1

    f_mask = d_ft_shift*mask
    return f_mask


def inv_FouTransf(image):

    f_ishift = np.fft.ifftshift(image)
    img_back = cv.idft(f_ishift)
    img_back = cv.magnitude(img_back[:, :, 0], img_back[:, :, 1])

    return img_back

def rms(sigma):
    rms = np.std(sigma)
    return rms

# Step 1: Import images
a = roi_image(path_a)
b = roi_image(path_b)

# Step 2: Convert the image to frequency domain
G_t0 = FouTransf(a)
G_t0_conj = G_t0.conj()
G_t1 = FouTransf(b)

# Step 3: Compute C(m, v)
C = G_t0_conj * G_t1

# Step 4: Convert the image to space domain to obtain Cov (p, q)
c_w = inv_FouTransf(C)

# Step 5: Compute Cross correlation
R_pq = c_w / (rms(a) * rms(b)) 

I'm a little confused because I've never use that technique. ¿The application es accurate?

HINT: eq (1) is : R(p,q) = Cov(p,q) / (sigma_t0 * sigma_t1). If more information is required the paper is: "An Automated Techinique or Obtaining Cloud Motion from Geostatiory Satellite Data Using Cross Correlation".

I found this source but I don't know if does something I'm trying.

Answer 1

If you are trying to do something similar to cv2.matchTemplate(), a working python implementation of the Normalized Cross-Correlation (NCC) method can be found in this repository:

########################################################################################
# Author: Ujash Joshi, University of Toronto, 2017                                     #
# Based on Octave implementation by: Benjamin Eltzner, 2014 <[email protected]>         #
# Octave/Matlab normxcorr2 implementation in python 3.5                                #
# Details:                                                                             #
# Normalized cross-correlation. Similiar results upto 3 significant digits.            #
# https://github.com/Sabrewarrior/normxcorr2-python/master/norxcorr2.py                #
# http://lordsabre.blogspot.ca/2017/09/matlab-normxcorr2-implemented-in-python.html    #
########################################################################################

import numpy as np
from scipy.signal import fftconvolve


def normxcorr2(template, image, mode="full"):
    """
    Input arrays should be floating point numbers.
    :param template: N-D array, of template or filter you are using for cross-correlation.
    Must be less or equal dimensions to image.
    Length of each dimension must be less than length of image.
    :param image: N-D array
    :param mode: Options, "full", "valid", "same"
    full (Default): The output of fftconvolve is the full discrete linear convolution of the inputs. 
    Output size will be image size + 1/2 template size in each dimension.
    valid: The output consists only of those elements that do not rely on the zero-padding.
    same: The output is the same size as image, centered with respect to the ‘full’ output.
    :return: N-D array of same dimensions as image. Size depends on mode parameter.
    """

    # If this happens, it is probably a mistake
    if np.ndim(template) > np.ndim(image) or \
            len([i for i in range(np.ndim(template)) if template.shape[i] > image.shape[i]]) > 0:
        print("normxcorr2: TEMPLATE larger than IMG. Arguments may be swapped.")

    template = template - np.mean(template)
    image = image - np.mean(image)

    a1 = np.ones(template.shape)
    # Faster to flip up down and left right then use fftconvolve instead of scipy's correlate
    ar = np.flipud(np.fliplr(template))
    out = fftconvolve(image, ar.conj(), mode=mode)

    image = fftconvolve(np.square(image), a1, mode=mode) - \
            np.square(fftconvolve(image, a1, mode=mode)) / (np.prod(template.shape))

    # Remove small machine precision errors after subtraction
    image[np.where(image < 0)] = 0

    template = np.sum(np.square(template))
    out = out / np.sqrt(image * template)

    # Remove any divisions by 0 or very close to 0
    out[np.where(np.logical_not(np.isfinite(out)))] = 0

    return out

The returned object from normxcorr2() is the cross correlation matrix.