探究网上的一个用MATLAB写的SIFT

这是第一个M文件 建立灰度图像的金字塔和高斯差分金字塔

function [pyr,imp] = build_pyramid(img,levels,scl)

img2 = img;

img2 = resample_bilinear(img2,1/2); %扩展成空间频域

%img2 = imresize(img2,2,'bilinear'); %expand to retain spatial frequencies

sigma=1.5; %拉普拉斯算子的放大系数

sigma2=1.5; %平滑系数

%levels=7；

%sig_delta = (1.6-.6)/levels;

scl=1.5;

%for i=1:levels

for i=1:7

if i==1

img3 = img2;

img2 = filter_gaussian(img2,7,.5); %轻度底层过滤器

end

imp{i}=img2;

A = filter_gaussian(img2,7,sigma); %计算高斯差分函数

B = filter_gaussian(A,7,sigma);

pyr{i} = A-B; %在单元序列中存储结果

if i==1

img2 = img3;

else

B = filter_gaussian(img2,7,sigma2); %向下采样

B = filter_gaussian(B,7,sigma2);

end

img2 = resample_bilinear(B,scl); %下一层进行采样

end

其中 filter_gaussian 和 resample_bilinear是另两个M文件 如下：

function image_out = filter_gaussian(img,order,sig)%输入为原始图像，输出为卷积后的图像尺度空间

img2 = img;

for i=1:floor(order/2) %补充图像边界使其能足够进行滤波处理

[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); %构造滤波系数矩阵

image_out = conv2(img2,f,'valid'); % 滤波处理

image_out = conv2(image_out,f','valid'); % 滤波处理

%/////////////////////////////////////////////////////////////////////////////////////////

function f = gauss1d(order,sig)%高斯函数：order为层数，sig为尺度因子

f=0;

i=0;

j=0;

%生成高斯系数

for x = -fix(order/2):1:fix(order/2)%fix取整

i = i + 1;

f(i) = 1/2/pi*exp(-((x^2)/(2*sig^2)));

end

f = f / sum(sum(f)); %归一化滤波

%/////////////////////////////////////////////////////////////////////////////////////////

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

function img2 = resample_bilinear(img, ratio)

img=double(img);

[h,w]=size(img); %获取图象大小

[y,x]=meshgrid( 1:ratio:h-1, 1:ratio:w-1 ); %为新图象建立X,Y值

[h2,w2] = size(x); %将图象设为维

x = x(:); %向量变换

y = y(:);

alpha = x - floor(x); %给每个点计算像素

beta = y - floor(y);

fx = floor(x); fy = floor(y);

inw = fy + (fx-1)*h; %给相邻点建立坐标值

ine = fy + fx*h;

isw = fy+1 + (fx-1)*h;

ise = fy+1 + fx*h;

img2 = (1-alpha).*(1-beta).*img(inw) + (1-alpha).*beta.*img(isw) + alpha.*(1-beta).*img(ine) + alpha.*beta.*img(ise);

img2 = reshape(img2,h2,w2); %转换回2维

由以上两个M文件 我在Command Window中输入下面的 得到DOG差分金字塔：

A=imread('F:\orl_zhifangtu\s3.jpg');

[pyr,imp]=build_pyramid(A,7,1.5);

[h,w]=size(pyr);

for i=1:w

figure

imagesc(pyr{i});

colormap gray;

end















顺序发乱了，从figure1到figure7分别就是第一层到第七层

下面的几个M函数是经过同尺度，上一个尺度，下一个尺度比较得到极大极小值点，从而得到特征点：

%/////////////////////////////////////////////////////////////////////////////////////////////

%

% find_features - scale space feature detector based upon difference of gaussian filters.

% selects features based upon their maximum response in scale space

%

% 关系式: maxima = find_features(pyr, img, thresh, radius, radius2, min_sep, edgeratio, disp_flag, img_flag)

%

% Parameters:

% pyr : 图像金字塔滤波后的单元序列 (built with build_pyramid)

% img : 原始图像 (only used for visualization)

% thresh : 用于搜索极值的阈值(minimum filter response considered)

% radius : 在当前尺度上，极值比较的范围

% radius2: 极值与邻近尺度进行比较的范围

% disp_flag: 1- 显示各个尺度上的离散图像. 0 - 不显示

% img_flag: 1 - 显示滤波反馈. 0 - 显示原始图像.

%

% 返回值:

%

% maxima - c一个由所选择的特征点组成的nX2 行列点阵

%

% Author:

% Scott Ettinger

% scott.m.ettinger@intel.com

%

% May 2002

%/////////////////////////////////////////////////////////////////////////////////////////////

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); %得到层数 7

