简介泊松融合是图像编辑中的经典算法源自 Pérez 等人 2003 年的 SIGGRAPH 论文常用于无缝拼接、纹理平滑与图像修复。这份 MATLAB 开源实现面向计算机视觉学习者、图像处理研究者及需要快速验证泊松方程效果的开发者完整对应论文方法进行复现重点解决区域融合后亮度不连续、边界痕迹明显等问题。压缩包约 213KB主要文件为 .m 脚本和示例图片核心包含 Poisson Fusion实现泊松融合与 Poisson Repair实现泊松修复两个 MATLAB 脚本并附有数张练习图片方便直接运行和对比原图、掩膜与目标结果。目前已有 710 人浏览学习。借助该代码可直观理解泊松方程在图像梯度场重构中的实际用法包括边界条件设置、像素迭代求解等细节无论是对照经典论文学习还是在此基础上调整参数开展实验都较为适合。1. 泊松融合到底解决什么问题做图像拼接、换脸、物体移植时最头疼的一件事明明抠图抠得很干净一贴到目标图上接缝处总有一圈肉眼可见的边界。颜色哪怕只差一点点边缘就像贴了道黑线或光晕。泊松融合Poisson Image Editing就是专门治这个问题的经典方案核心思路不是去调透明度、做羽化而是把问题转化成解一个泊松方程让融合区域的像素梯度去逼近源图像的梯度同时边界颜色严格等于目标图这样过渡区域就能做到理论上的无缝。这个方向从2003年提出来之后一直是图像编辑、材质合成、遥感影像拼接的标配算法。用MATLAB实现它适合做实验验证、毕业设计、算法对比也是理解梯度域图像处理最直接的一条路。本文按“原理 → 最小实现 → 验证 → 排错 → 进阶”推进每一步都给可运行的代码和参数说明。2. 从图像梯度到泊松方程把直觉翻译成矩阵2.1 为什么直接调像素颜色做不到无缝先看一个直观例子把源图里的一个苹果拖到目标图的桌面上。苹果区域内部颜色没问题问题出在苹果轮廓那一圈。最简单的方法是直接覆盖像素但源图和目标图在边界上的亮度通常有差异哪怕只差5个灰度级人眼也能在边缘感知到一条细线。羽化或者透明度混合只能把这条线变模糊不能真正消除它因为两者在边界处的梯度方向、大小本身就冲突。泊松融合换了一个思路不再直接拷贝颜色而是拷贝“变化量”。源图苹果内部苹果和背景的对比关系、高光的过渡、阴影的渐变这些信息可以浓缩成一个梯度场。只要融合区域内部的梯度场和源图一致那么人眼感知到的就是同一个苹果边界上颜色跟目标图硬性对齐就不会有跳变。这就像给一块拼图定了边框颜色内部怎么填由源图说了算。2.2 连续域方程和离散化拉普拉斯算子的角色在连续域里给定一个目标区域Ω边界∂Ω上的像素值已知等于目标图的原值要让融合结果 f 在Ω内部的梯度尽可能接近源图 g 的梯度等价于最小化下面这个积分E(f) ∬_Ω |∇f − ∇g|² dx dy且 f|∂Ω f_target|∂Ω把这个泛函做变分得到的就是泊松方程Δf Δg即 ∇·(∇f) ∇·(∇g)这里的Δ是拉普拉斯算子Δg其实就是源图像在Ω内部的二阶信息。也就是说算法要干的事就是把源图的拉普拉斯结果当作驱动项在目标区域里重解一遍这个方程同时边界值钉死在目标图上。离散化处理时图像上的拉普拉斯算子用一个3×3模板近似Δf(i,j) ≈ 4f(i,j) − f(i−1,j) − f(i1,j) − f(i,j−1) − f(i,j1)逐像素写出来就是一个线性方程组对Ω内的每一个像素列一条方程左边是5点差分模板右边是源图的离散拉普拉斯值。边界上那些邻居不在Ω内就把它移到等式右边值用目标图的像素。2.3 MATLAB中最终要解的线性系统长什么样把Ω内所有像素排成一个列向量x系数矩阵A就是每条方程按5点模板组出来的稀疏矩阵等号右边是b。MATLAB里求解就是一行命令 x A \ b。这个A不是全图那么大只覆盖融合区域。Ω内像素数记为n那么A就是一个n×n的稀疏矩阵每行只有5个非零值实际可能少于5因为边界邻居被移走了。系数矩阵A不需要手工el sparse循环填推荐先给Ω内每个像素编一个序号然后收集所有方程的下标三元组一次性sparse。常见做法是[h, w, ~] size(imTarget); maskIdx find(mask); % mask非零像素线性索引 im2var zeros(h, w); im2var(maskIdx) 1:numel(maskIdx); % 每个融合区域内像素分配到未知数编号 rows []; cols []; vals []; % 拉普拉斯模板系数 for k 1:numel(maskIdx) [i, j] ind2sub([h, w], maskIdx(k)); idx im2var(i, j); % 中心点系数4 rows(end1) idx; cols(end1) idx; vals(end1) 4; % 上邻居 if i 1 nid im2var(i-1, j); if nid 0, rows(end1)idx; cols(end1)nid; vals(end1)-1; else, b(k) b(k) imTarget(i-1, j, ch); end end % 下邻居、左邻居、右邻居同理... end A sparse(rows, cols, vals, n, n); x A \ b;这段代码的逻辑说明每个方程对应Ω内一个未知像素中心点占系数4四个方向若也在Ω内则系数−1若不在则把该方向的目标图像素值加到b上。注意b的初值要先存成源图的拉普拉斯值 Δg也就是 b -4g(i,j) g(i−1,j) g(i1,j) g(i,j−1) g(i,j1)这样方程才是4f_center − Σf_neighbor Δg。上面的代码片段只展示了矩阵构建实际写完整函数时b要先用源图的拉普拉斯填充再加边界修正项。2.4 三通道处理和矩阵规模预估彩色图需要按RGB三个通道分别求解。因为泊松方程是线性操作三通道互不干扰所以循环遍历通道或者把b拼成n×3矩阵一次求解都行。建议用后者MATLAB里A \ B对多列右端项很高效。矩阵规模预估假设融合区域是个500×400的矩形n大约20万A用稀疏存储只有约100万非零元MATLAB直接法求解大概需要几秒到十几秒取决于机器内存。如果Ω形状很细长不规则非零元结构稀疏但求解时间可能略涨。遇到内存不足时优先检查是否误用了full(A)。3. 最小可跑通的MATLAB实现从两张图到无缝结果3.1 输入准备源图、目标图、掩膜掩膜是关键输入它决定融合区域。掩膜为1的区域表示要把源图的梯度搬过来边界像素对齐目标图。掩膜可以由手动抠图生成也可以用ROI多边形工具。我做实验时习惯先用一个矩形区域把整个流程跑通再换真实掩膜。% 读图 imSrc im2double(imread(apple.png)); % 源图 imDst im2double(imread(desk.png)); % 目标图 % 掩膜生成: 以源图中一个矩形区域为例 mask false(size(imSrc,1), size(imSrc,2)); mask(50:250, 80:280) true; % 矩形融合区 % 掩膜大小必须与源图一致目标图可以更大 % 这里简化处理假设目标图已经裁剪到和源图一样大 mask imresize(mask, [size(imDst,1), size(imDst,2)], nearest);参数说明im2double把图像归一化到[0,1]避免求解时数值过大。掩膜用logical类型不要用uint8后面逻辑索引时省事。如果源图和目标图尺寸不一致常见做法是先对源图做仿射变换放到目标图上同时变换掩膜变换矩阵可以用fitgeotrans算。3.2 主函数梯度和拉普拉斯的MATLAB表达MATLAB里计算离散拉普拉斯最方便的是del2它对数组做拉普拉斯算子并除以4。也可以用imfilter配合[0 1 0; 1 -4 1; 0 1 0]这个核我趋向于用后者因为系数直观且不受del2那个1/4因子的混淆影响。function result poissonFusion(imSrc, imDst, mask) [h, w, c] size(imDst); mask mask 0; n sum(mask(:)); im2var zeros(h, w); im2var(mask) 1:n; % 给区域像素编号 lapKernel [0 1 0; 1 -4 1; 0 1 0]; % 标准五点拉普拉斯 % 对每个通道构建右端项 b b zeros(n, c); for ch 1:c lapSrc imfilter(imSrc(:, :, ch), lapKernel, replicate); b(:, ch) lapSrc(mask); % 源图拉普拉斯值作为驱动项 end % 构建系数矩阵 A rows zeros(5*n, 1); cols zeros(5*n, 1); vals zeros(5*n, 1); cnt 0; offsets [-1 0; 1 0; 0 -1; 0 1]; for idx 1:n [i, j] ind2sub([h, w], maskIdx(idx)); cnt cnt 1; rows(cnt) idx; cols(cnt) idx; vals(cnt) 4; for k 1:4 ni i offsets(k, 1); nj j offsets(k, 2); if ni 1 || ni h || nj 1 || nj w continue; end nid im2var(ni, nj); if nid 0 cnt cnt 1; rows(cnt) idx; cols(cnt) nid; vals(cnt) -1; else % 邻居在区域外从b中减去边界像素值 for ch 1:c b(idx, ch) b(idx, ch) imDst(ni, nj, ch); end end end end A sparse(rows(1:cnt), cols(1:cnt), vals(1:cnt), n, n); % 求解并写回目标图 x A \ b; result imDst; for ch 1:c valsCh zeros(h, w); valsCh(mask) x(:, ch); result(:, :, ch) imDst(:, :, ch); result(:, :, ch) result(:, :, ch) .* ~mask valsCh .* mask; end end逻辑说明先算源图的拉普拉斯存进b作为泊松方程右端项接着构建A矩阵5点模板中对区域外邻居的处理是把目标图像素值加到b相当于把该未知数消掉了最后直接求解线性系统。这里b的边界修正加的是imDst的像素值不是imSrc的。x是三列矩阵分别对应R、G、B三个通道的解最后一步把解覆盖到目标图的掩膜区域。3.3 参数调整和数值上的两个细节求解器直接用反斜杠背后走稀疏LU分解对n≤50万的系统稳定且快。如果n极大可以考虑改用pcg但要提供好的预条件子。数值上有一个容易翻车的点源图和目标图的数据类型必须都是double且范围一致。如果一张是uint8一张是double离散拉普拉斯会差255倍融合结果直接全黑或全白。另外imfilter的边界处理选replicate这对区域内部像素没有影响因为有边界修正逻辑replicate只影响掩膜边缘向外延伸时拉普拉斯模板里的值为0的情况。用symmetric也不会差太多但结果会有12个像素的细微差异。4. 结果验证怎么确认融合做对了4.1 重建测试把融合结果变回梯度场判断泊松融合代码是否正确的最简单方法是对已经融合的图再算一次拉普拉斯对比源图的拉普拉斯。在Ω内部这两个应该几乎相等差别只来自线性求解器的数值误差一般小于1e-6。lapResult imfilter(result, lapKernel, replicate); lapSrc imfilter(imSrc, lapKernel, replicate); diffMap zeros(h, w); for ch 1:c tmp abs(lapResult(:, :, ch) - lapSrc(:, :, ch)); diffMap max(diffMap, tmp); end diffRegion diffMap(mask); fprintf(拉普拉斯最大误差: %.6f\n, max(diffRegion(:)));这段代码是验证的核心重点看mask区域内部的误差。如果误差在1e−2量级说明A矩阵或b构建有问题。注意diffMap在mask边界上也包含了一个像素宽的过渡区因为那边有边界项的贡献这部分可以忽略。4.2 客观指标亮度一致性和接缝对比融合是否成功还可以用两个客观指标量化。第一是接缝处梯度变化率把融合结果沿mask边界求梯度模数值越小说明过渡越平滑。第二是内部结构保持计算融合结果与源图在mask区域内的梯度相关系数。% 接缝梯度: 对mask取边界 se strel(diamond, 1); boundaryMask imdilate(mask, se) ~mask; [gx, gy] gradient(result(:,:,1)); edgeGrad sqrt(gx.^2 gy.^2); seamValue mean(edgeGrad(boundaryMask));边界mask用膨胀减去原mask得到目标是在融合区域的轮廓线上计算平均梯度幅值。这个值如果明显大于融合区域内部平均值说明接缝还有可见痕迹。正常情况seamValue应该和区域内部平均梯度在同一数量级不会高出数倍。4.3 可视化检查的三个观察点指标靠谱但最终判断还是靠眼睛。把源图、目标图、融合结果、差值图拼在一张画布上重点看三个位置第一是融合区域的边界线那边不应该有深色或亮色轮廓第二是原本高梯度区域比如物体边缘那边细节不能糊掉第三是平坦区域比如纯色桌面那边不应该出现源图纹理的透传。差值图用源图减融合结果再取绝对值看结构残留。我用过一个经验阈值差值图在mask内95%的像素不超过0.0501范围基本可以交付。5. 泊松融合MATLAB常见问题排查现象、原因、解决5.1 融合结果出现整体色偏尤其是边界发黑发亮现象是融合区域颜色从边缘向内部扩散了一层暗色或亮色晕。原因基本是b的边界修正写错了常见的情况是把边界邻居的像素值直接加进b却没有和拉普拉斯符号对齐。拉普拉斯模板中心为−4、邻居为1时方程右边要把边界值移项也就是加的是正值的边界像素如果模板写反了中心为正那就要改成负的。解决方法是先把拉普拉斯模板统一写成[0 1 0; 1 -4 1; 0 1 0]b里用源图的lapSrc(mask)这样边界修正就是加目标图像素值不容易错。我踩过这个坑时排查办法是把n1的矩形区域打印出A和b手算一遍验证比盯着屏幕猜快得多。5.2 融合结果变成源图直接贴上去没有无缝过渡现象是mask内部完全是源图的像素边界一圈硬边还在。原因多半是A构建时条件不满足方程退化成了逐像素拷贝也就是对角线系数被设置成1而邻居系数为0。常见触发场景是imfilter用了conv模式且核没归一化或者掩膜和图像大小不一致导致mask(maskIdx)取到的值和源图像素错位。解决思路先验证mask和imSrc、imDst尺寸完全一致打印size确认再确认im2var(mask)1:n这行执行后区域内像素编号确实是按列优先顺序。最后用前面4.1的拉普拉斯误差测试如果误差接近0但边界仍硬那就是结果写回时公式不对检查最后一步赋值是否把融合结果覆盖到mask外部了。5.3 求解特别慢几分钟都出不来现象是整个函数跑完要几分钟。原因基本是A不小心变成了满矩阵。常见误操作是构建A之前用了full(A)调试或者sparse调用时把n算成了图像总像素。另一个高频原因是循环里反复拼接rows、cols、vals数组MATLAB会不断重新分配内存。解决方法是预分配数组如rows zeros(5*n, 1)这种写法或者把整个构建改成向量化思路。n超过50万时直接法也会变慢可以换pcg加ichol预条件子。我测试过一个20万像素的区域向量化构建之后再求解整体时间从80秒降到4秒瓶颈基本在构建而不是求解。5.4 融合边缘有一圈淡淡的轮廓线但区域内部没问题现象是内部颜色过渡都自然唯独融合区域最外一两像素有一条隐约的亮线或暗线。原因是掩膜边缘不够精确二值mask在边界产生了阶梯效应拉普拉斯模板在边界附近取的邻居分布不对称导致那一两行像素的方程系数和内部不一致。解决办法有两个一个是对mask做一次腐蚀或膨胀让边界落在源图纹理的低频区另一个是用抗锯齿化处理把mask边缘模糊0.51像素再做归一化。注意模糊之后mask不再是逻辑值系数矩阵构建要按mask0.5来判定b也要乘上权重。这个做法在边界处等效于软约束比硬切干净。5.5 彩色图像融合后色彩饱和度明显下降现象是融合区域颜色发灰像蒙了一层雾。原因是把RGB三个通道分开求解但b三个通道的拉普拉斯值差别很大而A完全一样解出来三个通道的数值相对关系被压缩了。严格说这个问题的根源不在泊松融合本身而是在于RGB色彩空间通道间相关性高梯度差异在突变区域相互抵消。解决方法是换成Lab或YUV空间只对亮度通道做泊松求解色度通道直接拷贝或者线性混合。实际做物体移植时这个改动对肤色、植物绿色这些高饱和区域效果改善非常明显而改动成本只是加一个rgb2lab和lab2rgb的转换。6. 进阶用法混合梯度和多区域并行融合混合梯度是泊松融合最常用的变体做法是把源图和目标图的梯度逐点比较取模更大的那个作为融合方程的右端项。好处是能同时保留源图的纹理和目标图的背景结构适合纹理替换、去除遮挡这类任务。在MATLAB里实现起来只需要改b的构建gradSrcX diff([imSrc(:,2:end,:) imSrc(:,1,:)], 1, 2); gradSrcY diff([imSrc(2:end,:,:); imSrc(1,:,:)], 1, 1); gradDstX diff([imDst(:,2:end,:) imDst(:,1,:)], 1, 2); gradDstY diff([imDst(2:end,:,:); imDst(1,:,:)], 1, 1); magSrc sqrt(gradSrcX.^2 gradSrcY.^2); magDst sqrt(gradDstX.^2 gradDstY.^2); useSrc magSrc magDst; gradMixX gradSrcX .* useSrc gradDstX .* ~useSrc; gradMixY gradSrcY .* useSrc gradDstY .* ~useSrc; lapMix div(gradMixX, gradMixY); % 需要用散度算子或从梯度重建拉普拉斯注意diff在MATLAB里的方向约定diff(X, 1, 2)算水平差分得到的矩阵会比原图少一列上面代码用循环填充的方式维持了尺寸。实际做这个变体时我更喜欢直接构建散度算子矩阵把梯度场转成散度项这样不会引入周期边界假设。散度算子是梯度算子的转置负值在离散网格上最终化简出来仍然是5点模板只不过每一项的系数是两边梯度的符号差。多区域融合的常见做法是如果源图里有多块不连通区域要分别融合不要整个mask一起求解。n太大时矩阵条件数变差且区域间如果有间隙泊松方程会把两块区域的边界也一起平滑这不是想要的。我的习惯是每个连通域单独调用一次主函数最后把结果合回一张图。如果多个区域位置相邻但语义不同比如把一只鸟和一架飞机同时贴到天空分开solve比一次solve在边界控制上更精确。验证多区域结果时只需要对每个连通域重复4.1的拉普拉斯误差测试。实践下来我发现一个经验mask区域边界外扩23个像素再截断融合结果会比原始mask边缘更干净因为边界约束的邻域更广方程对边界噪声更鲁棒。代价是融合区域比预期略微扩大如果对区域范围有严格限制就不能用这招。写到这里想起最早我在这套代码上翻车最狠的一次就是b的负号问题整个融合结果像被人糊了一层滤镜。后来养成一个习惯每次拿到新图片先跑矩形测试用例确认拉普拉斯误差在1e−6量级再去换真实掩膜。这个习惯帮我避开了至少七八次莫名其妙的色偏问题。泊松融合这个方向一旦理解成“梯度域的搬运工”代码量其实很少难的永远在边界处理和验证环节。希望帮到你。本文还有配套的精品资源点击获取
简介泊松融合是图像编辑中的经典算法源自 Pérez 等人 2003 年的 SIGGRAPH 论文常用于无缝拼接、纹理平滑与图像修复。这份 MATLAB 开源实现面向计算机视觉学习者、图像处理研究者及需要快速验证泊松方程效果的开发者完整对应论文方法进行复现重点解决区域融合后亮度不连续、边界痕迹明显等问题。压缩包约 213KB主要文件为 .m 脚本和示例图片核心包含 Poisson Fusion实现泊松融合与 Poisson Repair实现泊松修复两个 MATLAB 脚本并附有数张练习图片方便直接运行和对比原图、掩膜与目标结果。目前已有 710 人浏览学习。借助该代码可直观理解泊松方程在图像梯度场重构中的实际用法包括边界条件设置、像素迭代求解等细节无论是对照经典论文学习还是在此基础上调整参数开展实验都较为适合。1. 泊松融合到底解决什么问题做图像拼接、换脸、物体移植时最头疼的一件事明明抠图抠得很干净一贴到目标图上接缝处总有一圈肉眼可见的边界。颜色哪怕只差一点点边缘就像贴了道黑线或光晕。泊松融合Poisson Image Editing就是专门治这个问题的经典方案核心思路不是去调透明度、做羽化而是把问题转化成解一个泊松方程让融合区域的像素梯度去逼近源图像的梯度同时边界颜色严格等于目标图这样过渡区域就能做到理论上的无缝。这个方向从2003年提出来之后一直是图像编辑、材质合成、遥感影像拼接的标配算法。用MATLAB实现它适合做实验验证、毕业设计、算法对比也是理解梯度域图像处理最直接的一条路。本文按“原理 → 最小实现 → 验证 → 排错 → 进阶”推进每一步都给可运行的代码和参数说明。2. 从图像梯度到泊松方程把直觉翻译成矩阵2.1 为什么直接调像素颜色做不到无缝先看一个直观例子把源图里的一个苹果拖到目标图的桌面上。苹果区域内部颜色没问题问题出在苹果轮廓那一圈。最简单的方法是直接覆盖像素但源图和目标图在边界上的亮度通常有差异哪怕只差5个灰度级人眼也能在边缘感知到一条细线。羽化或者透明度混合只能把这条线变模糊不能真正消除它因为两者在边界处的梯度方向、大小本身就冲突。泊松融合换了一个思路不再直接拷贝颜色而是拷贝“变化量”。源图苹果内部苹果和背景的对比关系、高光的过渡、阴影的渐变这些信息可以浓缩成一个梯度场。只要融合区域内部的梯度场和源图一致那么人眼感知到的就是同一个苹果边界上颜色跟目标图硬性对齐就不会有跳变。这就像给一块拼图定了边框颜色内部怎么填由源图说了算。2.2 连续域方程和离散化拉普拉斯算子的角色在连续域里给定一个目标区域Ω边界∂Ω上的像素值已知等于目标图的原值要让融合结果 f 在Ω内部的梯度尽可能接近源图 g 的梯度等价于最小化下面这个积分E(f) ∬_Ω |∇f − ∇g|² dx dy且 f|∂Ω f_target|∂Ω把这个泛函做变分得到的就是泊松方程Δf Δg即 ∇·(∇f) ∇·(∇g)这里的Δ是拉普拉斯算子Δg其实就是源图像在Ω内部的二阶信息。也就是说算法要干的事就是把源图的拉普拉斯结果当作驱动项在目标区域里重解一遍这个方程同时边界值钉死在目标图上。离散化处理时图像上的拉普拉斯算子用一个3×3模板近似Δf(i,j) ≈ 4f(i,j) − f(i−1,j) − f(i1,j) − f(i,j−1) − f(i,j1)逐像素写出来就是一个线性方程组对Ω内的每一个像素列一条方程左边是5点差分模板右边是源图的离散拉普拉斯值。边界上那些邻居不在Ω内就把它移到等式右边值用目标图的像素。2.3 MATLAB中最终要解的线性系统长什么样把Ω内所有像素排成一个列向量x系数矩阵A就是每条方程按5点模板组出来的稀疏矩阵等号右边是b。MATLAB里求解就是一行命令 x A \ b。这个A不是全图那么大只覆盖融合区域。Ω内像素数记为n那么A就是一个n×n的稀疏矩阵每行只有5个非零值实际可能少于5因为边界邻居被移走了。系数矩阵A不需要手工el sparse循环填推荐先给Ω内每个像素编一个序号然后收集所有方程的下标三元组一次性sparse。常见做法是[h, w, ~] size(imTarget); maskIdx find(mask); % mask非零像素线性索引 im2var zeros(h, w); im2var(maskIdx) 1:numel(maskIdx); % 每个融合区域内像素分配到未知数编号 rows []; cols []; vals []; % 拉普拉斯模板系数 for k 1:numel(maskIdx) [i, j] ind2sub([h, w], maskIdx(k)); idx im2var(i, j); % 中心点系数4 rows(end1) idx; cols(end1) idx; vals(end1) 4; % 上邻居 if i 1 nid im2var(i-1, j); if nid 0, rows(end1)idx; cols(end1)nid; vals(end1)-1; else, b(k) b(k) imTarget(i-1, j, ch); end end % 下邻居、左邻居、右邻居同理... end A sparse(rows, cols, vals, n, n); x A \ b;这段代码的逻辑说明每个方程对应Ω内一个未知像素中心点占系数4四个方向若也在Ω内则系数−1若不在则把该方向的目标图像素值加到b上。注意b的初值要先存成源图的拉普拉斯值 Δg也就是 b -4g(i,j) g(i−1,j) g(i1,j) g(i,j−1) g(i,j1)这样方程才是4f_center − Σf_neighbor Δg。上面的代码片段只展示了矩阵构建实际写完整函数时b要先用源图的拉普拉斯填充再加边界修正项。2.4 三通道处理和矩阵规模预估彩色图需要按RGB三个通道分别求解。因为泊松方程是线性操作三通道互不干扰所以循环遍历通道或者把b拼成n×3矩阵一次求解都行。建议用后者MATLAB里A \ B对多列右端项很高效。矩阵规模预估假设融合区域是个500×400的矩形n大约20万A用稀疏存储只有约100万非零元MATLAB直接法求解大概需要几秒到十几秒取决于机器内存。如果Ω形状很细长不规则非零元结构稀疏但求解时间可能略涨。遇到内存不足时优先检查是否误用了full(A)。3. 最小可跑通的MATLAB实现从两张图到无缝结果3.1 输入准备源图、目标图、掩膜掩膜是关键输入它决定融合区域。掩膜为1的区域表示要把源图的梯度搬过来边界像素对齐目标图。掩膜可以由手动抠图生成也可以用ROI多边形工具。我做实验时习惯先用一个矩形区域把整个流程跑通再换真实掩膜。% 读图 imSrc im2double(imread(apple.png)); % 源图 imDst im2double(imread(desk.png)); % 目标图 % 掩膜生成: 以源图中一个矩形区域为例 mask false(size(imSrc,1), size(imSrc,2)); mask(50:250, 80:280) true; % 矩形融合区 % 掩膜大小必须与源图一致目标图可以更大 % 这里简化处理假设目标图已经裁剪到和源图一样大 mask imresize(mask, [size(imDst,1), size(imDst,2)], nearest);参数说明im2double把图像归一化到[0,1]避免求解时数值过大。掩膜用logical类型不要用uint8后面逻辑索引时省事。如果源图和目标图尺寸不一致常见做法是先对源图做仿射变换放到目标图上同时变换掩膜变换矩阵可以用fitgeotrans算。3.2 主函数梯度和拉普拉斯的MATLAB表达MATLAB里计算离散拉普拉斯最方便的是del2它对数组做拉普拉斯算子并除以4。也可以用imfilter配合[0 1 0; 1 -4 1; 0 1 0]这个核我趋向于用后者因为系数直观且不受del2那个1/4因子的混淆影响。function result poissonFusion(imSrc, imDst, mask) [h, w, c] size(imDst); mask mask 0; n sum(mask(:)); im2var zeros(h, w); im2var(mask) 1:n; % 给区域像素编号 lapKernel [0 1 0; 1 -4 1; 0 1 0]; % 标准五点拉普拉斯 % 对每个通道构建右端项 b b zeros(n, c); for ch 1:c lapSrc imfilter(imSrc(:, :, ch), lapKernel, replicate); b(:, ch) lapSrc(mask); % 源图拉普拉斯值作为驱动项 end % 构建系数矩阵 A rows zeros(5*n, 1); cols zeros(5*n, 1); vals zeros(5*n, 1); cnt 0; offsets [-1 0; 1 0; 0 -1; 0 1]; for idx 1:n [i, j] ind2sub([h, w], maskIdx(idx)); cnt cnt 1; rows(cnt) idx; cols(cnt) idx; vals(cnt) 4; for k 1:4 ni i offsets(k, 1); nj j offsets(k, 2); if ni 1 || ni h || nj 1 || nj w continue; end nid im2var(ni, nj); if nid 0 cnt cnt 1; rows(cnt) idx; cols(cnt) nid; vals(cnt) -1; else % 邻居在区域外从b中减去边界像素值 for ch 1:c b(idx, ch) b(idx, ch) imDst(ni, nj, ch); end end end end A sparse(rows(1:cnt), cols(1:cnt), vals(1:cnt), n, n); % 求解并写回目标图 x A \ b; result imDst; for ch 1:c valsCh zeros(h, w); valsCh(mask) x(:, ch); result(:, :, ch) imDst(:, :, ch); result(:, :, ch) result(:, :, ch) .* ~mask valsCh .* mask; end end逻辑说明先算源图的拉普拉斯存进b作为泊松方程右端项接着构建A矩阵5点模板中对区域外邻居的处理是把目标图像素值加到b相当于把该未知数消掉了最后直接求解线性系统。这里b的边界修正加的是imDst的像素值不是imSrc的。x是三列矩阵分别对应R、G、B三个通道的解最后一步把解覆盖到目标图的掩膜区域。3.3 参数调整和数值上的两个细节求解器直接用反斜杠背后走稀疏LU分解对n≤50万的系统稳定且快。如果n极大可以考虑改用pcg但要提供好的预条件子。数值上有一个容易翻车的点源图和目标图的数据类型必须都是double且范围一致。如果一张是uint8一张是double离散拉普拉斯会差255倍融合结果直接全黑或全白。另外imfilter的边界处理选replicate这对区域内部像素没有影响因为有边界修正逻辑replicate只影响掩膜边缘向外延伸时拉普拉斯模板里的值为0的情况。用symmetric也不会差太多但结果会有12个像素的细微差异。4. 结果验证怎么确认融合做对了4.1 重建测试把融合结果变回梯度场判断泊松融合代码是否正确的最简单方法是对已经融合的图再算一次拉普拉斯对比源图的拉普拉斯。在Ω内部这两个应该几乎相等差别只来自线性求解器的数值误差一般小于1e-6。lapResult imfilter(result, lapKernel, replicate); lapSrc imfilter(imSrc, lapKernel, replicate); diffMap zeros(h, w); for ch 1:c tmp abs(lapResult(:, :, ch) - lapSrc(:, :, ch)); diffMap max(diffMap, tmp); end diffRegion diffMap(mask); fprintf(拉普拉斯最大误差: %.6f\n, max(diffRegion(:)));这段代码是验证的核心重点看mask区域内部的误差。如果误差在1e−2量级说明A矩阵或b构建有问题。注意diffMap在mask边界上也包含了一个像素宽的过渡区因为那边有边界项的贡献这部分可以忽略。4.2 客观指标亮度一致性和接缝对比融合是否成功还可以用两个客观指标量化。第一是接缝处梯度变化率把融合结果沿mask边界求梯度模数值越小说明过渡越平滑。第二是内部结构保持计算融合结果与源图在mask区域内的梯度相关系数。% 接缝梯度: 对mask取边界 se strel(diamond, 1); boundaryMask imdilate(mask, se) ~mask; [gx, gy] gradient(result(:,:,1)); edgeGrad sqrt(gx.^2 gy.^2); seamValue mean(edgeGrad(boundaryMask));边界mask用膨胀减去原mask得到目标是在融合区域的轮廓线上计算平均梯度幅值。这个值如果明显大于融合区域内部平均值说明接缝还有可见痕迹。正常情况seamValue应该和区域内部平均梯度在同一数量级不会高出数倍。4.3 可视化检查的三个观察点指标靠谱但最终判断还是靠眼睛。把源图、目标图、融合结果、差值图拼在一张画布上重点看三个位置第一是融合区域的边界线那边不应该有深色或亮色轮廓第二是原本高梯度区域比如物体边缘那边细节不能糊掉第三是平坦区域比如纯色桌面那边不应该出现源图纹理的透传。差值图用源图减融合结果再取绝对值看结构残留。我用过一个经验阈值差值图在mask内95%的像素不超过0.0501范围基本可以交付。5. 泊松融合MATLAB常见问题排查现象、原因、解决5.1 融合结果出现整体色偏尤其是边界发黑发亮现象是融合区域颜色从边缘向内部扩散了一层暗色或亮色晕。原因基本是b的边界修正写错了常见的情况是把边界邻居的像素值直接加进b却没有和拉普拉斯符号对齐。拉普拉斯模板中心为−4、邻居为1时方程右边要把边界值移项也就是加的是正值的边界像素如果模板写反了中心为正那就要改成负的。解决方法是先把拉普拉斯模板统一写成[0 1 0; 1 -4 1; 0 1 0]b里用源图的lapSrc(mask)这样边界修正就是加目标图像素值不容易错。我踩过这个坑时排查办法是把n1的矩形区域打印出A和b手算一遍验证比盯着屏幕猜快得多。5.2 融合结果变成源图直接贴上去没有无缝过渡现象是mask内部完全是源图的像素边界一圈硬边还在。原因多半是A构建时条件不满足方程退化成了逐像素拷贝也就是对角线系数被设置成1而邻居系数为0。常见触发场景是imfilter用了conv模式且核没归一化或者掩膜和图像大小不一致导致mask(maskIdx)取到的值和源图像素错位。解决思路先验证mask和imSrc、imDst尺寸完全一致打印size确认再确认im2var(mask)1:n这行执行后区域内像素编号确实是按列优先顺序。最后用前面4.1的拉普拉斯误差测试如果误差接近0但边界仍硬那就是结果写回时公式不对检查最后一步赋值是否把融合结果覆盖到mask外部了。5.3 求解特别慢几分钟都出不来现象是整个函数跑完要几分钟。原因基本是A不小心变成了满矩阵。常见误操作是构建A之前用了full(A)调试或者sparse调用时把n算成了图像总像素。另一个高频原因是循环里反复拼接rows、cols、vals数组MATLAB会不断重新分配内存。解决方法是预分配数组如rows zeros(5*n, 1)这种写法或者把整个构建改成向量化思路。n超过50万时直接法也会变慢可以换pcg加ichol预条件子。我测试过一个20万像素的区域向量化构建之后再求解整体时间从80秒降到4秒瓶颈基本在构建而不是求解。5.4 融合边缘有一圈淡淡的轮廓线但区域内部没问题现象是内部颜色过渡都自然唯独融合区域最外一两像素有一条隐约的亮线或暗线。原因是掩膜边缘不够精确二值mask在边界产生了阶梯效应拉普拉斯模板在边界附近取的邻居分布不对称导致那一两行像素的方程系数和内部不一致。解决办法有两个一个是对mask做一次腐蚀或膨胀让边界落在源图纹理的低频区另一个是用抗锯齿化处理把mask边缘模糊0.51像素再做归一化。注意模糊之后mask不再是逻辑值系数矩阵构建要按mask0.5来判定b也要乘上权重。这个做法在边界处等效于软约束比硬切干净。5.5 彩色图像融合后色彩饱和度明显下降现象是融合区域颜色发灰像蒙了一层雾。原因是把RGB三个通道分开求解但b三个通道的拉普拉斯值差别很大而A完全一样解出来三个通道的数值相对关系被压缩了。严格说这个问题的根源不在泊松融合本身而是在于RGB色彩空间通道间相关性高梯度差异在突变区域相互抵消。解决方法是换成Lab或YUV空间只对亮度通道做泊松求解色度通道直接拷贝或者线性混合。实际做物体移植时这个改动对肤色、植物绿色这些高饱和区域效果改善非常明显而改动成本只是加一个rgb2lab和lab2rgb的转换。6. 进阶用法混合梯度和多区域并行融合混合梯度是泊松融合最常用的变体做法是把源图和目标图的梯度逐点比较取模更大的那个作为融合方程的右端项。好处是能同时保留源图的纹理和目标图的背景结构适合纹理替换、去除遮挡这类任务。在MATLAB里实现起来只需要改b的构建gradSrcX diff([imSrc(:,2:end,:) imSrc(:,1,:)], 1, 2); gradSrcY diff([imSrc(2:end,:,:); imSrc(1,:,:)], 1, 1); gradDstX diff([imDst(:,2:end,:) imDst(:,1,:)], 1, 2); gradDstY diff([imDst(2:end,:,:); imDst(1,:,:)], 1, 1); magSrc sqrt(gradSrcX.^2 gradSrcY.^2); magDst sqrt(gradDstX.^2 gradDstY.^2); useSrc magSrc magDst; gradMixX gradSrcX .* useSrc gradDstX .* ~useSrc; gradMixY gradSrcY .* useSrc gradDstY .* ~useSrc; lapMix div(gradMixX, gradMixY); % 需要用散度算子或从梯度重建拉普拉斯注意diff在MATLAB里的方向约定diff(X, 1, 2)算水平差分得到的矩阵会比原图少一列上面代码用循环填充的方式维持了尺寸。实际做这个变体时我更喜欢直接构建散度算子矩阵把梯度场转成散度项这样不会引入周期边界假设。散度算子是梯度算子的转置负值在离散网格上最终化简出来仍然是5点模板只不过每一项的系数是两边梯度的符号差。多区域融合的常见做法是如果源图里有多块不连通区域要分别融合不要整个mask一起求解。n太大时矩阵条件数变差且区域间如果有间隙泊松方程会把两块区域的边界也一起平滑这不是想要的。我的习惯是每个连通域单独调用一次主函数最后把结果合回一张图。如果多个区域位置相邻但语义不同比如把一只鸟和一架飞机同时贴到天空分开solve比一次solve在边界控制上更精确。验证多区域结果时只需要对每个连通域重复4.1的拉普拉斯误差测试。实践下来我发现一个经验mask区域边界外扩23个像素再截断融合结果会比原始mask边缘更干净因为边界约束的邻域更广方程对边界噪声更鲁棒。代价是融合区域比预期略微扩大如果对区域范围有严格限制就不能用这招。写到这里想起最早我在这套代码上翻车最狠的一次就是b的负号问题整个融合结果像被人糊了一层滤镜。后来养成一个习惯每次拿到新图片先跑矩形测试用例确认拉普拉斯误差在1e−6量级再去换真实掩膜。这个习惯帮我避开了至少七八次莫名其妙的色偏问题。泊松融合这个方向一旦理解成“梯度域的搬运工”代码量其实很少难的永远在边界处理和验证环节。希望帮到你。本文还有配套的精品资源点击获取
