matlab特征匹配部分程序整理.doc

Affine
function img2 = affine(img, r, sx, sy, a, b, xo, yo)
R = [cos(r) -sin(r) 0;
sin(r) cos(r) 0;
0 0 1];
S = [sx a 0;
b sy 0;
0 0 1];
T2 = [1 0 -xo;
0 1 -yo;
0 0 1];
T1 = [1 0 xo;
0 1 yo;
0 0 1];
A = T1R*S*T2;
[h,w]=size(img);
[x,y] = meshgrid(1:h,1:w);
x=x(:);
y=y(:);
XP = [x';y';ones(1,length(x))];
X = A*XP;
z = find(X(1,:)>w-1);
X(1,z) = w-1;
z = find(X(1,:)<=1);
X(1,z) = 1.5;
z= find(X(2,:)>h-1);
X(2,z) = h-1;
z= find(X(2,:)<1.5);
X(2,z) = 1.5;
x = X(1,:);
y = X(2,:);
alpha = x - floor(x); %calculate alphas and betas for each point
beta = y - floor(y);
fx = floor(x); fy = floor(y);
inw = fy + (fx-1)*h; %create index（索引） for neighboring pixels
ine = fy + fx*h;
isw = fy+1 + (fx-1)*h;
ise = fy+1 + fx*h;
img2 = (1-alpha).*(1-beta).*img(inw) + ... %interpolate同样）
(1-alpha).*beta.*img(isw) + ...
alpha.*(1-beta).*img(ine) + ...
alpha.*beta.*img(ise);
img2 = reshape(img2,h,w);
imagesc(img2);
%build_pyramid
function [pyr,imp] = build_pyramid(img,levels,scl)
img2 = img;
img2 = resample_bilinear(img2,1/2); %expand to retain spatial frequencies
%img2 = imresize(img2,2,'bilinear'); %expand to retain spatial frequencies
sigma=1.5; %variance for laplacian filter
sigma2=1.5; %variance for downsampling
sig_delta = (1.6-.6)/levels;
for i=1:levels
if i==1
img3 = img2;
img2 = filter_gaussian(img2,7,.5); %slightly filter bottom level
end
imp{i}=img2;
A = filter_gaussian(img2,7,sigma); %calculate difference of gaussians
B = filter_gaussian(A,7,sigma);
pyr{i} = A-B; %store result in cell array
if i==1
img2 img3;
else
B = filter_gaussian(img2,7,sigma2); %anti-alias for downsampling
B = filter_gaussian(B,7,sigma2);
end
img2 = resample_bilinear(B,scl); %downsample for next level
end
%show_pyramid(pyr) %show pyramid if desired
%///////////////////////////////////////////////////////////////////////////////
function show_pyramid(pyr)
close all
[h,w] = size(pyr);
for i=1:
figure
imagesc(pyr{i});
colormap gray;
end
3\
function key = construct_key(px, py, img, sz)
pct = .75;
[h,w] = size(img);
[yoff,xoff] = meshgrid(-1:1,-1:1);
yoff = yoff(:)*pct;
xoff = xoff(:)*pct;
for i = 1:size(yoff,1)
ctrx = px + xoff(i)*sz*2; %method using interpolated values
ctry = py + yoff(i)*sz*2;
[y,x] = meshgrid(ctry-sz:sz/3:ctry+sz,ctrx-sz:sz/3:ctrx+sz);
y=y(:);
x=x(:);
t = 0;
c = 0;
for k=1:size(y,1)
if x(k)1 & y(k)1
t = t + interp(img,x(k),y(k));
c=c+1;
end
end
if c==0
c
end
t t/c;
key(i) = t;
end
key = key/sum(key);
4\detect_features
function [features,pyr,imp,keys] = detect_features(img, scl, disp_flag, thresh, radius, radius2, radius3, min_sep, edgeratio)
if ~exist('scl')
scl = 1.5;
end
if ~exist('thresh')
thresh 3;
end
if ~exist('radius')
radius = 4;
end
if ~exist('radius2')
radius2 = 4;
end
if ~exist('radius3')
radius3 = 4;
end
if ~exist('min_sep')
min_sep = .04;
end
if ~exist('edgeratio')
edgeratio = 5;
end

