探究网上的一个用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

探究网上的一个用MATLAB写的SIFT相关推荐

用MATLAB写一个自动生成福利彩票双色球号码的程序
用MATLAB写一个自动生成福利彩票双色球号码的程序规则红色球:1-33号任选6个蓝色球:1-16号任选1个 red = randi([1,33],1,6); disp('红色球'); fpri ...
用Matlab写脚本求解线性方程组，让大家拥有一个线性方程组计算器
用Matlab写脚本求解线性方程组,让大家拥有一个线性方程组计算器有句话说的好:一个合格的程序员,不会写出"摧毁地球"的方法.他会写一个方法叫"摧毁行星",然 ...
用matlab写一个利用人工地震波定位掩埋物的程序
利用人工地震波定位掩埋物是地质勘探中常用的方法之一.下面是一个使用MATLAB编写的基本程序框架,可以用于实现这一目的. 程序主要分为以下几个步骤: 1.读取地震数据文件使用MATLAB中的load ...
在网上看到一个故事，我觉得挺感人
在网上看到一个故事,我觉得挺感人: 戴着一副黑框眼镜,脸上挂着笑脸,这就是吴福胜,来自轻工与食品工程学院096班的湖南小伙子,刚刚取得了2010年度的中信大锰奖学金一等奖,奖金一万元. 孝敬--面前的 ...
匿名内部类探究——它是一个实例
目录我的学习过程匿名内部类概述匿名内部类探究代码验证(匿名内部类是一个实例) 结论我的学习过程昨天想学习一下Java8新特性,看到Lambda表达式可以替代匿名内部类.我对匿名内部类不太理 ...
一个用js写的沙漏程序 hourglass
1.效果演示这是一个用js写的动态沙漏程序,在文本框中输入数字,回车,即可以看到一个相应层数的沙漏在不停的变化. 难度在于上半部分缺失的沙子形成的镂空图形.下半部分多出的沙子形 ...
原来，我连一个URL都写不对…
来自我的铁粉公号 @闹闹吃鱼 " 阅读本文需要 6 分钟" 作者 l 闹闹来源 l 闹闹吃鱼(ID:UpUp005) 又到了每周的分享时间了,今天分享一下在网络协议中,URI的相 ...
MATLAB写歌曲（艺术与科学的 MEET）
MATLAB写歌曲(艺术与科学的 MEET) 本人认为用 matlab 可以极大锻炼对于音乐和艺术的理性科学认识,我编的这个是孙燕姿的<遇见>,以表达这种艺术和科学的完美邂逅. 理解这样一 ...
python如何播放视频_在网上看到一个视频！怎么下载下来啊？
原标题:在网上看到一个视频!怎么下载下来啊? 大家早,我是云景,之前我给大家连续分享过三篇关于下载网页视频的技巧,但是有些同学觉得难!因为里面用到了Python的知识. 其实,我之前也分享过一个更加简 ...

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

探究网上的一个用MATLAB写的SIFT相关推荐

最新文章

热门文章