Matlab normxcorr2 implemented in python

For some reason there is no direct implementation of normxcorr2 from Matlab or Octave in numpy or scipy. This is a near exact normxcorr2 taken from Octave's implementation using numpy and scipy's fftconvolve.

# Based on Octave implementation by Benjamin Eltzner 

import numpy as np
from scipy.signal import fftconvolve


def normxcorr2(template, image, mode="full"):
    # 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)
    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)

    out[np.where(np.logical_not(np.isfinite(out)))] = 0
    return out

Using fftconvolve as it is significantly faster than convolve2d or correlate2d. Also allows for N-dimensions.
Python file available on GitHub: normxcorr2-python