if ~exist('disp_flag')
disp_flag = 0;
end
if size(img,3) > 1
img = rgb2gray(img);
end
% Computation of the maximum number of levels:
Lmax = floor(min(log(2*size(img)/12)/log(scl)));
%build image pyramid and difference of gaussians filter response pyramid
[pyr,imp] = build_pyramid(img,Lmax,scl);
%get the feature points
pts = find_features(pyr,img,scl,thresh,radius,radius2,disp_flag,1);
%classify points and create sub-pixel and sub-scale adjustments
[features,keys] = refine_features(img,pyr,scl,imp,pts,radius3,min_sep,edgeratio);
5\drawbox
function drawbox(ang,sz,xc,yc,clr)
ang=-ang;
b = [1 1; -1 1; -1 -1; 1 -1; 1 1] * sz;
t = [ sin(ang) cos(ang);cos(ang) -sin(ang)];
b = (t*b')';
b(:,2)=b(:,2)+xc;
b(:,1)=b(:,1)+yc;
plot(b(:,2),b(:,1),'color',clr);
6\ eliminate_edges
function [f2,k2] = eliminate_edges( features,keys);%消除边缘
f2 = features(find(features(:,5)>0),:);
k2 = keys(find(features(:,5)>0),:);
7\ f_class
function [cl, or]=f_class(r)
a = [r(length(r)) r r(1) r(2)];
b = (a(1:length(r)+1)-a(2:length(r)+2));
b = b./abs(b);
d = ones(1,length(r)).*(b(1:length(r))>b(2:length(r)+1));
c = find(d);
m = [(r(c))' c'];
m = sortrows(m,1);
or(1) = (m(1,2)-1)*pi/length(r)*2;
cl = 2;
dst = m(:,1)-m(1,1);
i=find(dst < .3*(max(r)-min(r)));
cl = length(i);
or = (m(i,2)-1)*pi/length(r)*2;
8/ filter_gaussian
function image_out = filter_gaussian(img,order,sig)

img2 = img;
for i=1:floor(order/2) %pad image borders with enough for filter order
[h,w] = size(img2);
img2 = [img2(1,1) img2(1,:) img2(1,w);
img2(:,1) img2 img2(:,w);
img2(h,1) img2(h,:) img2(h,w)];
end
f = gauss1d(order,sig); %create filter coefficient matrix%创建滤波器系数矩阵

image_out = conv2(img2,f,'valid'); % do the filtering
image_out = conv2(image_out,f','valid'); % do the filtering
function f = gauss1d(order,sig)
f=0;
i=0;
j=0;
%generate gaussian coefficients %产生高斯系数
for x = -fix(order/2):1:fix(order/2)
i = i 1;
f(i) = 1/2/pi*exp(-((x^2)/(2*sig^2)));
end
f = f / sum(sum(f)); %normalize filter%正常化滤波
function f = gauss2d(order,sig)
f=0;
i=0;
j=0;
%generate gaussian coefficients %产生高斯系数
for x = -fix(order/2):1:fix(order/2)
j=j+1;
i=0;
for y = -fix(order/2):1:fix(order/2)
i=i+1;
f(i,j) = 1/2/pi*exp(-((x^2+y^2)/(2*sig^2)));
end
end
f = f / sum(sum(f)); %normalize filter%正常化滤波
9/ filter_laplacian
function image_out = filter_gaussian(img,order,sig)
h1 = gauss1d(order,sig); %create filter coefficient matrix
h2 = conv2(h1, [.5 0 -.5]);
h3 = conv2(h2,h2);
h3 = h3/sum(abs(h3));
order = length(h3);
img2 = img;
for i=1:floor(order/2) %pad image borders with enough for filter order
[h,w] = size(img2);
img2 = [img2(1,1) img2(1,:) img2(1,w);
img2(:,1) img2 img2(:,w);
img2(h,1) img2(h,:) img2(h,w)];
end
image_out = conv2(img2,h3','valid'); % do the filtering
image_out = image_out(:,floor(order/2)+1:end-floor(order/2));
image_out2 = conv2(img2,h3,'valid'); % do the filtering
image_out2 = image_out2(floor(order/2)+1:end-floor(order/2),:);
image_out = -image_out-image_out2;
function f gauss1d(order,sig)
f=0;
i=0;
j=0;
%generate gaussian coefficients
for x = -fix(order/2):1:fix(order/2)
i = i + 1;
f(i) = 1/2/pi*exp(-((x^2)/(2*sig^2)));
end
f = f / sum(sum(f)); %normalize filter
function f = gauss2d(order,sig)
f=0;
i=0;
j=0;
%generate gaussian coefficients
for x = -fix(order/2):1:fix(order/2)
j=j+1;
i=0;
for y = -fix(order/2):1:fix(order/2)
i=i+1;
f(i,j) = 1/2/pi*exp(-((x^2+y^2)/(2*sig^2)));
end
end
f = f / sum(sum(f)); %normalize filter
10/ find_extrema
function [mx] = find_extrema(img,thresh,radius)
%img abs(img);
[h,w] = size(img);
% get interior image subtracting radius pixels from border
p = img(radius+1:h-radius, radius+1:w-radius);
%get pixels above threshold
[yp,xp] = find(p>thresh);
yp=yp+radius; xp=xp+radius;
pts = yp+(xp-1)*h;
%create offset list for immediate neighborhood

z=ones(3,3);
z(2,2)=0;
[y,x]=find(z);
y=y-2; x=x-2;
if size(pts,2)>size(pts,1)
pts = pts';
end
%test for max within immediate neighborhood
if size(pts,1)>0
maxima=ones(length(pts),1);
for i=1:length(x)
pts2 = pts + y(i) + x(i)*h;
maxima = maxima & img(pts)>img(pts2);
end
xp = xp(find(maxima)); %save maxima
yp = yp(find(maxima));
pts = yp+(xp-1)*h; %create new index list of good points
end
%create offset list for radius of pixels
[y,x]=meshgrid(-20:20,-20:20);
z = (x.^2+y.^2)<=radius^2 & (x.^2+y.^2)>(1.5)^2; %include points within radius without immediate neighborhood
[y,x]=find(z);
x=x-21; y=y-21;
%create offset list for ring of pixels
[y2,x2]=meshgrid(-20:20,-20:20);
z = (x2.^2+y2.^2)<=(radius)^2 & (x2.^2+y2.^2)>(radius-1)^2;
[y2,x2]=find(z);
x2=x2-21; y2=y2-21;
maxima = ones(length(pts),1);
%test within radius of pixels (done after first test for slight speed increase)
if size(pts,1)>0
for i = 1:length(x)
pts2 = pts + y(i) + x(i)*h;
maxima = maxima & img(pts)>img(pts2); %test points
end

xp = xp(find(maxima)); %save maxima from immediate neighborhood
yp = yp(find(maxima));
pts = yp+(xp-1)*h; %create new index list
mx = [yp xp];
else
mx = [];
end
11\ find_features
function maxima = find_features(pyr img, scl, thresh, radius, radius2, disp_flag, img_flag)
% pts = find_features(pyr,img,scl,thresh,radius,radius2,disp_flag,1);
levels = size(pyr);
levels = levels(2);
mcolor = [ 0 1 0; %color array for display of features at different scales
略
[himg,wimg] = size(img); %get size of images
[h,w] = size(pyr{2});
for i=2:levels-1
[h,w] = size(pyr{i});
[h2,w2] = size(pyr{i+1});
%find maxima
mx = find_extrema(pyr{i},thresh,radius); %find maxima at current scale level
mx2 = round((mx-1)/scl) + 1; %find coords in level above
mx_above = neighbor_max(pyr{i},pyr{i+1},mx,mx2,radius2); %do neighbor comparison in scale space above
if i>1
mx2 = round((mx-1)*scl) + 1; %find coords in level below
mx_below = neighbor_max(pyr{i},pyr{i-1},mx,mx2,radius2); %do comparison in scale below
maxima{i} = plist(mx, mx_below & mx_above); %get coord list for retained maxima and minima
else
maxima{i} = plist(mx, mx_above);
end
%find minima
%if i==11,
% keyboard;
%end;
mx = find_extrema(-pyr{i},thresh,radius); %find minima at current scale level
mx2 = round((mx-1)/scl) + 1; %find coords in level above
mx_above = neighbor_max(-pyr{i},-pyr{i+1},mx,mx2,radius2); %do neighbor comparison in scale space above
if i>1
mx2 = round((mx-1)*scl) + 1; %find coords in level below
mx_below = neighbor_max(-pyr{i},-pyr{i-1},mx,mx2,radius2); %do comparison in scale below
mxtemp = plist(mx, mx_below & mx_above); %get coord list for retained maxima and minima
else
mxtemp = plist(mx, mx_above);
end
maxima{i} = [maxima{i}; mxtemp]; %combine maxima and minima into list for return
%display results if desired
if disp_flag > 0
figure
if img_flag == 0
tmp=resample_bilinear(img,himg/h);
imagesc(tmp);
colormap gray;
show_plist(maxima{i},mcolor(mod(i-1,7)+1,:),'+');
else
imagesc(pyr{i});
colormap gray;
show_plist(maxima{i},mcolor(mod(i-1,7)+1,:),'+');
end
end
end
% Compare a vector of pixels with its neighbors in another scale
function v = neighbor_max(img1,img2,i,i2,radius) % i and i2 are column vectors of r,c coords
if (size(i2,1))==0 | size(img2,1)<11 | size(img2,2)<11
v=zeros(length(i),1);
else
[h,w] = size(img1);
[h2,w2] = size(img2);
[y,x]=meshgrid(-20:20,-20:20); %create set of offsets within radius
z = (x.^2+y.^2)<=radius^2;
[y,x]=find(z);
x=x-21; y=y-21;
radius=ceil(radius);
bound = ones(size(i2,1),2)*[h2-radius 0;0 w2-radius]; %create boundary listing
i2 = i2 - ((i2 > bound).*(i2-bound+1)); %test bounds to make all points within image
i2 = i2 + ((i2 < radius+1).*(radius-i2+1));
i2 = vec(i2,h2); %create indices from x,y coords
i = vec(i,h);
p = img1(i);
res = ones(length(i),1);
for j=1:length(x) %check against all points within radius
itest = i2 + x(j) + h2*y(j);
p2 = img2(itest);
res = res & (p>=p2);
end
v = res; %store results in binary vector
end
function v = vec(points,h)
y = points(:,1);
x = points(:,2);
v = y + (x-1)*h; %create index vectors
function p = plist(points, flags)
p = points(find(flags),:);
12\ fit_parabola
function [c]=fit_parabola(x1, x2, x3, y1, y2, y3)
z = [x1^2 x1 1; x2^2 x2 1; x3^2 x3 1];
c = inv(z)*[y1; y2; y3];
13/ fit_paraboloid
function [yctr,xctr] = fit_paraboloid(data)
%Create Solution Matrix
[yp,xp] = meshgrid(-1:1,-1:1);
yp=yp(:); xp=xp(:);
M=zeros(9,6);
for i=1:length(yp)
x=xp(i); y=yp(i);
M(i,:) = [x^2 y^2 x*y x y 1];
end
g = data';
g = g(:);
cf = inv(M'*M)*M'*g;
a=cf(1); b=cf(2); c=cf(3); d=cf(4); e=cf(5);
xctr = (-2*b*d+e*c)/(4*a*b-c*c);
yctr = (-2*a*e-d*c)/(4*a*b-c*c);
%A = [a c/2; c/2 b];
%[u,s,v] = svd(A);
14/ fmransac_test2
clear p im hc hcm hcmf hcmfi;
% Play with these parameters:
keytype = 2; % The type of feature key to use
gdorder = 4; % The highest order Gaussian derivative filters to use
gdsigma = 5; % The standard deviation of the Gaussian derivatives
steer = 0; % Steer filters in gradient direction for rotational invariance
% Create an image with lots of "inside" and "outside" corners in it

im1 = corners('prob', 0.5);
[ny, nx] = size(im1);
if 0
% Transform the image twice
T = maketform('projective', ...
[0 0; ny 0; ny nx; 0 nx], ...
[0 0; ny 0.3*nx; ny 0.7*nx; 0 nx]);
im2 = imtransform(im1, T, 'FillValues', 255);
else
im2 = affine(im1, pi/2, 1, 1, 0, 0, nx/2, ny/2);
end
% Run the Harris corner detector on both images
p1 = harris_pts(im1, 'smoothing', 1.0, 'nmsrad', 2, 'relthresh', 0.1);
% Project the corners in the first image through the transformation
p2 = fliplr(tformfwd(fliplr(p1), T));
% Compute the feature keys
switch keytype
case 0
% Test