%在不同尺度上的特征显示的颜色序列

mcolor = [ 0 1 0;0 1 0;1 0 0;.2 .5 0;0 0 1;1 0 1;0 1 1;1 .5 0;.5 1 0;0 1 .5;.5 1 .5];

[himg,wimg] = size(img); %获得图像大小

[h,w] = size(pyr{2});

for i=2:levels-1

[h,w] = size(pyr{i});

[h2,w2] = size(pyr{i+1});

%找到maxima点

mx = find_extrema(pyr{i},thresh,radius); %在当前尺度上找到极大值点

mx2 = round((mx-1)/scl) + 1; %找到上一层的坐标

mx_above = neighbor_max(pyr{i},pyr{i+1},mx,mx2,radius2); %在上一层尺度空间上进行极值比较

if i>1

mx2 = round((mx-1)*scl) + 1; %找到下一层的坐标

mx_below = neighbor_max(pyr{i},pyr{i-1},mx,mx2,radius2); %在下一层尺度空间上进行极值比较

maxima{i} = plist(mx, mx_below & mx_above); %所获得极值的坐标列表

else

maxima{i} = plist(mx, mx_above);

end

%find minima

%if i==11,

% keyboard;

%end;

mx = find_extrema(-pyr{i},thresh,radius); %在当前尺度上找到极小值点

mx2 = round((mx-1)/scl) + 1; %找到上一层的坐标

mx_above = neighbor_max(-pyr{i},-pyr{i+1},mx,mx2,radius2); %在上一层尺度空间上进行极值比较

if i>1

mx2 = round((mx-1)*scl) + 1; %找到下一层的坐标

mx_below = neighbor_max(-pyr{i},-pyr{i-1},mx,mx2,radius2); %在下一层尺度空间上进行极值比较

mxtemp = plist(mx, mx_below & mx_above); %所获得极值的坐标列表

else

mxtemp = plist(mx, mx_above);

end

maxima{i} = [maxima{i}; mxtemp]; %合并最大值和最小值，放入列表返回

%显示结果

% disp_flag=1;

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

%/////////////////////////////////////////////////////////////////////////////////////////////

%

% 找出极值 -

% 在一个灰度尺度图像上寻找局部极值。再一次检查在给出的半径内所有像素值来确认局部极值。在高于给出的阈值的像素值将会认为是一个有效的极值。

%

% 关系式: m = find_extrema(img,thresh,radius,min_separation)

%

% 参数:

% img : 图像矩阵

% thresh : 阈值

% radius : 像素区域

%

% Returns:

%

% m - 一个由所选择的特征点组成的nX2 行列点阵

% Author:

% Scott Ettinger

% scott.m.ettinger@intel.com

%

% May 2002

%/////////////////////////////////////////////////////////////////////////////////////////////

function [mx] = find_extrema(img,thresh,radius)%输入为图像矩阵，输出为用一个nX2行列点阵表示的极值点的坐标

%img = abs(img);

[h,w] = size(img);

% 获取内部图像与边界区域像素值的差

p = img(radius+1:h-radius, radius+1:w-radius);

% 获得极值点的像素值

