基于图像surf特征提取和深度信息的图像配准及三维重建matlab仿真

目录

1.图像surf特征提取

2.基于surf的图像配准

3.三维重建

4.MATLAB程序

5.仿真结果

---

1.图像surf特征提取

       SURF是一种高效的局部特征描述符，它结合了尺度空间理论和Haar小波特性，旨在实现对图像中尺度不变和旋转不变的关键点检测以及特征描述。

      尺度空间：为了检测在不同尺度下的图像特征点，SURF采用高斯金字塔构造多尺度空间。给定原图像I(x,y)，其在尺度 σ 下的高斯模糊图像表示为：

其中 G 是二维高斯核函数。

       Difference of Gaussian (DoG)：SURF利用DoG作为尺度空间极值检测器，它是相邻尺度的高斯差分：

其中σ1​ 和 σ2​ 是连续尺度上的高斯卷积核的标准差。

       极值点搜索：在尺度空间内，寻找那些在小邻域内对比度最大（或最小）的点，这些点被认为是潜在的关键点。

       非极大值抑制：对每一个初步检测到的极值点，进行非极大值抑制操作以消除边缘响应，保留真正的局部极值点。

      关键点方向赋值：对每个关键点周围的图像区域计算Haar小波响应的主方向，该方向被用于后续的特征描述符构建和旋转不变性。

       积分图像：首先，通过构建图像的积分图像加速运算，以便快速计算Haar小波特征。

       Haar-wavelet矩形特征：在关键点周围的一个圆形邻域内，选取若干固定大小和方向的小矩形区域（类似于Haar小波基），计算它们的灰度值差异。

      描述符构建：对选定的矩形区域集合，统计正负差异的累加值，形成一个向量，这就是SURF描述符。描述符通常会经过归一化处理以增强光照变化的鲁棒性。

例如，一个SURF描述符D 可以表达为：

其中i 表示描述符的维度索引，n 是矩形区域的数量，wj​ 是权重，Haar(Rj​) 是第j 个矩形区域的Haar小波响应。

在实际应用中，SURF描述符还会进一步二值化和量化，简化匹配过程并提高速度。

2.基于surf的图像配准

       图像配准是指将一幅图像（参考图像）与另一幅图像（浮动图像）进行精确的空间对应和变换，以实现两幅图像的准确叠加或融合。基于SURF特征的图像配准主要分为以下几个步骤：

      首先，分别在参考图像和浮动图像中应用SURF算法进行特征点检测和描述。通过寻找图像尺度空间的极值点（关键点）并生成方向向量，进而创建描述符：

• 关键点定位：在尺度空间中寻找对比度最大的点，并通过非极大值抑制(NMS)确定真正的关键点。
• 方向赋值：在关键点周围邻域内计算Haar小波响应，确定主方向。
• 描述符生成：在关键点周围邻域按照特定方向划分若干小区域，计算Haar-like特征并构成描述符。

       对参考图像和浮动图像中的SURF描述符进行匹配，常用方法有Brute-force匹配、KD树匹配或FLANN等。匹配过程中，通常使用汉明距离或L2距离衡量两个描述符之间的相似性：

      RANSAC（Random Sample Consensus）算法：通过迭代选择一组匹配对拟合一个基本几何模型（如单应性矩阵H），计算内点(inliers)数量，选出最有可能反映真实变换的模型。通过RANSAC或其他优化方法得到足够的内点后，可以更精细地估计最优变换参数，如使用最小二乘法对内点进行拟合。

3.三维重建

当拥有图像序列或多视图图像，并且每张图像都有对应的深度信息时，可以进行三维重建：

1. 深度映射： 将配准后的图像中的每个像素点与其深度值关联起来，形成深度图。

2. 点云构建： 将每张图像中同一物理点在不同视角下的投影点及其深度值组合，生成三维空间中的点云数据。

3. 三维表面重建： 利用点云数据通过体绘制、多边形网格化（如 Marching Cubes 算法）等方式构建三维表面模型。

       综上所述，基于SURF特征提取和深度信息的图像配准及三维重建涉及到了图像处理、特征检测、匹配、变换模型估计、深度信息融合等多个环节，最终目的是从二维图像中恢复出精确的三维场景模型。在实际应用中，这一系列过程往往更加复杂，涉及到更多的优化技术和误差处理机制。

