% MATCHBYMONOGENICPHASE - match image feature points using monogenic phase data
%
% Function generates putative matches between previously detected
% feature points in two images by looking for points that have minimal
% differences in monogenic phase data within windows surrounding each point.
% Only points that correlate most strongly with each other in *both*
% directions are returned. This is a simple-minded N^2 comparison.
%
% This matcher performs rather well relative to normalised greyscale
% correlation. Typically there are more putative matches found and fewer
% outliers. There is a greater computational cost in the pre-filtering stage
% but potentially the matching stage is much faster as each pixel is effectively
% encoded with only 3 bits. (Though this potential speed is not realized in this
% implementation)
%
% Usage: [m1,m2] = matchbymonogenicphase(im1, p1, im2, p2, w, dmax, ...
% nscale, minWaveLength, mult, sigmaOnf)
%
% Arguments:
% im1, im2 - Images containing points that we wish to match.
% p1, p2 - Coordinates of feature pointed detected in im1 and
% im2 respectively using a corner detector (say Harris
% or phasecong2). p1 and p2 are [2xnpts] arrays though
% p1 and p2 are not expected to have the same number
% of points. The first row of p1 and p2 gives the row
% coordinate of each feature point, the second row
% gives the column of each point.
% w - Window size (in pixels) over which the phase bit codes
% around each feature point are matched. This should
% be an odd number.
% dmax - Maximum search radius for matching points. Used to
% improve speed when there is little disparity between
% images. Even setting it to a generous value of 1/4 of
% the image size gives a useful speedup.
% nscale - Number of filter scales.
% minWaveLength - Wavelength of smallest scale filter.
% mult - Scaling factor between successive filters.
% sigmaOnf - Ratio of the standard deviation of the Gaussian
% describing the log Gabor filter's transfer function in
% the frequency domain to the filter center frequency.
%
%
% Returns:
% m1, m2 - Coordinates of points selected from p1 and p2
% respectively such that (putatively) m1(:,i) matches
% m2(:,i). m1 and m2 are [2xnpts] arrays defining the
% points in each of the images in the form [row;col].
%
%
% I have had good success with the folowing parameters:
%
% w = 11; Window size for correlation matching, 7 or greater
% seems fine.
% dmax = 50;
% nscale = 1; Just one scale can give very good results. Adding
% another scale doubles computation
% minWaveLength = 10;
% mult = 4; This is irrelevant if only one scale is used. If you do
% use more than one scale try values in the range 2-4.
% sigmaOnf = .2; This results in a *very* large bandwidth filter. A
% large bandwidth seems to be very important in the
% matching performance.
%
% See Also: MATCHBYCORRELATION, MONOFILT
% Copyright (c) 2005 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.
% May 2005 - Original version adapted from matchbycorrelation.m
function [m1,m2,cormat] = matchbymonogenicphase(im1, p1, im2, p2, w, dmax, ...
nscale, minWaveLength, mult, sigmaOnf)
orientWrap = 0;
[f1, h1f1, h2f1, A1] = ...
monofilt(im1, nscale, minWaveLength, mult, sigmaOnf, orientWrap);
[f2, h1f2, h2f2, A2] = ...
monofilt(im2, nscale, minWaveLength, mult, sigmaOnf, orientWrap);
% Normalise filter outputs to unit vectors (should also have masking for
% unreliable filter outputs)
for s = 1:nscale
% f1{s} = f1{s}./A1{s}; f2{s} = f2{s}./A2{s};
% h1f1{s} = h1f1{s}./A1{s}; h1f2{s} = h1f2{s}./A2{s};
% h2f1{s} = h2f1{s}./A1{s}; h2f2{s} = h2f2{s}./A2{s};
% Try quantizing oriented phase vector to 8 octants to see what
% effect this has (Performance seems to be reduced only slightly)
f1{s} = sign(f1{s}); f2{s} = sign(f2{s});
h1f1{s} = sign(h1f1{s}); h1f2{s} = sign(h1f2{s});
h2f1{s} = sign(h2f1{s}); h2f2{s} = sign(h2f2{s});
end
% Generate correlation matrix
cormat = correlationmatrix(f1, h1f1, h2f1, p1, ...
f2, h1f2, h2f2, p2, w, dmax);
[corrows,corcols] = size(cormat);
% Find max along rows give strongest match in p2 for each p1
[mp2forp1, colp2forp1] = max(cormat,[],2);
% Find max down cols give strongest match in p1 for each p2
[mp1forp2, rowp1forp2] = max(cormat,[],1);
% Now find matches that were consistent in both directions
p1ind = zeros(1,length(p1)); % Arrays for storing matched indices
p2ind = zeros(1,length(p2));
indcount = 0;
for n = 1:corrows
if rowp1forp2(colp2forp1(n)) == n % consistent both ways
indcount = indcount + 1;
p1ind(indcount) = n;
p2ind(indcount) = colp2forp1(n);
end
end
% Trim arrays of indices of matched points
p1ind = p1ind(1:indcount);
p2ind = p2ind(1:indcount);
% Extract matched points from original arrays
m1 = p1(:,p1ind);
m2 = p2(:,p2ind);
%-------------------------------------------------------------------------
% Function that does the work. This function builds a 'correlation' matrix
% that holds the correlation strength of every point relative to every other
% point. While this seems a bit wasteful we need all this data if we want
% to find pairs of points that correlate maximally in both directions.
function cormat = correlationmatrix(f1, h1f1, h2f1, p1, ...
f2, h1f2, h2f2, p2, w, dmax)
if mod(w, 2) == 0 | w < 3
error('Window size should be odd and >= 3');
end
r = (w-1)/2; % 'radius' of correlation window
[rows1, npts1] = size(p1);
[rows2, npts2] = size(p2);
if rows1 ~= 2 | rows2 ~= 2
error('Feature points must be specified in 2xN arrays');
end
% Reorganize monogenic phase data into a 4D matrices for convenience
[im1rows,im1cols] = size(f1{1});
[im2rows,im2cols] = size(f2{1});
nscale = length(f1);
phase1 = zeros(im1rows,im1cols,nscale,3);
phase2 = zeros(im2rows,im2cols,nscale,3);
for s = 1:nscale
phase1(:,:,s,1) = f1{s}; phase1(:,:,s,2) = h1f1{s}; phase1(:,:,s,3) = h2f1{s};
phase2(:,:,s,1) = f2{s}; phase2(:,:,s,2) = h1f2{s}; phase2(:,:,s,3) = h2f2{s};
end
% Initialize correlation matrix values to -infinity
cormat = repmat(-inf, npts1, npts2);
% For every feature point in the first image extract a window of data
% and correlate with a window corresponding to every feature point in
% the other image. Any feature point less than distance 'r' from the
% boundary of an image is not considered.
% Find indices of points that are distance 'r' or greater from
% boundary on image1 and image2;
n1ind = find(p1(1,:)>r & p1(1,:)r & p1(2,:)r & p2(1,:)r & p2(2,:)