简介：这是一套基于MATLAB的泊松融合与泊松修复实现代码，复现了Pérez等人于2003年提出的经典泊松图像编辑方法，适合图像处理学习者、算法复现者以及需要做无缝拼接或图像修复的开发者参考。压缩包内共两个MATLAB脚本及若干练习图片，整体仅213KB，分别实现泊松融合和泊松修复，可直接运行并对比效果。代码核心流程完整，涵盖梯度场计算、引导向量场构造、边界条件处理及泊松方程求解等关键步骤，结构清晰便于二次开发。项目附带的练习图片覆盖典型融合与修复场景，方便读者对照结果验证算法。目前已有711人学习下载。通过阅读和调试，读者能快速理解泊松方程在图像编辑中的作用，掌握从掩膜选取到融合结果输出的完整链路，也可基于此扩展实现更多图像处理算法。

1. 泊松融合是什么：让拼接痕迹凭空消失的“无缝换肤”黑匣子

把一张贴纸嵌进风景照片，边缘却看不到任何拼接痕迹，这不是磨皮滤镜的玄学，而是泊松融合（Poisson Fusion）的经典应用。这个名为MATLAB_Poisson_Fusion_HJT的项目，是一套用MATLAB实现的泊松方程求解器，输入源图、目标图与掩膜，就能输出平滑过渡的融合结果。它适合图像处理初学者理解梯度域编辑原理，也适合研究者快速跑通基线算法。核心思想不拷贝源图像素，而是拷贝源图像的“梯度”，让融合区域的纹理变化与源图一致，同时边界像素与目标图完全相等，从而抹掉生硬过渡。

2. 项目架构与技术选型：从梯度场到线性方程组的落地设计

2.1 技术选型：MATLAB为什么适合泊松方程求解

MATLAB的强项是矩阵运算，泊松融合最终要解一个大规模稀疏线性系统。在C++里你可能要接Eigen或者自己写稀疏存储，而在MATLAB里，sparse矩阵配合A\b左除运算符，一行代码就能完成求解。图像在MATLAB中天然就是二维数组，imread读进来可以直接参与数值计算，掩膜可以用逻辑矩阵表示，坐标映射用sub2ind和ind2sub处理，省去了大量格式转换代码。

因此这套代码把精力集中在“如何把泊松方程离散成矩阵”上，而不是重复造求解器的轮子。对于一张512×512的图像，掩膜区域可能包含几万到十几万个未知数，MATLAB的稀疏直接法在普通笔记本上几秒内就能完成。如果目标是验证算法效果，而不是部署到生产环境，MATLAB是这个场景下性价比最高的选择。

2.2 目录结构与函数职责：每个文件负责什么

一个合格的MATLAB泊松融合工程，通常会把“跑通演示”和“核心算法”分开。参考这份项目代码的组织方式，典型目录如下：

MATLAB_Poisson_Fusion_HJT/ ├── run_demo.m % 主入口脚本，演示完整融合流程 ├── poisson_fusion.m % 核心融合函数，接收源图/目标图/掩膜 ├── get_laplacian_matrix.m % 构建五点拉普拉斯稀疏矩阵 ├── compute_gradient_field.m % 计算源图像引导梯度场 ├── build_right_hand_side.m % 组装泊松方程右侧向量 ├── images/ │ ├── source.jpg % 源图像（要迁移的内容） │ ├── target.jpg % 目标图像（背景） │ └── mask.png % 掩膜，白色区域为融合区域 └── README.md

run_demo.m负责读图、调用核心函数、显示结果；poisson_fusion.m内部依次调用后面三个子函数，完成梯度计算、矩阵组装、方程求解和像素回写。这样的分层让二次开发变得直接：要改融合逻辑，只需要动poisson_fusion.m；要换求解法，改最后一步的A\b为pcg即可。

2.3 函数接口设计：输入输出与易用性

核心函数poisson_fusion的接口设计很干净，三个输入一个输出，不暴露中间变量。函数签名和参数说明如下：

function result = poisson_fusion(source, target, mask) % source: [rows, cols, 3] double，值域[0,1] % target: [rows, cols, 3] double % mask: [rows, cols] logical，1表示重建区域，0保持原样 % result: double类型图像，尺寸与source相同

用im2double把uint8图像转成double，是调用前的必备步骤。因为uint8做减法会溢出，比如10减20会折回245，梯度场彻底算错。mask必须是逻辑矩阵，如果读进来是uint8，要用mask = imread(...) > 0转换。接口层不检查这些，所以主脚本里最好自己加校验，避免把错误传给核心函数。

2.4 核心流程拆解：五步完成的梯度域编辑