4.MATLAB程序

clc;
clear;
close all;
warning off;
addpath(genpath(pwd));
rng('default')

global dpset;
global Iset;
global params;


load 'params.mat' %加载相机参数
params   = cam_params;
files   = 'test/';
Depth   = dir(strcat(files, '*.mat'));%加载深度信息
Image   = dir(strcat(files, '*.jpg'));%加载图片数据
Iset  = cell(length(Image),1);

dpset = cell(length(Depth), 1);
for i=1:length(Image)
    Iset{i} = strcat(files, Image(i).name);
    dpset{i} = strcat(files, Depth(i).name);
end

 
%三维重建 
Tmat = func_3Dbuilders(dpset, Iset, params);
 

%产生点云
% 读取图像总数
Num = numel(Iset);

% 准备局部变量，用于存储全部XYZ坐标和RGB颜色值
Psoy = [];
Isoy = [];

% 遍历每一张图像
for i = 1:Num

    % 加载RGB图像和深度图像
    im    = imread(Iset{i});                   % 读取RGB图像
    ydp   = load(dpset{i});           % 读取深度图像数据
    d_img = double(ydp.depth_array)/1000; % 将深度数据从毫米转换为米

    % 获取图像尺寸
    [Rs, Cs] = size(d_img);
    Cs2      = repmat((1:Rs)', 1, Cs);                % 创建列坐标矩阵
    Rs2      = ones(Rs, 1) * (1:Cs);                  % 创建行坐标矩阵

    % 获取相机外部参数矩阵
    M        = [params.R params.T];             % 组合旋转矩阵R和平移矩阵T

    % 从深度图像中计算XYZ坐标
    d        = d_img(:)';                              % 将深度图像展平
    d_w      = [d .* Rs2(:)'; d .* Cs2(:)'; d];    % 将深度值与行、列坐标组合
    d_pos    = params.Kdepth \ d_w;            % 使用深度相机内参矩阵Kdepth计算3D坐标

    % 将XYZ坐标转换至第一张图像参考框架
    HT       = [Tmat{i}.R Tmat{i}.T; [0 0 0 1]];
    tmp1     = HT * [d_pos; ones(1, size(d_pos, 2))];
    tmp2     = tmp1(1:3, :);                  % 提取XYZ坐标值

    % 根据3D坐标计算RGB坐标
    tmp3     = M * [d_pos; ones(1, size(d_pos, 2))];

    % 将RGB坐标转换至图像像素坐标
    I0       = params.Krgb * tmp3;
    Ipx      = [I0(1, :) ./ I0(3, :); I0(2, :) ./ I0(3, :)]; % 归一化像素坐标
    Ipx(isnan(Ipx)) = 1;                                         % 处理NaN值
    Ipx      = round(Ipx);                                           % 转换为整数像素坐标

    % 限制像素坐标范围
    row1     = Ipx(1, :); row1(row1 < 1) = 1; row1(row1 > Cs) = Cs;        % 限制行坐标范围
    row2     = Ipx(2, :); row2(row2 < 1) = 1; row2(row2 > Rs) = Rs;        % 限制列坐标范围
    Ipx      = [row1 ; row2];                                              % 更新像素坐标

    % 根据像素坐标获取RGB图像中的颜色值
    IR       = im(sub2ind(size(im), Ipx(2, :), Ipx(1, :), ones(1, size(Ipx, 2))));
    IG       = im(sub2ind(size(im), Ipx(2, :), Ipx(1, :), 2*ones(1, size(Ipx, 2))));
    IB       = im(sub2ind(size(im), Ipx(2, :), Ipx(1, :), 3*ones(1, size(Ipx, 2))));
    rgb      = [IR; IG; IB];                                        % 组合RGB颜色值

    % 将当前图像的XYZ坐标和RGB颜色值追加至累计列表
    Psoy     = [Psoy tmp2];
    Isoy     = [Isoy rgb];
end

% 将所有点的XYZ坐标和RGB颜色值组合为一个点云对象
xyz = pointCloud(Psoy', 'color', Isoy');  % 使用pointCloud函数创建点云对象
showPointCloud(xyz);
4070

5.仿真结果