[37] | 1 | % MATCHBYCORRELATION - match image feature points by correlation |
---|
| 2 | % |
---|
| 3 | % Function generates putative matches between previously detected |
---|
| 4 | % feature points in two images by looking for points that are maximally |
---|
| 5 | % correlated with each other within windows surrounding each point. |
---|
| 6 | % Only points that correlate most strongly with each other in *both* |
---|
| 7 | % directions are returned. |
---|
| 8 | % This is a simple-minded N^2 comparison. |
---|
| 9 | % |
---|
| 10 | % Usage: [m1,m2] = matchbycorrelation(im1, p1, im2, p2, w, dmax) |
---|
| 11 | % |
---|
| 12 | % Arguments: |
---|
| 13 | % im1, im2 - Images containing points that we wish to match. |
---|
| 14 | % p1, p2 - Coordinates of feature pointed detected in im1 and |
---|
| 15 | % im2 respectively using a corner detector (say Harris |
---|
| 16 | % or phasecong2). p1 and p2 are [2xnpts] arrays though |
---|
| 17 | % p1 and p2 are not expected to have the same number |
---|
| 18 | % of points. The first row of p1 and p2 gives the row |
---|
| 19 | % coordinate of each feature point, the second row |
---|
| 20 | % gives the column of each point. |
---|
| 21 | % w - Window size (in pixels) over which the correlation |
---|
| 22 | % around each feature point is performed. This should |
---|
| 23 | % be an odd number. |
---|
| 24 | % dmax - (Optional) Maximum search radius for matching |
---|
| 25 | % points. Used to improve speed when there is little |
---|
| 26 | % disparity between images. Even setting it to a generous |
---|
| 27 | % value of 1/4 of the image size gives a useful |
---|
| 28 | % speedup. If this parameter is omitted it defaults to Inf. |
---|
| 29 | % |
---|
| 30 | % |
---|
| 31 | % Returns: |
---|
| 32 | % m1, m2 - Coordinates of points selected from p1 and p2 |
---|
| 33 | % respectively such that (putatively) m1(:,i) matches |
---|
| 34 | % m2(:,i). m1 and m2 are [2xnpts] arrays defining the |
---|
| 35 | % points in each of the images in the form [row;col]. |
---|
| 36 | |
---|
| 37 | % Copyright (c) 2004 Peter Kovesi |
---|
| 38 | % School of Computer Science & Software Engineering |
---|
| 39 | % The University of Western Australia |
---|
| 40 | % http://www.csse.uwa.edu.au/ |
---|
| 41 | % |
---|
| 42 | % Permission is hereby granted, free of charge, to any person obtaining a copy |
---|
| 43 | % of this software and associated documentation files (the "Software"), to deal |
---|
| 44 | % in the Software without restriction, subject to the following conditions: |
---|
| 45 | % |
---|
| 46 | % The above copyright notice and this permission notice shall be included in |
---|
| 47 | % all copies or substantial portions of the Software. |
---|
| 48 | % |
---|
| 49 | % The Software is provided "as is", without warranty of any kind. |
---|
| 50 | |
---|
| 51 | % February 2004 - Original version |
---|
| 52 | % May 2004 - Speed improvements + constraint on search radius for |
---|
| 53 | % additional speed |
---|
| 54 | % August 2004 - Vectorized distance calculation for more speed |
---|
| 55 | % (thanks to Daniel Wedge) |
---|
| 56 | |
---|
| 57 | |
---|
| 58 | function [m1,m2,cormat] = matchbycorrelation(im1, p1, im2, p2, w, dmax) |
---|
| 59 | |
---|
| 60 | if nargin == 5 |
---|
| 61 | dmax = Inf; |
---|
| 62 | end |
---|
| 63 | |
---|
| 64 | im1 = double(im1); |
---|
| 65 | im2 = double(im2); |
---|
| 66 | |
---|
| 67 | % Subtract image smoothed with an averaging filter of size wXw from |
---|
| 68 | % each of the images. This compensates for brightness differences in |
---|
| 69 | % each image. Doing it now allows faster correlation calculation. |
---|
| 70 | im1 = im1 - filter2(fspecial('average',w),im1); |
---|
| 71 | im2 = im2 - filter2(fspecial('average',w),im2); |
---|
| 72 | |
---|
| 73 | % Generate correlation matrix |
---|
| 74 | cormat = correlatiomatrix(im1, p1, im2, p2, w, dmax); |
---|
| 75 | [corrows,corcols] = size(cormat); |
---|
| 76 | |
---|
| 77 | % Find max along rows give strongest match in p2 for each p1 |
---|
| 78 | [mp2forp1, colp2forp1] = max(cormat,[],2); |
---|
| 79 | |
---|
| 80 | % Find max down cols give strongest match in p1 for each p2 |
---|
| 81 | [mp1forp2, rowp1forp2] = max(cormat,[],1); |
---|
| 82 | |
---|
| 83 | % Now find matches that were consistent in both directions |
---|
| 84 | p1ind = zeros(1,length(p1)); % Arrays for storing matched indices |
---|
| 85 | p2ind = zeros(1,length(p2)); |
---|
| 86 | indcount = 0; |
---|
| 87 | for n = 1:corrows |
---|
| 88 | if rowp1forp2(colp2forp1(n)) == n % consistent both ways |
---|
| 89 | indcount = indcount + 1; |
---|
| 90 | p1ind(indcount) = n; |
---|
| 91 | p2ind(indcount) = colp2forp1(n); |
---|
| 92 | end |
---|
| 93 | end |
---|
| 94 | |
---|
| 95 | % Trim arrays of indices of matched points |
---|
| 96 | p1ind = p1ind(1:indcount); |
---|
| 97 | p2ind = p2ind(1:indcount); |
---|
| 98 | |
---|
| 99 | % Extract matched points from original arrays |
---|
| 100 | m1 = p1(:,p1ind); |
---|
| 101 | m2 = p2(:,p2ind); |
---|
| 102 | |
---|
| 103 | |
---|
| 104 | %------------------------------------------------------------------------- |
---|
| 105 | % Function that does the work. This function builds a correlation matrix |
---|
| 106 | % that holds the correlation strength of every point relative to every |
---|
| 107 | % other point. While this seems a bit wasteful we need all this data if |
---|
| 108 | % we want to find pairs of points that correlate maximally in both |
---|
| 109 | % directions. |
---|
| 110 | % |
---|
| 111 | % This code assumes im1 and im2 have zero mean. This speeds the |
---|
| 112 | % calculation of the normalised correlation measure. |
---|
| 113 | |
---|
| 114 | function cormat = correlatiomatrix(im1, p1, im2, p2, w, dmax) |
---|
| 115 | |
---|
| 116 | if mod(w, 2) == 0 |
---|
| 117 | error('Window size should be odd'); |
---|
| 118 | end |
---|
| 119 | |
---|
| 120 | [rows1, npts1] = size(p1); |
---|
| 121 | [rows2, npts2] = size(p2); |
---|
| 122 | |
---|
| 123 | % Initialize correlation matrix values to -infinty |
---|
| 124 | cormat = -ones(npts1,npts2)*Inf; |
---|
| 125 | |
---|
| 126 | if rows1 ~= 2 | rows2 ~= 2 |
---|
| 127 | error('Feature points must be specified in 2xN arrays'); |
---|
| 128 | end |
---|
| 129 | |
---|
| 130 | [im1rows, im1cols] = size(im1); |
---|
| 131 | [im2rows, im2cols] = size(im2); |
---|
| 132 | |
---|
| 133 | r = (w-1)/2; % 'radius' of correlation window |
---|
| 134 | |
---|
| 135 | % For every feature point in the first image extract a window of data |
---|
| 136 | % and correlate with a window corresponding to every feature point in |
---|
| 137 | % the other image. Any feature point less than distance 'r' from the |
---|
| 138 | % boundary of an image is not considered. |
---|
| 139 | |
---|
| 140 | % Find indices of points that are distance 'r' or greater from |
---|
| 141 | % boundary on image1 and image2; |
---|
| 142 | n1ind = find(p1(1,:)>r & p1(1,:)<im1rows+1-r & ... |
---|
| 143 | p1(2,:)>r & p1(2,:)<im1cols+1-r); |
---|
| 144 | |
---|
| 145 | n2ind = find(p2(1,:)>r & p2(1,:)<im2rows+1-r & ... |
---|
| 146 | p2(2,:)>r & p2(2,:)<im2cols+1-r); |
---|
| 147 | |
---|
| 148 | for n1 = n1ind |
---|
| 149 | % Generate window in 1st image |
---|
| 150 | w1 = im1(p1(1,n1)-r:p1(1,n1)+r, p1(2,n1)-r:p1(2,n1)+r); |
---|
| 151 | % Pre-normalise w1 to a unit vector. |
---|
| 152 | w1 = w1./sqrt(sum(sum(w1.*w1))); |
---|
| 153 | |
---|
| 154 | % Identify the indices of points in p2 that we need to consider. |
---|
| 155 | if dmax == inf |
---|
| 156 | n2indmod = n2ind; % We have to consider all of n2ind |
---|
| 157 | |
---|
| 158 | else % Compute distances from p1(:,n1) to all available p2. |
---|
| 159 | p1pad = repmat(p1(:,n1),1,length(n2ind)); |
---|
| 160 | dists2 = sum((p1pad-p2(:,n2ind)).^2); |
---|
| 161 | % Find indices of points in p2 that are within distance dmax of |
---|
| 162 | % p1(:,n1) |
---|
| 163 | n2indmod = n2ind(find(dists2 < dmax^2)); |
---|
| 164 | end |
---|
| 165 | |
---|
| 166 | % Calculate noralised correlation measure. Note this gives |
---|
| 167 | % significantly better matches than the unnormalised one. |
---|
| 168 | |
---|
| 169 | for n2 = n2indmod |
---|
| 170 | % Generate window in 2nd image |
---|
| 171 | w2 = im2(p2(1,n2)-r:p2(1,n2)+r, p2(2,n2)-r:p2(2,n2)+r); |
---|
| 172 | cormat(n1,n2) = sum(sum(w1.*w2))/sqrt(sum(sum(w2.*w2))); |
---|
| 173 | end |
---|
| 174 | end |
---|