网站建设平台排行榜,国外免费空间哪个好,企业网站建设方案书目录,世界营销大师排名第一次来请先看这篇文章#xff1a;【图像拼接#xff08;Image Stitching#xff09;】关于【图像拼接论文源码精读】专栏的相关说明#xff0c;包含专栏内文章结构说明、源码阅读顺序、培养代码能力、如何创新等#xff08;不定期更新#xff09;
【图像拼接论文源码精…第一次来请先看这篇文章【图像拼接Image Stitching】关于【图像拼接论文源码精读】专栏的相关说明包含专栏内文章结构说明、源码阅读顺序、培养代码能力、如何创新等不定期更新
【图像拼接论文源码精读】专栏文章目录
【源码精读】As-Projective-As-Possible Image Stitching with Moving DLTAPAP第一部分全局单应Global homography 文章目录 【图像拼接论文源码精读】专栏文章目录前言1. 跑通代码得到结果1.1 准备工作1.2 尝试运行1.3 得到拼接结果 2. 源码解读看懂原理2.1 参数与输入图像2.2 图像对齐2.2.1 总体流程2.2.2 siftMatch函数2.2.3 homoRANSAC函数2.2.4 calcHomo函数2.2.5 homographyAlign函数2.2.6 小结 2.3 接缝融合blendTexture函数2.3.1 找重叠区域边界并得到基于感知的接缝切割算法的【数据项】2.3.2 可视化重叠区域边界新增代码2.3.3 基于感知的接缝切割算法的【平滑项】2.3.4 图割梯度融合得到拼接结果2.3.5 可视化拼接结果带接缝2.3.6 沿着接缝找没对齐的地方进行修补论文算法核心代码对应论文中3.2和3.3部分2.3.7 论文算法核心代码对应论文中3.4和3.5部分 3. 总结思考试图创新 前言
论文题目Seam-guided local alignment and stitching for large parallax images——大视差图像的接缝引导局部对齐和拼接
论文地址Seam-guided local alignment and stitching for large parallax images
论文源码https://github.com/tlliao/Seam-guided-local-alignment
注matlab源码相关算法是C封装的。主要是对matlab的学习与理解。
论文精读【图像拼接】论文精读Seam-guided local alignment and stitching for large parallax images
配合【论文精读】专栏对应文章阅读效果更佳
注请重点关注代码段中的注释有一些讲解的东西直接写在代码段的注释中了同时多关注红色字体和绿色字体
1. 跑通代码得到结果
任何源码下载下来后第一件事就是先跑通。
无论你是否要在该工作的基础上创新总是需要得到拼接结果在论文的实验部分做对比。
所以请务必跑通得到拼接结果。
1.1 准备工作
本地需要matlab环境我是MATLAB R2018b选择一个适中的matlab版本即可。最新的matlab2023可以使用局部函数了类似jupyter notebook感兴趣的同学可以试试最新版。
源码下载下来后用matlab打开项目界面如下 左侧是文件目录结构中级是当前所选文件代码下面是命令行窗口右侧是工作区运行后显示相关变量的值
1.2 尝试运行
点击上面菜单栏【运行】出现如下报错 检查后发现当前目录下没有vlfeat-0.9.21文件夹 我们在其他论文的源码中将vlfeat-0.9.21文件夹整个复制到该目录下即可 并将报错位置的代码修改为
% run ../vlfeat-0.9.21/toolbox/vl_setup;
run vlfeat-0.9.21/toolbox/vl_setup;再次运行后 代码已经可以完整的跑通了。
1.3 得到拼接结果
我们发现源码并没有显示拼接结果。在main.m中末尾添加显示结果
f figure;
imshow(seam_cut,border,tight);
hold on;再次运行后得到1_l.jpg和1_r.jpg的拼接结果 原图为: 1_l.jpg 1_r.jpg
附上原作者README中的使用方法
2. 源码解读看懂原理
本节将按照main.m中的代码顺序进行【模块化】讲解包括与论文中的算法对应、matlab语法和函数学习、变量的类型和值、函数功能等方面。代码段中包含原作者的注释和我做的注释与讲解结合着阅读。
2.1 参数与输入图像
代码如下
%% Setup VLFeat toolbox.
%----------------------
addNeedingPaths;
% 刚下载下来的源码没有vlfeat-0.9.21文件夹
% run ../vlfeat-0.9.21/toolbox/vl_setup;
run vlfeat-0.9.21/toolbox/vl_setup;% setup parameters
% Parameters of SIFT detection
parameters.peakthresh 0;
parameters.edgethresh 500;% % Parameters of RANSAC via fundamental matrix
parameters.minPtNum 4; % minimal number for model fitting
parameters.iterNum 2000; % maximum number of trials
parameters.thDist 0.01; % distance threshold for inliersimgpath Imgs/;
img_format 4_*.jpg; % 以论文中图1为例
dir_folder dir(strcat(imgpath, img_format)); %连接路径列出文件夹中的内容path1 sprintf(%s%s,imgpath, dir_folder(1).name); % 路径转换成字符串
path2 sprintf(%s%s,imgpath, dir_folder(2).name); %
img1 im2double(imread(path1)); % target image 图像转成双精度后面SIFT的输入要求是双精度
img2 im2double(imread(path2)); % reference image主要包含SIFT算法参数、RANSAC算法参数、输入两张图像的路径、得到双精度的目标图和参照图。
2.2 图像对齐
代码如下
%% image alignment
fprintf( image alignment...);
tic; % 计时器
[pts1, pts2] siftMatch(img1, img2, parameters);
[matches_1, matches_2] homoRANSAC(pts1, pts2, parameters);
init_HcalcHomo(matches_1, matches_2);
[warped_img1, warped_img2] homographyAlign(img1, img2, init_H);
fprintf(done (%fs)\n, toc); % 输出时间2.2.1 总体流程
对两张输入图像应用SIFT算法得到两张图像的特征点pts1pts2通过RANSAC算法得到两张图像的匹配点即内点inliersmatches_1matcher_2通过匹配点得到初始的单应矩阵init_H根据单应矩阵init_H得到翘曲后的两张图像warped_img1,warped_img2
接下来我们具体解读上述步骤中的函数。
2.2.2 siftMatch函数
函数功能得到两张图的特征点输入两张图像和算法参数返回值为两张图像的特征点。 对应根目录下文件siftMatch.m代码如下
function [pts1, pts2] siftMatch( img1, img2, parameters )
%--------------------------------------
% SIFT keypoint detection and matching.
%--------------------------------------
peakthresh parameters.peakthresh;
edgethresh parameters.edgethresh;
%fprintf( Keypoint detection and matching...);tic;
[ kp1,ds1 ] vl_sift(single(rgb2gray(img1)),PeakThresh, peakthresh,edgethresh, edgethresh);
[ kp2,ds2 ] vl_sift(single(rgb2gray(img2)),PeakThresh, peakthresh,edgethresh, edgethresh);
matches vl_ubcmatch(ds1, ds2);
%fprintf(done (%fs)\n,toc);% extract match points position
pts1 kp1(1:2,matches(1,:));
pts2 kp2(1:2,matches(2,:)); end前文【源码精读】As-Projective-As-Possible Image Stitching with Moving DLTAPAP第一部分全局单应Global homography已经讲过包括每个函数、每个变量的含义非常详细。本节不再赘述。
2.2.3 homoRANSAC函数
函数功能使用MATLAB自带的RANSAC函数去除特征匹配点中的异常值得到内点索引进而得到所有的正常匹配点。
对应根目录下文件homoRANSAC.m代码如下
function [matches_1, matches_2] homoRANSAC(pts1, pts2, parameters)
% using fundamental matrix for robust fittingminPtNum parameters.minPtNum; % minimal number of points to estimate H and e
iterNum parameters.iterNum; % maximum iterations
thDist parameters.thDist; % distance threshold
% ptNum size(pts1, 2); % number of points%% perform coordinate normalization
[normalized_pts1, ~] normalise2dpts([pts1; ones(1,size(pts1, 2))]);
[normalized_pts2, ~] normalise2dpts([pts2; ones(1,size(pts2, 2))]);
points [normalized_pts1, normalized_pts2]; % n×6 [x1 y1 1 x2 y2 1]% (points)这部分定义了一个匿名函数在 MATLAB 中使用 符号表示创建一个函数。
% 括号中的 points 是函数的输入参数表示这个函数接受一个名为 points 的参数。
% 这么写是为了适配下面的ransac函数的输入必须是function_handle
fitmodelFcn (points)calcNormHomo(points); % fit function拟合函数 9×1
evalmodelFcn (homo, points)calcDistofHomo(homo, points); % 距离 1×nrng(0); % 设置随机种子让每次重启程序的RANSAC一致
% matlab自带的ransac算法其实我觉得用APAP那一套RANSAC更方便一些
[~, inlierIdx] ransac(points,fitmodelFcn,evalmodelFcn,minPtNum,thDist,MaxNumTrials,iterNum);inliers1 pts1(:, inlierIdx);
inliers2 pts2(:, inlierIdx);matches_1 inliers1;
matches_2 inliers2;% delete duplicate feature match
[~, ind1] unique(matches_1, rows);
[~, ind2] unique(matches_2, rows);
ind intersect(ind1, ind2);
matches_1 matches_1(:, ind);
matches_2 matches_2(:, ind);endfunction [ homo ] calcNormHomo(points) % estimate H_inf and e via DLT
% 求解单应矩阵的过程先随便计算一个矩阵H得到模型。RANSAC第一步npts1 points(:, 1:3);
npts2 points(:, 4:6);%% calculation the initial H0 and e0
Equation_matrix zeros(2*size(npts1, 2), 9);
for i1:size(npts1, 2)xi npts1(1,i); yi npts1(2,i);xi_ npts2(1,i); yi_ npts2(2,i);tmp_coeff1 [xi, yi, 1, 0, 0, 0, -xi*xi_, -yi*xi_, -xi_];tmp_coeff2 [0, 0, 0, xi, yi, 1, -xi*yi_, -yi*yi_, -yi_];Equation_matrix(2*i-1:2*i, :) [tmp_coeff1; tmp_coeff2];
end[~,~,v] svd(Equation_matrix, 0);
norm_homo reshape(v(1:9, end), 3, 3);
homo norm_homo(:); % 按列排序排成1列endfunction dist calcDistofHomo(homo, points) % calculate the projective error
% 计算投影误差看手写部分RANSAC第二步
% RANSAC最后一步就是用自带的ransac函数迭代找比这个误差小的作为内点更新H返回模型和索引pts1 points(:, 1:3);
pts2 points(:, 4:6);H reshape(homo(1:9),3,3); % 和上面norm_homo一样tmp1 (H(1,:)*pts1)./pts1(3,:);
tmp2 (H(2,:)*pts1)./pts1(3,:);
tmp3 (H(3,:)*pts1)./pts1(3,:);
mapped_pts2(1,:) tmp1./tmp3;
mapped_pts2(2,:) tmp2./tmp3;
dist sqrt(sum((mapped_pts2-pts2(1:2,:)).^2, 1));end
RANSAC过程以及代码对应
随机抽出mm4个样本数据计算出单应性矩阵H’记为模型M使用剩余的点对计算投影误差检测误差小于设置阈值的点的个数如果个数多于之前记录的最优值则替换H’并记录此次的点数为最优值是否达到设置的迭代次数若是则输出H’并结束否则重新迭代回到步骤1
步骤1对应calcNormHomo函数步骤2对应calcDistofHomo函数步骤3对应homoRANSAC中的ransac函数。
其中calcNormHomo函数是计算单应矩阵的过程详解见【源码精读】As-Projective-As-Possible Image Stitching with Moving DLTAPAP第一部分全局单应Global homography。有手写的公式推导。
重点讲一下投影误差计算calcDistofHomo函数。
计算出H’对应函数中的H后使用如下公式进行投影误差计算 E p r o j ∑ i 0 n ( ( x i ′ − h 11 x i h 12 y i h 13 h 31 x i h 32 y i h 33 ) 2 ( y i ′ − h 21 x i h 22 y i h 23 h 31 x i h 32 y i h 33 ) 2 ) E_{proj} \sqrt{\sum_{i 0}^n((x_i - \frac{h_{11}x_ih_{12}y_ih_{13}}{h_{31}x_ih_{32}y_ih_{33}})^2(y_i - \frac{h_{21}x_ih_{22}y_ih_{23}}{h_{31}x_ih_{32}y_ih_{33}})^2)} Eproji0∑n((xi′−h31xih32yih33h11xih12yih13)2(yi′−h31xih32yih33h21xih22yih23)2)
代码中的tmp1为 h 11 x i h 12 y i h 13 h_{11}x_ih_{12}y_ih_{13} h11xih12yih13tmp2为 h 21 x i h 22 y i h 23 h_{21}x_ih_{22}y_ih_{23} h21xih22yih23tpm3为 h 31 x i h 32 y i h 33 h_{31}x_ih_{32}y_ih_{33} h31xih32yih33。过程就是计算上面的公式得到投影误差。
最后经过matlab自带的ransac函数得到正常匹配点索引进而得到匹配点。
2.2.4 calcHomo函数
函数功能计算DLT直接线性变换得到全局单应矩阵。 对应根目录下文件calcHomo.m代码如下
function H calcHomo(pts1,pts2)
%% use Direct linear tranformation (DLT) to calculate homography
% approxmation: H*[pts1; ones(1,size(pts1,2))] [pts2; ones(1,size(pts2,2))]
% Normalise point distribution.
data_pts [ pts1; ones(1,size(pts1,2)) ; pts2; ones(1,size(pts2,2)) ];
[ dat_norm_pts1,T1 ] normalise2dpts(data_pts(1:3,:));
[ dat_norm_pts2,T2 ] normalise2dpts(data_pts(4:6,:));
data_norm [ dat_norm_pts1 ; dat_norm_pts2 ]; % 6×n同APAP中全局单应部分%-----------------------
% Global homography (H).
%-----------------------
%fprintf(DLT (projective transform) on inliers\n);
% Refine homography using DLT on inliers.
%fprintf( Refining homography (H) using DLT...);
[ h,~,~,~ ] homography_fit(data_norm); % 计算DLT得到9×1
H T2\(reshape(h,3,3)*T1); % 得到初始的全局单应矩阵end其中homography_fit函数为
function [P A C1 C2] homography_fit(X,A,W,C1,C2)
% 得到9×1的单应矩阵P
x1 X(1:3,:); % 两个图的点3×n[x y 1]
x2 X(4:6,:);if nargin 5H vgg_H_from_x_lin(x1,x2,A,W,C1,C2);
else[H A C1 C2] vgg_H_from_x_lin(x1,x2); % 只传一个参数会进入这里
endP H(:); % 9 ×1end再进一步vgg_H_from_x_lin函数为
function [H, A, C1, C2] x(xs1,xs2,A,W,C1,C2)% 函数功能计算DLT% 返回值含义% H3×3单应矩阵h33为1% AAh 0 的A最小二乘计算先新变化% C1、C2归一化矩阵% H vgg_H_from_x_lin(xs1,xs2)%% Compute H using linear method (see Hartley Zisserman Alg 3.2 page 92 in% 1st edition, Alg 4.2 page 109 in 2nd edition). % Point preconditioning is inside the function.%% The format of the xs [p1 p2 p3 ... pn], where each p is a 2 or 3% element column vector.[r,c] size(xs1); % r是行数此处使用为3c是点的个数if (size(xs1) ~ size(xs2)) % 如果两个匹配点数量不同报错error (Input point sets are different sizes!)endif (size(xs1,1) 2) % 如果是非齐次的则齐次化xs1 [xs1 ; ones(1,size(xs1,2))];xs2 [xs2 ; ones(1,size(xs2,2))];end% condition points只传了两组特征点坐标进入这里if nargin 2 C1 vgg_conditioner_from_pts(xs1); % 两个点集的调节矩阵C2 vgg_conditioner_from_pts(xs2);xs1 vgg_condition_2d(xs1,C1); % 得到归一化后的坐标xs2 vgg_condition_2d(xs2,C2);endif nargin 6B A;B(1:2:end,:)W*A(1:2:end,:);B(2:2:end,:)W*A(2:2:end,:);% Extract nullspace[u,s,v] svd(B, 0); s diag(s);elseA [];ooo zeros(1,3);for k1:c % 取每个点p1 xs1(:,k); % 得到齐次坐标p2 xs2(:,k);A [ A;p1*p2(3) ooo -p1*p2(1)ooo p1*p2(3) -p1*p2(2)]; % 得到2n × 9Ah0 的A: [ x y 1 0 0 0 -xx -yx -x]% [ 0 0 0 x y 1 -xy -yy -y]% 反正是解方程等式右边是0矩阵-A和A一样end% Extract nullspace[u,s,v] svd(A, 0); %svd分解s diag(s);endnullspace_dimension sum(s eps * s(1) * 1e3);if nullspace_dimension 1fprintf(Nullspace is a bit roomy...);endh v(:,9); %要v的最右列向量9×1H reshape(h,3,3); % 变成 3×3%deconditionH C2\H*C1; % 调节回去去归一化H H/H(3,3); % h33变成1
endvgg_condition_2d函数
function pc z(p,C)
% 函数功能使用调节矩阵调节点集得到归一化后的坐标点
% function pc vgg_condition_2d(p,C);
%
% Condition a set of 2D homogeneous or nonhomogeneous points using conditioner C[r,c] size(p);
if r 2pc vgg_get_nonhomg(C * vgg_get_homg(p));
elseif r 3pc C * p; % 3维点每个点左乘调节矩阵
elseerror (rows ! 2 or 3);
endend上面函数跟着阅读注释即可。重点理解vgg_H_from_x_lin函数在文件夹modespecific中都是一些通用方法。
2.2.5 homographyAlign函数
函数功能得到参照图和目标图翘曲后的图像 对应根目录下文件homographyAlign.m代码如下
function [ homo1, homo2] homographyAlign( img1,img2,init_H )
%input: target image and reference image, saliency map of the two images
%output: homography-warped target and reference, with their corresponding
%saliency maps
tform projective2d(init_H); % 根据全局单应得到2d投影对象与imwarp同时使用
img1mask imwarp(true(size(img1,1),size(img1,2)), tform, nearest); % 得到img1的掩码img1To2 imwarp(img1, tform); % 得到变换后的图像
img1To2 cat(3,img1mask, img1mask, img1mask).*img1To2; % 用掩码再处理一下
%注可视化上面两个目标图几乎没有区别用掩码在三个通道分别处理一下可能是为了修正图像pt [1, 1, size(img1,2), size(img1,2);1, size(img1,1), 1, size(img1,1);1, 1, 1, 1]; % 图像四个角点列向量
H_pt init_H*pt; % 得到变换后的点
H_pt H_pt(1:2,:)./repmat(H_pt(3,:),2,1); %归一化齐次转非齐次% calculate the convas
% 计算画布大小适配两张图像
off ceil([ 1 - min([1 H_pt(1,:)]) 1 ; 1 - min([1 H_pt(2,:)]) 1 ]);
cw max(size(img1To2,2)max(1,floor(min(H_pt(1,:))))-1, size(img2,2)off(1)-1);
ch max(size(img1To2,1)max(1,floor(min(H_pt(2,:))))-1, size(img2,1)off(2)-1);% 目标图翘曲后的画布肯定能容纳下参照图所以就用目标图的画布大小
homo1 zeros(ch,cw,3); homo2 zeros(ch,cw,3); % % 填充两张图像到画布中
homo1(floor(min(H_pt(2,:)))off(2)-1:floor(min(H_pt(2,:)))off(2)-2size(img1To2,1),...floor(min(H_pt(1,:)))off(1)-1:floor(min(H_pt(1,:)))off(1)-2size(img1To2,2),:) img1To2;
homo2(off(2):(off(2)size(img2,1)-1),off(1):(off(1)size(img2,2)-1),:) img2;end注计算画布那里原作者的注释写错了是canvas不是convas。
在主函数中添加可视化代码
figure;
imshow(warped_img1);
title(warpedimg1);
figure;
imshow(warped_img2);
title(warpedimg2);% 平均融合除2非重叠区域暗重叠区域原亮度不除2则是非重叠区域原亮度重叠区域更亮
% 原论文图1(b)的是TFA得到的而不是此处的单应变换
output_canvas(:,:,1) (warped_img1(:,:,1)warped_img2(:,:,1))/2;
output_canvas(:,:,2) (warped_img1(:,:,2)warped_img2(:,:,2))/2;
output_canvas(:,:,3) (warped_img1(:,:,3)warped_img2(:,:,3))/2;figure;
imshow(output_canvas);
title(output_canvas);翘曲后的两张图像和拼接结果如下使用论文图1中的数据
2.2.6 小结
本节详细讲解了两张图像对齐的流程首先通过SIFT得到特征点然后经过RANSAC筛选得到正常匹配点最后将参照图和目标图翘曲。添加了可视化代码得到拼接结果实现了论文中图1的重叠区域保持原亮度非重叠区域暗一些的效果。 读代码的方法 1. 紧跟变量流包括特征点流、图像流等 2. 把函数单拿出来测试matlab的工作区可以显示变量的结构和值 比如SIFT和RANSAC的部分就是特征点流一系列的处理都是围绕特征匹配点展开的重点关注特征点的向量如何变化包括特征点的结构、取值、行列变换、齐次化等。
前面都是一些铺垫和准备工作是其他论文中已经实现的预对齐步骤。下一节是本论文的核心部分我将重点讲解并与论文中的部分对应。
2.3 接缝融合blendTexture函数
2.3.1和2.3.3对应论文的3.1部分Conventional Seam-cutting和算法1中的步骤1是Perception-based-seam-cutting论文中基于感知的接缝算法得到的初始接缝通过数据项和平滑项经过图割和梯度融合得到拼接结果。
2.3.6对应论文的3.2 SSIM-based Seam Evaluation和3.3 Components Extraction部分算法1中步骤2345。SSIM误差是论文Quality evaluation-based iterative seam estimation for image stitching中提出的根据SSIM误差和Ostu算法得到的阈值比较得到未对齐的接缝像素点和块。
2.3.7对应论文的3.4 Patch Alignment和3.5 Seam Merging部分算法1中步骤6789101112。将未对齐的块分离开使用SIFTflow算法获得流向量用公式4中的sigmoid函数修正流向量进而对齐块计算新的接缝用新的接缝块替换原图中的接缝块。
代码如下
%% image composition
fprintf( seam cutting...);tic;
[seam_cut] blendTexture(warped_img1, warped_img2);
fprintf(done (%fs)\n, toc);本节我们将一段一段讲解blendTexture.m文件下对应的代码也是本论文的核心。
2.3.1 找重叠区域边界并得到基于感知的接缝切割算法的【数据项】
seam-cutting预处理
%% pre-process of seam-cutting
w1 imfill(imbinarize(rgb2gray(warped_img1), 0),holes); % 翘曲图灰度化、二值化、黑色变白色
w2 imfill(imbinarize(rgb2gray(warped_img2), 0),holes);
A w1; B w2;
C A B; % mask of overlapping region重叠区域像素值为1其余为0
[ sz1, sz2 ]size(C);
ind find(C); % index of overlapping region重叠区域索引
nNodes size(ind,1); % 重叠区域像素个数
revindC zeros(sz1*sz2,1); % 反向索引初始全是0整张图
revindC(C) 1:length(ind); % 重叠区域中的索引顺序得到重叠区域边界以及数据项作为图割算法输入参数
%% terminalWeights, choose source and sink nodes
% 这段代码是在找重叠区域边界border_B findBorder(B); %参照图边界
border_C findBorder(C); %重叠区域边界imgseedA border_B border_C; % 参照图与重叠区域相交的两条边
imgseedB ~imgseedA border_C; % 目标图翘曲后与重叠区域相交的两条边% data term
twzeros(nNodes,2); % 初始化数据项权重
tw(revindC(imgseedA),1)inf; % 第一列对应于标签1的权重设置为正无穷
tw(revindC(imgseedB),2)inf; % 第二列对应于标签2的权重设置为正无穷terminalWeightstw; % data term得到图割算法中的数据项% 上面代码的目的是为了图割算法中通过最小化数据项权重来实现图割算法将边界设置为无穷大方便后续计算。其中找边界用的是错位差分法。对应的函数文件为findBorder.m
function [ border_img ] findBorder( mask_img )
% give a mask image, find its border, (boundary points)
% border of mask image
[sz1, sz2] size(mask_img);
mask_R(mask_img-[mask_img(:,2:end) false(sz1,1)])0; %右边界
mask_L(mask_img-[false(sz1,1) mask_img(:,1:end-1)])0; %左边界
mask_D(mask_img-[mask_img(2:end,:);false(1,sz2)])0; % 下边界
mask_U(mask_img-[false(1,sz2);mask_img(1:end-1,:)])0; %上边界
border_img mask_R | mask_L | mask_D | mask_U;% 差分操作列方向往右平移一个像素然后通过差分判断是否大于0从而得到边界
% 最后合并边界信息边界像素值为1其他地方为0
end2.3.2 可视化重叠区域边界新增代码
我们添加可视化代码显示重叠区域边界
%% 可视化重叠区域边界[m,~] find(C); % 找最大最小行
minrow min(m);
maxrow max(m);mask~(imgseedB | imgseedA); %掩码 重叠区域边界为0
% 将两个翘曲图像对齐并添加权重跟前面显示图像拼接结果一样0.5是平均融合亮度与原图一致
imglap 0.5.*imadd(warped_img1.*cat(3,maskC,maskC,maskC), warped_img2.*cat(3,maskC,maskC,maskC));% free_seed中没有最小最大行索引对应的行
freeseed imgseedA;
freeseed(minrow,:) 0;
freeseed(maxrow,:) 0;
free_seed imgseedA (~freeseed);% 种子区域权重设置为0
tw(revindC(free_seed),1) 0;
tw(revindC(free_seed),2) 0;
imgseed warped_img1.*cat(3,A-C,A-C,A-C) warped_img2.*cat(3,B-C,B-C,B-C) imglapcat(3,freeseed,free_seed,imgseedB);
% imgseed warped_img1.*cat(3,A-C,A-C,A-C) warped_img2.*cat(3,B-C,B-C,B-C) imglap% cat(3,A-C,A-C,A-C)目标图重叠区域边界为红色R通道高亮
% cat(3, B - C, B - C, B - C)参照图重叠区域边界标记为蓝色B通道高亮
% cat(3,freeseed,free_seed,imgseedB): 会影响部分边界颜色变为绿色G通道高亮不加这一项边界都是黑色
% 三色组合显示重叠区域边界% 可视化
figure,imshow(imgseed);
title(seed image on two warped images);% 如果要保存的话
% pngout sprintf(Overlapping_result.png);
% imwrite(imgseed,pngout);可视化结果 注这里代码与Perception-based-seam-cutting中的类似将找边界封装成函数了。
2.3.3 基于感知的接缝切割算法的【平滑项】
计算平滑项
%% calculate edgeWeights
% 四个边界为0的掩码
CL1C[C(:,2:end) false(sz1,1)]; % 与C比最右一列为0
CL2[false(sz1,1) CL1(:,1:end-1)]; % 最左一列为0
CU1C[C(2:end,:);false(1,sz2)]; % 最下一行
CU2[false(1,sz2);CU1(1:end-1,:)]; % 最上一行% edgeWeights: sigmoid-metric difference mapsigmoid差异图
[imgdif_sig, ~] calcSigmoidDiff(warped_img1, warped_img2, C);% sigmoid method用于衡量图像中相邻区域的差异即图像变化的平滑度
DL (imgdif_sig(CL1)imgdif_sig(CL2))./2; % 水平方向上平均值
DU (imgdif_sig(CU1)imgdif_sig(CU2))./2; % 垂直方向上平均值% smoothness term
edgeWeights[revindC(CL1) revindC(CL2) DL1e-8 DL1e-8;revindC(CU1) revindC(CU2) DU1e-8 DU1e-8];% 每一行一个边界包含各边界索引和对应的平滑项值小偏移量1e-8防止分母为0通过sigmoid函数计算差异calcSigmoidDiff.m
unction [ imgdif_sig, imgdif ] calcSigmoidDiff(img1, img2, C)
% 函数功能计算两张图像差异返回一个sigmoid函数差异图% edgeWeights: Euclidean-weighted norm
% 两张图像三通道分量
ang_1img1(:,:,1); sat_1img1(:,:,2); val_1img1(:,:,3);
ang_2img2(:,:,1); sat_2img2(:,:,2); val_2img2(:,:,3);% baseline difference map两个图像每个像素位置上RGB通道的差异
imgdif sqrt( ( (ang_1.*C-ang_2.*C).^2 (sat_1.*C-sat_2.*C).^2 (val_1.*C-val_2.*C).^2 )./3 ); % sigmoid-metric difference map
a_rgb 0.06; % bin of histogram
beta4/a_rgb; % beta
gammaexp(1); % base number
para_alpha histOstu(imgdif(C), a_rgb); % parameter:tau
imgdif_sig 1./(1power(gamma,beta*(-imgdifpara_alpha))); % difference map with logistic function
imgdif_sig imgdif_sig.*C; % difference to compute the smoothness term % 使用 Sigmoid 函数将 imgdif 映射到 (0, 1) 范围内得到 imgdif_sig。
% 然后通过与掩码 C 相乘仅保留 C 中为真的区域的值这样可以计算平滑项smoothness term。
% 对应论文中计算平滑项的公式end上面数据项与平滑项的理论来自于文章Perception-based energy functions in seam-cutting。对应的论文精读为
2.3.4 图割梯度融合得到拼接结果
图割得到接缝的标签矩阵值为1是接缝的像素然后梯度融合得到拼接结果
%% graph-cut labeling
[~, labels] graphCutMex(terminalWeights, edgeWeights); AsA;
BsB;
As(ind(labels1))false; % mask of target seam
Bs(ind(labels0))false; % mask of reference seam
imgout gradient_blend(warped_img1, As, warped_img2); SE_seam strel(diamond, 1);
As_seam imdilate(As, SE_seam) A;
Cs_seam As_seam Bs; % 重叠区域最终的接缝像素值为1其余为02.3.5 可视化拼接结果带接缝
我们添加可视化代码显示拼接缝
%% 可视化接缝
imgseam imgout.*cat(3,(A|B)-Cs_seam,(A|B)-Cs_seam,(A|B)-Cs_seam) cat(3,Cs_seam,zeros(sz1,sz2),zeros(sz1,sz2));figure,imshow(imgseam);
title(final stitching seam); SE_seam中的strel第二个参数条件接缝粗细。
2.3.6 沿着接缝找没对齐的地方进行修补论文算法核心代码对应论文中3.2和3.3部分
找到接缝上所有像素点的坐标给接缝上所有像素点定义一个20×20的块每个块计算SSIM误差如果所得的SSIM误差最大值都在范围内则是合理的接缝否则根据Ostu算法得到阈值tou找到接缝上比该阈值大的块和像素点
该模块代码如下
%% find potential artifacts along the seam for patch mark
% 沿着拼接缝找没对齐的地方然后修补% extract pixels on the seam and evaluate the patch error
% 提取接缝处的像素并评估补丁误差seam_pts contourTracingofSeam(Cs_seam); % 接缝上所有点的坐标[行列]是整张图的索引
[ssim_error, ~, patch_coor] evalSSIMofSeam(warped_img1, warped_img2, C, seam_pts, patch_size); %每个像素点的SSIM误差和块范围% 对应论文3.3部分
% mark misaligned local regions
% max(Q) kmean(Q)
if max(ssim_error)1.5*mean(ssim_error) % 接缝所有的像素点都在误差允许范围内则接缝合理直接就是结果seam_cutimgout;return;
endT graythresh(ssim_error); % Ostu的阈值tou, 以图4为例T0.1882
artifacts_pixels seam_pts(ssim_errorT,:); % 找到所有比阈值大的像素坐标
artifacts_patchs patch_coor(ssim_errorT,:); % 比阈值大的块的范围% 错位像素、和错位像素的块的掩码
artifacts_masks false(sz1,sz2);
mask_pixels false(sz1,sz2);
for i1:size(artifacts_patchs,1)artifacts_masks(artifacts_patchs(i,1):artifacts_patchs(i,2),artifacts_patchs(i,3):artifacts_patchs(i,4))1;mask_pixels(artifacts_pixels(i,1),artifacts_pixels(i,2))1;
end% add modification to artifacts_masks: connect neighboring patches if they
% are close enough
% 如果相邻的块足够近则连上
artifacts_masks imclose(artifacts_masks, strel(square,10));找到接缝上所有像素点的坐标contourTracingofSeam.m:
function [ BoundaryPts ] contourTracingofSeam( mask_seam )
% tracing the image seam of the stitched image
% B_seam: binary image with only seam mask
% BoundaryPts: contour coordinates [rows, cols]
% the contour points are in white region% 函数功能拼接图的接缝轮廓跟踪
% 输入mask_seam逻辑二值接缝
% 返回值BoundaryPts论文坐标[行列]Movement [0, 1;1, 1;1, 0;1,-1;0,-1;-1,-1;-1, 0;-1, 1];% eight directions1--E, 2--SE, 3--S, 4--SW, 5--W, 6--NW, 7--N, 8--NE% 坐标系8个方向[sz1,sz2] size(mask_seam);conv_kernel ones(3,3)./9;conv_mask imfilter(double(mask_seam), conv_kernel); % 对mask_seam卷积conv_mask conv_mask.*mask_seam;[se_row, se_col] ind2sub([sz1,sz2], find(conv_mask2/9)); % 0.2222是起点ind2sub找与条件相同的索引start_pts [se_row(1), se_col(1)]; % 接缝的起点坐标
% end_pts [se_row(2), se_col(2)];max_num 2*sum(mask_seam(:)); % 因为接缝的索引是1和就是接缝像素的个数×2就是行列坐标BoundaryPts zeros(max_num,2); % 存接缝所在的行列BoundaryPtsNO 1; % 序号BoundaryPts(BoundaryPtsNO,1)start_pts(1); % 初始化为起点BoundaryPts(BoundaryPtsNO,2)start_pts(2);EndFlag false;% 从起点开始八个方向上搜索for i1:8tmpi start_pts(1) Movement(i,1);tmpj start_pts(2) Movement(i,2);if tmpi1 tmpj1 tmpisz1 tmpjsz2 mask_seam(tmpi, tmpj)0ClockDireaction i; break;endend% 绕起点1圈第一个为0的像素记为ClockDireaction% 0 0 0% 0 0.2222 0.3333% 0 0 0% 八个方向为% 6 7 8% 5 起点 1% 4 3 2%% current version needs revision _by ltl 2022 4/18BoundaryPtsNO BoundaryPtsNO 1;while (~EndFlag) for k0:1:7 % 顺时针找点tmpiBoundaryPts(BoundaryPtsNO-1,1) Movement(mod(kClockDireaction-1,8)1,1); tmpjBoundaryPts(BoundaryPtsNO-1,2) Movement(mod(kClockDireaction-1,8)1,2);if (tmpi1 || tmpj1 || tmpisz1 || tmpjsz2)continue;end if mask_seam(tmpi,tmpj)1 %find the first white point in clockwise in the 8-neighborhoodbreak; % 找到了就跳出循环endendif ismember([tmpi,tmpj],BoundaryPts,rows) break; % 如果点在集合里说明找完了endBoundaryPts(BoundaryPtsNO,1) tmpi; % 将找到的新的点添加到点集中BoundaryPts(BoundaryPtsNO,2) tmpj;BoundaryPtsNO BoundaryPtsNO 1; % 更新序号更新方向ClockDireaction mod(kClockDireaction4,8)1;
% if tmpiend_pts(1) tmpjend_pts(2)
% EndFlag true;
% endif BoundaryPtsNOmax_numfprintf( Warning! searching number exceeds the max_num in contour tracing, please find the BUG!\n);EndFlag true;endendBoundaryPts BoundaryPts(1:BoundaryPtsNO-1,:); % 截取接缝所在的点的有效范围去掉没有存储的部分% fprintf(Contour tracing finished! total %d pixels traced.\n, BoundaryPtsNO-1);end计算SSIM误差evalSSIMofSeam.m
function [ denoised_signal, eval_signal, patch_coor ] evalSSIMofSeam(img1, img2, C_lap, seam_pts, patchsize)
% evaluate the seam according to patch difference between input images (img1,img2)
% 函数功能通过论文3.2中的公式2评估接缝对应算法流程中的步骤2
% 具体思路通过遍历接缝上的像素点seam_pts 中的坐标在图像中提取相应的图像块通过 patch_coor 记录的坐标范围并计算这些图像块之间的 SSIM 误差。
% img1,img2: 翘曲后的参照图和目标图
% C_lap: 重叠区域mask
% seam_pts: 接缝上所有的像素点集
% patchsize块大小为21即每个像素点的块大小是20×20
% denoised_signal去噪后的SSIM误差n×1n为接缝上像素点的个数值就是误差
% eval_signal: 去噪前的SSIM误差
% patch_coor: 每个像素点的块范围n×4[y1 y2 x1 x2]
% 这部分是Quality evaluation-based iterative seam estimation for image stitching中写到的bound_num size(seam_pts,1); % 接缝像素点个数
eval_signal zeros(bound_num,1);
patch_coor zeros(bound_num, 4); % 记录坐标范围for i1:bound_numi_bound seam_pts(i,1);j_bound seam_pts(i,2);% 每个接缝上的点所在的块的范围20×20中心是接缝像素点y1 max(i_bound-(patchsize-1)/2, 1);y2 min(i_bound(patchsize-1)/2, size(img1,1));x1 max(j_bound-(patchsize-1)/2, 1);x2 min(j_bound(patchsize-1)/2, size(img1,2));patch_coor(i,:) [y1, y2, x1, x2]; % 图像块patch_mask C_lap(y1:y2,x1:x2); % 掩码img1_crop img1(y1:y2,x1:x2,:).*cat(3,patch_mask,patch_mask,patch_mask);img2_crop img2(y1:y2,x1:x2,:).*cat(3,patch_mask,patch_mask,patch_mask);ssim_error1 ssim(img1_crop(:,:,1), img2_crop(:,:,1));ssim_error2 ssim(img1_crop(:,:,2), img2_crop(:,:,2));ssim_error3 ssim(img1_crop(:,:,3), img2_crop(:,:,3));ssim_error (ssim_error1 ssim_error2 ssim_error3)/3; %三个通道平均SSIM误差eval_signal(i) (1-ssim_error)/2; % 记录所有接缝像素点的SSIM误差去噪前
enddenoised_signal signalDenoise(eval_signal); % 去噪end添加可视化代码显示接缝上错位的像素对应的块
%% 可视化错位的块
imgseam imgout.*cat(3,(A|B)-artifacts_masks,(A|B)-artifacts_masks,(A|B)-artifacts_masks) cat(3,artifacts_masks,zeros(sz1,sz2),zeros(sz1,sz2));figure,imshow(imgseam);
title(final stitching seam);说明初始的接缝在红色部分的对齐效果不好。
2.3.7 论文算法核心代码对应论文中3.4和3.5部分
该模块是将上一步得到的未对齐的接缝块重新修补计算得到新的接缝。具体步骤是
得到接缝上未对齐的块的部分用SIFTflow算法和sigmoid函数调整得到流向量重新对齐未对齐的块计算接缝
代码如下
%% delete photometric misaligned patches, preserve geometric misaligned patches for correspondences insertion
% 删除光度错位斑块保留几何错位斑块以便插入对应点
% 纹理融合和修复[L,n] bwlabel(artifacts_masks); % 连通区域标签化L中每个连通区域的像素值为1n个连通区域
As2 As; % 目标图和参照图掩码
Bs2 Bs;
for i1:ntmp_L Li; % 逻辑矩阵每个连通区域的逻辑标签矩阵% 找到连通区域坐标范围s_y, e_y, s_x, e_x[tmpm, tmpn]ind2sub([sz1,sz2],find(tmp_L)); s_y min(tmpm); e_y max(tmpm);s_x min(tmpn); e_x max(tmpn);% 提取连通区域在原图像和掩码中的范围crop_img1 warped_img1(s_y:e_y,s_x:e_x,:);crop_img2 warped_img2(s_y:e_y,s_x:e_x,:);s_c_img1 As(s_y:e_y,s_x:e_x);s_c_img2 Bs(s_y:e_y,s_x:e_x);% SIFTflow的向量流调整未对齐的块得到调整后的翘曲图[w_c_img1, w_c_img2]realignmentviaSIFTflow(crop_img1, crop_img2, s_c_img1);% 未对齐的部分重新计算接缝[seam_As, seam_Bs] blendTexture_clean(w_c_img1, w_c_img2, s_c_img1, s_c_img2);% 未对齐的部分替换新的接缝As2(s_y:e_y,s_x:e_x)seam_As;Bs2(s_y:e_y,s_x:e_x)seam_Bs;warped_img1(s_y:e_y,s_x:e_x,:)w_c_img1; warped_img2(s_y:e_y,s_x:e_x,:)w_c_img2;
endseam_cut gradient_blend(warped_img1, As2, warped_img2); SIFTflow算法与sigmoid平滑流向量realignmentviaSIFTflow.m
function [warpI1, warpI2 ] realignmentviaSIFTflow(im1, im2, mask_p)
% 函数功能使用SIFT flow重新对齐两张图像对应论文的3.4部分
% im1、im2两张输入图像
% mask_p: 目标图掩码
% warpl1、warpl2: 重新对齐后的两张图像
% 该方法的论文——SIFT Flow: Dense Correspondence across Scenes and Its Applications
%% pre-process
% SIFTflow的默认参数块大小3步长1
cellsize3;
gridspacing1;
SIFTflowpara.alpha2*255;
SIFTflowpara.d40*255;
SIFTflowpara.gamma0.005*255;
SIFTflowpara.nlevels4;
SIFTflowpara.wsize2;
SIFTflowpara.topwsize10;
SIFTflowpara.nTopIterations 60;
SIFTflowpara.nIterations 30;%% sift flow
% 算法1中步骤8
sift1 mexDenseSIFT(im1,cellsize,gridspacing);
sift2 mexDenseSIFT(im2,cellsize,gridspacing);
[vx,vy,~]SIFTflowc2f(sift2,sift1,SIFTflowpara);% vx、vy是水平位移和垂直位移
[h_im1,w_im1,nchannels]size(im1);
[h_vx, w_vx]size(vx);
[py, px] ind2sub([h_im1,w_im1],find(mask_p));
seam_pts [px, py]; % 目标图中需要处理的像素位置%% smoothly realignment
% 论文3.4 公式4算法1中步骤9% 创建网格
[xx1,yy1]meshgrid(1:w_im1,1:h_im1);
[XX,YY]meshgrid(1:w_vx,1:h_vx);
vec_XY [XX(:), YY(:)];% 通过检查Seam的末端像素来确定Seam的方向即Seam是从右到左、从上到下还是从下到上
orth_v [1,0];%[m_vy, -m_vx];
if sum(mask_p(:,end))h_im1 % if seam is right-leftorth_v [-1,0];%[m_vy, -m_vx];
end
if sum(mask_p(1,:))w_im1 % if seam is up-downorth_v [0,1];%[m_vy, -m_vx];
end
if sum(mask_p(end,:))w_im1 % if seam is down-uporth_v [0,-1];%[m_vy, -m_vx];
end% 计算投影变换引入sigmoid以便在接缝块附加产生光滑效果
corner_x orth_v*[0, 0, w_im1-1, w_im1-1; 0, h_im1-1, 0, h_im1-1];
max_x max(corner_x);
min_x min(corner_x);
proj_x (sum(repmat(orth_v,length(vec_XY),1).*(vec_XY-1),2)-min_x)/(max_x-min_x);
proj_y 1./(1exp(-8.*(proj_x-0.5))); % 论文中公式4图4可视化了beta148。实验中为8这里写死了
smooth_v reshape(proj_y, [h_vx, w_vx]);
smooth_vx vx.*smooth_v; % 用光流得到光滑的投影偏移用于调整图像块的位置
smooth_vy vy.*smooth_v;%% vector flow calculation
% 计算向量流
XX1XXsmooth_vx;
YY1YYsmooth_vy;
XX1min(max(XX1,1),w_im1); YY1min(max(YY1,1),h_im1);%% patch re-alignment
% 双三次插值根据重新对齐的坐标合并
warpI1 zeros(h_vx,w_vx,nchannels);
warpI2 im2;
for i1:nchannelsfoo1interp2(xx1,yy1,im1(:,:,i),XX1,YY1,bicubic);warpI1(:,:,i)foo1;
endend重新计算接缝blendTexture_clean.m:
略。过程与blendTexture.m基本一致。
添加可视化代码显示修正后块中的接缝
%% 可视化未对齐的块的接缝
imgout gradient_blend(warped_img1, As, warped_img2); SE_seam strel(diamond, 1);
As_seam imdilate(As, SE_seam) C;
Cs_seam As_seam Bs;
imgseam imgout.*cat(3,(C|C)-Cs_seam,(C|C)-Cs_seam,(C|C)-Cs_seam) cat(3,Cs_seam,zeros(sz1,sz2),zeros(sz1,sz2));figure,imshow(imgseam);
title(final stitching seam); 至此本文的代码已解读完毕。
3. 总结思考试图创新
回到论文的最开始作者阐述的创新点
解决基于接缝感知的算法得到的接缝有错位的问题。解决方法SIFTflowsigmoid平滑。可以将这套体系加入到任何接缝切割方法的框架中。
阅读完代码后相信你已经对论文有了更深入的理解。那么读懂了之后你自己写论文的创新点不就来了吗把本文的图像融合框架加入到你自己方法的框架中。以后的论文涉及到接缝融合那么基本上就都是这篇了而不是之前的基于感知的接缝融合算法了。将其他方法应用到你的框架里是最好的创新方法。Stable Linear Structures and Seam Measurements for Parallax Image Stitching这篇文章不知道大家看没看过TCSVT的基本就是没啥创新把LPC改写了一遍LPC中没细写的比如Perception-based-seam-cutting给丰富到论文中了。其实LPC等传统方法的接缝融合都用了Perception-based-seam-cutting。我当时就是看了这篇论文都能中TCSVT我就投的TCSVT直接被拒了现在想想水还是太深。但TCSVT初审是真的快如果你的创新点自己觉得还不错但又没有那么好我第一个推荐投TCSVT。
再回顾一下算法
用Perception-based-seam-cutting计算一个初始接缝用Quality evaluation-based iterative seam estimation for image stitching中提到的SSIM误差评估一下接缝如果接缝中有错位的像素则修改怎么判别呢就是大于Ostu阈值的就是错位的像素。将错位的像素和块提取出来用SIFTflow计算两个错位的像素块得到流向量用sigmoid函数平滑流向量根据平滑后的流向量重新对齐之前错位的块对上面的块重新计算接缝将修正过的块以及接缝替换原图中的错位的块。
复现一下beta取148时的不同接缝效果 beta1 beta4 beta8 初始接缝Perception-based-seam-cutting vs 本文算法修正后
初始接缝Perception-based-seam-cutting
本文算法修正后 可以看到优化后的接缝明显优于初始的接缝。因为初始的接缝破坏了重叠区域原有的结构接缝穿过了马路牙子导致错位而用本文的算法修正后效果会好很多。
其他的结果留给同学们自己复现吧。
其他想法经过测试算法还是有点慢的考虑能不能提升速度。如果速度能有优势的话又可以丰富论文内容了。
感谢同学们阅读本文如果对你有所帮助点个赞点个收藏吧我们下一篇论文源码精读再会