整个融合流程可以拆成五步：第一步，提取源图的引导梯度场，通常用前向差分；第二步，根据掩膜找出所有需要参与求解的像素，建立从像素线性索引到方程未知量编号的映射；第三步，在掩膜区域内离散化泊松方程，得到稀疏系数矩阵和右侧向量；第四步，对RGB三个通道分别求解稀疏线性方程组；第五步，把解写回目标图对应像素位置。

用一句话记忆这个链路：它不是直接拷贝源图的像素值，而是拷贝源图的“梯度”，再通过求解拉普拉斯方程把梯度场还原成像素值。边界上的像素因为被硬编码为目标图的值，所以融合区域和目标背景之间不存在任何人为的权重混合，边缘自然连续。

3. 核心代码拆解：稀疏矩阵、边界条件与求解的完整链路

3.1 从掩膜到未知量编号：find与map的配对

图像有几十万像素，但我们只对掩膜内的子集求解。第一步是把掩膜内像素的位置和线性索引取出来，并为它们建立稀疏编号。

% 输入: mask - 逻辑矩阵 (rows*cols)，1表示融合区域 % 输出: map - 全图长度向量，把像素线性索引映射到方程未知量编号 function map = buildMapping(mask) rows = size(mask,1); cols = size(mask,2); pixIdx = find(mask); % 掩膜内所有像素的线性索引 n = numel(pixIdx); % 未知量总数 map = zeros(rows*cols,1); % 全图映射 map(pixIdx) = 1:n; % 掩膜内编号为1..n，其他为0 end

find返回的是单索引，也就是“从左上角开始逐列扫描”的编号。map的作用是把像素坐标与方程未知量对应起来，后面构建稀疏矩阵时，需要知道某个像素的邻居是否也在掩膜内，在的话就在它的系数行里加1。这一步看着简单，但漏掉map就会在组装矩阵时把索引搞乱。

3.2 引导场计算：用前向差分代替连续微分

泊松方程右侧的散度项来自源图像的梯度场。在离散域里，我们用像素差代替导数，常见做法是前向差分：

function [gradX, gradY] = computeGradient(img) % img: double图像，尺寸 [rows, cols] % gradX: 水平梯度（沿列方向） % gradY: 垂直梯度（沿行方向） gradX = diff(img, 1, 2); % 每行内后一列减前一列 gradY = diff(img, 1, 1); % 每列内后一行减前一行 % 为了与图像尺寸一致，在右侧和下侧补零 gradX = cat(2, gradX, zeros(size(img,1),1)); gradY = cat(1, gradY, zeros(1,size(img,2))); end

diff(img,1,2)是水平方向差分，结果比原图少一列，所以补零。补零相当于假设最右侧和最下侧像素的梯度为0，这对融合效果几乎没有影响，因为那些像素通常不在掩膜内。实际工程里，也可以把gradX和gradY先限定在掩膜内部再计算，减少无效计算。

3.3 泊松方程离散化：五点拉普拉斯模板与稀疏矩阵组装

二维拉普拉斯算子的离散形式是：当前像素乘以-4，上下左右四个邻居各乘以1。这一模板在掩膜内逐点应用，就可以得到线性方程组。下面的代码用直观的循环实现，虽然效率不是最高，但每个系数怎么来的看得清清楚楚。

function A = buildLaplacianMatrix(mask, map) rows = size(mask,1); cols = size(mask,2); n = sum(mask(:)); A = sparse(n, n); % 预分配稀疏矩阵 % 邻居偏移：左右上下 neighbors = [0 -1; 0 1; -1 0; 1 0]; for i = 1:rows for j = 1:cols if ~mask(i,j) continue; % 跳过非融合区域 end varIdx = map(sub2ind([rows, cols], i, j)); % 当前未知量编号 A(varIdx, varIdx) = -4; % 对角线系数 for k = 1:4 ni = i + neighbors(k,1); nj = j + neighbors(k,2); if ni < 1 || ni > rows || nj < 1 || nj > cols continue; % 越界邻居，跳过（由边界条件处理） end if mask(ni,nj) % 邻居也是未知量，系数为1 nbIdx = map(sub2ind([rows, cols], ni, nj)); A(varIdx, nbIdx) = 1; end % 如果邻居不在掩膜内，它的值属于右侧常数项，这里不处理 end end end end

neighbors数组定义四个方向偏移。sub2ind把行列坐标转成线性索引。这个双重循环在演示项目中完全够用，但如果你要处理4K图像，一定要改成向量化构建，否则单是组装矩阵就要好几分钟。向量化思路是把所有掩膜内像素的坐标批量取出来，用repmat构造三元组，再一次性交给sparse。

3.4 右侧向量：源图散度加上边界像素值

方程组的形式是所有掩膜内像素的系数乘以未知数等于右侧项。右侧项来自两部分：一是源图像梯度的散度，二是当邻居在掩膜外时，那个邻居的目标图像值作为已知常数移到右边。

