上QQ阅读APP看书，第一时间看更新
1.5 SIFT特征提取与描述
人类在识别一个物体时，不管这个物体是远还是近，都能对它进行正确的辨识，这就是所谓的“尺度不变性”。同样，当这个物体发生旋转时，我们照样可以正确地识别它，这就是所谓的“旋转不变性”……那么，如何让机器与人类一样具有这种能力呢？这就是图像局部不变性特征要解决的问题。对于图像不变性特征提取的方法，核心是“不变性”三个字。
提到图像局部不变性特征，有两个人不得不提及，一个是Lindeberg，另一个是Lowe。如果将局部不变性特征方法比作一个孩子，那么，Lindeberg就是这个孩子的父亲，Lowe是母亲。Lindeberg奠定了局部不变性特征方法的理论基础，播下了局部不变性特征方法的种子，而Lowe则将这颗种子孕育成一种能具体实现的方法。由于在Lindeberg的尺度空间理论中，各种高斯微分算子与哺乳动物的视网膜和视觉皮层的感受域剖面有着高度的相似性，所以，尺度空间理论经常与生物视觉相关联，有人也称图像局部不变性特征方法为基于生物视觉的不变性方法。
1.5.1 SIFT算法
Lowe总结了现有的特征提取方法，并提出了一种新型高效的特征检测描述方法——尺度不变特征变换（Scale Invariant Feature Transform, SIFT）。SIFT的主要思路是：首先建立图像的尺度空间表示，然后在尺度空间中搜索图像的极值点，由极值点再建立特征描述向量，最后用特征描述向量进行相似度匹配。采用SIFT方法提取的图像特征具有放缩不变性、旋转不变性、仿射不变性，还有一定的抗光照变化和抗视点变换性能。该特征还具有高度的可区分性，能够在一个具有大量特征数据的数据库中精确的匹配。SIFT还使用金字塔算法大大缩小了提取特征时的运算量。计算方法分为四个步骤：尺度空间极值提取、特征点定位、特征方向赋值和提取特征点描述。
1．尺度空间极值检测
尺度空间理论是检测不变性特征的基础。Witkin（1983）提出了尺度空间理论，他主要讨论了一维信号平滑处理的问题。Koenderink（1984）把这种理论扩展到二维图像，并证明高斯卷积核是实现尺度变换的唯一变换核。
一幅二维图像在不同尺度下的尺度空间表示可由图像与高斯核卷积得到：
其中，G（x,y,σ）为高斯核函数：
其中，（x,y）为图像点的像素坐标，I（x,y）为图像数据。σ称为尺度空间因子，它是高斯正态分布的方差，反映了图像被平滑的程度，其值越小表征图像被平滑程度越小，相应尺度也越小。L（x,y,σ）代表了图像的尺度空间。
为高效地在尺度空间内检测出稳定的特征点，Lowe使用尺度空间中差分高斯（Difference of Gaussian, DoG）极值作为判断依据。DoG算子定义为两个不同尺度的高斯核的差分，是归一化（Laplacian of Gaussian, LoG）算子的近似。设k为两个相邻尺度间的比例因子，则DoG算子定义如下：
D（x,y,σ）构造方式如图1.5.1所示，其建立高斯图像（图1.5.2中左列）与DoG（图1.5.2中右列）两个金字塔。高斯图像金字塔分为多组，每组间又分为多层。一组中的多层间不同的是尺度，相邻层间尺度相差一个比例因子k。在S个尺度间隔内变化尺度因子，如使加倍，则k应为
。为了在整个金字塔内获取DoG极值，应在高斯金字塔中生成S+3层高斯平滑图像。下一组图像的最底层由上一组中尺度为2σ的图像进行因子为2的降采样得到，其中σ为上一组中最底层图像的尺度因子。DoG金字塔由相邻的高斯图像金字塔相减得到。所生成的高斯图像金字塔与DoG金字塔如图1.5.2所示。
图1.5.1 高斯图像金字塔（S＝2）与DoG金字塔
图1.5.2 所生成的高斯图像金字塔与DoG金字塔
金字塔中每个高斯图像的σ为：
其中，σ0为基础尺度因子；o、s分别为图像所在的图像组坐标、组内层坐标，o∈omin+［0, …,O－1］,s∈［0, …,S－1］;omin是第一个金字塔组的坐标，通常omin取0或者－1，当设为－1的时候，则图像在计算高斯尺度空间前先扩大一倍。在Lowe的算法实现中，以上参数的取值为：σ0＝1.6×, omin＝－1,S＝3。
此外，空间坐标x是组坐标o 的函数。设x0为第0组内的空间坐标，则有：x＝2ox0,x0∈［0, …,N0－1］×［0, …,M0－1］，其中（N0,M0）是第0组中图像的分辨率。设（N0,M0）为第0组中的图像分辨率，则其他组的分辨率为：
金字塔构造完后开始检测DoG局部极值。其中每个像素需要跟同一尺度的周围邻域8个像素和相邻尺度对应位置的周围邻域9×2个像素总共26个像素进行比较，如图1.5.3所示。仅当被检测点的DoG值大于此26个像素点或小于此26个像素点时才将该点判定为极值点并保存以进行后续计算。
图1.5.3 DoG空间局部极值检测
2．确定关键点位置及尺度
通过拟合三维二次函数以精确确定关键点的位置和尺度，得到原图像SIFT候选特征点集合，如图1.5.4所示。
图1.5.4 原图像的SIFT候选特征点集合
在得到原图像的SIFT候选特征点集合X0后，要从中筛选稳定的点作为该图像的SIFT特征点，其组成的集合用X表示。由于X0中对比度较低的点对噪声较敏感，而位于边缘上的点难以准确定位，为确保SIFT特征点的稳定性，必须将这两种点剔除。
1）剔除对比度低的点
将候选特征点x的偏移量定义为Δx，其对比度为D（x）的绝对值D（x），对x的DoG函数表达式（1.5.3）进行泰勒级数展开：
式中，由于x为DoG函数的极值点，所以
＝0，解方程得Δx＝
，通过多次迭代得到最终候选点的精确位置及尺度
，将其代入式（1.5.6）求得D（
），求其绝对值可得
。设对比度阈值为Tc，则低对比度点的剔除公式为：
2）剔除边缘点
由于边缘梯度方向上主曲率值较大，而边缘方向上曲率较小，在边缘上得到的DoG函数的极值点与非边缘区域的点相比，其主曲率比值较大，因此可以将主曲率比值大于一定阈值的点视为位于边缘上的点将其剔除。
由于候选点的DoG函数D（x）的主曲率与2×2的Hessian矩阵H的特征值成正比：
式（1.5.8）中，Dxx、Dxy、Dyy为候选点邻域对应位置的像素差分。令α为H 的最大特征值，β为H的最小特征值，令γ＝
，则D（x）的主曲率比值与γ成正比。由H的迹和行列式的值可得：
只与两特征值之比有关，与特征值自身大小无关，当两特征值相等时最小，且随着γ的增大而增大。设主曲率比值阈值为Tγ，则边缘点的剔除公式为：
剔除低对比度和边缘点后所剩的特征点如图1.5.5所示。
3）关键点方向确定
通过确定关键点的方向，可以使特征描述符以与方向相关的方式构造，从而使算子具有旋转不变性。关键点的方向利用其邻域像素的梯度分布特性确定。对每个高斯图像，每个点L（x,y）的梯度的模m（x,y）与方向θ（x,y）可以通过下式计算得到：
其中，L（x,y）所用尺度为关键点所在尺度。
图1.5.5 剔除低对比度和边缘点后所剩的特征点
1.5.2 SIFT特征描述
对每个关键点，在以其为中心的邻域窗口内利用直方图的方式统计邻域像素的梯度分布。此直方图有36个柱，每柱10°，共360°。每个加入直方图的邻域像素样本的权重由该像素的梯度模与高斯权重确定，此高斯窗的σ为关键点的尺度的1.5倍。加入高斯窗的目的是增强离关键点近的邻域点对关键点的影响。
直方图的峰值反映了关键点所处邻域梯度的主方向。完成直方图统计后，找到直方图的最高峰值以确定关键点的方向。关键点的方向可以由离最高峰值最近的三个柱值通过抛物线插值精确得到。
在梯度方向直方图中，当存在一个大于等于主峰值80%能量的峰值时，则添加一个新关键点，此关键点的坐标、尺度与当前关键点相同，但方向为由此峰值确定的方向。因此一个关键点可能产生多个坐标、尺度相同，方向不同的关键点。这样做的目的是增强匹配的鲁棒性。
至此，特征点检测完毕，特征描述前的准备工作已经完成。每个关键点含有三个信息：坐标、尺度、方向。
首先将坐标轴旋转为关键点的方向，以确保旋转不变性。
接下来以关键点为中心取8×8的窗口。图1.5.6左部分的中心为当前关键点的位置，每个小格代表关键点邻域所在尺度空间的一个像素，箭头方向代表该像素的梯度方向，箭头长度代表梯度模值，图中的圈代表高斯加权的范围（越靠近关键点的像素梯度方向信息贡献越大）。然后在每4×4的小块上计算8个方向的梯度方向直方图，绘制每个梯度方向的累加值，即可形成一个种子点，如图1.5.6右部分所示。此图中一个关键点由2×2共4个种子点组成，每个种子点有8个方向向量信息。这种邻域方向性信息联合的思想增强了算法抗噪声的能力，同时对于含有定位误差的特征匹配也提供了较好的容错性。
实际计算过程中，为了增强匹配的稳健性，Lowe建议对每个关键点使用4×4共16个种子点来描述，这样对于一个关键点就可以产生128个数据，即最终形成128维的SIFT特征向量。此时SIFT特征向量已经去除了尺度变化、旋转等几何变形因素的影响，再继续将特征向量的长度归一化，则可以进一步去除光照变化的影响。
图1.5.6 特征点的特征向量构造
综上所述，提取图像SIFT特征点及其特征向量的流程如图1.5.7所示。
图1.5.7 提取SIFT特征点及其特征向量的流程
1.5.3 实例精讲
【例1.5.1】 SIFT算法实现的MATLAB源程序，读者可以根据程序及相关的注释对SIFT作进一步的理解。
********************************************************************
function[pos, scale, orient, desc]＝ SIFT(im, octaves, intervals, object＿mask, contrast
＿threshold, curvature＿threshold, interactive)
% 功能：提取灰度图像的尺度不变特征(SIFT特征)
% 输入：
% im——灰度图像，该图像的灰度值在0～1范围内
% octaves——金字塔的组数：octaves(默认值为4)
% intervals——该输入参数决定每组金字塔的层数(默认值为 2)
% object＿mask——确定图像中尺度不变特征点的搜索区域
% contrast＿threshold——对比度阈值(默认值为0.03)
% curvature＿threshold——曲率阈值(默认值为10.0)
% interactive——函数运行显示标志，将其设定为1，则显示算法运行时间和过程的相关信息
% 如果将其设定为2，则仅显示最终运行结果(default＝ 1)
% 输出：
% pos——Nx2 矩阵，每一行包括尺度不变特征点的坐标(x, y)
% scale——Nx3 矩阵，每一行包括尺度不变特征点的尺度信息
% orient——Nx1 向量，每个元素是特征点的主方向，其范围为[－pi, pi)
% desc——Nx128 矩阵，每一行包含特征点的特征向量
% 参考文献：
%[1]David G.Lowe, ＂Distinctive Image Features from Sacle－Invariant Keypoints＂,
% accepted for publication in the International Journal of Computer
% Vision,2004.
%[2]David G.Lowe, ＂Object Recognition from Local Scale－Invariant Features＂,
% Proc.of the International Conference on Computer Vision, Corfu,
% September 1999.
% 设定输入量的默认值
if ～exist(＇octaves＇)
octaves ＝ 4;
end
if ～exist(＇intervals＇)
intervals ＝ 2;
end
if～exist(＇object＿mask＇)
object＿mask＝ ones(size(im));
end
if size(object＿mask)～＝ size(im)
object＿mask＝ ones(size(im));
end
if～exist(＇contrast＿threshold＇)
contrast＿threshold＝ 0.02;
end
if～exist(＇curvature＿threshold＇)
curvature＿threshold＝ 10.0;
end
if ～exist(＇interactive＇)
interactive ＝ 1;
end
% 检验输入灰度图像的像素灰度值是否已归一化到[0,1]
if((min(im(:)) 0)｜(max(im(:)) 1))
fprintf(2, ＇Warning:image not normalized to[0,1].＼n＇);
end
% 将输入图像经过高斯平滑处理，采用双线性差值将其扩大一倍
if interactive ＝ 1
fprintf(2, ＇Doubling image size for first octave...＼n＇);
end
tic;
antialias＿sigma＝ 0.5;
if antialias＿sigma＝＝ 0
signal ＝ im;
else
g＝ gaussian＿filter(antialias＿sigma);
if exist(＇corrsep＇)＝＝ 3
signal ＝ corrsep(g, g, im);
else
signal ＝ conv2(g, g, im, ＇same＇);
end
end
signal ＝ im;
[X Y]＝ meshgrid(1:0.5:size(signal,2),1:0.5:size(signal,1));
signal ＝ interp2(signal, X, Y, ＇*linear＇);
% 降采样率；
subsample ＝[0.5];
%%%%%%下一步是生成高斯和差分高斯(DoG)金字塔%%%%%%
%%%%%%高斯金字塔含有s+3层，差分高斯金字塔含有s+2层%%%%%%
if interactive ＝ 1
fprintf(2, ＇Prebluring image...＼n＇);
end
preblur＿sigma＝ sqrt(sqrt(2)＾2 －(2*antialias＿sigma)＾2);
if preblur＿sigma＝＝ 0
gauss＿pyr{1,1}＝ signal;
else
g＝ gaussian＿filter(preblur＿sigma);
if exist(＇corrsep＇)＝＝ 3
gauss＿pyr{1,1}＝ corrsep(g, g, signal);
else
gauss＿pyr{1,1}＝ conv2(g, g, signal, ＇same＇);
end
end
clear signal
pre＿time＝ toc;
if interactive ＝ 1
fprintf(2, ＇Preprocessing time %.2f seconds.＼n＇, pre＿time);
end
% 第一组第一层的sigma
initial＿sigma＝ sqrt((2*antialias＿sigma)＾2 + preblur＿sigma ＾2);
% 记录每一层和每一个尺度的sigma
absolute＿sigma＝ zeros(octaves, intervals+3);
absolute＿sigma(1,1)＝ initial＿sigma* subsample(1);
% 记录产生金字塔的滤波器的尺寸和标准差
filter＿size＝ zeros(octaves, intervals+3);
filter＿sigma＝ zeros(octaves, intervals+3);
% 生成高斯和差分高斯金字塔
if interactive ＝ 1
fprintf(2, ＇Expanding the Gaussian and DOG pyramids...＼n＇);
end
tic;
for octave ＝ 1:octaves
if interactive ＝ 1
fprintf(2, ＼＇ tProcessing octave %d:image size %d x %d subsample %.1f＼n＇, octave,
size(gauss＿pyr{octave,1},2), size(gauss＿pyr{octave,1},1), subsample(octave));
fprintf(2, ＇＼t＼tInterval 1 sigma %f＼n＇, absolute＿sigma(octave,1));
end
sigma＝ initial＿sigma;
g＝ gaussian＿filter(sigma);
filter＿size(octave,1)＝ length(g);
filter＿sigma(octave,1)＝ sigma;
DOG＿pyr{octave}＝ zeros(size(gauss＿pyr{octave,1},1), size(gauss＿pyr{octave,1},2),
intervals+2);
for interval ＝ 2:(intervals+3)
% 计算生成下一层几何采样金字塔的标准差
% 其中，sigma＿i+1 ＝ k*sigma.
% sigma＿i+1 ＾2 ＝ sigma＿f, i ＾2 + sigma＿i ＾2
% (k*sigma＿i)＾2 ＝ sigma＿f, i ＾2 + sigma＿i ＾2
% 因此
% sigma＿f, i＝ sqrt(k ∧2 － 1)sigma＿i
% 对于扩展的组(span the octave), k＝ 2 ＾(1/intervals)
% 所以
% sigma＿f, i＝ sqrt(2 ＾(2/intervals)－ 1)sigma＿i
sigma＿f＝ sqrt(2 ＾(2/intervals)－ 1)*sigma;
g＝ gaussian＿filter(sigma＿f);
sigma＝(2 ＾(1/intervals))*sigma;
% 记录sigma的值
absolute＿sigma(octave, interval)＝ sigma* subsample(octave);
% 记录滤波器的尺寸和标准差
filter＿size(octave, interval)＝ length(g);
filter＿sigma(octave, interval)＝ sigma;
if exist(＇corrsep＇)＝＝ 3
gauss＿pyr{octave, interval}＝ corrsep(g, g, gauss＿pyr{octave, interval－1});
else
gauss＿pyr{octave, interval}＝ conv2(g, g, gauss＿pyr{octave, interval－1}, ＇same＇);
end
DOG＿pyr{octave}(:, :, interval－1) ＝ gauss＿pyr{octave, interval} － gauss＿pyr
{octave, interval－1};
if interactive ＝ 1
fprintf(2, ＇＼t＼tInterval % d sigma % f＼n ＇, interval, absolute＿sigma(octave,
interval));
end
end
if octave octaves
sz ＝ size(gauss＿pyr{octave, intervals+1});
[X Y]＝ meshgrid(1:2:sz(2),1:2:sz(1));
gauss＿pyr{octave+1,1}＝ interp2(gauss＿pyr{octave, intervals+1}, X, Y, ＇*nearest＇);
absolute＿sigma(octave+1,1)＝ absolute＿sigma(octave, intervals+1);
subsample ＝[subsample subsample(end)*2];
end
end
pyr＿time＝ toc;
if interactive ＝ 1
fprintf(2, ＇Pryamid processing time %.2f seconds.＼n＇, pyr＿time);
end
% 在交互模式下显示高斯金字塔
if interactive ＝ 2
sz ＝ zeros(1,2);
sz(2)＝(intervals+3)*size(gauss＿pyr{1,1},2);
for octave ＝ 1:octaves
sz(1)＝ sz(1)+ size(gauss＿pyr{octave,1},1);
end
pic ＝ zeros(sz);
y＝ 1;
for octave ＝ 1:octaves
x ＝ 1;
sz ＝ size(gauss＿pyr{octave,1});
for interval ＝ 1:(intervals + 3)
pic(y:(y+sz(1)－1), x:(x+sz(2)－1))＝ gauss＿pyr{octave, interval};
x ＝ x + sz(2);
end
y＝ y+ sz(1);
end
fig ＝ figure;
clf;
imshow(pic);
resizeImageFig(fig, size(pic),0.25);
fprintf(2, ＇The gaussian pyramid(0.25 scale).＼nPress any key to continue.＼n＇);
pause;
close(fig)
end
% 在交互模式下显示差分高斯金字塔
if interactive ＝ 2
sz ＝ zeros(1,2);
sz(2)＝(intervals+2)*size(DOG＿pyr{1}(:, :,1),2);
for octave ＝ 1:octaves
sz(1)＝ sz(1)+ size(DOG＿pyr{octave}(:, :,1),1);
end
pic ＝ zeros(sz);
y＝ 1;
for octave ＝ 1:octaves
x ＝ 1;
sz ＝ size(DOG＿pyr{octave}(:, :,1));
for interval ＝ 1:(intervals + 2)
pic(y:(y+sz(1)－1), x:(x+sz(2)－1))＝ DOG＿pyr{octave}(:, :, interval);
x ＝ x + sz(2);
end
y＝ y+ sz(1);
end
fig ＝ figure;
clf;
imagesc(pic);
resizeImageFig(fig, size(pic),0.25);
fprintf(2, ＇The DOG pyramid(0.25 scale).＼nPress any key to continue.＼n＇);
pause;
close(fig)
end
% 下一步是查找差分高斯金字塔中的局部极值，并通过曲率和照度进行检验
curvature＿threshold＝((curvature＿threshold+ 1)＾2)/curvature＿threshold;
% 二阶微分核
xx ＝[1 －2 1];
yy＝ xx＇;
xy＝[1 0 －1;0 0 0; －1 0 1]/4;
raw＿keypoints＝[];
contrast＿keypoints＝[];
curve＿keypoints＝[];
% 在高斯金字塔中查找局部极值
if interactive ＝ 1
fprintf(2, ＇Locating keypoints...＼n＇);
end
tic;
loc＝ cell(size(DOG＿pyr));
for octave ＝ 1:octaves
if interactive ＝ 1
fprintf(2, ＇＼tProcessing octave %d＼n＇, octave);
end
for interval ＝ 2:(intervals+1)
keypoint＿count＝ 0;
contrast＿mask＝ abs(DOG＿pyr{octave}(:, :, interval)) ＝ contrast＿threshold;
loc{octave, interval}＝ zeros(size(DOG＿pyr{octave}(:, :, interval)));
if exist(＇corrsep＇)＝＝ 3
edge ＝ 1;
else
edge＝ ceil(filter＿size(octave, interval)/2);
end
for y＝(1+edge):(size(DOG＿pyr{octave}(:, :, interval),1)－edge)
for x＝(1+edge):(size(DOG＿pyr{octave}(:, :, interval),2)－edge)
if object＿mask(round(y*subsample(octave)), round(x*subsample(octave)))＝＝ 1
if((interactive ＝ 2)｜(contrast＿mask(y, x)＝＝ 1))
% 通过空间核尺度检测最大值和最小值
DOG＿pyr{octave}((y－1):(y+1), (x－1):(x+1), (interval－1):(interval+1));
pt＿val ＝ tmp(2,2,2);
if((pt＿val ＝＝ min(tmp(:)))｜(pt＿val ＝＝ max(tmp(:))))
% 存储对灰度大于对比度阈值的点的坐标
raw＿keypoints＝[raw＿keypoints; x*subsample(octave)y*
subsample(octave)];
if abs(DOG＿pyr{octave}(y, x, interval)) ＝ contrast＿threshold
contrast＿keypoints＝[contrast＿keypoints; raw＿keypoints(end, :)];
% 计算局部极值的Hessian矩阵
Dxx＝ sum(DOG＿pyr{octave}(y, x－1:x+1, interval).* xx);
Dyy＝ sum(DOG＿pyr{octave}(y－1:y+1, x, interval).* yy);
Dxy＝ sum(sum(DOG＿pyr{octave}(y－1:y+1, x－1:x+1, interval).* xy));
% 计算Hessian矩阵的直迹和行列式
Tr＿H＝ Dxx+ Dyy;
Det＿H＝ Dxx*Dyy－ Dxy ＾2;
% 计算主曲率
curvature＿ratio＝(Tr＿H ＾2)/Det＿H;
if((Det＿H ＝ 0)&(curvature＿ratio curvature＿threshold))
% 存储主曲率小于阈值的极值点的坐标(非边缘点)
curve＿keypoints＝[curve＿keypoints; raw＿keypoints(end, :)];
% 将该点的位置的坐标设为 1，并计算点的数量
loc{octave, interval}(y, x)＝ 1;
keypoint＿count＝ keypoint＿count+ 1;
end
end
end
end
end
end
end
if interactive ＝ 1
fprintf(2, ＇＼t＼t%d keypoints found on interval %d＼n＇, keypoint＿count, interval);
end
end
end
keypoint＿time＝ toc;
if interactive ＝ 1
fprintf(2, ＇Keypoint location time %.2f seconds.＼n＇, keypoint＿time);
end
% 在交互模式下显示特征点检测的结果
if interactive ＝ 2
fig ＝ figure;
clf;
imshow(im);
hold on;
plot(raw＿keypoints(:,1), raw＿keypoints(:,2), ＇y+ ＇);
resizeImageFig(fig, size(im),2);
fprintf(2, ＇DOG extrema(2x scale).＼nPress any key to continue.＼n＇);
pause;
close(fig);
fig ＝ figure;
clf;
imshow(im);
hold on;
plot(contrast＿keypoints(:,1), contrast＿keypoints(:,2), ＇y+ ＇);
resizeImageFig(fig, size(im),2);
fprintf(2, ＇Keypoints after removing low contrast extrema(2x scale).＼nPress any key to
continue.＼n＇);
pause;
close(fig);
fig ＝ figure;
clf;
imshow(im);
hold on;
plot(curve＿keypoints(:,1), curve＿keypoints(:,2), ＇y+ ＇);
resizeImageFig(fig, size(im),2);
fprintf(2, ＇Keypoints after removing edge points using principal curvature filtering(2x
scale).＼nPress any key to continue.＼n＇);
pause;
close(fig);
end
clear raw＿keypoints contrast＿keypoints curve＿keypoints
%%%%%% 下一步是计算特征点的主方向%%%%%%
% 在特征点的一个区域内计算其梯度直方图
g＝ gaussian＿filter(1.5 * absolute＿sigma(1, intervals+3)/subsample(1));
zero＿pad＝ ceil(length(g)/2);
% 计算高斯金字塔图像的梯度方向和幅值
if interactive ＝ 1
fprintf(2, ＇Computing gradient magnitude and orientation...＼n＇);
end
tic;
mag＿thresh＝ zeros(size(gauss＿pyr));
mag＿pyr＝ cell(size(gauss＿pyr));
grad＿pyr＝ cell(size(gauss＿pyr));
for octave ＝ 1:octaves
for interval ＝ 2:(intervals+1)
% 计算x, y的差分
diff＿x ＝ 0.5*(gauss＿pyr{octave, interval}(2:(end－1),3:(end))－gauss＿pyr
{octave, interval}(2:(end－1),1:(end－2)));
diff＿y ＝ 0.5*(gauss＿pyr{octave, interval}(3:(end),2:(end－1))－gauss＿pyr
{octave, interval}(1:(end－2),2:(end－1)));
% 计算梯度幅值
mag＝ zeros(size(gauss＿pyr{octave, interval}));
mag(2:(end－1),2:(end－1))＝ sqrt(diff＿x.＾ 2 + diff＿y.＾ 2);
% 存储高斯金字塔梯度幅值
mag＿pyr{octave, interval}＝ zeros(size(mag)+2*zero＿pad);
mag＿pyr{octave, interval}((zero＿pad+1):(end－zero＿pad), (zero＿pad+1):(end－zero＿pad))
＝ mag;
% 计算梯度主方向
grad＝ zeros(size(gauss＿pyr{octave, interval}));
grad(2:(end－1),2:(end－1))＝ atan2(diff＿y, diff＿x);
grad(find(grad ＝＝ pi))＝ －pi;
% 存储高斯金字塔梯度主方向
grad＿pyr{octave, interval}＝ zeros(size(grad)+2*zero＿pad);
grad＿pyr{octave, interval}((zero＿pad+1):(end－zero＿pad), (zero＿pad+1):(end－zero＿pad))
＝ grad;
end
end
clear mag grad
grad＿time＝ toc;
if interactive ＝ 1
fprintf(2, ＇Gradient calculation time %.2f seconds.＼n＇, grad＿time);
end
%%%%%%下一步是确定特征点的主方向%%%%%%
%%%%%% 方法：通过寻找每个关键点的子区域内梯度直方图的峰值来解决%%%%%%
% 将灰度直方图分为36等分，每隔 10度一份
num＿bins ＝ 36;
hist＿step＝ 2*pi/num＿bins;
hist＿orient＝[－pi:hist＿step:(pi－hist＿step)];
% 初始化关键点的位置、方向和尺度信息
pos ＝[];
orient ＝[];
scale ＝[];
% 给关键点确定主方向
if interactive ＝ 1
fprintf(2, ＇Assigining keypoint orientations...＼n＇);
end
tic;
for octave ＝ 1:octaves
if interactive ＝ 1
fprintf(2, ＇＼tProcessing octave %d＼n＇, octave);
end
for interval ＝ 2:(intervals + 1)
if interactive ＝ 1
fprintf(2, ＇＼t＼tProcessing interval %d ＇, interval);
end
keypoint＿count＝ 0;
% 构造高斯加权掩模
g＝ gaussian＿filter(1.5 * absolute＿sigma(octave, interval)/subsample(octave));
hf＿sz ＝ floor(length(g)/2);
g ＝ g＇*g;
loc＿pad＝ zeros(size(loc{octave, interval})+2*zero＿pad);
loc＿pad((zero＿pad+1):(end－zero＿pad), (zero＿pad+1):(end－zero＿pad))＝
loc{octave, interval};
[iy ix]＝find(loc＿pad＝＝1);
for k ＝ 1:length(iy)
x ＝ ix(k);
y＝ iy(k);
wght＝g.*mag＿pyr{octave, interval}((y－hf＿sz):(y+hf＿sz), (x－hf＿sz):(x+hf＿sz));
grad＿window＝grad＿pyr{octave, interval}((y－hf＿sz):(y+hf＿sz), (x－hf＿sz):(x+
hf＿sz));
orient＿hist＝zeros(length(hist＿orient),1);
for bin＝1:length(hist＿orient)
diff＝ mod(grad＿window－ hist＿orient(bin)+ pi,2*pi)－ pi;
orient＿hist(bin)＝orient＿hist(bin)+sum(sum(wght.*max(1 －abs(diff)/hist＿
step,0)));
end
% 运用非极大抑制法查找主方向直方图的峰值
peaks＝ orient＿hist;
rot＿right＝[peaks(end); peaks(1:end－1)];
rot＿left＝[peaks(2:end); peaks(1)];
peaks(find(peaks rot＿right))＝ 0;
peaks(find(peaks rot＿left))＝ 0;
% 提取最大峰值的值和其索引位置
[max＿peak＿val ipeak]＝ max(peaks);
% 将大于等于最大峰值80%的直方图也确定为特征点的主方向
peak＿val ＝ max＿peak＿val;
while(peak＿val 0.8*max＿peak＿val)
% 最高峰值最近的三个柱值通过抛物线插值精确得到
A ＝[];
b ＝[];
for j ＝ －1:1
A＝[A; (hist＿orient(ipeak)+hist＿step*j).＾2(hist＿orient(ipeak)+
hist＿step*j)1];
bin＝ mod(ipeak+ j + num＿bins－ 1, num＿bins)+ 1;
b＝[b; orient＿hist(bin)];
end
c ＝ pinv(A)*b;
max＿orient＝ －c(2)/(2*c(1));
while(max＿orient －pi)
max＿orient＝ max＿orient+ 2*pi;
end
while(max＿orient ＝ pi)
max＿orient＝ max＿orient－ 2*pi;
end
% 存储关键点的位置、主方向和尺度信息
pos＝[pos; [(x－zero＿pad)(y－zero＿pad)]*subsample(octave)];
orient＝[orient; max＿orient];
scale＝[scale; octave interval absolute＿sigma(octave, interval)];
keypoint＿count＝ keypoint＿count+ 1;
peaks(ipeak)＝ 0;
[peak＿val ipeak]＝ max(peaks);
end
end
if interactive ＝ 1
fprintf(2, ＇(%d keypoints)＼n＇, keypoint＿count);
end
end
end
clear loc loc＿pad
orient＿time＝ toc;
if interactive ＝ 1
fprintf(2, ＇Orientation assignment time %.2f seconds.＼n＇, orient＿time);
end
% 在交互模式下显示关键点的尺度和主方向信息
if interactive ＝ 2
fig ＝ figure;
clf;
imshow(im);
hold on;
display＿keypoints(pos, scale(:,3), orient, ＇y＇);
resizeImageFig(fig, size(im),2);
fprintf(2, ＇Final keypoints with scale and orientation(2x scale).＼nPress any key to
continue.＼n＇);
pause;
close(fig);
end
%%%%%% SIFT算法的最后一步是特征向量生成%%%%%%
orient＿bin＿spacing＝ pi/4;
orient＿angles＝[－pi:orient＿bin＿spacing:(pi－orient＿bin＿spacing)];
grid＿spacing＝ 4;
[x＿coords y＿coords]＝ meshgrid([－6:grid＿spacing:6]);
feat＿grid＝[x＿coords(:)y＿coords(:)]＇;
[x＿coords y＿coords]＝ meshgrid([－(2*grid＿spacing－0.5):(2*grid＿spacing－0.5)]);
feat＿samples＝[x＿coords(:)y＿coords(:)]＇;
feat＿window＝ 2*grid＿spacing;
desc ＝[];
if interactive ＝ 1
fprintf(2, ＇Computing keypoint feature descriptors for %d keypoints＇, size(pos,1));
end
for k ＝ 1:size(pos,1)
x ＝ pos(k,1)/subsample(scale(k,1));
y＝ pos(k,2)/subsample(scale(k,1));
% 将坐标轴旋转为关键点的方向，以确保旋转不变性
M＝[cos(orient(k))－sin(orient(k)); sin(orient(k))cos(orient(k))];
feat＿rot＿grid＝M*feat＿grid+ repmat([x; y],1, size(feat＿grid,2));
feat＿rot＿samples＝M*feat＿samples+ repmat([x; y],1, size(feat＿samples,2));
% 初始化特征向量
feat＿desc＝ zeros(1,128);
for s＝ 1:size(feat＿rot＿samples,2)
x＿sample＝ feat＿rot＿samples(1, s);
y＿sample＝ feat＿rot＿samples(2, s);
% 在采样位置进行梯度插值
[X Y]＝ meshgrid((x＿sample－1):(x＿sample+1), (y＿sample－1):(y＿sample+1));
G＝ interp2(gauss＿pyr{scale(k,1), scale(k,2)}, X, Y, ＇*linear＇);
G(find(isnan(G)))＝ 0;
diff＿x＝ 0.5*(G(2,3)－ G(2,1));
diff＿y＝ 0.5*(G(3,2)－ G(1,2));
mag＿sample＝ sqrt(diff＿x ＾2 + diff＿y ＾2);
grad＿sample＝ atan2(diff＿y, diff＿x);
if grad＿sample＝＝ pi
grad＿sample＝ －pi;
end
% 计算x、y方向上的权重
x＿wght＝ max(1 －(abs(feat＿rot＿grid(1, :)－ x＿sample)/grid＿spacing),0);
y＿wght＝ max(1 －(abs(feat＿rot＿grid(2, :)－ y＿sample)/grid＿spacing),0);
pos＿wght＝ reshape(repmat(x＿wght.*y＿wght,8,1),1,128);
diff＝ mod(grad＿sample－ orient(k)－ orient＿angles+ pi,2*pi)－ pi;
orient＿wght＝ max(1 － abs(diff)/orient＿bin＿spacing,0);
orient＿wght＝ repmat(orient＿wght,1,16);
% 计算高斯权重
g＝ exp(－((x＿sample－x)＾2+(y＿sample－y)＾2)/(2*feat＿window ＾2))/(2*pi*feat＿
window ＾2);
feat＿desc＝ feat＿desc+ pos＿wght.*orient＿wght*g*mag＿sample;
end
% 将特征向量的长度归一化，则可以进一步去除光照变化的影响
feat＿desc＝ feat＿desc/norm(feat＿desc);
feat＿desc(find(feat＿desc 0.2))＝ 0.2;
feat＿desc＝ feat＿desc/norm(feat＿desc);
% 存储特征向量
desc＝[desc; feat＿desc];
if(interactive ＝ 1)&(mod(k,25)＝＝ 0)
fprintf(2, ＇.＇);
end
end
desc＿time＝ toc;
% 调整采样偏差
sample＿offset＝ －(subsample－ 1);
for k ＝ 1:size(pos,1)
pos(k, :)＝ pos(k, :)+ sample＿offset(scale(k,1));
end
if size(pos,1) 0
scale ＝ scale(:,3);
end
% 在交互模式下显示运行过程耗时
if interactive ＝ 1
fprintf(2, ＇＼nDescriptor processing time %.2f seconds.＼n＇, desc＿time);
fprintf(2, ＇Processing time summary:＼n＇);
fprintf(2, ＇＼tPreprocessing:＼t%.2f s＼n＇, pre＿time);
fprintf(2, ＇＼tPyramid:＼t%.2f s＼n＇, pyr＿time);
fprintf(2, ＇＼tKeypoints:＼t%.2f s＼n＇, keypoint＿time);
fprintf(2, ＇＼tGradient:＼t%.2f s＼n＇, grad＿time);
fprintf(2, ＇＼tOrientation:＼t%.2f s＼n＇, orient＿time);
fprintf(2, ＇＼tDescriptor:＼t%.2f s＼n＇, desc＿time);
fprintf(2, ＇Total processing time %.2f seconds.＼n＇, pre＿time+ pyr＿time+ keypoint＿
time+ grad＿time+ orient＿time+ desc＿time);
end
********************************************************************
在上述实现SIFT算法的MATLAB程序中，调用了其他一些函数如下：
********************************************************************
function g＝ gaussian＿filter(sigma)
% 功能：生成一维高斯滤波器
% 输入：
% sigma——高斯滤波器的标准差
% 输出：
% g——高斯滤波器
sample ＝ 7.0/2.0;
n ＝ 2*round(sample * sigma)+1;
x＝1:n;
x＝x－ceil(n/2);
g＝ exp(－(x.＾2)/(2*sigma ＾2))/(sigma*sqrt(2*pi));
********************************************************************
********************************************************************
function im ＝ pgmread(fname);
% 功能：读入．pgm图像
global SIZE＿LIMIT;
if SIZE＿LIMIT
fprintf(1, ＇SIZE＿LIMIT recognized by pgmread＼n＇);
end
[fid, msg]＝ fopen(fname, ＇r＇);
if(fid ＝＝ －1)
error(msg);
end
[pars type]＝ pnmReadHeader(fid);
if(pars＝＝－1)
fclose(fid);
error([fname ＇:cannot parse pgm header＇]);
end
if ～(type ＝＝ ＇P2 ＇｜type ＝＝ ＇P5 ＇)
fclose(fid);
error([fname ＇:Not of type P2 or P5.＇]);
end
xdim ＝ pars(1);
ydim ＝ pars(2);
maxval ＝ pars(3);
sz ＝ xdim * ydim;
if SIZE＿LIMIT
if sz ＝ 16384
ydim ＝ floor(16384/xdim);
sz ＝ xdim * ydim;
fprintf(1, ＇truncated image size:cols %d rows %d＼n＇, xdim, ydim)
end
end
if(type ＝＝ ＇P2 ＇)
[im, count] ＝ fscanf(fid, ＇%d＇, sz);
elseif(type ＝＝ ＇P5 ＇)
count ＝ 0;
im ＝[];
stat ＝ fseek(fid, －sz, ＇eof＇);
if ～stat
[im, count] ＝ fread(fid, sz, ＇uchar＇);
end
if(count ～＝ sz)
fprintf(1, ＇Warning:File ended early! %s＼n＇, fname);
fprintf(1, ＇...Padding with %d zeros.＼n＇, sz－count);
im ＝[im ; zeros(sz－count,1)];
end
else
fclose(fid);
error([fname ＇:Not of type P2 or P5.＇]);
end
im ＝ reshape(im, xdim, ydim)＇;
fclose(fid);
********************************************************************
********************************************************************
function[pars, type]＝ pnmReadHeader(fid)
% 功能：读入并解析一个pnm格式图像的头文件
pars ＝ －1;
type ＝ ＇Unknown＇;
TheLine ＝ fgetl(fid);
szLine ＝ size(TheLine);
endLine ＝ szLine(2);
if(endLine 2)
fprintf(1, [＇Unrecognized PNMfile type＼n＇]);
pars ＝ －1;
return;
end
type ＝ TheLine(1:2);
ok ＝ 0;
if(type(1)＝＝ ＇P＇)
if(type(2)＝＝ ＇1 ＇｜type(2)＝＝ ＇2 ＇｜...
type(2)＝＝ ＇3 ＇｜type(2)＝＝ ＇4 ＇｜...
type(2)＝＝ ＇5 ＇｜type(2)＝＝ ＇6 ＇｜...
type(2)＝＝ ＇B＇｜type(2)＝＝ ＇L＇)
ok ＝ 1;
else
fprintf(1, [＇Unrecognized PNMfile type:＇type ＇＼n＇]);
end
end
if(type(1)＝＝ ＇F＇)
if(type(2)＝＝ ＇P＇｜type(2)＝＝ ＇U＇)
ok ＝ 1;
else
fprintf(1, [＇Unrecognized PNMfile type:＇type ＇＼n＇]);
end
end
if ～ok
pars ＝ －1;
return;
end
current ＝ 3;
parIndex＝1;
while(parIndex 4)
while(current endLine)
TheLine ＝ fgetl(fid);
if(TheLine ＝＝ －1)
fprintf(1, ＇Unexpected EOF＼n＇);
pars ＝ －1;
return;
end
szLine ＝ size(TheLine);
endLine ＝ szLine(2);
current ＝ 1;
end
[token, count, errmsg, nextindex]＝ ...
sscanf(TheLine(current:endLine), ＇%s＇,1);
nextindex ＝ nextindex+current－1;
if(count＝＝0)
if(nextindex endLine)
current ＝ nextindex;
else
pars ＝ －1;
fprintf(1, ＇Unexpected EOF＼n＇);
return;
end
else
if token(1)＝＝ ＇# ＇
current ＝ endLine+1;
else
[pars(parIndex), count, errmsg, nextindex]＝ ...
sscanf(TheLine(current:endLine), ＇%d＇,1);
if ～(count＝＝1)
fprintf(1, ＇Confused reading pgm header＼n＇);
pars＝－1;
return;
end
parIndex ＝ parIndex+1;
current ＝ current+nextindex－1;
end
end
end
********************************************************************
********************************************************************
function resizeImageFig(h, sz, frac)
% function resizeImageFig(h, sz, frac)
% 功能：根据句柄重获图像的大小
% sz ＝ 图像尺寸
% frac(默认值＝ 1)
if(nargin 3)
frac ＝ 1;
end
pos ＝ get(h, ＇Position＇);
set(h, ＇Units＇, ＇pixels＇, ＇Position＇, ...
[pos(1), pos(2)+pos(4)－frac*sz(1), ...
frac*sz(2), frac*sz(1)]);
set(gca, ＇Position＇, [0 0 1 1], ＇Visible＇, ＇off＇);
********************************************************************
********************************************************************
function hh＝ display＿keypoints(pos, scale, orient, varargin)
% 功能：在原始图像上显示特征点
% 输入：
% pos——特征点的位置矩阵
% scale——特征点的尺度矩阵
% orient——特征点的主方向向量
% 输出：
% hh——返回向量的线句柄
hold on;
alpha ＝ 0.33;
beta ＝ 0.33;
autoscale ＝ 1.5;
plotarrows ＝ 1;
sym ＝ ＇ ＇;
filled ＝ 0;
ls ＝ ＇－ ＇;
ms ＝ ＇ ＇;
col ＝ ＇ ＇;
varin ＝ nargin － 3;
while(varin 0)&isstr(varargin{varin}),
vv ＝ varargin{varin};
if ～isempty(vv)&strcmp(lower(vv(1)), ＇f＇)
filled ＝ 1;
nin ＝ nin－1;
else
[l, c, m, msg]＝ colstyle(vv);
if ～isempty(msg),
error(sprintf(＇Unknown option＂%s＂.＇, vv));
end
if ～isempty(l), ls ＝ l; end
if ～isempty(c), col ＝ c; end
if ～isempty(m), ms ＝ m; plotarrows ＝ 0; end
if isequal(m, ＇.＇), ms ＝ ＇ ＇; end % Don＇t plot ＇.＇
varin ＝ varin－1;
end
end
if varin 0
autoscale ＝ varargin{varin};
end
x ＝ pos(:,1);
y＝ pos(:,2);
u ＝ scale.*cos(orient);
v ＝ scale.*sin(orient);
if prod(size(u))＝＝1, u ＝ u(ones(size(x))); end
if prod(size(v))＝＝1, v ＝ v(ones(size(u))); end
if autoscale,
u ＝ u*autoscale; v ＝ v*autoscale;
end
ax ＝ newplot;
next ＝ lower(get(ax, ＇NextPlot＇));
hold＿state＝ ishold;
x ＝ x(:).＇; y＝ y(:).＇;
u ＝ u(:).＇; v ＝ v(:).＇;
uu ＝[x; x+u; repmat(NaN, size(u))];
vv ＝[y; y+v; repmat(NaN, size(u))];
h1 ＝ plot(uu(:), vv(:), [col ls]);
if plotarrows,
hu ＝[x+u－alpha*(u+beta*(v+eps)); x+u; ...
x+u－alpha*(u－beta*(v+eps)); repmat(NaN, size(u))];
hv ＝[y+v－alpha*(v－beta*(u+eps)); y+v; ...
y+v－alpha*(v+beta*(u+eps)); repmat(NaN, size(v))];
hold on
h2 ＝ plot(hu(:), hv(:), [col ls]);
else
h2 ＝[];
end
if ～isempty(ms),
hu ＝ x; hv ＝ y;
hold on
h3 ＝ plot(hu(:), hv(:), [col ms]);
if filled, set(h3, ＇markerfacecolor＇, get(h1, ＇color＇)); end
else
h3 ＝[];
end
if～hold＿state, hold off, view(2); set(ax, ＇NextPlot＇, next); end
if nargout 0, hh＝[h1; h2; h3]; end
********************************************************************