[yp,xp] = find(p>thresh);

yp=yp+radius; xp=xp+radius;

pts = yp+(xp-1)*h;

%给最近邻点构造反馈列表

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

%在最近邻里找最大值

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)); %保存最大极值 这里不应该是[xp,yp]=find(maxima); ???

yp = yp(find(maxima));

pts = yp+(xp-1)*h; %构造选好的特征点的索引列表

end

%为区域像素构造反馈列表

[y,x]=meshgrid(-20:20,-20:20);

z = (x.^2+y.^2)<=radius^2 & (x.^2+y.^2)>(1.5)^2; %包括区域内的点，但不包括最近邻的点

[y,x]=find(z);

x=x-21;

y=y-21;

%为环形内的像素构造反馈列表

[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);

%检测区域内的像素点 (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); %测试点

end

xp = xp(find(maxima)); %从最近邻获取的极值点

yp = yp(find(maxima));

pts = yp+(xp-1)*h; %创建新的索引

mx = [yp xp];

else

mx = [];

end

%//////////////////////////////////////////////////////////////////////////////////////////////

%

% 在不同的尺度上将像素点的向量与它的近邻进行比较

%

%//////////////////////////////////////////////////////////////////////////////////////////////

function v = neighbor_max(img1,img2,i,i2,radius) % i和i2是坐标的列向量

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); %创建区域内点的反馈集合

z = (x.^2+y.^2)<=radius^2;

[y,x]=find(z);

x=x-21; y=y-21;

radius=ceil(radius); %得到大于等于radius的整数

bound = ones(size(i2,1),2)*[h2-radius 0;0 w2-radius]; %创建分界线

i2 = i2 - ((i2 > bound).*(i2-bound+1)); %测试边界使每个点在图像里

i2 = i2 + ((i2 < radius+1).*(radius-i2+1));

i2 = vec(i2,h2); %创建索引

i = vec(i,h);

p = img1(i);

res = ones(length(i),1);

for j=1:length(x) %再一次检查点在区域内

itest = i2 + x(j) + h2*y(j);

p2 = img2(itest);

res = res & (p>=p2);

end

v = res; %在二值向量中存储结果

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),:);

我在Command Window中输入下面的 得到极值点 即最初的特征点显示出来是这样

A=imread('F:\orl_zhifangtu\s3.jpg');

[pyr,imp]=build_pyramid(A,7,1.5);

pts = find_features(pyr,A,1.5,3,4,4,1,1);

结果









顺序发反了 可以看到第一层金字塔上特征点最多

如果看下面两个M文件：

%/////////////////////////////////////////////////////////////////////////////////////////////

%

% detect_features - scale space feature detector based upon difference of gaussian filters.

% selects features based upon their maximum response in scale space

% 特征检测 - 基于高斯差分滤波器的尺度空间特征检测。在尺度极值空间中选择特征。

% Usage关系式: [features,pyr,imp,keys] = detect_features(img, scl, disp_flag, thresh, radius, radius2, radius3, min_sep, edgeratio)

%

% Parameters:

% 参数：

% img : 原始图像

% scl : 图像金字塔层之间的尺度比例因子

% thresh : 极值搜索的阈值(minimum filter response considered)

% radius : radius for maxima comparison within current scale

% 与当前尺度进行比较的最大值半径

% radius2: radius for maxima comparison between neighboring scales

% 近邻尺度比较得出的最大值半径

% radius3: radius for edge rejection test 边缘抑制测试半径

% min_sep : minimum separation for maxima selection.最大值选择上的最小分割点

% edgeratio: maximum ratio of eigenvalues of feature curvature for edge rejection.

% 对于边缘抑制，特征曲率与特征值的最大比例

% disp_flag: 1- display each scale level on separate figure. 0 - no display

% 1表示：显示各尺度的离散图像 0表示：不显示

% Returns:

% 返回值：

