Phase-Only Correlation in Python

In a meeting at work the other day, the topic of image misalignment in the X and Y directions came up, and it reminded me that I had written a Phase-Only Correlation (POC) program back in graduate school.

In image processing, it is common to need an algorithm that detects the positional shift between two images, and many different approaches have been developed for that purpose.

At the time, I was working on super-resolution, and I used this algorithm because of its accuracy. I found the old code, so I am posting it here. Unfortunately, programs like this are not as easy to find online as you might expect.

This version does not estimate rotation. If you also need rotation, you can apply a log-polar transform before alignment. I had code for that too, but I was not sure whether it still worked, so I left it out this time.

#! /usr/bin/env python
# -*- coding: utf-8
import numpy
import scipy, scipy.fftpack
from numpy import pi, sin, cos
from scipy.optimize import leastsq

def zero_padding(src, dstshape, pos = (0, 0)):
    y, x = pos
    dst = numpy.zeros(dstshape)
    dst[y:src.shape[0] + y, x:src.shape[1] + x] = src
    return dst

def pocfunc_model(alpha, delta1, delta2, r, u):
    N1, N2 = r.shape
    V1, V2 = map(lambda x: 2 * x + 1, u)
    return lambda n1, n2: alpha / (N1 * N2) * sin((n1 + delta1) * V1 / N1 * pi) * sin((n2 + delta2) * V2 / N2 * pi)\
                                            / (sin((n1 + delta1) * pi / N1) * sin((n2 + delta2) * pi / N2))


def pocfunc(f, g, windowfunc = numpy.hanning, withlpf = True):
    m = numpy.floor(map(lambda x: x / 2.0, f.shape))
    u = map(lambda x: x / 2.0, m)

    # hanning window
    hy = windowfunc(f.shape[0])
    hx = windowfunc(f.shape[1])
    hw = hy.reshape(hy.shape[0], 1) * hx
    f = f * hw
    g = g * hw

    # compute 2d fft
    F = scipy.fftpack.fft2(f)
    G = scipy.fftpack.fft2(g)
    G_ = numpy.conj(G)
    R = F * G_ / numpy.abs(F * G_)

    if withlpf == True:
        R = scipy.fftpack.fftshift(R)
        lpf = numpy.ones(map(lambda x: x + 1.0, m))
        lpf = zero_padding(lpf, f.shape, u)
        R = R * lpf
        R = scipy.fftpack.fftshift(R)

    return scipy.fftpack.fftshift(numpy.real(scipy.fftpack.ifft2(R)))


def poc(f, g, fitting_shape = (9, 9)):
    # compute phase-only correlation
    center = map(lambda x: x / 2.0, f.shape)
    m = numpy.floor(map(lambda x: x / 2.0, f.shape))
    u = map(lambda x: x / 2.0, m)

    r = pocfunc(f, g)

    # least-square fitting
    max_pos = numpy.argmax(r)
    peak = (max_pos / f.shape[1], max_pos % f.shape[1])
    max_peak = r[peak[0], peak[1]]

    mf = numpy.floor(map(lambda x: x / 2.0, fitting_shape))
    fitting_area = r[peak[0] - mf[0] : peak[0] + mf[0] + 1,\
                     peak[1] - mf[1] : peak[1] + mf[1] + 1]

    p0 = [0.5, -(peak[0] - m[0]) - 0.02, -(peak[1] - m[1]) - 0.02]
    y, x = numpy.mgrid[-mf[0]:mf[0] + 1, -mf[1]:mf[1] + 1]
    y = y + peak[0] - m[0]
    x = x + peak[1] - m[1]
    errorfunction = lambda p: numpy.ravel(pocfunc_model(p[0], p[1], p[2], r, u)(y, x) - fitting_area)
    plsq = leastsq(errorfunction, p0)
    return (plsq[0][0], plsq[0][1], plsq[0][2])

Usage

I recommend using OpenCV. I tested it by taking the image on the left below and shifting it in GIMP by 59 pixels in X and 16 pixels in Y. The result looks correct. The first value in the output is the correlation score between the two images, the second is the offset in Y, and the third is the offset in X.

import cv2
from poc import *
img1 = cv2.imread("lena.jpg", 0)
img2 = cv2.imread("lena_shift_y16x59.jpg", 0)
poc(img1, img2)
# result (0.97930538205507089, 15.997296364421981, 58.996715087791408)

I also uploaded the code to GitHub: daisukekobayashi/pocpy