function SPHP_stitching(in_name, edge_list, ref, tar, warp_type, zeroR_ON)
% 输入参数:
% in_name: 一个cell数组，包含要拼接的图像的文件名
% edge_list: 一个矩阵，每一行包含两个图像的索引，表示这两个图像是相邻的
% ref: 参考图像的索引
% tar: 目标图像的索引（通常是要拼接到的图像）
% warp_type, zeroR_ON:
img_n = size(in_name, 1); % 图像的数量
edge_n = size(edge_list, 1); % 相邻图像对的数量
% 检查edge_list中的索引对是否有自环（即i不等于j）
for ei = 1 : edge_n
i = edge_list(ei, 1);
j = edge_list(ei, 2);
if i == j
fprintf('Index pair error.\n'); % 打印错误信息
pause; % 暂停程序
end
end
% 加载图像
I = cell(img_n, 1); % 初始化一个cell数组来存储图像
for i = 1 : img_n
I{i} = imread(in_name{i}); % 读取图像
end
% 预处理
% T{i}是一个坐标变换矩阵，用于将图像坐标（原点在图像的左上角，x向右，y向下）
% 转换为标准坐标（原点在图像中心，x向右，y向上）
max_size = 1000 * 1000; % 最大像素数量限制
imgw = zeros(img_n, 1); % 初始化图像的宽度数组
imgh = zeros(img_n, 1); % 初始化图像的高度数组
resized_in_name = cell(img_n, 1); % 初始化一个cell数组来存储调整大小后的图像的文件名
T = cell(img_n, 1); % 初始化一个cell数组来存储坐标变换矩阵
for i = 1 : img_n
% 如果图像的像素数量超过了最大限制，则对其进行下采样
if numel(I{i}(:, :, 1)) > max_size
I{i} = imresize(I{i}, sqrt(max_size / numel(I{i}(:, :, 1))));
end
resized_in_name{i} = sprintf('%s_resized%s', in_name{i}(1:end-4), in_name{i}(end-3:end)); % 生成新的文件名
imwrite(I{i}, resized_in_name{i}); % 保存调整大小后的图像
imgw(i) = size(I{i}, 2); % 获取图像的宽度
imgh(i) = size(I{i}, 1); % 获取图像的高度
% 计算坐标变换矩阵
T{i} = [1 0 -(imgw(i)+1)/2; ...
0 -1 (imgh(i)+1)/2; ...
0 0 1];
end
% 将图像转换为双精度浮点数格式
for i = 1 : img_n
I{i} = im2double(I{i});
end
% 计算图像对之间的单应性矩阵
H = cell(img_n, img_n); % H{i,j}是从I{j}到I{i}的变换矩阵（使用标准坐标）
for i = 1 : img_n
H{i, i} = eye(3); % 对角线上的元素是单位矩阵
end
for ei = 1 : edge_n
i = edge_list(ei, 1);
j = edge_list(ei, 2);
tmpH = sift_mosaic(I{i}, I{j}); % 使用SIFT算法计算从I{i}到I{j}的单应性矩阵（使用图像坐标）
disp(tmpH)
% 将坐标从图像坐标转换为标准坐标
H{j, i} = T{j} * tmpH * inv(T{i});
H{i, j} = inv(H{j, i}); % 计算逆变换
end
mdltpara = [];
mdlt_results = [];
% 对于每张图像I，计算从I到参考图像的单应性矩阵，以及从I到目标图像的单应性矩阵
% 构造图
G = sparse([edge_list(:, 1); edge_list(:, 2)], ...
[edge_list(:, 2); edge_list(:, 1)], ...
[ones(edge_n, 1); ones(edge_n, 1)], img_n, img_n); % 使用稀疏矩阵表示图
disp(edge_list(:, 1))
% 计算从第i张图像到参考图像的单应性矩阵
for i = 1 : img_n
if i == ref
continue;
end
[~, path, ~] = graphshortestpath(G, i, ref); % 计算从第i张图像到参考图像的最短路径
tmpH = eye(3);
for ppi = 1 : numel(path)-1
tmpH = H{path(ppi+1), path(ppi)} * tmpH; % 累积单应性矩阵
end
H{ref, i} = tmpH;
H{i, ref} = inv(H{ref, i}); % 计算逆变换
end
% 计算从第i张图像到目标图像的单应性矩阵
for i = 1 : img_n
if i == tar
continue;
end
[~, path, ~] = graphshortestpath(G, i, tar); % 计算从第i张图像到目标图像的最短路径
tmpH = eye(3);
for ppi = 1 : numel(path)-1
tmpH = H{path(ppi+1), path(ppi)} * tmpH; % 累积单应性矩阵
end
H{tar, i} = tmpH;
H{i, tar} = inv(H{tar, i}); % 计算逆变换
end
% 规范化所有的单应性矩阵，使得h9 = 1
for hi = 1 : numel(H)
if H{hi}
H{hi} = H{hi} / H{hi}(3, 3);
end
end
% 提取从参考图像到目标图像的单应性矩阵的参数，以便计算我们的warp
H0 = H{tar, ref};
A = H0(1:2, 1:2);
t = H0(1:2, 3);
h31 = H0(3, 1);
h32 = H0(3, 2);
B = A - t*[h31 h32];
c = sqrt(h31^2+h32^2);
theta = atan2(-h32, -h31); % 计算ub1和ub2所需的角度
disp(H0)
disp(A)
disp(t)
disp(h31)
disp(h32)
return;
%% compute ub1 and ub2
% 这部分代码的目的是计算图像变换后的边界值ub1和ub2，以及为优化过程准备一个成本表totalcost_table。
% 初始化一个用于存储计算结果的数组cu，其大小与图像数量img_n相同。
% cu数组将用于存储每张图像经过变换后在参考图像中的x坐标。
cu = zeros(img_n, 1);
% 对每一张图像进行遍历。
for i = 1 : img_n
% 使用apply_transform函数将第i张图像的中心点(0, 0)变换到参考图像的坐标系中。
% H{ref, i}是一个单应性矩阵，描述了从第i张图像到参考图像的变换。
% 返回的tmpx和tmpy是经过变换后的x和y坐标。
[tmpx, tmpy] = apply_transform(0, 0, H{ref, i});
% 再次使用apply_transform函数，这次是为了应用一个旋转变换。
% theta是旋转角度，这里没有给出其定义，可能是在代码的其他部分定义的。
% 这个旋转变换将上一步得到的tmpx和tmpy进一步变换。
% 注意这里的~[...]表示我们不关心这个输出值，即用~忽略了第二个输出。
[cu(i), ~] = apply_transform(tmpx, tmpy, [cos(theta), sin(theta), 0; ...
-sin(theta), cos(theta), 0; ...
0, 0, 1]);
end
% 从cu数组中找到最小和最大的x坐标值，这些值将用于定义优化的边界。
oriub1 = min(cu);
oriub2 = max(cu);
% 定义采样间隔s_itv，这里设置为20。
s_itv = 20; % sample interval
% 使用meshgrid函数创建两个二维网格矩阵offset_table1和offset_table2。
% 这些矩阵将用于在优化过程中尝试不同的偏移量组合。
% meshgrid的参数定义了网格的范围和步长，这里offset_table1的范围是-300到600，步长为10；
% offset_table2的范围是-300到200，步长也为10。
[offset_table1, offset_table2] = meshgrid(-300:10:600, -300:10:200);
% 初始化totalcost_table为零矩阵，其大小与offset_table1相同。
% 这个矩阵将用于存储在优化过程中计算的总成本。
totalcost_table = zeros(size(offset_table1));
% 遍历offset_table1和offset_table2的所有元素，这两个表格定义了可能的偏移量组合
for oi = 1 : size(offset_table1, 1)
for oj = 1 : size(offset_table1, 2)
% 根据当前的偏移量组合更新ub1和ub2的值
ub1 = oriub1 + offset_table1(oi, oj);
ub2 = oriub2 + offset_table2(oi, oj);
% 如果ub2和ub1的差小于120，则将该组合的成本设置为NaN，并继续下一个组合
if ub2 - ub1 < 120
totalcost_table(oi, oj) = nan;
continue;
end
% 计算用于图像变换的系数c1para
c1para = compute_c1_warp_coeff(H0, t, c, theta, ub1, ub2, zeroR_ON);
% 定义旋转变换的逆矩阵
invR = [cos(theta), sin(theta), 0; ...
-sin(theta), cos(theta), 0; ...
0, 0, 1];
% 初始化成本列表，用于存储每个图像的成本
cost_list = zeros(1, img_n);
% 遍历每个图像
for i = 1 : img_n
% 在图像上创建网格，用于后续的坐标变换
x = linspace(1, imgw(i), ceil(imgw(i)/s_itv));
y = linspace(1, imgh(i), ceil(imgh(i)/s_itv));
[x, y] = meshgrid(x, y);
% 应用变换T{i}来改变坐标
[x, y] = apply_transform(x, y, T{i});
% 计算从第i个图像到参考图像的变换矩阵tmpH
tmpH = invR * H{ref, i};
% 应用变换tmpH得到在参考图像中的坐标
[u, v] = apply_transform(x, y, tmpH);
% 根据u的值定义三个区域
reg1_mask = u < ub1;
reg2_mask = u > ub1 & u < ub2;
reg3_mask = u > ub2;
% 初始化用于存储Jacobian矩阵元素的映射
J11_map = zeros(size(x));
J12_map = zeros(size(x));
J21_map = zeros(size(x));
J22_map = zeros(size(x));
% 遍历每个区域
for reg = 1 : 3
% 根据区域设置相应的掩码
if reg == 1
reg_mask = u < ub1;
elseif reg == 3
reg_mask = u > ub2;
else
reg_mask = ~(u < ub1 | u > ub2);
end
% 如果没有有效的点，则继续下一个区域
if nnz(reg_mask) == 0
continue;
end
% 提取满足当前区域条件的坐标
x1 = x(reg_mask);
y1 = y(reg_mask);
u1 = u(reg_mask);
v1 = v(reg_mask);
% 根据区域设置变换矩阵A
if reg == 1
A = c1para.H * H{ref, i};
elseif reg == 2
A = tmpH;
elseif reg == 3
A = c1para.S * H{ref, i};
end
% 归一化变换矩阵A
A = A / A(3, 3);
% 提取变换矩阵A的元素
h1 = A(1, 1); h2 = A(1, 2); h3 = A(1, 3); h4 = A(2, 1); h5 = A(2, 2); h6 = A(2, 3); h7 = A(3, 1); h8 = A(3, 2);
% 计算Jacobian矩阵的元素
J11 = h1./(h7*x1 + h8*y1 + 1) - (h7*(h3 + h1*x1 + h2*y1))./(h7*x1 + h8*y1 + 1).^2;
J12 = h2./(h7*x1 + h8*y1 + 1) - (h8*(h3 + h1*x1 + h2*y1))./(h7*x1 + h8*y1 + 1).^2;
J21 = h4./(h7*x1 + h8*y1 + 1) - (h7*(h6 + h4*x1 + h5*y1))./(h7*x1 + h8*y1 + 1).^2;
J22 = h5./(h7*x1 + h8*y1 + 1) - (h8*(h6 + h4*x1 + h5*y1))./(h7*x1 + h8*y1 + 1).^2;
% 如果是第二个区域，还需要进一步计算Jacobian矩阵
if reg == 2
if c1para.zeroR_ON == 0
JT11 = (2*c1para.a(1)*u1 + c1para.a(2)) .* v1 + (2*c1para.b(1)*u1 + c1para.b(2));
JT12 = c1para.a(1)*u1.^2 + c1para.a(2)*u1 + c1para.a(3);
JT21 = (2*c1para.e(1)*u1 + c1para.e(2)) .* v1 + (2*c1para.f(1)*u1 + c1para.f(2));
JT22 = c1para.e(1)*u1.^2 + c1para.e(2)*u1 + c1para.e(3);
else
JT11 = (3*c1para.a(1)*u1.^2 + 2*c1para.a(2)*u1 + c1para.a(3)) .* v1 + (2*c1para.b(1)*u1 + c1para.b(2));
JT12 = c1para.a(1)*u1.^3 + c1para.a(2)*u1.^2 + c1para.a(3)*u1 + c1para.a(4);
JT21 = (3*c1para.e(1)*u1.^2 + 2*c1para.e(2)*u1 + c1para.e(3)) .* v1 + (2*c1para.f(1)*u1 + c1para.f(2));
JT22 = c1para.e(1)*u1.^3 + c1para.e(2)*u1.^2 + c1para.e(3)*u1 + c1para.e(4);
end
tmp11 = JT11.*J11 + JT12.*J21;
tmp12 = JT11.*J12 + JT12.*J22;
tmp21 = JT21.*J11 + JT22.*J21;
tmp22 = JT21.*J12 + JT22.*J22;
J11 = tmp11; J12 = tmp12; J21 = tmp21; J22 = tmp22;
end
% 更新Jacobian矩阵映射
J11_map(reg_mask) = J11;
J12_map(reg_mask) = J12;
J21_map(reg_mask) = J21;
J22_map(reg_mask) = J22;
end
% 计算平均的alpha和beta值
% 计算平均的alpha值
% 通过将J11_map和J22_map中的所有元素相加，然后除以x的元素数量的两倍来得到
avg_alpha = (sum(J11_map(:)) + sum(J22_map(:))) / (numel(x)*2);
% 计算平均的beta值
% 通过将J21_map中的所有元素相加，然后减去J12_map中的所有元素之和，
% 最后再除以x的元素数量的两倍来得到
avg_beta = (sum(J21_map(:)) + (-sum(J12_map(:)))) / (numel(x)*2);
% 计算cost
% cost是通过计算每个map与其对应的平均值的差的平方和来得到的
cost = sum( (J11_map(:)-avg_alpha).^2 + (J12_map(:)-(-avg_beta)).^2 + ...
(J21_map(:)-avg_beta).^2 + (J22_map(:)-avg_alpha).^2 );
% 将计算得到的cost保存到cost_list中的第i个位置
cost_list(i) = cost;
end
% 计算cost_list中所有cost的总和
totalcost = sum(cost_list);
% 将计算得到的总cost保存到totalcost_table中的(oi, oj)位置
totalcost_table(oi, oj) = totalcost;
end
end
% 在totalcost_table中找到最小值的位置索引
% idx1和idx2分别表示该最小值在totalcost_table中的行索引和列索引
[idx1, idx2] = find(totalcost_table == min(totalcost_table(:)));
% 使用找到的索引idx1和idx2从offset_table1中获取对应的偏移量，并更新ub1的值
ub1 = oriub1 + offset_table1(idx1, idx2);
% 使用找到的索引idx1和idx2从offset_table2中获取对应的偏移量，并更新ub2的值
ub2 = oriub2 + offset_table2(idx1, idx2);
%% Compute parameters/coefficients of our warp
% 调用compute_c1_warp_coeff函数来计算warp的c1系数
% 输入参数包括：H0, t, c, theta, ub1, ub2, 和 zeroR_ON
% 返回的c1para将包含计算得到的warp系数
c1para = compute_c1_warp_coeff(H0, t, c, theta, ub1, ub2, zeroR_ON);
%% c1 warp for each image
% 这一部分代码将对每个图像应用c1类型的warp（扭曲/变形）
img_grid_size = 10; % for warping
% 设置用于warp的网格大小为10
[c1out, c1omask, Td] = texture_mapping(resized_in_name, imgw, imgh, img_grid_size, T, warp_type, H, ref, tar, c1para, mdltpara, 0, 'white');
% 调用texture_mapping函数来对resized_in_name中的图像进行纹理映射（可能是warp的一种形式）
% 输入参数包括：图像名称、图像宽度、图像高度、网格大小、变换矩阵T、warp类型、单应性矩阵H、参考图像、目标图像、c1warp系数、模型参数、填充选项和填充颜色
% 函数返回：warp后的图像c1out、对应的遮罩c1omask以及可能的变换数据Td
for i = 1 : img_n
% 遍历所有图像
%figure; imshow(c1out{i}); % the warped result of each image
% （可选）为每个warp后的图像创建一个新窗口并显示它
% 注意：这行代码被注释掉了，因此不会实际执行显示图像的操作
end
for i = 1 : img_n
eval(['delete ' resized_in_name{i}]);
% 遍历所有图像并删除对应的内存中的变量
% 使用eval函数执行删除操作，这可能是因为resized_in_name{i}包含了需要删除的变量的名称
% 注意：通常不推荐使用eval，因为它可能带来安全风险和执行效率问题
% 一个更好的做法可能是直接将图像数据存储在一个数组中，然后通过索引来管理它们
end
%% composite by linear blending
% 通过线性混合来合成图像的代码段
for i = 1 : img_n
% 遍历所有图像
if i == 1
% 对于第一张图像，直接将其设置为输出
out = c1out{1};
out_mask = c1omask{1};
% 计算输出图像的中心
[r, c] = find(c1omask{1}(:, :, 1));
out_center = [mean(r) mean(c)];
else % blend out and c1out{i}
% 对于其他图像，将其与之前的输出进行混合
% 计算当前图像的中心
[r, c] = find(c1omask{i}(:, :, 1));
out_i_center = [mean(r) mean(c)];
% 计算加权遮罩
vec = out_i_center - out_center; % 从out_center到out_i_center的向量
intsct_mask = c1omask{i}(:, :, 1) & out_mask(:, :, 1); % 交集遮罩，单通道
out_only_mask = out_mask(:, :, 1) - intsct_mask;
out_i_only_mask = c1omask{i}(:, :, 1) - intsct_mask;
[r, c] = find(intsct_mask(:, :, 1));
idx = sub2ind(size(c1omask{i}(:, :, 1)), r, c);
out_wmask = zeros(size(c1omask{i}(:, :, 1)));
proj_val = (r - out_center(1))*vec(1) + (c- out_center(2))*vec(2); % 向量的内积
% 计算加权遮罩的值（重叠区域的权重图）
out_wmask(idx) = (proj_val - (min(proj_val)+(1e-3))) / ...
((max(proj_val)-(1e-3)) - (min(proj_val)+(1e-3)));
% 混合图像
mask1 = out_mask(:, :, 1)&(out_wmask==0);
mask2 = out_wmask;
mask3 = c1omask{i}(:, :, 1)&(out_wmask==0);
mask1 = cat(3, mask1, mask1, mask1); mask2 = cat(3, mask2, mask2, mask2); mask3 = cat(3, mask3, mask3, mask3);
out = out.*(mask1+(1-mask2).*(mask2~=0)) + c1out{i}.*(mask2+mask3);
% 更新遮罩和中心
out_mask = out_mask | c1omask{i};
out_center = out_i_center; % 用当前图像的中心更新输出中心
end
end
c1all = out;
% 计算分母，用于处理背景色
denom = zeros(size(c1out{1}));
for i = 1 : img_n
denom = denom + c1omask{i};
end
bgcolor = 'white'; % 设置背景色
% 根据背景色设置未覆盖区域的像素值
if strcmp(bgcolor, 'white')
c1all(denom==0) = 1; % 如果是白色背景，则设置为1（白色）
else
c1all(denom==0) = 0; % 否则设置为0（黑色）
end
figure; imshow(c1all); % 显示最终结果
% for i = 1 : img_n
% imwrite(c1out{i}, sprintf('%s_out%d.png', warp_type, i));
% end
% imwrite(c1all, [warp_type '_outall.png']);
SPHP_stitching算法解读（一）
于 2024-02-02 17:28:11 首次发布
3508
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