% features - matrix with one row for each feature consisting of the following:

% [x position, y position, scale(sub-level), size of feature on image, edge flag,

% edge orientation, curvature of response through scale space ]

% 特征：一个矩阵，矩阵中的每一行代表一个特征，组成：x位置，y位置，子级别的尺度，图像上特征的大小，边缘标志，边缘方向，尺度空

% 间的曲率响应。

% pyr, imp - filter response and image pyramids

% 滤波响应和图像金字塔

% keys - key values generated for each feature by construct_key.m

% 由construct_key.m为每个特征产生的关键值。

% Notes: 推荐参数值

% recommended parameter values are:

% scl = 1.5; thresh = 3; radius = 4; radius2 = 4; radius3 = 4; min_sep = .04; edgeratio = 5;

%

% Author:

% Scott Ettinger

% scott.m.ettinger@intel.com

%

% May 2002

%//////////////////////////////////////////////////////////////////////////

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

% img=imread('ima1.jpg');

if size(img,3) > 1

img = rgb2gray(img);

end

% 计算层数的极值

Lmax = floor(min(log(2*size(img)/12)/log(scl)));

% 建立图像金字塔和高斯差分滤波器反馈金字塔

[pyr,imp] = build_pyramid(img,Lmax,scl);

% 获得特征点

pts = find_features(pyr,img,scl,thresh,radius,radius2,disp_flag,1);

%对点进行分类并创造基本像素点和基尺度调节器 这是什么？??

[features,keys] = refine_features(img,pyr,scl,imp,pts,radius3,min_sep,edgeratio);

其中最后一个函数是：

%/////////////////////////////////////////////////////////////////////////////////////////////

%

% refine_features - scale space feature detector based upon difference of gaussian filters.

% selects features based upon their maximum response in scale space

%

% Usage: features = refine_features(img, pyr, scl, imp, pts, radius, min_separation, edgeratio)

%

% Parameters:

%

% img : original image

% pyr : cell array of filtered image pyramid

% scl : scaling factor between levels of the image pyramid

% imp : image pyramid cell array

% pts : cell array of selected points on each pyramid level

% radius : radius for edge rejection test

% min_separation : minimum separation distance for maxima rejection.

% edgeratio: maximum ratio of eigenvalues of feature curvature for edge rejection.

%

% Returns:

%

% features - matrix with one row for each feature consisting of the following:

% [x loc, y loc, scale value, size, edge flag, edge orientation, scale space curvature]

%

% where:

% x loc and y loc are the x and y positions on the original image

% scale value is the sub level adjusted scale value

% size is the size of the feature in pixels on the original image

% edge flag is zero if the feature is classified as an edge

% edge orientation is the angle made by the edge through the feature point

% scale space curvature is a rough confidence measure of feature prominence

%

% Author:

% Scott Ettinger

% scott.m.ettinger@intel.com

%

% May 2002

%/////////////////////////////////////////////////////////////////////////////////////////////

function [features,keys] = refine_features(img, pyr, scl, imp,pts, radius, min_separation, edgeratio)

[ho,wo]=size(img);

[h2,w2]=size(imp{2});

%给环上的像素点创建反馈

[ry2,rx2]=meshgrid(-20:20,-20:20);

z=(rx2.^2+ry2.^2)<=(radius)^2 & (rx2.^2+ry2.^2)>(radius-1)^2;

[ry2,rx2]=find(z);

rx2=rx2-21;

ry2=ry2-21;

ndx = 1;

%在金字塔的每层上进行循环

for j=2:length(imp)-1

p=pts{j}; %得到当前层的滤波反馈

img=imp{j}; %得到当前层的图像

[dh,dw]=size(img);

np = 1;

min_sep = min_separation*max(max(pyr{j})); %对于有效的极大值计算极小值分割

for i=1:size(p,1)

ptx = p(i,2);

pty = p(i,1);

if p(i,1) < radius+3 %确保近邻在图像区域内部

