1 | % MATCHBYMONOGENICPHASE - match image feature points using monogenic phase data |
---|
2 | % |
---|
3 | % Function generates putative matches between previously detected |
---|
4 | % feature points in two images by looking for points that have minimal |
---|
5 | % differences in monogenic phase data within windows surrounding each point. |
---|
6 | % Only points that correlate most strongly with each other in *both* |
---|
7 | % directions are returned. This is a simple-minded N^2 comparison. |
---|
8 | % |
---|
9 | % This matcher performs rather well relative to normalised greyscale |
---|
10 | % correlation. Typically there are more putative matches found and fewer |
---|
11 | % outliers. There is a greater computational cost in the pre-filtering stage |
---|
12 | % but potentially the matching stage is much faster as each pixel is effectively |
---|
13 | % encoded with only 3 bits. (Though this potential speed is not realized in this |
---|
14 | % implementation) |
---|
15 | % |
---|
16 | % Usage: [m1,m2] = matchbymonogenicphase(im1, p1, im2, p2, w, dmax, ... |
---|
17 | % nscale, minWaveLength, mult, sigmaOnf) |
---|
18 | % |
---|
19 | % Arguments: |
---|
20 | % im1, im2 - Images containing points that we wish to match. |
---|
21 | % p1, p2 - Coordinates of feature pointed detected in im1 and |
---|
22 | % im2 respectively using a corner detector (say Harris |
---|
23 | % or phasecong2). p1 and p2 are [2xnpts] arrays though |
---|
24 | % p1 and p2 are not expected to have the same number |
---|
25 | % of points. The first row of p1 and p2 gives the row |
---|
26 | % coordinate of each feature point, the second row |
---|
27 | % gives the column of each point. |
---|
28 | % w - Window size (in pixels) over which the phase bit codes |
---|
29 | % around each feature point are matched. This should |
---|
30 | % be an odd number. |
---|
31 | % dmax - Maximum search radius for matching points. Used to |
---|
32 | % improve speed when there is little disparity between |
---|
33 | % images. Even setting it to a generous value of 1/4 of |
---|
34 | % the image size gives a useful speedup. |
---|
35 | % nscale - Number of filter scales. |
---|
36 | % minWaveLength - Wavelength of smallest scale filter. |
---|
37 | % mult - Scaling factor between successive filters. |
---|
38 | % sigmaOnf - Ratio of the standard deviation of the Gaussian |
---|
39 | % describing the log Gabor filter's transfer function in |
---|
40 | % the frequency domain to the filter center frequency. |
---|
41 | % |
---|
42 | % |
---|
43 | % Returns: |
---|
44 | % m1, m2 - Coordinates of points selected from p1 and p2 |
---|
45 | % respectively such that (putatively) m1(:,i) matches |
---|
46 | % m2(:,i). m1 and m2 are [2xnpts] arrays defining the |
---|
47 | % points in each of the images in the form [row;col]. |
---|
48 | % |
---|
49 | % |
---|
50 | % I have had good success with the folowing parameters: |
---|
51 | % |
---|
52 | % w = 11; Window size for correlation matching, 7 or greater |
---|
53 | % seems fine. |
---|
54 | % dmax = 50; |
---|
55 | % nscale = 1; Just one scale can give very good results. Adding |
---|
56 | % another scale doubles computation |
---|
57 | % minWaveLength = 10; |
---|
58 | % mult = 4; This is irrelevant if only one scale is used. If you do |
---|
59 | % use more than one scale try values in the range 2-4. |
---|
60 | % sigmaOnf = .2; This results in a *very* large bandwidth filter. A |
---|
61 | % large bandwidth seems to be very important in the |
---|
62 | % matching performance. |
---|
63 | % |
---|
64 | % See Also: MATCHBYCORRELATION, MONOFILT |
---|
65 | |
---|
66 | % Copyright (c) 2005 Peter Kovesi |
---|
67 | % School of Computer Science & Software Engineering |
---|
68 | % The University of Western Australia |
---|
69 | % http://www.csse.uwa.edu.au/ |
---|
70 | % |
---|
71 | % Permission is hereby granted, free of charge, to any person obtaining a copy |
---|
72 | % of this software and associated documentation files (the "Software"), to deal |
---|
73 | % in the Software without restriction, subject to the following conditions: |
---|
74 | % |
---|
75 | % The above copyright notice and this permission notice shall be included in |
---|
76 | % all copies or substantial portions of the Software. |
---|
77 | % |
---|
78 | % The Software is provided "as is", without warranty of any kind. |
---|
79 | |
---|
80 | % May 2005 - Original version adapted from matchbycorrelation.m |
---|
81 | |
---|
82 | |
---|
83 | function [m1,m2,cormat] = matchbymonogenicphase(im1, p1, im2, p2, w, dmax, ... |
---|
84 | nscale, minWaveLength, mult, sigmaOnf) |
---|
85 | |
---|
86 | orientWrap = 0; |
---|
87 | [f1, h1f1, h2f1, A1] = ... |
---|
88 | monofilt(im1, nscale, minWaveLength, mult, sigmaOnf, orientWrap); |
---|
89 | |
---|
90 | [f2, h1f2, h2f2, A2] = ... |
---|
91 | monofilt(im2, nscale, minWaveLength, mult, sigmaOnf, orientWrap); |
---|
92 | |
---|
93 | % Normalise filter outputs to unit vectors (should also have masking for |
---|
94 | % unreliable filter outputs) |
---|
95 | for s = 1:nscale |
---|
96 | % f1{s} = f1{s}./A1{s}; f2{s} = f2{s}./A2{s}; |
---|
97 | % h1f1{s} = h1f1{s}./A1{s}; h1f2{s} = h1f2{s}./A2{s}; |
---|
98 | % h2f1{s} = h2f1{s}./A1{s}; h2f2{s} = h2f2{s}./A2{s}; |
---|
99 | |
---|
100 | % Try quantizing oriented phase vector to 8 octants to see what |
---|
101 | % effect this has (Performance seems to be reduced only slightly) |
---|
102 | f1{s} = sign(f1{s}); f2{s} = sign(f2{s}); |
---|
103 | h1f1{s} = sign(h1f1{s}); h1f2{s} = sign(h1f2{s}); |
---|
104 | h2f1{s} = sign(h2f1{s}); h2f2{s} = sign(h2f2{s}); |
---|
105 | end |
---|
106 | |
---|
107 | % Generate correlation matrix |
---|
108 | cormat = correlationmatrix(f1, h1f1, h2f1, p1, ... |
---|
109 | f2, h1f2, h2f2, p2, w, dmax); |
---|
110 | |
---|
111 | [corrows,corcols] = size(cormat); |
---|
112 | |
---|
113 | % Find max along rows give strongest match in p2 for each p1 |
---|
114 | [mp2forp1, colp2forp1] = max(cormat,[],2); |
---|
115 | |
---|
116 | % Find max down cols give strongest match in p1 for each p2 |
---|
117 | [mp1forp2, rowp1forp2] = max(cormat,[],1); |
---|
118 | |
---|
119 | % Now find matches that were consistent in both directions |
---|
120 | p1ind = zeros(1,length(p1)); % Arrays for storing matched indices |
---|
121 | p2ind = zeros(1,length(p2)); |
---|
122 | indcount = 0; |
---|
123 | for n = 1:corrows |
---|
124 | if rowp1forp2(colp2forp1(n)) == n % consistent both ways |
---|
125 | indcount = indcount + 1; |
---|
126 | p1ind(indcount) = n; |
---|
127 | p2ind(indcount) = colp2forp1(n); |
---|
128 | end |
---|
129 | end |
---|
130 | |
---|
131 | % Trim arrays of indices of matched points |
---|
132 | p1ind = p1ind(1:indcount); |
---|
133 | p2ind = p2ind(1:indcount); |
---|
134 | |
---|
135 | % Extract matched points from original arrays |
---|
136 | m1 = p1(:,p1ind); |
---|
137 | m2 = p2(:,p2ind); |
---|
138 | |
---|
139 | |
---|
140 | %------------------------------------------------------------------------- |
---|
141 | % Function that does the work. This function builds a 'correlation' matrix |
---|
142 | % that holds the correlation strength of every point relative to every other |
---|
143 | % point. While this seems a bit wasteful we need all this data if we want |
---|
144 | % to find pairs of points that correlate maximally in both directions. |
---|
145 | |
---|
146 | function cormat = correlationmatrix(f1, h1f1, h2f1, p1, ... |
---|
147 | f2, h1f2, h2f2, p2, w, dmax) |
---|
148 | |
---|
149 | if mod(w, 2) == 0 | w < 3 |
---|
150 | error('Window size should be odd and >= 3'); |
---|
151 | end |
---|
152 | |
---|
153 | r = (w-1)/2; % 'radius' of correlation window |
---|
154 | |
---|
155 | [rows1, npts1] = size(p1); |
---|
156 | [rows2, npts2] = size(p2); |
---|
157 | |
---|
158 | if rows1 ~= 2 | rows2 ~= 2 |
---|
159 | error('Feature points must be specified in 2xN arrays'); |
---|
160 | end |
---|
161 | |
---|
162 | % Reorganize monogenic phase data into a 4D matrices for convenience |
---|
163 | |
---|
164 | [im1rows,im1cols] = size(f1{1}); |
---|
165 | [im2rows,im2cols] = size(f2{1}); |
---|
166 | nscale = length(f1); |
---|
167 | phase1 = zeros(im1rows,im1cols,nscale,3); |
---|
168 | phase2 = zeros(im2rows,im2cols,nscale,3); |
---|
169 | |
---|
170 | for s = 1:nscale |
---|
171 | phase1(:,:,s,1) = f1{s}; phase1(:,:,s,2) = h1f1{s}; phase1(:,:,s,3) = h2f1{s}; |
---|
172 | phase2(:,:,s,1) = f2{s}; phase2(:,:,s,2) = h1f2{s}; phase2(:,:,s,3) = h2f2{s}; |
---|
173 | end |
---|
174 | |
---|
175 | % Initialize correlation matrix values to -infinity |
---|
176 | cormat = repmat(-inf, npts1, npts2); |
---|
177 | |
---|
178 | % For every feature point in the first image extract a window of data |
---|
179 | % and correlate with a window corresponding to every feature point in |
---|
180 | % the other image. Any feature point less than distance 'r' from the |
---|
181 | % boundary of an image is not considered. |
---|
182 | |
---|
183 | % Find indices of points that are distance 'r' or greater from |
---|
184 | % boundary on image1 and image2; |
---|
185 | n1ind = find(p1(1,:)>r & p1(1,:)<im1rows+1-r & ... |
---|
186 | p1(2,:)>r & p1(2,:)<im1cols+1-r); |
---|
187 | |
---|
188 | n2ind = find(p2(1,:)>r & p2(1,:)<im2rows+1-r & ... |
---|
189 | p2(2,:)>r & p2(2,:)<im2cols+1-r); |
---|
190 | |
---|
191 | for n1 = n1ind |
---|
192 | |
---|
193 | % Identify the indices of points in p2 that we need to consider. |
---|
194 | if dmax == inf |
---|
195 | n2indmod = n2ind; % We have to consider all of n2ind |
---|
196 | |
---|
197 | else % Compute distances from p1(:,n1) to all available p2. |
---|
198 | p1pad = repmat(p1(:,n1),1,length(n2ind)); |
---|
199 | dists2 = sum((p1pad-p2(:,n2ind)).^2); |
---|
200 | % Find indices of points in p2 that are within distance dmax of |
---|
201 | % p1(:,n1) |
---|
202 | n2indmod = n2ind(find(dists2 < dmax^2)); |
---|
203 | end |
---|
204 | |
---|
205 | % Generate window in 1st image |
---|
206 | w1 = phase1(p1(1,n1)-r:p1(1,n1)+r, p1(2,n1)-r:p1(2,n1)+r, :, :); |
---|
207 | |
---|
208 | for n2 = n2indmod |
---|
209 | % Generate window in 2nd image |
---|
210 | w2 = phase2(p2(1,n2)-r:p2(1,n2)+r, p2(2,n2)-r:p2(2,n2)+r, :, :); |
---|
211 | % Compute dot product as correlation measure |
---|
212 | cormat(n1,n2) = w1(:)'*w2(:); |
---|
213 | |
---|
214 | % *** Need to add mask stuff |
---|
215 | end |
---|
216 | end |
---|
217 | |
---|