1. 项目概述重温2013年数学建模B题第一问十多年前的数学建模赛题现在看起来可能有些“复古”但其中蕴含的建模思想、数据处理方法和编程技巧至今依然极具价值。2013年高教社杯全国大学生数学建模竞赛B题“碎纸片的拼接复原”其第一问是一个经典的图像处理和优化问题。题目给出了来自同一页印刷文字文件的碎纸片图像仅纵切要求我们建立模型和算法将这些碎片重新拼接复原。这听起来像是一个有趣的拼图游戏但用Matlab来实现就需要我们将感性的图像匹配转化为严谨的数学计算和逻辑判断。当年参赛时我和队友们为这个问题绞尽脑汁试过多种方案。如今回顾结合更成熟的工具和理解我们可以更优雅、更系统地解决它。这个问题的核心在于如何量化两张碎纸片边缘的“匹配程度”。对于仅纵切的情况我们只需要关注碎片的左边缘和右边缘。最直观的想法就是比较像素将碎片右边缘的像素列与另一碎片左边缘的像素列进行对比相似度最高的就认为是最可能的邻居。但“相似度”如何定义是简单的像素差值求和还是更高级的相关系数如何处理纸张的纹理、墨迹的深浅、扫描产生的噪声这些都是建模时需要仔细考虑的。本文将带你从头开始用Matlab复现解决这个问题的完整流程。我们不仅会写出可运行的代码更会深入探讨每一步背后的“为什么”为什么选择这种相似度度量为什么预处理要这样做当匹配出现歧义时又该如何决策我会分享当年我们踩过的坑以及后来总结出的更稳健的技巧。无论你是正在学习数学建模的学生还是希望提升Matlab图像处理能力的爱好者这篇详尽的“复盘”都能提供一条清晰的路径和许多可直接复用的代码片段。2. 问题拆解与核心思路设计面对一叠扫描出来的、顺序混乱的纵切碎纸片图像我们的目标是恢复其原始顺序。这本质上是一个排序问题但排序的依据不是数字的大小而是碎片边缘内容的连续性。我们需要设计一个算法自动找出哪张碎片应该放在另一张碎片的右边。2.1 核心思路基于边缘匹配的贪婪算法最主流且有效的思路是基于边缘匹配的贪婪算法。其核心步骤可以分解为数据读取与预处理将所有碎片图像读入Matlab统一尺寸并可能进行二值化、去噪等操作将图像数据转化为便于计算的矩阵。特征提取针对每一张碎片我们关心它的两个竖直边缘最左边的一列像素和最右边的一列像素。这些列向量就是后续匹配的“特征”。相似度计算定义一个函数用于计算任意两张碎片之间“左边缘-右边缘”的匹配程度。这个函数的值越大或越小表示这两张碎片是相邻的可能性越高。全局排序与拼接从某一张碎片开始比如最左端的碎片其左边缘应为全白根据相似度计算找到它最匹配的右邻居然后将该邻居作为新的当前碎片继续寻找下一个右邻居直到所有碎片被串联起来。这就是“贪婪”的由来每一步都只选择当前最优的匹配。这个思路清晰直接但难点和精髓全在细节里如何定义“相似度”才能抗噪声、效果好如何确定起点如果贪婪算法中途匹配错误导致“死胡同”怎么办2.2 相似度度量的选择与比较这是整个模型的心脏。我们至少需要比较一个列向量碎片A的右边缘与另一个列向量碎片B的左边缘。假设我们处理的是二值图像黑色文字为1白色背景为0列向量是0和1的序列。方法一绝对差值和Sum of Absolute Differences, SAD这是最直观的方法。将两个列向量逐像素相减取绝对值然后求和。值越小说明两个边缘越相似。相似度_SAD sum(abs(edge_right_A - edge_left_B))优点计算简单快速物理意义明确。缺点对噪声敏感。如果边缘有少量噪点会导致差值和不稳定。更重要的是它没有考虑整体的相关性。方法二归一化互相关Normalized Cross-Correlation, NCC这种方法在图像匹配中非常常用。它计算两个向量的相关系数值域在[-1, 1]之间。1表示完全正相关-1表示完全负相关0表示不相关。相似度_NCC corrcoef(edge_right_A, edge_left_B); similarity 相似度_NCC(1,2)优点对图像的亮度线性变化不敏感抗噪能力比SAD强能更好地捕捉形状的相似性。缺点计算量稍大于SAD。方法三基于边界黑白像素比例的方法考虑到这是文字文档我们可以利用先验知识文字行的边界其黑白像素的分布应该具有连续性。例如我们可以计算右边缘列中黑色像素值为1的比例与左边缘列中黑色像素的比例进行比较。但单纯比较比例丢失了空间位置信息效果通常不如前两种。实操心得在2013年我们实际解题和后续的多次复现中归一化互相关NCC是稳健性和效果最好的选择。它能够有效应对扫描图像常见的亮度不均、轻微污渍等问题。SAD虽然快但在遇到复杂情况时容易出错。因此本文后续将主要基于NCC来构建模型。2.3 确定拼接起点与终点在纵切情况下原文档最左侧碎片的左边缘和最右侧碎片的右边缘理论上应该是全白的即像素值全部为背景色。我们可以利用这一特性。寻找最左碎片遍历所有碎片的左边缘列计算该列中非背景色如黑色像素的数量。数量最少 ideally 为0的那个碎片极有可能是最左边的碎片。我们可以将它作为拼接的起点。寻找最右碎片同理遍历所有碎片的右边缘列找到非背景色像素数量最少的碎片作为拼接的终点或验证终点。验证有时可能因为噪声或装订线痕迹没有绝对全白的边缘。我们可以设定一个阈值比如非白色像素数少于总像素的1%或者同时结合左、右边缘的信息来综合判断。3. 完整Matlab实现流程与代码详解下面我们按照操作流程一步步实现算法。假设我们已经将19张纵切碎片图像pic1.png,pic2.png, ...,pic19.png保存在名为fragments的文件夹中。3.1 环境准备与数据读取首先我们需要将图像读入Matlab并存储为统一的数据结构。clear; clc; close all; % 1. 设置路径和参数 fragment_dir ./fragments/; % 碎片图像所在文件夹 file_list dir(fullfile(fragment_dir, *.png)); % 获取所有png文件 num_fragments length(file_list); % 碎片总数应为19 % 2. 初始化存储结构 fragments cell(num_fragments, 1); % 用元胞数组存储每个碎片的图像矩阵 fragment_names cell(num_fragments, 1); % 3. 循环读取图像 for i 1:num_fragments file_path fullfile(fragment_dir, file_list(i).name); img imread(file_path); % 读入图像得到三维矩阵 (高度, 宽度, 3[RGB]) fragments{i} img; % 存储原始彩色图像后续根据需要处理 fragment_names{i} file_list(i).name; % 可选显示一下看看 % subplot(4,5,i); imshow(img); title(sprintf(Frag %d, i)); end disp([成功读取 , num2str(num_fragments), 张碎片图像。]);3.2 图像预处理灰度化与二值化彩色图像RGB包含三个通道计算量大且对于黑白文档不必要。我们首先将其转换为灰度图像然后二值化将文字和背景分离。% 4. 预处理转换为灰度图并二值化 bin_fragments cell(num_fragments, 1); % 存储二值化图像 gray_fragments cell(num_fragments, 1); left_edges cell(num_fragments, 1); % 存储左边缘列向量 right_edges cell(num_fragments, 1); % 存储右边缘列向量 for i 1:num_fragments % 4.1 灰度化 gray_img rgb2gray(fragments{i}); gray_fragments{i} gray_img; % 4.2 二值化 - 这是关键步骤 % 使用全局阈值 Otsus method 自动确定阈值 level graythresh(gray_img); % level 是一个介于0-1之间的归一化阈值 bin_img imbinarize(gray_img, level); % 二值化背景为0(黑)文字为1(白)注意 % 通常imshow显示时0是黑1是白。但为了计算方便我们通常希望背景是1白文字是0黑。 % 所以我们需要反转图像。 bin_img ~bin_img; % 逻辑非操作反转图像。现在背景1(白)文字0(黑)。 bin_fragments{i} bin_img; % 4.3 提取边缘特征 [height, width] size(bin_img); left_edges{i} bin_img(:, 1); % 最左边一列尺寸为 [height x 1] right_edges{i} bin_img(:, end); % 最右边一列尺寸为 [height x 1] % 可选可视化检查预处理效果 % figure; % subplot(1,3,1); imshow(fragments{i}); title(原始彩色); % subplot(1,3,2); imshow(gray_img); title(灰度); % subplot(1,3,3); imshow(bin_img); title(二值化(背景白文字黑)); % pause(0.5); end注意事项二值化是预处理的核心阈值选择不当会导致文字断裂或背景噪点过多。graythresh函数使用的Otsu方法在大多数文档图像上效果很好。如果效果不佳可以考虑手动调整阈值或使用自适应阈值法imbinarize(gray_img, ‘adaptive’)。图像反转~bin_img是为了让后续计算更符合直觉匹配时我们希望相同的部分都是背景或都是文字贡献正相关。3.3 构建相似度矩阵接下来我们需要计算任意两个碎片i作为左碎片和j作为右碎片之间的匹配度。我们使用归一化互相关NCC。% 5. 计算归一化互相关NCC相似度矩阵 % 矩阵 S 的维度为 [num_fragments x num_fragments] % S(i, j) 表示将碎片 i 放在左边碎片 j 放在右边时它们边缘的匹配程度。 S zeros(num_fragments, num_fragments); for i 1:num_fragments right_edge_i double(right_edges{i}); % 将逻辑型转换为双精度浮点型便于corrcoef计算 for j 1:num_fragments if i j S(i, j) -inf; % 自己和自己不匹配设为负无穷 continue; end left_edge_j double(left_edges{j}); % 计算相关系数矩阵 corr_matrix corrcoef(right_edge_i, left_edge_j); % corrcoef返回一个2x2矩阵[1, r; r, 1]。我们需要的是r即(1,2)或(2,1)元素。 ncc_value corr_matrix(1, 2); % NCC值在[-1,1]之间。越接近1匹配越好。 % 为了后续排序方便我们直接存储NCC值。 S(i, j) ncc_value; end end % 可视化相似度矩阵热图 figure; imagesc(S); colorbar; title(碎片间NCC相似度矩阵 (行i右边缘 - 列j左边缘)); xlabel(潜在右邻居碎片 j); ylabel(当前碎片 i);这个矩阵S是我们的“地图”。S(i,j)的值越高说明碎片i的右边缘和碎片j的左边缘越相似j越有可能是i的右邻居。3.4 确定起点与贪婪拼接现在我们利用相似度矩阵S和“最左碎片左边缘应全白”的先验知识开始拼接。% 6. 确定起点寻找最可能是最左边的碎片左边缘最“白” % 在我们的二值图像中背景白 1文字黑 0。 % 因此左边缘列向量中所有像素值越接近1说明越白。 left_whiteness zeros(num_fragments, 1); for i 1:num_fragments left_whiteness(i) mean(left_edges{i}); % 计算左边缘列的平均值。背景全白1全黑0。 end [~, start_idx] max(left_whiteness); % 平均值最大的即最白的作为起点。 disp([推测的起始碎片编号: , num2str(start_idx), (文件名: , fragment_names{start_idx}, )]); % 7. 贪婪算法进行拼接 unused true(num_fragments, 1); % 标记碎片是否已被使用 current_idx start_idx; unused(current_idx) false; % 标记起点已使用 order current_idx; % 记录拼接顺序 for step 1:num_fragments-1 % 总共需要找 n-1 次邻居 % 从相似度矩阵 S 的当前行中找出未使用碎片里相似度最高的那个 candidate_scores S(current_idx, :); candidate_scores(~unused) -inf; % 将已使用碎片的得分设为负无穷排除它们 [best_score, next_idx] max(candidate_scores); % 判断匹配质量。如果最佳匹配得分太低可能意味着匹配错误或已到终点。 if best_score 0.3 % 这是一个经验阈值可以根据实际情况调整 warning(在步骤 %d当前碎片 %d 的最佳匹配得分(%.3f)过低。可能已到达最右端或出现错误。, ... step, current_idx, best_score); break; end % 更新状态 current_idx next_idx; unused(current_idx) false; order [order, current_idx]; % 将新碎片编号加入顺序列表 fprintf(步骤 %2d: 碎片 %2d - 碎片 %2d (匹配得分: %.4f)\n, ... step, order(end-1), order(end), best_score); end % 检查是否所有碎片都已使用 if sum(unused) 0 disp(恭喜所有碎片拼接完成。); else disp([警告仍有 , num2str(sum(unused)), 张碎片未被使用。拼接可能不完整。]); unused_indices find(unused); disp(未使用的碎片编号:); disp(unused_indices); end disp(最终的碎片顺序 (从左到右):); disp(order);3.5 可视化拼接结果得到顺序后我们需要将按照这个顺序排列的碎片图像横向拼接起来直观地查看复原效果。% 8. 根据顺序可视化拼接结果 % 获取单张碎片图像的尺寸所有碎片尺寸应相同 sample_bin bin_fragments{1}; [img_height, img_width] size(sample_bin); % 创建一张空白的大画布用于横向拼接 % 宽度 单张碎片宽度 * 碎片数量 % 高度 单张碎片高度 reconstructed_img zeros(img_height, img_width * length(order), logical); % 使用逻辑矩阵存储二值图像 for i 1:length(order) frag_idx order(i); % 计算当前碎片在大图中的列范围 col_start (i-1) * img_width 1; col_end i * img_width; % 将碎片图像贴到对应位置 reconstructed_img(:, col_start:col_end) bin_fragments{frag_idx}; end % 显示拼接结果 figure; imshow(reconstructed_img); title(基于贪婪算法(NCC)的拼接复原结果); xlabel(像素位置); ylabel(像素位置); % 为了更清晰我们也可以显示灰度图拼接结果更接近原貌 reconstructed_gray zeros(img_height, img_width * length(order), uint8); for i 1:length(order) frag_idx order(i); col_start (i-1) * img_width 1; col_end i * img_width; reconstructed_gray(:, col_start:col_end) gray_fragments{frag_idx}; end figure; imshow(reconstructed_gray); title(基于贪婪算法(NCC)的拼接复原结果 (灰度图));4. 算法优化与鲁棒性提升上述基础贪婪算法在理想情况下碎片清晰、切割整齐、无噪声可以工作得很好。但实际扫描图像总会存在各种问题直接使用上述算法可能会在某个步骤做出错误匹配导致后续全盘皆错。我们需要增加算法的鲁棒性。4.1 处理匹配歧义多候选与回溯基础贪婪算法只选择“当前最优”这是一种“目光短浅”的策略。当最佳匹配和第二佳匹配得分非常接近时盲目选择最佳的可能不是全局最优解。优化策略K-最佳候选与回溯记录多候选在每一步寻找邻居时不只看第一名而是保留前K个例如K3得分最高的候选碎片。深度优先搜索与回溯尝试沿着每一条候选路径进行拼接。如果某条路径走到后面发现无法继续所有候选匹配得分都很低或者用完了所有碎片但顺序不合理则回溯到上一个决策点尝试下一个候选。评估完整路径为每一条完整的拼接路径定义一个“总得分”例如路径上所有匹配对的NCC得分之和。选择总得分最高的路径作为最终结果。这种方法计算量会增大K越大计算量指数增长但对于解决歧义、提高准确率非常有效。对于19个碎片K2或3的回溯搜索是可行的。% 这是一个简化的多候选回溯思路框架并非完整可执行代码用于说明逻辑。 function best_order backtrack_reconstruction(S, start_idx, K) num_frags size(S, 1); best_order []; best_total_score -inf; % 定义一个递归搜索函数 function dfs(current_order, used, current_total_score) if length(current_order) num_frags % 找到一条完整路径 if current_total_score best_total_score best_total_score current_total_score; best_order current_order; end return; end current_idx current_order(end); % 获取当前碎片的所有可能右邻居得分并排除已使用的 scores S(current_idx, :); scores(used) -inf; % 找出前K个最佳候选 [sorted_scores, sorted_indices] sort(scores, descend); topK_indices sorted_indices(1:min(K, sum(~isinf(sorted_scores)))); topK_scores sorted_scores(1:min(K, sum(~isinf(sorted_scores)))); for k 1:length(topK_indices) next_idx topK_indices(k); next_score topK_scores(k); if next_score 0.2 % 如果候选得分太低剪枝不再探索这条路径 continue; end new_used used; new_used(next_idx) true; dfs([current_order, next_idx], new_used, current_total_score next_score); end end initial_used false(1, num_frags); initial_used(start_idx) true; dfs([start_idx], initial_used, 0); end4.2 利用双边缘信息进行全局优化贪婪算法只用了“右边缘匹配左边缘”的单向信息。实际上我们可以利用更全局的信息。我们可以将问题转化为一个**旅行商问题TSP**的变种城市每个碎片是一个“城市”。距离从城市i到城市j的“距离”定义为1 - NCC(i, j)因为NCC越大表示越相似距离应越小。目标找到一条访问所有城市一次且仅一次的路径使得总距离最小。这条路径就是碎片的顺序。这样我们就将局部匹配问题转化为了一个全局优化问题。可以使用TSP的经典算法来求解如动态规划对于19个碎片状态空间是18!仍然巨大但比暴力枚举好、模拟退火、遗传算法等。Matlab的全局优化工具箱可以提供帮助。% 思路示例将相似度矩阵转换为距离矩阵并使用近似算法如最近邻法求解TSP % 注意这是一个启发式方法不一定得到最优解但通常比基础贪婪法更稳健。 D 1 - S; % 将相似度转换为相异度距离 for i 1:num_fragments D(i,i) inf; % 自己到自己的距离设为无穷大 end % 使用最近邻法Nearest Neighbor作为TSP的启发式解法 tsp_order zeros(1, num_fragments); tsp_order(1) start_idx; visited false(1, num_fragments); visited(start_idx) true; for i 2:num_fragments last_city tsp_order(i-1); [~, next_city] min(D(last_city, ~visited)); % 在未访问城市中找距离最小的 visited(next_city) true; tsp_order(i) next_city; end disp(基于TSP最近邻法得到的顺序:); disp(tsp_order); % 然后可以比较tsp_order和基础贪婪法的order选择总距离更小的那个。4.3 预处理增强边界扩展与噪声过滤有时碎片边缘在扫描时会有少量像素的缺失或污渍。我们可以考虑在提取边缘特征时不只取最边上一列而是取边缘的若干列例如2-3列来计算一个更宽泛的“边缘区域”的特征比如计算这几列像素的平均向量再进行匹配。这相当于给匹配增加了一个小的容错空间。另外在二值化后可以使用形态学操作如bwareaopen去除面积过小的孤立噪点避免这些噪点对边缘匹配产生干扰。% 示例提取多列作为边缘特征并去噪 width_of_edge_region 3; % 使用边缘的3列像素 for i 1:num_fragments bin_img bin_fragments{i}; % 去噪移除面积小于50像素的连通区域 bin_img_cleaned bwareaopen(bin_img, 50); [height, width] size(bin_img_cleaned); % 提取左边缘区域前3列的均值向量 left_region bin_img_cleaned(:, 1:min(width_of_edge_region, width)); left_edges{i} mean(left_region, 2); % 对行求平均得到 height x 1 的列向量 % 提取右边缘区域后3列的均值向量 right_region bin_img_cleaned(:, max(1, width-width_of_edge_region1):end); right_edges{i} mean(right_region, 2); end % 注意使用均值向量后特征值变为0~1之间的浮点数计算NCC仍然适用。5. 常见问题排查与调试技巧在实际运行代码时你可能会遇到各种问题。以下是一些常见情况及解决方法。5.1 拼接结果出现明显错位或乱码症状复原后的文字行在碎片交界处断开、错位无法阅读。可能原因与排查二值化效果差这是最常见的原因。检查bin_fragments中的图像文字是否清晰连贯背景是否干净尝试调整二值化阈值。可以将level graythresh(gray_img)改为手动设置例如level 0.5或使用自适应二值化bin_img imbinarize(gray_img, ‘adaptive’);。起点选择错误算法从一个错误的碎片开始拼接。检查left_whiteness的值确认被选为起点的碎片其左边缘是否确实最白。可以可视化显示所有碎片的左边缘列figure; for i1:19; subplot(4,5,i); plot(left_edges{i}); title(i); end。真正的左边缘应该是一条接近1的水平线。相似度度量失效如果图像噪声极大或文字非常稀疏NCC可能不稳定。可以尝试换用其他度量如余弦相似度cosine_sim dot(a,b)/(norm(a)*norm(b))或者结合SAD和NCC。贪婪算法陷入局部最优这就是为什么需要4.1和4.2中的优化策略。尝试运行回溯算法或TSP方法看结果是否改善。5.2 算法未能使用所有碎片症状最终order中的碎片数量少于19控制台提示有碎片未使用。可能原因与排查匹配阈值过高在贪婪算法的循环中我们设置了if best_score 0.3的跳出条件。如果某个碎片的真实邻居匹配得分恰好低于0.3可能因为边缘损坏算法会提前终止。尝试降低这个阈值到0.2或0.15或者直接注释掉这个break语句让算法强制遍历但可能接错。存在重复或极度相似的碎片本题应不存在如果相似度矩阵中某个碎片与其他多个碎片的匹配得分都极高且接近可能导致算法混乱。检查相似度矩阵S看是否有异常高的非对角线元素。回溯或TSP方法使用4.1或4.2的全局方法它们的目标是使用所有碎片通常能避免此问题。5.3 运行速度慢特别是对于回溯算法症状当碎片数量增多如本题附件2、3有更多碎片或K值设得较大时回溯算法运行时间很长。优化策略剪枝在回溯搜索中如果当前路径的局部匹配得分已经很低或者当前部分路径的总得分远低于已知的较优解可以提前终止该分支的搜索如上文代码中的if next_score 0.2就是一种剪枝。降低K值对于大多数情况保留前2个候选K2已经能很大程度改善贪婪法的不足同时计算量可控。使用启发式TSP求解器如模拟退火、遗传算法等它们能在可接受时间内为大规模TSP问题找到近似最优解。Matlab的Global Optimization Toolbox提供了simulannealbnd和ga函数。5.4 如何评估拼接结果的正确性对于本题由于我们有原文档的参考尽管在解题时不知道可以人工目视检查复原后的文字是否可读、行文是否连贯。在无法目视判断时可以定义一些量化指标平均匹配得分最终顺序中相邻碎片对的NCC得分的平均值。越高越好。平滑度指标计算拼接后在接缝处相邻列的像素差异总和。越低越好。文字行连续性高级使用OCR工具识别拼接前后的文字计算编辑距离或语言模型的困惑度perplexity连续性好的文本这些指标会更优。最后分享一个我调试时常用的小技巧可视化每一步的匹配对。在贪婪算法循环中不仅打印文本日志还可以实时显示当前碎片和它最佳匹配碎片的边缘像素对比图直观感受匹配质量。% 在贪婪算法循环内部添加可视化调试代码 if best_score 0.3 % 只在匹配质量尚可时显示 figure(100); clf; subplot(2,1,1); plot(right_edges{current_idx}, b-); hold on; plot(left_edges{next_idx}, r-); legend(当前碎片右边缘, 候选碎片左边缘); title(sprintf(匹配对比: 碎片%d - 碎片%d (得分: %.3f), current_idx, next_idx, best_score)); subplot(2,1,2); imshow([bin_fragments{current_idx}, bin_fragments{next_idx}]); title(二值图像拼接预览); drawnow; pause(0.5); % 暂停半秒观察 end通过这种“人机交互”式的调试你能迅速定位是哪个匹配环节出了问题从而有针对性地调整预处理、相似度度量或算法逻辑。数学建模的魅力就在于将一个现实问题抽象为数学模型和算法再通过编程实现和反复调试最终让计算机替你完成复杂的识别与推理。2013年B题第一问正是这样一个经典的入门案例。