function b = buildRightHandSide(source, target, mask, map, gradX, gradY) rows = size(mask,1); cols = size(mask,2); n = sum(mask(:)); b = zeros(n, 1); for i = 1:rows for j = 1:cols if ~mask(i,j) continue; end varIdx = map(sub2ind([rows, cols], i, j)); % 离散散度 = gradX(i,j) - gradX(i,j-1) + gradY(i,j) - gradY(i-1,j) % 这里用前向差分近似，边界位置用条件乘法保护 div = gradX(i,j) - (gradX(i,j-1) * (j>1)) ... + gradY(i,j) - (gradY(i-1,j) * (i>1)); b(varIdx) = div; % 边界像素处理：若邻居在掩膜外，把目标图像值移入右侧 neighbors = [0 -1; 0 1; -1 0; 1 0]; for k = 1:4 ni = i + neighbors(k,1); nj = j + neighbors(k,2); if ni < 1 || ni > rows || nj < 1 || nj > cols || mask(ni,nj) continue; end b(varIdx) = b(varIdx) + target(ni, nj); % 系数为1的边界贡献 end end end end

这里容易踩坑的是j=1时gradX(i,j-1)越界，所以乘了(j>1)作为保护。实际更稳妥的写法是先把gradX在右边补零，再用if显式判断，避免条件乘法让阅读者费解。散度符号必须与矩阵对角线符号保持统一，否则解出来的是“反扩散”，图像会变成噪声。

3.5 求解与回写：MATLAB左除与双精度图像

核心方程组装完毕后，求解只在最后一行。A\b是MATLAB里最典型的稀疏线性求解，对于对称正定的泊松矩阵，它内部会走Cholesky分解。如果矩阵规模很大，也可以替换为x = pcg(A, b, 1e-6, 1000)。

% 融合函数的核心片段 function result = poisson_fusion(source, target, mask) % source/target: double类型图像，值范围[0,1] % mask: 逻辑矩阵 result = target; % 先复制目标图，融合区域会被覆盖 for ch = 1:3 % 取单通道 src = source(:,:,ch); tgt = target(:,:,ch); [gradX, gradY] = computeGradient(src); map = buildMapping(mask); A = buildLaplacianMatrix(mask, map); b = buildRightHandSide(src, tgt, mask, map, gradX, gradY); x = A \ b; % 核心求解：MATLAB自动选择稀疏直接法 % 把解写回目标图对应像素 pixIdx = find(mask); result(:,:,ch) = target(:,:,ch); result(pixIdx) = x; end end

回写时注意result(pixIdx)是线性索引，MATLAB会按列填充，正好与求解时的编号顺序一致。由于三个通道独立求解，如果源图和目标图整体色温差异大，融合区域可能偏色；这个问题的后续处理我会在第五章讲到。

4. 跑通项目：从加载图像到验收融合效果

4.1 环境准备：版本、工具箱与依赖

这份代码依赖MATLAB的图像处理工具箱（Image Processing Toolbox），因为imread、imshow、roipoly都来自那里。实际上sub2ind是MATLAB基础函数，但图像输出用了imshow，所以工具箱还是需要的。建议使用R2018a或更高版本，低版本在稀疏矩阵运算上性能差异不会太明显，但新版对sparse的并行处理更好。

4.2 最小复现：一段脚本完成“源图+掩膜+目标图”融合

有了核心函数，跑通一个最小示例只要几行。这里给出可以复制的完整脚本：

% run_demo.m % 读取输入图像，注意统一尺寸 source = im2double(imread('images/source.jpg')); target = im2double(imread('images/target.jpg')); mask = imread('images/mask.png') > 0; % 转逻辑型 % 确保掩膜与图像尺寸相同，否则会报错 if ~isequal(size(mask), size(source(:,:,1))) error('掩膜尺寸与源图不一致'); end % 执行泊松融合 result = poisson_fusion(source, target, mask); % 展示并保存 subplot(1,3,1); imshow(source); title('源图'); subplot(1,3,2); imshow(target); title('目标图'); subplot(1,3,3); imshow(result); title('泊松融合'); imwrite(result, 'images/result.png');

im2double把uint8图像转到0到1的double范围，这一步很关键。如果直接用uint8参与减法，溢出时会折返，比如10减20会变成245，整个梯度场就完全错了。所以一定记得先转double。如果源图和目标图尺寸不一致，也需要提前用imresize统一。

4.3 效果对比：直接拷贝、羽化混合与本项目的差别

拿到结果后，要验证它是否有效。你可以用同样的源图和掩膜，分别用下面三种方式处理：直接拷贝，将源图掩膜内像素直接覆盖到目标图，边界会有明显的台阶；线性羽化，对掩膜边界做高斯模糊权重混合，边界柔和但内部纹理可能发虚；泊松融合，保持源图梯度场，边界无缝，内部细节清晰。