p(i,1) = radius+3;

end

if p(i,1) > dh-radius-3

p(i,1) = dh-radius-3;

end

if p(i,2) < radius+3

p(i,2) = radius+3;

end

if p(i,2) > dw-radius-3

p(i,2) = dw-radius-3;

end

%adjust to sub pixel maxima location

r = pyr{j}(pty-1:pty+1,ptx-1:ptx+1); %获得3X3的相邻像素

[pcy,pcx] = fit_paraboloid(r); %找到抛物线中心点所对应的点

if abs(pcy)>1 %由于抛物线拟合中的奇异性，忽略极值的偏移

pcy=0;

pcx=0;

end

if abs(pcx)>1

pcx=0;

pcy=0;

end

p(i,1) = p(i,1) + pcy; %调整中心点

p(i,2) = p(i,2) + pcx;

ptx = p(i,2);

pty = p(i,1);

px=(pts{j}(i,2)+pcx - 1)*scl^(j-2) + 1; %计算出点在金字塔第二层中的位置

py=(pts{j}(i,1)+pcy - 1)*scl^(j-2) + 1;

%calculate Sub-Scale level adjustment

y1 = interp(pyr{j-1},(p(i,2)-1)*scl+1, (p(i,1)-1)*scl+1); %利用插值在周围尺度层上获取反馈

y3 = interp(pyr{j+1},(p(i,2)-1)/scl+1, (p(i,1)-1)/scl+1);

y2 = interp(pyr{j},p(i,2),p(i,1));

coef = fit_parabola(0,1,2,y1,y2,y3); % 拟合相邻3个尺度的点到抛物线上

scale_ctr = -coef(2)/2/coef(1); %在尺度空间上找最大值

if abs(scale_ctr-1)>1 %由于抛物线拟合的奇异性忽略极值

scale_ctr=0;

end

%消除边缘点且执行最小值分割

rad2 = radius * scl^(scale_ctr-1); %调整半径大小来计算新的尺度值

o=0:pi/8:2*pi-pi/8; %在半径周围的测试点创建环形点

rx = (rad2)*cos(o);

ry = (rad2)*sin(o);

rmax = 1e-9; %初始化最大值和最小值

rmin = 1e9;

gp_flag = 1;

pval = interp(pyr{j},ptx,pty); %在特征点中心获得反馈

rtst = [];

%在特征点的环形周围检测点

for k=1:length(rx)

rtst(k) = interp(pyr{j},ptx+rx(k),pty+ry(k)); %用二次线性插值来得到环形点值

if pval> 0 %对环形中的每个点计算它们到特征点的距离

rtst(k) = pval - rtst(k);

else

rtst(k) = rtst(k) - pval;

end

gp_flag = gp_flag * rtst(k)>min_sep; %在有背景噪音的情况下测试有效极大值

if rtst(k)>rmax %找到极小和极大值

rmax = rtst(k);

end

if rtst(k)<rmin

rmin = rtst(k);

end

end

fac = scl/(wo*2/size(pyr{2},2)); %由于采样的边缘效应计算大小偏移

[cl,or] = f_class(rtst); %特征分类以及获取方向

if rmax/rmin > edgeratio %测试边界

gp_flag=0;

