clear;close all
%% read pictures
directory = ['6/'];
files = dir(directory);
files = files(3:end);
%% shuffle
randIndex = randperm(numel(files));
files = files(randIndex);
%% preprocess
N = numel(files);
dataset = {};
cnt = 1;
for i = 1:N
if files(i).name(1) ~= '.'
im = imread(strcat(directory,files(i).name));
im = double(imrotate(imresize(im, [480, 640]), 0))/255;
dataset{cnt} = im;
cnt = cnt + 1;
end
end
%% load dataset and run
order_list = get_order_list(dataset);
panorama = mymosaic(dataset, order_list);
imshow(panorama);
imwrite(panorama,'panorama.png');
%% 质量评估
if size(panorama, 3) == 3
grayPanorama = rgb2gray(panorama);
else
grayPanorama = panorama;
end
meanValue = mean2(grayPanorama);
stdValue = std2(grayPanorama);
imageEntropy = entropy(grayPanorama);
[Gx, Gy] = imgradientxy(grayPanorama, 'prewitt');
gradientMagnitude = mean2(sqrt(double(Gx).^2 + double(Gy).^2));
brisqueScore = brisque(panorama);
niqeScore = niqe(panorama);
piqeScore = piqe(grayPanorama);
fprintf('\n===== 图像质量评估结果 =====\n');
fprintf('基础指标:\n');
fprintf(' 图像均值 (亮度): %.2f\n', meanValue);
fprintf(' 图像标准差 (对比度): %.2f\n', stdValue);
fprintf(' 图像熵 (信息量): %.4f\n', imageEntropy);
fprintf(' 梯度幅值 (清晰度): %.2f\n', gradientMagnitude);
fprintf('\n=====无参考质量指标=====\n');
fprintf(' BRISQUE: %.2f (越低越好)\n', brisqueScore);
fprintf(' NIQE: %.2f (越低越好)\n', niqeScore);
fprintf(' PIQE: %.2f (越低越好)\n', piqeScore);
%% 特征提取函数
function [features, points] = get_features(img)
if ndims(img) == 3
gray_img = rgb2gray(img);
else
gray_img = img;
end
points = detectSURFFeatures(gray_img);
[features, points] = extractFeatures(gray_img, points);
end
%% 生成拼接顺序
function [order_list] = get_order_list(dataset)
level_matrix = get_level_matrix(dataset);
N = numel(dataset);
for i = 1:N
[~, index] = sort(level_matrix(i,:), 'descend');
rank_matrix(i,:) = index;
end
base = find_best(rank_matrix);
order_list(1) = base;
j = base;
while length(order_list) < N
k = 1;
while ismember(rank_matrix(j,k), order_list)
k = k + 1;
end
order_list(end+1) = rank_matrix(j,k);
j = order_list(end);
end
end
%% 选择基准图像
function [base] = find_best(rank_matrix)
N = length(rank_matrix(:,1));
rank_list = zeros(1, N);
freq_matrix = zeros(N, N);
for i = 1:N
table = tabulate(rank_matrix(:, i));
j = length(table(:, 2)) + 1;
for k = j:N
table(k, :) = 0;
end
freq_matrix(:, i) = table(:, 2);
end
for i = 1:N
for j = 1:N-1
rank_list(i) = rank_list(i) + (N - j) * freq_matrix(i, j);
end
end
[~, base] = max(rank_list);
end
%% 生成相似度矩阵
function [level_matrix] = get_level_matrix(dataset)
N = numel(dataset);
level_matrix = zeros(N, N);
[feature_map, ~] = get_feature_map(dataset);
for i = 1:N
for j = 1:N
if i ~= j
level_matrix(i, j) = kd_match(feature_map(i).cluster', feature_map(j).cluster');
end
end
end
end
%% KD树匹配
function [match_points_num] = kd_match(descs1, descs2)
n1 = size(descs1,2);
match = zeros(n1, 1);
kdtree = KDTreeSearcher(descs2');
for i = 1:size(descs1,2)
desc = descs1(:, i);
[idx, ~] = knnsearch(kdtree, desc', 'K', 2);
nn_1 = descs2(:,idx(1));
nn_2 = descs2(:,idx(2));
if sum((desc - nn_1).^2)/sum((desc - nn_2).^2) < 0.6
match(i) = idx(1);
else
match(i) = 0;
end
end
match_points_num = sum(match ~= 0);
end
%% 提取特征
function [feature_map, points_map] = get_feature_map(dataset)
N = numel(dataset);
for i = 1:N
[feature, points] = get_features(dataset{i});
feature_map(i)=struct('cluster',feature);
points_map(i)=struct('cluster',points);
end
end
%% 普通融合：替换基于边界距离的融合
function [img_b] = map_pairs(H, img_i, img_b)
hei_i = size(img_i,1);
wid_i = size(img_i,2);
hei_b = size(img_b,1);
wid_b = size(img_b,2);
% 计算图像i映射到图像b后的四角坐标（用于确定扩展范围）
ul = H*[1 1 1]'; ul = ul/ul(end); % 左上
ur = H*[wid_i, 1, 1]'; ur = ur/ur(end); % 右上
bl = H*[1, hei_i, 1]'; bl = bl/bl(end); % 左下
br = H*[wid_i, hei_i, 1]'; br = br/br(end); % 右下
% 扩展图像b，确保能容纳映射后的图像i
pad_up = 0; pad_left = 0; pad_down = 0; pad_right = 0;
if max(br(1),ur(1)) > wid_b
pad_right = round(max(br(1),ur(1)) - wid_b + 30);
img_b = padarray(img_b, [0, pad_right], 'post');
end
if max(br(2), bl(2)) > hei_b
pad_down = round(max(br(2), bl(2)) - hei_b + 30);
img_b = padarray(img_b, [pad_down, 0], 'post');
end
if min(ul(1), bl(1)) <= 0
pad_left = round(-min(ul(1), bl(1)) + 30);
img_b = padarray(img_b, [0, pad_left], 'pre');
end
if min(ul(2), ur(2)) <= 0
pad_up = round(-min(ul(2), ur(2)) + 30);
img_b = padarray(img_b, [pad_up, 0], 'pre');
end
H_inv = inv(H); % 逆单应矩阵
% 生成图像b的坐标网格并映射到图像i
[y_b, x_b] = meshgrid(...
round(pad_up + min(ul(2), ur(2))):round(pad_up + max(bl(2), br(2))), ...
round(pad_left + min(ul(1), bl(1))):round(pad_left + max(br(1), ur(1))));
y_b = y_b(:); x_b = x_b(:);
xy = H_inv * [x_b - pad_left, y_b - pad_up, ones(size(x_b,1),1)]';
x_i = int64(xy(1,:)' ./ xy(3,:)'); % 图像i的x坐标
y_i = int64(xy(2,:)' ./ xy(3,:)'); % 图像i的y坐标
% 普通融合：重叠区域取平均，非重叠区域直接覆盖
% 1. 重叠区域（同时在img_i和img_b范围内）：平均融合
indices = x_i > 0 & x_i <= wid_i & y_i > 0 & y_i <= hei_i ...
& x_b > 0 & x_b <= size(img_b,2) & y_b > 0 & y_b <= size(img_b,1);
% 对RGB三通道取平均
img_b(sub2ind(size(img_b), y_b(indices), x_b(indices))) = ...
(img_i(sub2ind(size(img_i), y_i(indices), x_i(indices))) + ...
img_b(sub2ind(size(img_b), y_b(indices), x_b(indices)))) / 2;
if ndims(img_b) == 3
% 确保所有索引向量具有相同的大小
y_idx = y_b(indices);
x_idx = x_b(indices);
len = length(indices);
% 绿色通道 - 修正索引维度
green_indices = sub2ind(size(img_b), y_idx, x_idx, 2*ones(len, 1));
img_b(green_indices) = ...
(img_i(sub2ind(size(img_i), y_i(indices), x_i(indices), 2*ones(len, 1))) + ...
img_b(green_indices)) / 2;
% 蓝色通道 - 修正索引维度
blue_indices = sub2ind(size(img_b), y_idx, x_idx, 3*ones(len, 1));
img_b(blue_indices) = ...
(img_i(sub2ind(size(img_i), y_i(indices), x_i(indices), 3*ones(len, 1))) + ...
img_b(blue_indices)) / 2;
end
% 2. 非重叠区域（仅在img_i范围内）：直接复制img_i的像素
indices = x_i > 0 & x_i <= wid_i & y_i > 0 & y_i <= hei_i ...
& ~(x_b > 0 & x_b <= size(img_b,2) & y_b > 0 & y_b <= size(img_b,1));
% 确保所有索引向量具有相同的大小
y_idx = y_b(indices);
x_idx = x_b(indices);
len = length(indices);
img_b(sub2ind(size(img_b), y_idx, x_idx)) = ...
img_i(sub2ind(size(img_i), y_i(indices), x_i(indices)));
if ndims(img_b) == 3
green_indices = sub2ind(size(img_b), y_idx, x_idx, 2*ones(len, 1));
img_b(green_indices) = img_i(sub2ind(size(img_i), y_i(indices), x_i(indices), 2*ones(len, 1)));
blue_indices = sub2ind(size(img_b), y_idx, x_idx, 3*ones(len, 1));
img_b(blue_indices) = img_i(sub2ind(size(img_i), y_i(indices), x_i(indices), 3*ones(len, 1)));
end
end
%% 全景拼接主函数
function [img_b] = mymosaic(img_input, order_list)
N = numel(img_input);
fprintf('Image %d is initial...\n',order_list(1));
img_b = img_input{order_list(1)};
for i = 2:N
fprintf('Blending image %d...\n',order_list(i));
img_i = img_input{order_list(i)};
img_b = stitch(img_i, img_b);
end
fprintf('Done!\n');
end
%% 单图拼接函数
function [img_b] = stitch(img_i, img_b)
[features_i, points_i] = get_features(img_i);
[features_b, points_b] = get_features(img_b);
pairs = matchFeatures(features_i, features_b);
matchedBoxPoints = points_i(pairs(:, 1), :);
matchedScenePoints = points_b(pairs(:, 2), :);
x_i = matchedBoxPoints.Location(:,1);
y_i = matchedBoxPoints.Location(:,2);
x_b = matchedScenePoints.Location(:,1);
y_b = matchedScenePoints.Location(:,2);
[H, ~] = ransac_est_homography(x_i, y_i, x_b, y_b, 10);
img_b = map_pairs(H, img_i, img_b);
end
%% RANSAC估计单应矩阵
function [H, inlier_ind] = ransac_est_homography(x1, y1, x2, y2, thresh)
N = numel(x1);
max_inliers = 0;
H = eye(3);
ssd = @(x,y) sum((x-y).^2);
for t = 1:1000
r_idx = randi(N,4,1);
H_t = est_homography(x2(r_idx),y2(r_idx),x1(r_idx),y1(r_idx));
inliers = 0;
for i = 1:N
t_xy = H_t*[x1(i), y1(i), 1]';
t_xy = t_xy/t_xy(end);
if ssd([x2(i), y2(i), 1]', t_xy) < thresh
inliers = inliers + 1;
end
end
if inliers > max_inliers
max_inliers = inliers;
H = H_t;
end
end
inlier_ind = find(arrayfun(@(i) ssd(H*[x1(i);y1(i);1], [x2(i);y2(i);1]) < thresh, 1:N));
end
%% 计算单应矩阵
function H = est_homography(X,Y,x,y)
A = zeros(2*length(x),9);
for i = 1:length(x)
A(2*i-1,:) = [x(i), y(i), 1, 0, 0, 0, -X(i)*x(i), -X(i)*y(i), -X(i)];
A(2*i,:) = [0, 0, 0, x(i), y(i), 1, -Y(i)*x(i), -Y(i)*y(i), -Y(i)];
end
[~,~,V] = svd(A);
H = reshape(V(:,9),3,3)';
H = H/H(3,3);
end错误使用 sub2ind (第 61 行)
下标参量必须为标量或大小相同的数组。
出错 pingjunjiaquanronghe>map_pairs (第 213 行)
green_indices = sub2ind(size(img_b), y_idx, x_idx, 2*ones(len, 1));
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
出错 pingjunjiaquanronghe>stitch (第 271 行)
img_b = map_pairs(H, img_i, img_b);
^^^^^^^^^^^^^^^^^^^^^^^^^^
出错 pingjunjiaquanronghe>mymosaic (第 254 行)
img_b = stitch(img_i, img_b);
^^^^^^^^^^^^^^^^^^^^
出错 pingjunjiaquanronghe (第 26 行)
panorama = mymosaic(dataset, order_list);
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
代码修正，然后输出全套代码
最新发布
博客提到在LC646解题思路中看到方法但不懂，还给出相似方法，即A<B按字符顺序为升序，A>B为降序，并给出了LC646的地址。
扫一扫
2418
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
