Phase Correlation

This is an implementation of the phase correlation for aligning images.  This implementation also works on images which are of different size.  I.e. stitching of overlapping images is possible as well as searching a small template in a large image.  Note that you may have to add a windowing function to the implementation to avoid boundary effects.

See also

• Normalised Cross-Correlation
• http://en.wikipedia.org/wiki/Phase_correlation

#!/usr/bin/env ruby
# Image-registration using phase correlation.
# Thanks to Denis Tessier for his initial contribution.
# http://en.wikipedia.org/wiki/Phase_correlation
require 'hornetseye'
include Hornetseye
class MultiArray
  def zeropad( *s )
    retval = MultiArray.new( typecode, *s ).fill!
    retval[ 0...shape[0], 0...shape[1] ] = self
    retval
  end
  def phasecorr( other )
    nominator = rfft * other.rfft.conj
    denominator = nominator.abs
    mask = denominator.abs > 1.0e-5
    ( nominator.mask( mask ) / denominator.mask( mask ) ).unmask( mask ).irfft
  end
  def zphasecorr( other )
    xshape = [ shape[0] + other.shape[0], shape[1] + other.shape[1] ]
    zeropad( *xshape ).phasecorr( other.zeropad( *xshape ) )
  end
end
syntax = <<END_OF_STRING
Align images using phase correlation
Syntax : phasecorrelation.rb <image1> <image2>
Example: phasecorrelation.rb image1.jpg image2.jpg
END_OF_STRING
if ARGV.size != 2
  puts syntax
  raise "Wrong number of command-line arguments"
end
image = MultiArray.load_grey8( ARGV[0] )
template = MultiArray.load_grey8( ARGV[1] )
shift = image.zphasecorr( template )
idx = MultiArray.lint( *shift.shape ).indgen!.mask( shift >= shift.max )[0]
shiftx = idx % shift.shape[0]
shifty = idx / shift.shape[0]
shiftx = shiftx - image.shape[0] - template.shape[0] if shiftx > image.shape[0]
shifty = shifty - image.shape[1] - template.shape[1] if shifty > image.shape[1]
minx = [ 0, shiftx ].min
miny = [ 0, shifty ].min
maxx = [ image.shape[0], template.shape[0] + shiftx ].max - 1
maxy = [ image.shape[1], template.shape[1] + shifty ].max - 1
offsetx = -minx
offsety = -miny
resultwidth  = maxx + 1 - minx
resultheight = maxy + 1 - miny
result1 = MultiArray.ubyte( resultwidth, resultheight ).fill!
result1[ offsetx...( offsetx + image.shape[0] ),
         offsety...( offsety + image.shape[1] ) ] = image / 2
result2 = MultiArray.ubyte( resultwidth, resultheight ).fill!
result2[ ( shiftx + offsetx )...( shiftx + offsetx + template.shape[0] ),
         ( shifty + offsety )...( shifty + offsety + template.shape[1] ) ] = template / 2
( result1 + result2 ).show