在实际跑这个项目时，最直观的感受是：泊松融合的结果不是两者像素的加权平均，而是重新“画”了一个区域，像素值来自求解方程，所以参数选择不当时效果可能显得“过于平滑”，但算法本身确实能抹平边界。

4.4 参数调整：掩膜尺寸、边界过渡带与多通道一致性

掩膜的大小直接影响未知量数量。掩膜越大，矩阵越大，求解越慢，但融合内容更完整；掩膜边缘如果紧贴内容边界，容易出现半透明感。另外，RGB三个通道是独立求解的，如果源图和目标图的整体色温差异大，融合区域可能会偏色。常见做法是保留一个全局颜色校正项，也就是在泊松方程里加入一个额外的常数修正，让融合区域的平均颜色更接近源图。

我一般会在poisson_fusion外面套一层颜色调整，先计算源图和目标图在掩膜区域的平均亮度差值，然后对源图整体偏移该差值，再做融合。这样可以在不改变梯度结构的前提下，减少色差导致的偏色问题。

5. 避坑与常见问题排查：MATLAB泊松融合的血泪经验

5.1 报错“矩阵维度不一致”或“索引超出范围”

这是初学者最常翻车的地方。现象是运行run_demo.m时，mask与source尺寸不一致，或者sub2ind坐标越界，MATLAB直接抛异常。原因通常是掩膜是手动绘制的，没有和图像对齐。解决方法是：在调用前显式检查尺寸，并且把掩膜从uint8转成logical，最稳妥的写法是：

mask = imresize(mask, [rows, cols]) > 0;

5.2 融合结果出现黑色网格或条纹

现象是结果图像上有一层规则的暗线或亮线，像是“网格纸”。原因基本可以锁定在稀疏矩阵组装时系数符号不对，或者散度计算用了不匹配的差分方向。我调试时吃过一次亏：把buildLaplacianMatrix里对角线写成4、邻居写成-1，结果整个系统变成“反扩散”，解出来的图像充满噪声。

解决办法是做一个零测试：让source和target都等于一张灰度渐变图，掩膜任意，理想结果应该与目标图一致，因为梯度场没变化。如果结果有网格，检查buildLaplacianMatrix中的正负号。

5.3 求解速度慢到无法接受，内存耗尽

现象是300×300的掩膜都要算很久，再大直接内存不足。原因是双重循环逐像素组装矩阵，且使用sparse(n,n)动态填充，效率极低。解决思路有两条：第一，向量化构建矩阵。用repmat构造三元组，再一次性交给sparse，速度能提升两个数量级。第二，用pcg迭代求解器替代直接法。

opts.tol = 1e-8; opts.maxit = 20000; x = pcg(A, b, opts.tol, opts.maxit);

5.4 融合结果整体偏灰或亮度偏移

现象是融合区域看起来“褪色”，像蒙了一层灰。原因通常是RGB通道独立求解后，像素值被饱和到[0,1]之外，或者源图与目标图的整体亮度差异较大。解决方法是在回写时对result做一次clamp，并且可以在求解前把源图减掉其掩膜内平均亮度，与目标图掩膜内平均亮度对齐。

在实际项目里，很多开源实现会在梯度场上乘一个缩放系数，用经验值校准。如果你发现偏灰，试着把源图整体亮度拉高10%再做融合，往往能改善视觉效果。

6. 二次开发思路：把泊松融合改造成你自己的工具箱

6.1 给融合函数增加通道并行

原函数用for ch = 1:3逐通道求解，简单但低效。你可以用parfor替换外层循环，前提是每个通道独立，不存在数据竞争。注意parfor里不要共用result变量，可以用临时cell数组收集三通道结果，最后合并。这样在多核CPU上能获得接近两倍的提速。

6.2 用roipoly交互式生成掩膜

把固定掩膜文件换成用户交互绘制，能极大提升可用性。MATLAB的roipoly函数支持在图像窗口画多边形：

figure; imshow(target); mask = roipoly(); % 用户在图上手动圈选区域

得到掩膜后，直接传给poisson_fusion即可。结合imrect加createMask，还可以做矩形选择，适合快速抠图工具。

6.3 一个验证数值正确性的技巧

在你改造矩阵组装逻辑之后，强烈建议做一次解析解校验：构造一个已知的线性渐变图，比如source = repmat(linspace(0,1,cols), rows, 1)，然后让目标图是常数0，掩膜是整个图。理论上泊松融合解出的结果应该等于源图，因为梯度场完全一致。如果数值偏差小于1e-6，说明你的矩阵构建和右侧向量都是对的。我用这个方法排除过很多“看起来正常其实坐标错位”的bug，从那以后每次改动我都会先跑一遍最小示例，再改矩阵和求解器，最后用渐变图验数值。这套流程让我少走了好几次弯路，希望帮到你。

本文还有配套的精品资源，点击获取