if cl ~= 2 %保留所有的交叉 (# ridges > 2)

gp_flag = 1;

else

ang = min(abs([(or(1)-or(2)) (or(1)-(or(2)-2*pi))]));

if ang < 6.5*pi/8; %保留角度大于145度的边缘

gp_flag = 1;

end

end

end

%save info

features(ndx,1) = (px-1)*wo/w2*fac +1; %保存X,Y值 (sub pixel adjusted)

features(ndx,2) = (py-1)*wo/w2*fac +1;

features(ndx,3) = j+scale_ctr-1; %保存尺度值 (sub scale adjusted in units of pyramid level)

features(ndx,4) = ((7-4)*scl^(j-2+scale_ctr-1))*wo/w2*fac; %保存原始图象的大小

features(ndx,5) = gp_flag; %保存边界标记

features(ndx,6) = or(1); %保存边界方向角

features(ndx,7) = coef(1); %保存尺度空间响应的曲率

%if features(ndx,1) > wo

%px 因为不知道这里是什么 就暂时把它们全做注释看

%end

keys(ndx,:) = construct_key(ptx,pty,imp{j},3 * scl^(scale_ctr-1));

ndx = ndx + 1;

end

end

其中的M函数

function v = interp(img,xc,yc) %点之间的双线性插值

px = floor(xc);

py = floor(yc);

alpha = xc-px;

beta = yc-py;

nw = img(py,px);

ne = img(py,px+1);

sw = img(py+1,px);

se = img(py+1,px+1);

v = (1-alpha)*(1-beta)*nw + ... %插值

(1-alpha)*beta*sw + ...

alpha*(1-beta)*ne + ...

alpha*beta*se;

其中另一个M函数

function key = construct_key(px, py, img, sz)

pct = .75;

%img=imread('ima1.jpg');

[h,w] = size(img);

[yoff,xoff] = meshgrid(-1:1,-1:1);

yoff = yoff(:)*pct;

xoff = xoff(:)*pct;

%px=0;

%py=0;

for i = 1:size(yoff,1)

ctrx = px + xoff(i)*sz*2; %采用插值方法

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)<w-1 && x(k)>1 && y(k)<h-1 && y(k)>1

t = t + interp(img,x(k),y(k));

c=c+1;

end

end

if c==0

continue; %It is right???It is "continue"?

else t = t/c;

end

key(i) = t;

end

key = key/sum(key);

在Command Window中这样输入：

A=imread('F:\orl_zhifangtu\s3.jpg');

[pyr,imp]=build_pyramid(A,7,1.5);

[features,pyr,imp,keys] = detect_features(A, 1.5, 1, 3, 4, 4, 4, 0.04, 5);

结果如下：











在Command Window中这样输入时：

imgo1=imread('F:\orl_zhifangtu\s1.jpg');

img1=imresize(imgo1,2,'bilinear');

imgo2=imread('F:\orl_zhifangtu\s3.jpg');

img2=imresize(imgo2,2,'bilinear');

>> [f1,pyr,imp,k1] = detect_features(img1);

[f1,k1] = eliminate_edges(f1,k1);

figure(1);

showfeatures(f1,img1,0)

[f2,pyr,imp,k2] = detect_features(img2);

[f2,k2] = eliminate_edges(f2,k2);

figure(2);

showfeatures(f2,img2,0)

结果：




如果我没理解错
这两个显示出来的是特征点吧 可是为什么是这样

其中包括的函数有：

function [f2,k2] = eliminate_edges( features,keys);%对由detect_features检测出的特征，消除其中的边缘特征，得到非边缘的特征及其关键值

f2 = features(find(features(:,5)>0),:);

k2 = keys(find(features(:,5)>0),:);

还有一个M文件是：

%/////////////////////////////////////////////////////////////////////////////////////////////

%

% showfeatures - visualize features generated by detect_features

% Usage: showfeatures(features,img)

%

% Parameters:

%

% img : original image

% features: matrix generated by detect_features

%

% Returns: nothing, generates figure

%

% Author:

% Scott Ettinger

% scott.m.ettinger@intel.com

%

% May 2002

%//////////////////////////////////////////////////////////////////////////

function [] = showfeatures(features,img,num_flag)

if ~exist('num_flag')

num_flag = 0;

end

figure(gcf);

imagesc(img)

hold on

colormap gray

for i=1:size(features,1)

x = features(i,1);

y = features(i,2);

sz = features(i,4);

if x>size(img,2)

x %这里错误的 我也不知道作者这里是要写什么？

end

