[37] | 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 | |
---|