2013年全国大学生数学建模竞赛 碎纸片拼接问题按照以下思路来写。其实这个问题是第二问，做过第一问的应该都知道这两问的难度根本不在一个档次上。第一问的时候，每个碎纸片都可以有一个最佳的匹配，但是在第二问就不行了。
因此，我们考虑对于所有的拼好的“图”，求其总不吻合度的最小值。
每两块碎纸片拼在一起，上下(或左右)边缘的灰度值之差的绝对值之和称为这两块碎纸片的不吻合度。
总不吻合度，就是一张图片中所有的不吻合度加在一起，这个很好理解。
现在我们把问题规范一下，以便于使用模拟退火算法：求解域11*19阶矩阵的集合，并且矩阵中的元素恰好是 1到209的一个排列
目标函数f(x):假设矩阵为 X,那么把x(i,j)作为拼图的第(i,j)块。可以得到一张图。这张图的总不吻合度就是f(x)我们设计的新解生成函数：在温度较高时，以更高的概率交换两行、两列的位置；在温度较低时，以更高的概率交换不吻合度最高的两块，或者随机交换两块。对模拟退火算法的修改建议感觉可以改进的地方：
1.可以用蒙特卡洛先做出一个较好的解，再从那个开始拼。
2.不难看出最后图片变成了许多大块，可以考虑在新解生成函数中也设计这种交换大块的生成方式。
3.可以尝试先拼出四个角和边，然后再往内层一层一层地拼，不知道会不会好一点。
请根据以上的思路，写出问题2的完整代码，数据放在附件3，bmp格式，000-208
%% 碎纸片拼接完整解决方案
% 2013年全国大学生数学建模竞赛B题第二问
% 双面碎纸片自动拼接系统
%% 步骤1: 数据读取与预处理
% 从txt文件读取209个72×72灰度图像碎片
fprintf('正在读取碎纸片数据...\n');
pieces = zeros(72, 72, 209); % 存储碎片数据
for k = 0:208
filename = sprintf('%03d.bmp', k); % 文件名格式为000.txt~208.txt
if exist(filename, 'file')
pieces(:,:,k+1) = dlmread(filename); % 读取灰度数据
else
error('文件缺失: %s', filename);
end
end
fprintf('成功读取%d个碎纸片数据\n', size(pieces,3));
%% 步骤3: 蒙特卡洛初始解生成（贪心算法）
fprintf('生成初始解...\n');
% 初始化参数
rows = 11; cols = 19;
solved = zeros(rows, cols); % 解矩阵
used = false(1, 209); % 标记已使用碎片
% 随机选择左上角碎片
start_idx = randi(209);
solved(1,1) = start_idx;
used(start_idx) = true;
% 第一行贪心填充
for j = 2:cols
min_diff = inf;
best_k = 0;
current_piece = pieces(:,:,solved(1,j-1));
% 寻找最佳匹配
for k = find(~used)
diff_val = calc_diff(current_piece, pieces(:,:,k), 'left_right');
if diff_val < min_diff
min_diff = diff_val;
best_k = k;
end
end
solved(1,j) = best_k;
used(best_k) = true;
end
% 逐列贪心填充
for j = 1:cols
for i = 2:rows
min_diff = inf;
best_k = 0;
current_piece = pieces(:,:,solved(i-1,j));
% 寻找最佳匹配
for k = find(~used)
diff_val = calc_diff(current_piece, pieces(:,:,k), 'top_bottom');
if diff_val < min_diff
min_diff = diff_val;
best_k = k;
end
end
solved(i,j) = best_k;
used(best_k) = true;
end
end
% 验证所有碎片是否使用
if any(~used)
error('初始解生成失败：有碎片未被使用');
end
%% 步骤4: 模拟退火优化
fprintf('开始模拟退火优化...\n');
% 初始化参数
T0 = 1000; % 初始温度
Tf = 1e-3; % 终止温度
alpha = 0.95; % 降温系数
L = 100; % 每个温度迭代次数
current_sol = solved; % 使用贪心解作为起点
current_cost = total_cost(current_sol, pieces);
best_sol = current_sol;
best_cost = current_cost;
T = T0;
iteration = 0;
cost_history = zeros(1, 1000); % 记录成本历史
temp_history = zeros(1, 1000); % 记录温度历史
while T > Tf
iteration = iteration + 1;
for iter = 1:L
% 温度依赖的扰动策略
p = (T - Tf)/(T0 - Tf)*0.5 + 0.3; % 高温大扰动概率高
if rand < p % 大扰动：行列交换
if rand < 0.5 % 行交换
i1 = randi(11); i2 = randi(11);
new_sol = current_sol;
new_sol([i1,i2], :) = new_sol([i2,i1], :);
else % 列交换
j1 = randi(19); j2 = randi(19);
new_sol = current_sol;
new_sol(:, [j1,j2]) = new_sol(:, [j2,j1]);
end
else % 小扰动：碎片交换
new_sol = current_sol;
pos1 = [randi(11), randi(19)];
pos2 = [randi(11), randi(19)];
temp = new_sol(pos1(1), pos1(2));
new_sol(pos1(1), pos1(2)) = new_sol(pos2(1), pos2(2));
new_sol(pos2(1), pos2(2)) = temp;
end
% 计算新解成本
new_cost = total_cost(new_sol, pieces);
delta = new_cost - current_cost;
% Metropolis接受准则
if delta < 0 || rand < exp(-delta/T)
current_sol = new_sol;
current_cost = new_cost;
if new_cost < best_cost
best_sol = new_sol;
best_cost = new_cost;
fprintf('迭代%d: 发现更优解 %.2f (T=%.4f)\n', iteration, best_cost, T);
end
end
end
% 记录历史数据
if iteration <= 1000
cost_history(iteration) = best_cost;
temp_history(iteration) = T;
end
T = T * alpha; % 降温
end
fprintf('优化完成! 最优不吻合度: %.2f\n', best_cost);
%% 步骤5: 结果可视化
% 拼接完整图像
fprintf('拼接完整图像...\n');
result_image = assemble_image(best_sol, pieces);
% 显示拼接结果
figure('Name', '碎纸片拼接结果', 'Position', [100, 100, 1200, 800]);
subplot(2,2,[1,3]);
imshow(result_image, []);
title(sprintf('拼接结果 (不吻合度: %.2f)', best_cost));
% 添加网格线显示碎片边界
hold on;
for i = 1:10 % 行边界
y = i*72 + 0.5;
plot([0.5, 19*72+0.5], [y, y], 'r-', 'LineWidth', 0.5);
end
for j = 1:18 % 列边界
x = j*72 + 0.5;
plot([x, x], [0.5, 11*72+0.5], 'r-', 'LineWidth', 0.5);
end
hold off;
% 显示碎片编号分布
subplot(2,2,2);
imagesc(best_sol);
colormap(jet(209)); % 使用209种颜色的色谱
colorbar;
title('碎片位置编号分布图');
axis equal tight;
% 显示优化过程
subplot(2,2,4);
yyaxis left;
plot(1:iteration, cost_history(1:iteration), 'b-', 'LineWidth', 1.5);
ylabel('不吻合度');
yyaxis right;
semilogy(1:iteration, temp_history(1:iteration), 'r-', 'LineWidth', 1.5);
ylabel('温度(对数)');
title('优化过程');
xlabel('迭代次数');
grid on;
legend('不吻合度', '温度');
% 保存结果
imwrite(result_image, 'reconstructed_image.png');
fprintf('拼接图像已保存为 reconstructed_image.png\n');
%% 步骤6: 交互式碎片查看（可选）
fprintf('添加交互式查看功能...\n');
figure('Name', '碎纸片拼接结果(交互式)', 'Position', [100, 100, 1200, 800]);
imshow(result_image, []);
title('点击任意碎片查看详情');
hold on;
% 绘制网格
for i = 1:10
y = i*72 + 0.5;
plot([0.5, 19*72+0.5], [y, y], 'r-', 'LineWidth', 0.5);
end
for j = 1:18
x = j*72 + 0.5;
plot([x, x], [0.5, 11*72+0.5], 'r-', 'LineWidth', 0.5);
end
hold off;
% 设置点击回调函数
set(gcf, 'WindowButtonDownFcn', @imageClickCallback);
% 回调函数定义
function imageClickCallback(~, event)
pos = round(event.IntersectionPoint(1:-1:2)); % 获取点击位置
row = ceil(pos(1)/72);
col = ceil(pos(2)/72);
if row >=1 && row <=11 && col>=1 && col<=19
% 获取碎片信息
idx = best_sol(row, col);
fprintf('碎片位置: (%d, %d)\n碎片编号: %d\n', row, col, idx);
% 显示碎片详情
fig = figure('Name', '碎片详情', 'Position', [200, 200, 800, 400]);
% 显示碎片图像
subplot(1,2,1);
imshow(pieces(:,:,idx), []);
title(sprintf('碎片%d (%d,%d)', idx, row, col));
% 显示边缘特征
subplot(1,2,2);
hold on;
% 左侧边缘
if col > 1
left_diff = calc_diff(pieces(:,:,best_sol(row, col-1)), pieces(:,:,idx), 'left_right');
plot(1:72, pieces(:,1,idx), 'b-', 'LineWidth', 1.5);
end
% 右侧边缘
if col < 19
right_diff = calc_diff(pieces(:,:,idx), pieces(:,:,best_sol(row, col+1)), 'left_right');
plot(1:72, pieces(:,end,idx), 'r-', 'LineWidth', 1.5);
end
% 上侧边缘
if row > 1
top_diff = calc_diff(pieces(:,:,best_sol(row-1, col)), pieces(:,:,idx), 'top_bottom');
plot(1:72, pieces(1,:,idx), 'g-', 'LineWidth', 1.5);
end
% 下侧边缘
if row < 11
bottom_diff = calc_diff(pieces(:,:,idx), pieces(:,:,best_sol(row+1, col)), 'top_bottom');
plot(1:72, pieces(end,:,idx), 'm-', 'LineWidth', 1.5);
end
hold off;
title('边缘灰度特征');
legend('左侧', '右侧', '上侧', '下侧');
xlabel('像素位置');
ylabel('灰度值');
grid on;
end
end
%% 步骤2: 辅助函数定义
% 计算两个碎片在指定方向的不吻合度
function diff = calc_diff(piece1, piece2, direction)
if strcmp(direction, 'left_right')
edge1 = piece1(:, end); % 右边缘
edge2 = piece2(:, 1); % 左边缘
elseif strcmp(direction, 'top_bottom')
edge1 = piece1(end, :); % 下边缘
edge2 = piece2(1, :); % 上边缘
end
diff = sum(abs(edge1 - edge2), 'all');
end
% 计算整个排列的总不吻合度
function cost = total_cost(X, pieces)
[rows, cols] = size(X);
cost = 0;
% 水平相邻计算
for i = 1:rows
for j = 1:cols-1
idx1 = X(i, j);
idx2 = X(i, j+1);
cost = cost + calc_diff(pieces(:,:,idx1), pieces(:,:,idx2), 'left_right');
end
end
% 垂直相邻计算
for j = 1:cols
for i = 1:rows-1
idx1 = X(i, j);
idx2 = X(i+1, j);
cost = cost + calc_diff(pieces(:,:,idx1), pieces(:,:,idx2), 'top_bottom');
end
end
end
% 图像拼接函数
function full_image = assemble_image(solution, pieces)
[rows, cols] = size(solution);
[h, w] = size(pieces(:,:,1));
% 创建空白画布
full_image = zeros(rows * h, cols * w, 'uint8');
% 按排列矩阵填充碎片
for i = 1:rows
for j = 1:cols
% 计算当前碎片在画布中的位置
row_start = (i-1)*h + 1;
row_end = i*h;
col_start = (j-1)*w + 1;
col_end = j*w;
% 放置碎片
full_image(row_start:row_end, col_start:col_end) = ...
pieces(:, :, solution(i, j));
end
end
end
最新发布
扫一扫
5320
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
