本文作者分享了在学习和实现SIFT算法过程中的心得，包括使用《SIFT算法原理与OpenCV源码解读》作为主要参考，通过阅读书籍、博客和GitHub代码逐步掌握SIFT特征检测、极值点定位、主方向计算和描述符生成。详尽的MATLAB代码附带注释，助你快速理解和实践SIFT技术。
最近在上数字摄影测量，老师给我们布置了一个大作业——使用MATLAB实现SIFT算法，我对这种编程作业一向比较感兴趣，就想着自己要动手把这个东西做出来。刚开始真的无从下手，看了很多博客，但是对于SIFT算法原理的理解仍旧十分浅显，而且理解原理和代码实现完全不是一个难度，有很多实现细节一般博客是不会介绍的。于是我又开始找书，最后我找到了一本叫《Sift算法原理与OpenCV源码解读》的书，其对OpenCV源码解读比较详细，对我帮助颇大。于是我以这本书为主，再精选了几篇博客，边理解边自己动手写代码，期间也参照了Github上其他人写的代码，最终成功复现了SIFT算法。这个过程真的不容易，资料的良莠不齐，原理的艰难晦涩，真的让我一度想要放弃，不过还好我坚持下来了，收获也是颇多的。
这篇博客我想分享我在学习过程中使用的资料以及实现的代码，代码上附带有详细注释，相信通过这篇博客和我的代码能够帮你快速理解和实现Sift算法。
推荐书本：
《Sift算法原理与OpenCV源码解读》：
推荐博客：
SIFT算法原理详解：
https://www.cnblogs.com/Alliswell-WP/p/SIFT.html
https://www.cnblogs.com/Alliswell-WP/p/SIFT.htmlSIFT特征详解：
https://www.cnblogs.com/wangguchangqing/p/4853263.htmlhttps://www.cnblogs.com/wangguchangqing/p/4853263.html经典图像特征点提取SIFT算法详解：
MATLAB代码：
特征提取代码：
function [pic_orig, fvec, keypoints_dir] = sift(img_file)
% sift算法提取特征
%% 图像预处理
disp("读入图像和图像预处理...");
pic_orig = imread(img_file); % 读入原始图像
% 判断原始是否为彩色图像，是彩色则转换为灰度图像
if(ndims(pic_orig) == 3) % nidms函数用于获取数组维数，彩色图像为3维，灰色图像为2维
pic_gray = rgb2gray(pic_orig);
end
%% 构建高斯金字塔和高斯差分金字塔
disp("正在构建高斯金字塔和高斯差分金字塔...");
tic
[M, N] = size(pic_gray);
sigma = 1.6;
nOctaves = floor(log2(min(M, N))) - 3; % 高斯金字塔的组数(floor函数用来取整，round函数用来四舍五入)
nOctaveLayers = 3; % % 能在每组检测出S个尺度的极值点，nOctaveLayers+3为高斯金字塔每组的图像层数，nOctaveLayers+2为高斯差分金字塔每组的图像层数
k = 2 ^ (1 / nOctaveLayers); % k为尺度空间间隔，即相邻两层尺度之间相差的比例因子
% 使用参考文档中的方法二得到高斯金字塔的第一组第一层图像
pic_gray = imresize(pic_gray, 2, 'bicubic'); % 插值放大图像两倍，'bicubic'为双三次插值，比双线性插值精度更高
sigma_init = sqrt(sigma^2 - (2*0.5)^2);
pic_gray = imgaussfilt(pic_gray, sigma_init); % 对插值后的图像进行高斯滤波
% 将像素值转化为double便于之后处理，这一步必须要做，MATLAB对图像的值使用uint8编码，对其值进行操作到小于0时会让值自动等于0
pic_gray = double(pic_gray);
% 构建高斯金字塔和高斯差分金字塔
gauss_pyr_img = cell(nOctaves, 1);
DoG_pyr_img = cell(nOctaves, 1);
pic_gray = imgaussfilt(pic_gray, sigma);
img_i = pic_gray; % 初始化迭代变量
for oct_i = 1 : nOctaves
gauss_oct_i = zeros([size(img_i) nOctaveLayers+3]); % 预分配空间，此为某一个Octave组中所有高斯模糊图像，是一个三维数组
gauss_oct_i(:, :, 1) = img_i; % 每个Octave组里第一张图片不用进行模糊处理，除第一张外其余是通过对前一组的倒数第三张进行下采样得到
for interval = 2 : nOctaveLayers + 3 % 对一个Octave内的除第一张外其余图像进行高斯滤波处理
sigma_prev = k^(interval - 2) * sigma;
sigma_diff = sqrt((k*sigma_prev)^2 - sigma_prev^2); % 需叠加的高斯模糊参数值
img_i = imgaussfilt(img_i, sigma_diff);
gauss_oct_i(:, :, interval) = img_i;
end
gauss_pyr_img{oct_i} = gauss_oct_i;
DoG_pyr_img{oct_i} = gauss_oct_i(:, :, 1 : end - 1) - gauss_oct_i(:, :, 2 : end); % 求每一个Octave的DoG图像
img_i = imresize(gauss_oct_i(:, :, end - 2), 0.5, 'bicubic'); % 对前一个Octave的倒数第三张图像进行降采样生成下一Octave的第一张高斯模糊图像
end
toc
%% 进行极值点探测
disp('正在检测极值点...')
tic
extrema = []; % 存放极值点容器
num = 0; % 极值总个数
contrastThreshold = 0.04; % 关键点排查阈值
threshold = 0.5 * contrastThreshold * 255 / nOctaveLayers; % 用于筛除像素值的绝对值过低的点
% 寻找极值(包括极大值和极小值)
for oct_i = 1 : size(DoG_pyr_img, 1)
DoG = DoG_pyr_img{oct_i};
[DoG_r, DoG_c, interval] = size(DoG);
for r = 6 : DoG_r - 5
for c = 6 : DoG_c - 5 % 长宽上下5个像素内不用来检测关键点
for interval_i = 2 : interval - 1
DoG_neighborhood = zeros(3, 3, 3);
DoG_neighborhood(:,:,1) = DoG((r-1):(r+1), (c-1):(c+1), interval_i-1);
DoG_neighborhood(:,:,2) = DoG((r-1):(r+1), (c-1):(c+1), interval_i);
DoG_neighborhood(:,:,3) = DoG((r-1):(r+1), (c-1):(c+1), interval_i+1);
val = DoG_neighborhood(2, 2, 2);
if (abs(val) > threshold && (val == max(DoG_neighborhood(:)) || val == min(DoG_neighborhood(:))))
num = num + 1;
sigma_i = k^(interval_i - 1) * sigma; % 计算该极值点尺度sigma
extrema(num, :) = [oct_i, interval_i, r, c, sigma_i, 0]; % 最后一位是预存的主方向角度
end
end
end
end
end
toc
%% 筛选极值点
disp('正在筛选极值点...')
keypoints = [];
tic
for kpt_i = 1 : num
% 获取当前极值点的信息
kpt = extrema(kpt_i, :);
oct = kpt(1);
interval = kpt(2);
kpt_r = kpt(3);
kpt_c = kpt(4);
should_delete = true; % 是否删除标识
% 极值点的精确定位
for try_i = 1 : 6 % 迭代次数超过5次（即当try_i == 6时），则判定为未找到准确特征点
% 取出相关DoG图像，interval可能会变化，所以每次需要重新取DoG图像
DoG = DoG_pyr_img{oct}(:, :, interval) ./ 255;
DoG_prev = DoG_pyr_img{oct}(:, :, interval - 1) ./ 255;
DoG_next = DoG_pyr_img{oct}(:, :, interval + 1) ./ 255;
% 计算一阶导
dD = [DoG(kpt_r, kpt_c + 1) - DoG(kpt_r, kpt_c - 1);
DoG(kpt_r + 1, kpt_c) - DoG(kpt_r - 1, kpt_c);
DoG_prev(kpt_r, kpt_c) - DoG_next(kpt_r, kpt_c)] ./ 2; % 使用有限差分方法求近似偏导，参考文档公式（5.11）
% 计算二阶导(Hessian矩阵)
dxx = (DoG(kpt_r, kpt_c + 1) + DoG(kpt_r, kpt_c - 1) - 2*DoG(kpt_r, kpt_c)) / 1; % 注意分母为1
dyy = (DoG(kpt_r + 1, kpt_c) + DoG(kpt_r - 1, kpt_c) - 2*DoG(kpt_r, kpt_c)) / 1;
dss = (DoG_next(kpt_r, kpt_c) + DoG_prev(kpt_r, kpt_c) - 2*DoG(kpt_r, kpt_c)) / 1;
% 混合导 注意分母为4
dxy = (DoG(kpt_r + 1, kpt_c + 1) - DoG(kpt_r + 1, kpt_c - 1) ...
- DoG(kpt_r - 1, kpt_c + 1) + DoG(kpt_r - 1, kpt_c - 1)) / 4;
dxs = (DoG_next(kpt_r, kpt_c + 1) - DoG_next(kpt_r, kpt_c - 1) ...
- DoG_prev(kpt_r, kpt_c + 1) + DoG_prev(kpt_r, kpt_c - 1)) / 4;
dys = (DoG_next(kpt_r + 1, kpt_c) - DoG_next(kpt_r - 1, kpt_c) ...
- DoG_prev(kpt_r + 1, kpt_c) + DoG_prev(kpt_r - 1, kpt_c)) / 4;
% 由二阶导合成Hessian矩阵
H = [dxx, dxy, dxs;
dxy, dyy, dys;
dxs, dys, dss];
% 求解偏差值
% 如果Hessian矩阵不可逆 去除极值点
if (det(H) == 0)
break;
end
x_hat = -H \ dD;
% 判断偏移程度，如果三个方向的偏移量均小于0.5，说明已经定位到精确点了
if (all(abs(x_hat) < 0.5)) % 阈值为 0.5
should_delete = false;
break;
end
% 调整极值点坐标 以便进行下一次插值计算
kpt_c = kpt_c + round(x_hat(1));
kpt_r = kpt_r + round(x_hat(2));
interval = interval + round(x_hat(3));
% 判断下调整后的坐标是否超过边界（（行边界，列边界，层边界））
if (kpt_r <= 1 || kpt_r >= size(DoG, 1) || ...
kpt_c <= 1 || kpt_c >= size(DoG, 2) || ...
interval < 2 || interval > nOctaveLayers + 1) % interval 不能取第一张与最后一张
break;
end
end
% 迭代次数超过5次，移除该极值点
if(try_i == 6)
continue;
end
% 如果该点已删除 处理下一个极值点
if (should_delete)
continue;
end
% 进一步筛选
% 排除低响应的极值点，低响应点容易受噪声干扰而不稳定
D_hat = DoG(kpt_r, kpt_c) + dD' * x_hat / 2; % 根据文档公式（5.19）计算极值点处的响应值
if (abs(D_hat) * nOctaveLayers < contrastThreshold)
continue;
end
% 消除DoG图像的边缘响应，并非所有的边缘点都能作为好的特征关键点，我们更喜欢角点作为关键点
% 计算Hessian矩阵的迹与行列式
trH = dxx + dyy;
detH = dxx * dyy - dxy * dxy;
edgeThreshold = 10;
if (detH <= 0 || trH^2 * edgeThreshold >= (edgeThreshold + 1)^2 * detH) % 根据文档公式（3.14）可以计算
continue;
end
%筛选完毕，保留筛选后的极值点
sigma_temp = k^(interval - 1) * sigma; % 计算该极值点尺度sigma
kpt_temp = [oct, interval, kpt_r, kpt_c, sigma_temp];
keypoints = [keypoints; kpt_temp];
end
toc
%% 求特征点主方向 kpts
disp("正在提取特征点主方向...");
tic
% 这一部分主要参照博客：https://www.cnblogs.com/Alliswell-WP/p/SIFT.html
num = size(keypoints, 1);
keypoints_dir = []; % 存储具有主方向的特征点
for kpt_i = 1 : num
% 取出特征点信息
kpt = keypoints(kpt_i, :);
oct_i = kpt(1); % 当前特征点所在组
interval = kpt(2); % 当前特征点所在层
kpt_r = kpt(3); % 当前特征点所在行
kpt_c = kpt(4); % 当前特征点所在列
scale = kpt(5); % 当前特征点所在图像的尺度。在将样本点添加到直方图时，需要使用该尺度对其梯度幅值进行高斯加权
radius = round(3 * 1.5 * scale); % 参与计算像素区域半径
gauss_pyr_i = gauss_pyr_img{oct_i}(:, :, interval); % 当前处理的高斯图像
% 遍历参与计算像素区域 (radius*2+1)*(radius*2+1)，生成梯度直方图
img_r = 0; img_c = 0; % 遍历到像素的绝对坐标 (img_r, img_c)
hist_temp = zeros(1, 36 + 4); % 两边各多出的2个空位便于插值运算 1 2 3 ~ 38 39 40
for i = -radius : radius % 从左往右，列数
img_c = kpt_c + i;
if (img_c <= 1 || img_c >= size(gauss_pyr_i, 2)) % 检查坐标的范围
continue;
end
for j = -radius : radius % 从上到下，行数
img_r = kpt_r + j;
if (img_r <= 1 || img_r >= size(gauss_pyr_i, 1))
continue;
end
% 计算梯度
dx = gauss_pyr_i(img_r, img_c + 1) - gauss_pyr_i(img_r, img_c - 1);
dy = gauss_pyr_i(img_r - 1, img_c) - gauss_pyr_i(img_r + 1, img_c);
% 计算梯度 幅值 与 幅角
mag = sqrt(dx^2 + dy^2); % 参照文档（3.17）求幅值
% 由于atan2的特性 计算出的度数为 -pi ~ pi 故需要经过一个转化
if (dx >= 0 && dy >= 0 || dx <= 0 && dy >= 0)
ang = atan2(dy, dx); % 参照文档公式（3.18）求幅角，注意：列数为横坐标x，行数为纵坐标y
else
ang = atan2(dy, dx) + 2*pi;
end
ang = ang * 360/(2 * pi); % 转化为角度制
% 计算落入的直方图柱数
bin = round(36 / 360 * ang); % 根据梯度的角度求出该落入直方图的哪一柱
if( bin >= 36 )
bin = bin - 36;
elseif( bin < 0 )
bin = bin + 36;
end
% 计算高斯加权项
w_G = exp(- (i^2 + j^2) / (2 * (1.5 * scale)^2));
% 存入histtemp 1 2 3 ~ 38 39 40
hist_temp(bin + 3) = hist_temp(bin + 3) + mag * w_G;
end
end
% 为了防止某个梯度方向角度因受到噪声的干扰而突变，需要对梯度方向直方图进行平滑处理
% 填充histtemp的空位 1 2 3 ~ 38 39 40
hist = zeros(1, 36); % 存储平滑处理后的梯度直方图
hist_temp(1 : 2) = hist_temp(37 : 38);
hist_temp(39 : 40) = hist_temp(3 : 4);
for k = 3 : 38 % 加权移动平均 长度为5
hist(k - 2) = (hist_temp(k - 2) + hist_temp(k + 2)) * (1/16) + ...
(hist_temp(k - 1) + hist_temp(k + 1)) * (4/16) + ...
hist_temp(k) * (6/16); % 参照文档公式（5.22）
end
% 遍历直方图求主方向与辅方向
hist_threshold = 0.8 * max(hist); % 计算幅度阈值
for k = 1 : 36
if (k == 1) % kl为第k柱左边一柱的索引
kl = 36;
else
kl = k - 1;
end
if (k == 36)
kr = 1; % kr为第k柱左边一柱的索引
else
kr = k + 1;
end
% 这里不仅要求直方图的柱高超过阈值 还要求比左右相邻的柱都高
if (hist(k) > hist(kl) && hist(k) > hist(kr) && hist(k) >= hist_threshold)
% 通过抛物线插值计算精确幅角公式
bin = k + 0.5 * (hist(kl) - hist(kr)) / (hist(kl) - 2*hist(k) + hist(kr)); % 根据参考文档公式（5.20）
if (bin < 0) % bin越界处理
bin = bin + 36;
elseif (bin >= 36)
bin = bin - 36;
end
% 计算精确幅角
ang = bin * (360/36); % ang为角度制
% 存储特征点的主方向。具有多个方向（主方向和辅方向）的特征点复制成多份，然后将方向值分别赋值给复制后的特征点
kpt_temp = [kpt(1 : 5), ang];
keypoints_dir = [keypoints_dir; kpt_temp];
end
end
end
toc
%% 计算生成特征描述向量
disp("正在计算生成特征描述向量...");
tic
% 此部分主要参考文档中的opencv源码
d = 4; % 关键点附近的区域划分为d*d个子区域，即描述子宽度
n = 8; % 即8个方向
fvec = []; % 即保存描述子的数组
% 遍历关键点，对每个特征点计算128维特征向量
for kpt_i = 1 : size(keypoints_dir, 1)
kpt = keypoints_dir(kpt_i, :); % 待处理的特征点信息
oct_i = kpt(1); % 当前特征点所在组
interval = kpt(2); % 当前特征点所在层
kpt_r = kpt(3); % 当前特征点所在行
kpt_c = kpt(4); % 当前特征点所在列
scale = kpt(5); % 当前特征点所在图像的尺度
ori = kpt(6); % 当前特征点的主方向
gauss_pyr_i = gauss_pyr_img{oct_i}(:, :, interval);
[M, N] = size(gauss_pyr_i); % 当前高斯模糊图像的行列信息
% 坐标轴旋转至主方向
cos_t = cos(ori * pi / 180); % cos接受的是弧度制
sin_t = sin(ori * pi / 180);
hist_width = 3 * scale; % 每个子区域边长
radius = round(hist_width * sqrt(2) * (d + 1) * 0.5); % 实际计算所需要的图像区域半径
radius = min(radius, floor(sqrt(M^2 + N^2))); % 为了避免半径过大，使其小于图像斜对角线长度
% 要把坐标压缩到d*d区域里，参考文档公式（5.23）
cos_t = cos_t / hist_width;
sin_t = sin_t / hist_width;
% 初始化一个存放直方图的容器，容量为(d+2)*(d+2)(n+2)
% 三个维度都多留了2个空位，在之后差值时会超过4*4*8
hist = zeros(d+2, d+2, n+2); % 子区域坐标行数，子区域坐标列数，直方图柱数
for i = -radius : radius % 从上到下 行数
for j = -radius : radius % 从左往右扫描 列数
% 计算旋转以后得坐标位置，原位置为(j, i)
c_rot = j * cos_t - i * sin_t;
r_rot = j * sin_t + i * cos_t;
% 将邻域的原点坐标从中心位置移到该区域左下角（加d/2），然后再减0.5是为了方便三线性插值的计算
rbin = r_rot + d/2 - 0.5;
cbin = c_rot + d/2 - 0.5;
% 计算像素图像坐标绝对(img_r, img_c)
img_r = kpt_r + i;
img_c = kpt_c + j;
% 确定邻域像素是否在d*d邻域正方形内，以及是否超过了图像边界
if (-1 < rbin && rbin < d && -1 < cbin && cbin < d && ...
1 < img_r && img_r < M && 1 < img_c && img_c < N)
% 计算梯度 注意 dy
dx = gauss_pyr_i(img_r, img_c + 1) - gauss_pyr_i(img_r, img_c - 1);
dy = gauss_pyr_i(img_r - 1, img_c) - gauss_pyr_i(img_r + 1, img_c);
% 计算梯度 幅值 与 幅角
mag = sqrt(dx^2 + dy^2);
% 由于atan2的特性 计算出的度数为 -pi ~ pi 故需要经过一个转化
if (dx >= 0 && dy >= 0 || dx <= 0 && dy >= 0)
ang = atan2(dy, dx);
else
ang = atan2(dy, dx) + 2*pi;
end
ang = ang * 360/(2 * pi); % 转化为角度制
% 判断幅角落在哪个直方图的柱内
obin = (ang - ori) * (n / 360); % 取值 0 ~ 7... |8
% 计算高斯加权后的幅值
w_G = exp(- (r_rot^2 + c_rot^2) / (0.5 * d^2)); % 计算高斯加权项 -1/((d/2)^2 * 2) -> -1/(d^2 * 0.5)
mag = mag * w_G;
% 计算三维坐标的整数部分，表示在d*d*d区域中属于哪个正方体
r0 = floor(rbin); % -1 0 1 2 3
c0 = floor(cbin);
o0 = floor(obin);
% 计算三维坐标的小数部分，相当于立方体内点坐标
rbin = rbin - r0;
cbin = cbin - c0;
obin = obin - o0;
% 柱数o0越界处理
if (o0 < 0)
o0 = o0 + n;
elseif (o0 >= n)
o0 = o0 - n;
end
% 根据三线性插值法，计算该像素对正方体的8个顶点的贡献大小
v_rco111 = rbin * cbin * obin;
v_rco110 = rbin * cbin * (1 - obin);
v_rco101 = rbin * (1 - cbin) * obin;
v_rco100 = rbin * (1 - cbin) * (1 - obin);
v_rco011 = (1 - rbin) * cbin * obin;
v_rco010 = (1 - rbin) * cbin * (1 - obin);
v_rco001 = (1 - rbin) * (1 - cbin) * obin;
v_rco000 = (1 - rbin) * (1 - cbin) * (1 - obin);
% 存入直方图中
% rbin 1~6 r0 -1~3
% cbin 1~6 c0 -1~3
% obin 1~8|9 10 插值会到9 o0 0~7
temp = zeros(2, 2, 2);
temp(:, :, 1) = [v_rco000, v_rco010; v_rco100, v_rco110];
temp(:, :, 2) = [v_rco001, v_rco011; v_rco101, v_rco111];
hist(r0+2 : r0+3, c0+2 : c0+3, o0+1 : o0+2) = hist(r0+2 : r0+3, c0+2 : c0+3, o0+1 : o0+2) + mag .* temp;
end
end
end
% 遍历子区域 从hist直方图中导出特征向量
fvec_i = []; % 用于存储kpt_i的特征描述子，一个128维的列向量
for i = 1 : d
for j = 1 : d
% 对幅角小于0或大于360的值进行调整，最终取中间4*4*8的值
hist(i + 1, j + 1, 1) = hist(i + 1, j + 1, 1) + hist(i + 1, j + 1, 9);
fvec_i = [fvec_i; reshape(hist(i + 1, j + 1, 1 : 8), 8, 1)];
end
end
% 为了使归一化后特征矢量中大于0.2的元素设置为0.2，但是为了避免累积误差，此处先进行反归一化处理，为后续归一化处理做准备
fvec_norm = norm(fvec_i, 2);
fvec_threshold = 0.2 * fvec_norm; % 反归一化阈值
fvec_i(fvec_i > fvec_threshold) = fvec_threshold; % 将向量中大于阈值thr的元素替换为thr，为后面归一化做准备
% 归一化特征向量
fvec_norm = norm(fvec_i, 2);
fvec_i = fvec_i ./ fvec_norm;
fvec = [fvec, fvec_i];
end
toc
%% 转换坐标尺度
% 转换坐标尺度
for kpt_i = 1 : size(keypoints_dir, 1)
keypoints_dir(kpt_i, 3 : 4) = keypoints_dir(kpt_i, 3 : 4) .* 2^(keypoints_dir(kpt_i, 1) - 2);
end
%% 显示带主方向关键点
figure
imshow(pic_orig)
hold on
for kpt_i = 1 : size(keypoints_dir, 1)
kpt = keypoints_dir(kpt_i, :);
oct_i = kpt(1);
kpt_r = kpt(3);
kpt_c = kpt(4);
scale = kpt(5);
ori = kpt(6);
radius = scale*2^(oct_i - 1);
color = [rand(), rand(), rand()];
rectangle('position', [kpt_c - radius, kpt_r - radius, radius*2, radius*2], 'curvature', [1, 1], ...
'EdgeColor', color, 'LineWidth', 1);
kpt_c2 = kpt_c + radius * cos(ori * pi / 180);
kpt_r2 = kpt_r - radius * sin(ori * pi / 180);
plot([kpt_c, kpt_c2], [kpt_r, kpt_r2], 'Color', color, 'LineWidth', 1);
% plot(kpt_c, kpt_r, 'r.', 'MarkerSize', 20, 'LineWidth', 1);
end
hold off特征提取结果图：
特征匹配代码：
%%% 对提取出来的特征点进行匹配
%% 预处理
clear
clc
[img1, fvec1, kpts1] = sift('D:\XFF\001.jpg');
[img2, fvec2, kpts2] = sift('D:\XFF\002.jpg');
%% 对两张图片进行垂直拼接
dif = abs(size(img1, 1) - size(img2, 1)); % 计算两张图片高度之间的绝对差值
if (size(img1, 1) > size(img2, 2)) % 给行数较少的图片底部填充白色像素，使两图片高度相同
img2 = [img2; 255*ones(dif, size(img2, 2), size(img2, 3))];
else
img1 = [img1; 255*ones(dif, size(img1, 2), size(img1, 3))];
end
img3 = [img1 img2]; % 拼接图像
%% 特征匹配
march = []; % 用于存储特征点的匹配对
d = []; % 用于存储匹配点对应的距离
threshold = 0.08; % 距离阈值
for kpt_i = 1 : size(fvec1, 2)
for kpt_j = 1 : size(fvec2, 2)
if (norm(fvec1(:, kpt_i) - fvec2(:, kpt_j), 2) < threshold) % 计算欧氏距离
d = [d, norm(fvec1(:, kpt_i) - fvec2(:, kpt_j))];
march = [march; kpt_i, kpt_j];
end
end
end
%% 绘图
figure
imshow(img3)
hold on
% 绘制左图特征点
for kpt_i = 1 : size(kpts1, 1)
kpt = kpts1(kpt_i, :);
plot(kpt(4), kpt(3), 'r.', 'MarkerSize', 10);
end
% 绘制右图特征点
for kpt_i = 1 : size(kpts2, 1)
kpt = kpts2(kpt_i, :);
plot(kpt(4) + size(img1, 2), kpt(3), 'r.', 'MarkerSize', 10);
end
% 绘制特征点连线
for i = 1 : size(march, 1)
kpt1 = kpts1(march(i, 1), :);
kpt2 = kpts2(march(i, 2), :);
line([kpt1(4), kpt2(4) + size(img1, 2)], [kpt1(3), kpt2(3)], 'LineWidth', 1, 'Color', [rand(),rand(),rand()])
end
hold off特征匹配结果：
最后说两句
通过这篇博客推荐的资料完全可以帮助你实现SIFT算法，这个过程很花时间和精力，但是成功没有捷径，唯有静下心来慢慢学，加油吧！
扫一扫
14万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
