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 |
---|