简介面向图像融合与泊松方程学习者的MATLAB复现资源对应Pérez等人2003年经典论文《Poisson image editing》中的核心思路。资源围绕泊松融合与泊松修复两条主线包含两个可直接运行的MATLAB脚本分别实现图像无缝融合与破损区域修复并附带多张练习图片便于观察融合前后效果、掩膜作用以及不同输入下的结果差异。整个压缩包大小约213KB上游未标注精确文件总数与类型清单从描述看以.m脚本和示例图像为主。截至目前已有710人学习下载适合正在学习图像处理、计算摄影或复现经典论文的初学者与研究人员。通过运行示例可直观理解泊松方程在梯度域编辑中的离散化构造、边界条件设置与线性方程组求解流程也能为后续实现更大规模图像融合算法提供可复用的脚本框架。 关于泊松融合很多刚开始接触图像合成的人都会有一个共同困惑明明两张图都选好了位置也摆正了直接用target(rows, cols) source(rows, cols)贴过去的一瞬间接缝处那条光晕和断层是怎么都藏不住。我最早做这个项目时也被这个问题卡了很久后来发现真正的解法根本不是像素级抠图而是去求一个泊松方程。这个思路一旦转过来很多合成问题都会豁然开朗。这篇文章我会直接围绕我在 MATLAB 里的实现来写包含泊松融合的原理推导怎么落到离散网格上、稀疏矩阵怎么组装、实际调试中反复踩到的坑以及混合梯度这种进阶玩法。适合已经会用 MATLAB 处理图像、想摆脱硬边缘拼图这个阶段的人也适合刚接触泊松方程、想通过实例理解它到底在干什么的读者。1. 直接拼图为什么永远有接缝感我最初的困惑我第一次做图像合成实验时用的是最直白的方式把源图像某块区域直接复制到目标图像上再手动做一点羽化。结果问题非常明显——物体边缘哪怕像素对齐了周围的光照、色调、明暗过渡依然是断开的。羽化只是把硬边变成软边本质上还是在两个不同的光照场之间做了线性插值所以这块区域看起来要么发灰要么像蒙了一层雾。后来我稍微进阶了一点学着用alpha混合给源区域中心分配更高的权重、边缘降低权重。这个方法对纹理简单、光照接近的场景有一定效果但只要目标背景和源图像的光照方向不一致过渡区域还是会出现明显的脏带因为 alpha 混合只是在 RGB 三个通道上做加权平均它并不知道图像内部的结构信息比如边缘方向、纹理所对应的空间变化率。真正改变我思路的是一篇关于泊松融合的文章。它提出的目标不是让两张图的像素值相加而是让融合后的区域其梯度结构和源图像尽量一致。换句话说拼接后那块内容只要变化趋势、边缘走向和源图一致整体明暗可以被目标背景重新决定。这正好解释了为什么人眼看接缝会敏感人类视觉系统对图像里的局部亮度变化率也就是梯度比对绝对亮度值更敏感。只要梯度连续亮度绝对值和周围差个零点几反而是能被大脑接受的。泊松融合还有一个非常实用的特点它天然满足边界连续性。合成区域边缘的像素值被强制约束为目标图像原本的值中间部分则完全由源图像的梯度场导演。于是结果常常是内容来自源图光照氛围来自目标图边界处没有过渡带也没有羽化痕迹。2. 泊松方程的离散化从物理场到线性方程组在连续数学里泊松方程长这样Δf div v。这里的 f 是待求解的像素强度场v 是引导场通常取源图像的梯度Δ 是拉普拉斯算子div 是散度。你不需要把它想得太神秘可以类比成热量传导区域边界温度固定内部热量梯度由源场决定最终温度分布在满足边界的条件下让内部梯度最接近引导场。图像是离散网格所以这个方程必须转成差分形式。对网格上某个像素 (i, j)二阶导用差分近似后拉普拉斯算子变成f(i-1, j) f(i1, j) f(i, j-1) f(i, j1) - 4f(i, j)。这是图像处理里最常用的离散拉普拉斯模板。引导场的散度 div v 也需要离散化。如果引导场是两个分量 gx 和 gy分别表示 x 方向和 y 方向的梯度那么散度大约等于 gx 在 x 方向上的差分加上 gy 在 y 方向上的差分。于是对 mask 内部每个像素我都能写出一个方程4f(i,j) - f(i-1,j) - f(i1,j) - f(i,j-1) - f(i,j1) divG(i,j)这里有个关键点如果某个邻居像素在 mask 外它的 f 值是已知的就是目标图像在这个位置的像素值所以要把这一项移到等式右边变成4f(i,j) - f(内部邻居...) divG(i,j) 已知边界邻居值把所有未知像素的方程堆在一起就是一个标准的稀疏线性方程组 A * x b。A 的每一行对应一个未知像素中心系数为 4内部邻居系数为 -1右端项由散度和边界值构成。求出 x 后把解写回 mask 区域就完成了泊松融合。这里我不推荐把整张图的像素都放进方程组哪怕非 mask 区域在矩阵里全设成单位行也没必要。正确做法是只对 mask 内的像素建立紧凑索引求解规模就是 mask 内的像素数量速度和解的稳定性都会好很多。3. MATLAB 实现从拉普拉斯到稀疏矩阵我实际写 MATLAB 版本的时候发现最顺手的路径是先算引导场再组装矩阵最后直接求解。以单通道为例我会先准备三张尺寸完全一致的图target、source、mask。function out poisson_fusion_single(target, source, mask) % 输入三张图需要保证大小一致 target im2double(target); source im2double(source); mask logical(mask); [h, w] size(target); % 1. 计算源图像的梯度场 [sgx, sgy] gradient(source); % 2. 计算引导场的散度 [gxx, ~] gradient(sgx); [~, gyy] gradient(sgy); divG gxx gyy; % 3. 对 mask 内部像素建立紧凑索引 maskIdx find(mask); % mask 内像素的线性索引 M numel(maskIdx); pixMap zeros(h, w); pixMap(maskIdx) 1:M; % 从原始索引映射到 1..M % 4. 组装稀疏矩阵 A 和右端项 b A sparse(M, M); b zeros(M, 1); for k 1:M [i, j] ind2sub([h, w], maskIdx(k)); A(k, k) 4; % 中心系数 nb [i-1, j; i1, j; i, j-1; i, j1]; for t 1:4 ni nb(t, 1); nj nb(t, 2); if ni 1 ni h nj 1 nj w nbLinear sub2ind([h, w], ni, nj); if mask(nbLinear) A(k, pixMap(nbLinear)) -1; else b(k) b(k) target(ni, nj); end end end b(k) b(k) divG(i, j); end % 5. 求解并写回 sol A \ b; out target; out(mask) sol; end这段代码是标准的教学写法用来理解逻辑非常清楚。求解 A \ b 的时候我建议先确认 mask 区域内像素数量级。如果只有几万像素直接反斜杠求解没有问题如果 mask 是几十万像素的大区域A 会达到几十万阶这时候速度会很慢内存占用也高建议换成迭代求解器比如pcg(A, b, 1e-6, 200)或者提前用ichol做预处理。我实测下来中等尺寸的 mask 用A \ b通常在几百毫秒到一两秒内完成真正需要注意的反而是 mask 不要做得太满下面会展开说。关于彩色图像我看到的很多教程会把 RGB 图像直接转成灰度再融合这是不对的。正确做法是对 R、G、B 三个通道分别计算梯度场、分别组装右端项、分别求解最后再把三个通道合回彩色。原理上讲三个通道共享同一个 A 矩阵所以可以只组一次系数阵只是 b 变成 M×3 的矩阵这样效率更高。这也是为什么很多 MATLAB 代码会把 A 提出循环、只对不同通道换 b 的原因。4. 实际操作中反复踩到的坑mask、边界、颜色漂移这部分我想认真说说我在调试里真实遇到的问题因为如果只是照着书上的公式写代码很可能跑出看起来像融合、其实颜色完全不对劲的结果。第一个坑是 mask 区域贴到了图像边界。泊松融合依赖 Dirichlet 边界条件也就是 mask 边缘的像素值需要用目标图像对应位置的值来约束。如果 mask 直接延伸到图像边缘那么边缘像素的四邻域会落在图像外代码里这些邻域就会被跳过中心系数仍然是 4导致方程在该像素处缺邻居却不减系数整个解会偏向异常值。我的处理习惯是让 mask 离图像边缘至少留出几个像素的空白误差会小很多。第二个坑是源图、目标图、mask 尺寸不一致。这是我早期最常犯的错误用一个 resize 过的 source 配上原始尺寸的 mask结果gradient函数计算的引导场和 mask 对不上。调试手段不复杂第一行就检查size(target) size(source)和size(mask) size(target)不一致就直接报错不要试图在循环里靠min兜底那样只会掩盖问题。第三个坑是最容易忽略的散度的符号搞反。拉普拉斯算子有两种查分符号惯例有的资料写 Δf f_up f_down f_left f_right - 4f有的写 Δf 4f - ...如果把右端项的 divG 按另一种惯例代入结果会变成反向扩散图像会出现类似浮雕的怪异效果。解决方法是构造一个很小的测试把 source 和 target 设置成同一张图mask 选一块区域此时理想结果应该和原图几乎一致。如果结果出现明显变暗或变灰大概率就是散度符号反了。第四个坑是关于稀疏矩阵中非 mask 像素被留成零行。我在早期版本里图省事给整张图的每一个像素都建立顺序编号然后给所有像素都写方程非 mask 像素直接放一个单位行。这样 N 会等于整个图像素数量如果图片是 1000×1000A 就是百万阶稀疏矩阵虽然稀疏但还是白白浪费内存和计算时间。改成只对 mask 内部像素编号后效率提升非常明显代码也更干净。颜色漂移问题也常在这个时候出现如果只用灰度图计算引导场再应用到 RGB三个通道共享一个散度场色彩细节会被磨掉所以彩色图一定要逐通道处理。第五个坑是关于 mask 边缘的过度锐利。泊松融合本身擅长保持边界自然但前提是 mask 边缘覆盖的区域确实是真正需要迁移的对象边缘。如果 mask 边缘选在一个渐变区域比如天空或皮肤高光中间那么即便数学上完全正确融合结果也会在 mask 边缘产生微小的亮度跳变。实际操作中我会用imdilate把 mask 向外膨胀一两个像素让融合区域稍微覆盖一点背景这样对话框边缘和外部的过渡会更有余地。5. 从经典泊松融合到混合梯度一个更稳的进阶方向经典的泊松融合直接把源图像的梯度作为引导场这在很多场景下效果很好但一旦源图像里有透明物体、遮挡或者高频纹理迁移需求就会出现源图纹理被完全保留包括不该保留的高光或反射的问题。这时可以用混合梯度mixing gradients来改进。混合梯度的核心思想很简单在 mask 内部逐像素比较源图梯度和目标图像梯度的模长谁更大就用谁。公式表达是如果源图梯度模长大于目标图梯度模长取源图梯度反之取目标图梯度。这样做的好处是源图中的强纹理仍然会被保留比如文字的笔画、物体的轮廓但源图中相对平坦、而目标图中原本有微弱渐变的地方会沿用目标图的渐变避免把目标背景里自然的明暗过渡覆盖掉。在 MATLAB 里实现这个改动非常顺手只需要在准备引导场时做一次逐像素比较[tgx, tgy] gradient(target); useSource (abs(sgx) abs(sgy)) (abs(tgx) abs(tgy)); gx useSource .* sgx ~useSource .* tgx; gy useSource .* sgy ~useSource .* tgy;这个改动虽然只多了几行代码但效果差别非常大。我做过一个小实验把一张含有反光物体的照片迁移到另一张背景图上经典泊松融合会把反光里的白色光斑也一起迁移过去导致结果里那块就像贴了一块白斑换成混合梯度后反光区域因为目标背景的梯度更强就自动选择了背景的渐变整体观感自然很多。混合梯度也有它需要注意的地方如果源图和目标图在同一个像素位置的梯度大小接近切换会变得很敏感可能造成微小抖动。比较合理的做法是在比较模长前做一次高斯平滑或者干脆给 source 梯度设一个小阈值偏移避免在平坦区域频繁切换。具体阈值要根据图像噪声水平调整我通常用 0.01 到 0.05 之间。除了混合梯度泊松融合在项目中还能扩展到图像修复、多图拼接去接缝、纹理迁移这些方向。实现思路上其实是同一个框架算梯度、组装稀疏矩阵、求解泊松方程。你只需要把引导场从源图像梯度换成被打乱后的梯度或者多张图像加权后的梯度很多看起来高大上的效果就都能在 MATLAB 里复现。6. 给新手的实验建议做一次最小闭环我建议第一次做泊松融合实验时不要直接拿大图开始先准备一张 200×200 左右的小图目标图是一张渐变背景源图是另一张有明显边缘图形的小图mask 直接画一个矩形。这样无论是调试散度符号还是观察边界过渡成本都极低跑一次只需要零点几秒。调试时可以把 A 矩阵的稀疏结构用spy(A)画出来看一眼正常情况应该是一条窄带因为每个像素只和上下左右四个邻居相连矩阵带状结构很清晰。如果看到一堆杂乱的非零块说明索引映射写错了。最后一个小经验mask 区域如果很大与其让 A 矩阵直接参与A \ b不如考虑用decomposition(A)或预处理共轭梯度法。因为泊松方程在大尺寸 mask 下本质上是对一个近似拉普拉斯矩阵求逆慢一点不奇怪但没必要让代码卡死在求解这一步。先控制 mask 尺寸、先跑通逻辑再谈优化这是最稳的顺序。本文还有配套的精品资源点击获取