if features(i,5) > 0 %检测边缘标记

if num_flag ~= 1

plot(x,y,'g+'); %在真正的特征点周围画方格

else

text(x,y,sprintf('%d',i),'color','m');

end

if abs(features(i,7))>1.8

drawbox(0,sz,x,y,[0 0 1]);

else

drawbox(0,sz,x,y,[0 .9 .2]);

end

else %画边界

ang = features(i,6);

px = [x-sz*cos(ang) x+sz*cos(ang)];

py = [y-sz*sin(ang) y+sz*sin(ang)];

plot(px,py,'r');

end

end

hold off;

如果在Command Window中这样输入：

imgo1=imread('F:\orl_zhifangtu\s1.jpg');

img=imresize(imgo1,2,'bilinear');

>> test

则结果：


我真的不知道这出来的又是什么鬼

其中test函数是：

%example file for feature detector

%close all

testpat=0; %set to zero to use current contents of img variable

if testpat==1

img=zeros(400,400);

img(150:200,100:150)=200;

img(100:110,100:110)=200;

img(100:118,140:158)=200;

img(100:120,185:205)=200;

img(100:125,235:260)=200;

img(100:130,290:320)=200;

img(150:185,40:75)=200;

img(190:290,190:290)=200;

img(300:302,100:102)=200;

img=filter_gaussian(img,7,sqrt(2));%输入为原始图像，高斯金字塔层数7以及尺度因子sigma为sqrt（2）；输出为图像尺度空间

img=img+randn(400,400)*5;

elseif testpat==2

img=zeros(400,400);

img(200:400,190:210) = 200;

img=filter_gaussian(img,7,sqrt(2));

img=img+randn(400,400)*5;

elseif testpat==3

img=zeros(400,400);

img(1:400,190:210) = 200;

img=filter_gaussian(img,7,sqrt(2));

img=img+randn(400,400)*5;

elseif testpat==4

img=zeros(400,400);

img(200:400,190:210) = 200;

img(200:220,200:400) = 200;

img=filter_gaussian(img,7,sqrt(2));

img=img+randn(400,400)*0;

elseif testpat==5

img=zeros(400,400);

img(200:400,200:210) = 200;

img(200:210,1:400) = 200;

img=filter_gaussian(img,7,sqrt(2));

img=img+randn(400,400)*5;

elseif testpat==6

img=zeros(400,400);

img(1:400,200:210) = 200;

img(200:210,1:400) = 200;

img=filter_gaussian(img,7,sqrt(2));

img=img+randn(400,400)*5;

img=-img;

end

close all

if size(img,3)>1

img = rgb2gray(img);

end

%imagesc(img);

threshold=3; %Threshold value for rejecting maxima/minima

disp_flag = 0; %change to a zero for a combined view of all scales

img_flag = 1; %change to a zero to see features plotted on original image

radius = 4;

radius2 = 4;

radius3 = 4;

min_sep = .04;

edgeratio = 5;

scl = 1.5;

%img = imread(file_name);

%[pyr,imp] = build_pyramid(img,12,scl);

%pts = find_features(pyr,img,scl,threshold,radius,radius2,min_sep,edgeratio,disp_flag,img_flag);

%[features] = getpts(img,pyr,scl,imp,pts,6,radius3,min_sep,edgeratio);

[features,pyr,imp,keys] = detect_features(img,scl,disp_flag,threshold,radius,radius2,radius3, min_sep,edgeratio);

showfeatures(features,img);

axis equal;

%enjoy...

后面还有几十个M函数 我真的理不清了 感觉后面有的和SIFT根本没关系 什么读取每一帧 还有计算时间什么的 哎呀好烦啊 这个SIFT！！！所有的M文件都在这里了，要看的话自己下载吧 我是已经绝望了啊啊啊啊啊 我已经无能为力了 SIFT我该拿你肿么办？
http://download.csdn.net/detail/wd1603926823/8916599