简介面向图像处理、计算机视觉方向的学生与工程师这份MATLAB代码以系统开源方式复现了Pérez等人的经典泊松图像编辑算法用于消除图像拼接与融合中的边界痕迹并支持局部区域的无缝修复。压缩包整体约213KB保持轻量内含Poisson Fusion与Poisson Repair两个核心脚本分别对应图像融合和图像修复两种典型任务另有若干示例图片便于直接运行并验证算法效果。代码按论文思路组织涉及泊松方程构造、边界条件设置、稀疏线性系统求解等关键环节注释与步骤清晰能帮助读者将数学公式与MATLAB实现互相印证在此基础上还能学习结构化矩阵组装与数值求解的编程技巧。已有711人学习下载适合具备基础MATLAB操作和图像处理概念的学习者也可作为课程设计、论文复现或科研入门的起始模板。1. 泊松融合是什么拖一块像素过去让接缝消失做图像合成时最头疼的事就是边缘。直接把源图的一块区域裁下来贴到目标图上无论选点选得多准放大看总有一圈生硬的轮廓。泊松融合Poisson Fusion的思路完全不同它不搬像素搬梯度。先提取源图区域的梯度场再把它作为约束条件去求解目标区域内的泊松方程最终重建出像素值。这个过程会让合成区域与周围环境在梯度域上保持一致接缝在视觉上彻底消失。这套 MATLAB 实现把整个管线拆成了可单独调用的函数从梯度计算、散度构建到稀疏矩阵求解每一步都有对应代码。适合两类人一是课程作业或论文复现需要无缝克隆、物体移除效果的学生二是想理解离散泊松方程如何落地到图像网格上、准备把这套逻辑移植到自研引擎里的开发者。2. 工程结构与数据流源码包拿到手先看懂五个文件的职责2.1 目录结构不是一坨代码而是按处理阶段分好的模块下载解压后典型的结构是这样的。核心逻辑全部在 src 目录下demo 数据单独放一份方便直接运行验证效果。MATLAB_Poisson_Fusion_HJT/ ├── main_demo.m # 入口脚本读图、设 mask、调融合、保存结果 ├── src/ │ ├── poisson_fusion.m # 融合主函数串联整个求解流程 │ ├── get_gradient.m # 用中心差分计算源图的 x/y 方向梯度 │ ├── build_A.m # 构建稀疏系数矩阵 A │ ├── build_b.m # 构建右侧向量 b合并梯度散度和边界条件 │ └── solve_poisson.m # 调 MATLAB 反斜杠求解器解 Ax b ├── data/ │ ├── source.png # 源图要抠出来的前景所在图 │ ├── target.png # 目标图合成背景 │ └── mask.png # 二值掩码白色区域为融合范围 └── README.md # 运行说明和参数解释拿到手的第一件事不是急着跑而是先确认 MATLAB 版本。这个实现用到了sparse、diff、imfilter和\求解器全是 MATLAB 基础功能不需要额外工具箱。理论上 R2016a 之后的版本都能跑我实际在 R2020b 和 R2023a 上都验证过没遇到兼容性问题。唯一要注意的是imfilter的边界处理参数老版本和新版本的默认行为有细微差别后面避坑章节会专门说。2.2 数据流动向从源图梯度到合成像素的六步管道整个融合过程可以拆成一条清晰的流水线。第一步读取源图、目标图和 mask做类型归一化所有图像统一转成 double 范围在 [0,1] 的矩阵。第二步对源图计算梯度场水平方向用[1, -1]滤波器垂直方向用[1; -1]得到两个梯度分量。第三步对梯度场求散度用相邻像素梯度差值的和来近似。第四步根据 mask 的大小构建稀疏矩阵 A内部区域节点对应拉普拉斯模板边界节点对应单位矩阵行。第五步构建右侧向量 b内部节点取散度值边界节点取目标图在该位置的像素值。第六步求解线性方程组把得到的解加回到目标图的基线灰度上得到最终融合结果。这六步的划分是有讲究的。它把「数据准备」和「求解」彻底分离get_gradient只负责梯度build_A只负责系数矩阵build_b只负责右端项。这样一旦你发现结果不对可以单独检查每一步的输出而不需要把整个流程重新捋一遍。我调这个项目时就是这样定位问题的——先看梯度图是否合理再看散度是否有异常尖峰最后才看求解结果。2.3 与常见图像编辑工具的区别为什么不用 photoshop 的混合选项可能你会问PS 里做个图层混合不也能过渡吗确实能但那是像素域的透明度或混合模式操作属于经验性的近似。泊松融合做的是梯度域的重构它在数学上保证融合区域内部的梯度场与源图梯度场尽可能一致同时边界像素值严格等于目标图原有值。这带来一个直接好处融合后的光照变化是连续的不会出现局部过曝或偏灰的色块。渐变过渡是模拟的梯度重构是本构的这是本质差异。这个仓库的价值就在于把后者抽象成了可复现的 MATLAB 代码。3. 核心模块拆解从梯度场到稀疏矩阵求解的完整链路3.1 梯度计算为什么用中心差分而不是 sobel泊松方程的核心输入是源图梯度。常见做法是用中心差分而非 sobel 算子原因是泊松方程在离散化时天然建立在相邻像素差值之上中心差分算子恰好与五邻点拉普拉斯模板形成一组自洽的离散对偶。代码实现通常长这样function [gx, gy] get_gradient(img) % 输入 img 为 double 型二维灰度矩阵 h_gauss fspecial(gaussian, [3, 3], 0.8); % 先做轻微平滑防噪 img_sm imfilter(img, h_gauss, replicate); gx diff(img_sm, 1, 2); % 水平梯度sizeM x (N-1) gy diff(img_sm, 1, 1); % 垂直梯度size(M-1) x N % 对齐尺寸最后一列/行补零方便后续散度计算 gx(:, end1) 0; gy(end1, :) 0; end这里有一个值得注意的设计diff的输出比原图少一列或一行所以代码里做了对齐补零。gx和gy的尺寸最后都和原图一致这样在计算散度时不需要处理边界特例。高斯平滑的核大小取 3×3sigma 取 0.8是抗噪和保细节之间的折中。如果 sigma 太大源图的边缘会被过度模糊融合结果的细节发闷如果完全不平滑纯色区域的传感器噪声会被当成真实梯度导致融合区域出现颗粒感。3.2 系数矩阵构建五邻点拉普拉斯模板的离散化泊松方程离散化之后是一个标准稀疏线性系统。对融合区域内部的每个像素方程形式是4 * f(i,j) - f(i-1,j) - f(i1,j) - f(i,j-1) - f(i,j1) b(i,j)右边 b(i,j) 的取值分为两种情况如果该像素在 mask 内部取源图该位置的梯度散度如果该像素在 mask 边界上则方程变成了f(i,j) target(i,j)直接把目标图的像素值作为硬约束。function A build_A(mask) [rows, cols] size(mask); N rows * cols; idx reshape(1:N, rows, cols); % 给每个像素一个全局编号 A sparse(N, N); % 务必用 sparse全矩阵会内存爆炸 for i 1:rows for j 1:cols id idx(i, j); if mask(i, j) A(id, id) 4; if i 1, A(id, idx(i-1, j)) -1; end if i rows, A(id, idx(i1, j)) -1; end if j 1, A(id, idx(i, j-1)) -1; end if j cols, A(id, idx(i, j1)) -1; end else A(id, id) 1; % 非融合区域直接约束为原值 end end end end按这个写法A 是一个 N×N 的方阵N 等于整张图的像素总数。对融合区域以外的像素系数矩阵置为单位行方程退化为f(i,j) b(i,j)。这种「全图统一编号 单位行填充」的策略非常实用因为 MATLAB 的\求解器天然要求所有变量都参与求解你要么把非融合区域也纳入方程组要么在求解后手动合并。纳入方程组虽然让稀疏矩阵的行数多了不少但省去了手工合并索引的步骤代码更简洁不容易出错。唯一的前提是你必须用sparse存储全图的 N 如果是 1024×768大约是 78 万个未知数全矩阵会直接卡死。顺带一提如果你要追求反查求解顺序的直观性把 A 中非零元的行、列号打印出来能看到清晰的对角带结构。3.3 右端向量构建边界像素决定结果下限内部像素提供梯度趋势右侧向量 b 的构建比 A 更微妙因为处理「边界像素」的标准在不同实现里并不统一。核心原则是融合区域内的非边界像素其 b 值取源图梯度的散度融合区域内的边界像素其 b 值取目标图在该坐标的原始像素值。这里的关键点是「覆盖」优先级——如果某个像素同时满足「在 mask 内」和「在 mask 的边界上」它的方程对应的是边界条件而非散度条件。function b build_b(gx, gy, target, mask) [rows, cols] size(mask); N rows * cols; b zeros(N, 1); idx reshape(1:N, rows, cols); % 计算源图的散度散度 gx(i,j) - gx(i-1,j) gy(i,j) - gy(i,j-1) div zeros(rows, cols); div(2:end, :) gy(2:end, :) - gy(1:end-1, :); div(:, 2:end) div(:, 2:end) gx(:, 2:end) - gx(:, 1:end-1); % 判断 mask 内部像素是否在边界上 se strel(disk, 1); mask_eroded imerode(mask, se); interior mask mask_eroded; boundary mask ~interior; for i 1:rows for j 1:cols id idx(i, j); if boundary(i, j) b(id) target(i, j); % 边界处直接取目标像素值 elseif interior(i, j) b(id) div(i, j); % 内部取源图梯度散度 else b(id) target(i, j); % mask 外也取目标值 end end end end这里「边界」的判定用了imerode——把 mask 向内腐蚀一个像素后mask 内但不属于腐蚀结果的就是边界环。这种判定方式比逐像素判断上下左右是否越界更稳健也更能保证边界是闭合连通的一圈。如果不用腐蚀而改用邻域判断边界像素可能落在非连通区结果会在某个角落出现撕裂感。实际调试中你可以把interior和boundary各画成一张二值图肉眼确认边界是否闭合这一步能省掉很多后端的折腾。3.4 主函数串联求解、重组、像素值恢复function result poisson_fusion(source, target, mask) % 统一到 double [0,1]避免 uint8 溢出问题 src im2double(source); tgt im2double(target); msk logical(mask); % 如果 mask 尺寸比图像小自动扩展 if ~isequal(size(msk), size(src(:,:,1))) msk imresize(msk, [size(src,1), size(src,2)], nearest); end % 逐通道处理 result zeros(size(src)); for c 1:size(src, 3) [gx, gy] get_gradient(src(:, :, c)); A build_A(msk); b build_b(gx, gy, tgt(:, :, c), msk); x A \ b; result(:, :, c) reshape(x, size(msk)); end end主函数的几个细节值得解释一下。首先是逐通道处理泊松方程本身是在灰度域上定义的彩色图拆成 R、G、B 三个通道分别求解再拼回来这是最常见也是最稳定的做法。想在 Lab 空间做可以只对 L 通道用泊松融合ab 通道直接拷贝源图区域效果在抗色偏上更好但那是二次开发的事。其次是A \ b这一步MATLAB 的稀疏求解器会自动选择方法对五邻点拉普拉斯矩阵这种对称正定结构走的是 Cholesky 分解路线速度和内存都表现不错。如果你发现求解耗时特别长先去查 A 有没有意外被转成 full 矩阵这是我排过最多的坑。最后是结果直接 reshape 回二维因为数据本身是浮点且在 [0,1] 范围保存前再用imwrite或转成 uint8 即可。4. 运行指南与避坑记录把 demo 跑通再处理三类翻车现场4.1 运行步骤从解压到看到融合图像只要四步第一步把仓库整个目录加入 MATLAB 路径。第二步在 MATLAB 命令行执行main_demo脚本会自动加载data目录里的三张图。第三步观察弹出的两张对比图左侧是直接裁剪拼接的朴素结果右侧是泊松融合结果可以放大检查边缘过渡。第四步程序自动把融合结果保存到result.png。cd MATLAB_Poisson_Fusion_HJT # 在 MATLAB 命令行里执行 # addpath(genpath(pwd)); # main_demo如果你没有现成的测试图也可以用任意两张图自己生成 mask。常见做法是先用imrect交互框选一块矩形区域作为 mask再保存成二值图。但矩形 mask 有个问题矩形的四条边都是硬边界融合区域边缘是全图最显眼的地方泊松融合对矩形边界的处理效果一般。建议用手绘套索或者从网上下载带 alpha 通道的 PNG 前景图来转 mask过渡会自然很多。4.2 避坑一mask 尺寸与图像不一致报错维度不匹配现象直接拿自己准备的 mask 跑poisson_fusionMATLAB 抛出Matrix dimensions must agree或者是下标越界。原因mask 的宽高和源图/目标图不一致build_A里的idx矩阵尺寸对不上索引就错位了。解决在主函数入口强制对齐尺寸按上面的代码用imresize把 mask 缩放到和源图一致。注意插值方法选nearest如果用双线性插值mask 边缘会变成灰色过渡logical转换之后产生锯齿空洞。Here is the corrected text:4.2 避坑二融合结果出现灰白色雾状覆盖现象融合后图像整体发灰像是蒙了一层半透明的雾气。原因数据类型没统一。如果源图是 uint8计算梯度时diff会得到负数uint8 会截断为 0梯度场惨不忍睹。解决一进入poisson_fusion就把图像转成 double 并归一化到 [0,1]全程用浮点运算。另一个被忽略的点是get_gradient中的高斯滤波如果输入图是 uint8imfilter的输出会被四舍五入到整数值散度会失真。跑完再转回 uint8 存储即可。4.3 避坑三稀疏矩阵被隐式转成满矩阵程序卡死或爆内存现象build_A执行完毕A 的维度看着正常但求解时内存占用飙升耗时暴涨甚至直接卡死。原因某处代码对 A 做了非稀疏操作——包括A(id, id)4这种单元素赋值如果 A 之前是 sparse单元素赋值会触发重排但仍然保留稀疏结构。真正的坑是调试时把 A 打印到命令行或者用full(A)检查这会把整个矩阵展开。解决永远不要对 A 用full查看结构用spy(A)查看单行元素用A(:, idx)会返回稀疏列也尽量避免。判断是否稀疏用issparse(A)。4.4 避坑四融合后边缘出现细小黑色或白色描边现象融合区域的轮廓线特别清晰像被人用画笔勾了一圈。原因边界条件处理不当。build_b中 boundary 判定用了imerode如果 mask 本身有孤立的小噪点腐蚀后会被拆分成细颈甚至断环边界像素的索引错位泊松方程的边界约束失效。解决在调用poisson_fusion之前先对 mask 做形态学闭运算再填充孔洞。常见做法是先imclose(mask, strel(disk, 3))再imfill(mask, holes)确保 mask 是干净闭合的一块连续区域。另一个可能原因是源图和目标图尺寸不一致导致求解结果在对齐边界时产生一圈残差。5. 二次开发把泊松融合改造成无缝克隆、纹理压平与局部调色5.1 改造一无缝克隆到任意位置通过平移变换实现最直接的改造是把融合区域从「原坐标」改为「任意坐标」也就是把源图前景抠出来贴到目标图任意位置。这里的思路很简单把 mask 平移到一个新位置再把源图的梯度场也平移到对应位置。需要修改的是主函数里梯度计算的区域索引而不是整个矩阵结构。具体做法是先对 mask 做平移变换imtranslate然后重算源图在新坐标区域的梯度再调用融合主流程。要注意平移后的 mask 边缘会引入插值误差所以平移后最好再做一次imbinarize来消除灰色边缘。验证方法也很直观在目标图上画一个矩形框把源图前景水平移动 50 像素后融合观察边缘是否仍然无缝衔接。梯度域操作的好处是平移不会改变梯度的相对关系融合后的光照仍然均匀。5.2 改造二纹理压平抑制细节突出大面积色块做皮肤磨皮或背景纹理清理时可以把源图梯度乘一个小于 1 的系数 α让融合区域内部的细节梯度减弱得到一个相对平坦但颜色连续的填充区域。α 取 0.2 时细节保留很少α 取 0.8 时细节基本保留。修改点只在一处build_b里计算散度时把div乘上alpha其余逻辑不用动。这个技巧也叫梯度缩放在图像修复inpainting场景中很实用比直接用高斯模糊糊上去要自然因为它保留了低频的颜色渐变趋势只是压平了高频噪声。5.3 改造三局部调色把源图颜色趋势注入目标图还有些场景做的是颜色迁移而不是结构融合。此时只需要修改梯度场的来源——把源图的梯度场替换为目标图本身的梯度场但把 b 的散度项改为源图与目标图梯度的混合加权。常见做法是div_mixed alpha * div_source (1 - alpha) * div_targetalpha 控制调色强度。这样一来融合区域的边缘仍然与目标图吻合内部色调整体向源图偏移。这是我个人最喜欢的一个改法因为只需要在build_b里加一行混合运算五分钟就能实验出不同效果。5.4 我坚持的验证习惯每次改完代码我都要做两个验证第一跑通 demo 看结果这一步只保证链路通第二把求解向量 x 从reshape的结果里提取出来手动计算一遍目标区域的残差norm(A*x - b)如果残差量级在 1e-10 以上说明 A 或 b 某处拼错了。从那以后我每次改涉及系数矩阵的代码都要先跑一遍残差校验再去看图否则边缘的怪象会浪费半天时间去猜。希望帮到你这份 MATLAB 泊松融合仓库是个极佳的起点跑通它再逐步扩展自己的图像编辑算法路径会顺畅很多。本文还有配套的精品资源点击获取