简介：Perona-Malik 扩散模型的修改版 Matlab 代码，面向图像处理、计算机视觉方向的大学生与科研人员，适用于课程设计、期末大作业和毕业设计中的去噪与边缘保持实验。代码基于偏微分方程热扩散思想，在原始模型上调整了扩散函数与迭代策略，便于对比不同参数下的平滑和细节保留效果。压缩包共24个文件，以11个m脚本（含主程序 demo.m、modifiedPM.m 及 PSNR/SSIM 评估函数）、多张 tif/png 测试图（含原始图与噪声图）和1个说明 txt 为主，整体约2.55MB，结构紧凑、可直接运行。目前已有77人学习浏览，适合希望快速上手并深入理解 Perona-Malik 改进思路的 Matlab 用户。代码采用参数化编程，注释清晰，附赠案例数据，改动参数即可观察去噪与边缘保持的变化，是开展图像处理实验的实用工具。

1. Perona-Malik 扩散模型：为什么修改它的 Matlab 代码值得花时间

当高斯滤波把边缘一起抹平、中值滤波把纹理当成噪声干掉的时候，Perona-Malik 扩散模型是图像去噪里少有的“保边”方案：它在梯度小的地方使劲扩散，在梯度大的地方踩刹车。你拿到的这包“修改Matlab代码”，大概率就是有人把原版 PM 的两个扩散系数、梯度估计方式和停止条件改过一遍的产物。注意这里的“扩散模型”指的是偏微分方程意义上的各向异性扩散，跟 AI 绘画里那套 latent diffusion 不是一回事。读这篇文章的人，要么在赶 Matlab 图像处理大作业或毕设，要么手上有真实含噪图想找一个能保住边缘的去噪算法。我会从原版原理讲到修改方向，再给你一套能直接跑的 pm_denoise 函数，最后把参数和踩坑点都摊开讲，值不值得投入，你跑完就有数了。

2. 原版 PM 的数学行为与三个参数：K、λ、迭代次数分别控制什么

2.1 从热扩散到各向异性扩散：为什么均匀模糊会毁掉边缘

Perona-Malik 模型脱胎于经典热扩散方程：

$$\frac{\partial I}{\partial t} = \mathrm{div}(c \cdot \nabla I)$$

当扩散系数 $c$ 是常数时，这个方程的解等价于不断对图像做高斯模糊，而且随着 $t$ 增大，所有结构都被均匀抹掉。边缘和噪声在数学上都是高频信息，线性扩散分不清谁是谁，所以高斯滤波的代价必然是边缘模糊。

PM 的改动只有一处：把常数扩散系数换成梯度模长的函数 $g(|\nabla I|)$。

$$\frac{\partial I}{\partial t} = \mathrm{div}(g(|\nabla I|) \cdot \nabla I)$$

这样图像内部平坦区域 $|\nabla I|$ 小，扩散系数趋近于 1，去噪效率高；边缘处 $|\nabla I|$ 大，扩散系数趋近于 0，边界被“冻结”。这个思想非常直觉，也是 PM 能成为经典的原因。但数学上的优雅不等于数值实现容易，后面要调的东西还很多。

2.2 两个经典扩散系数函数与 K 值的含义

原版论文给出了两个扩散系数函数，现在几乎所有 Matlab 实现里都能看到：

$$g_1(x) = \exp\left(-\left(\frac{x}{K}\right)^2\right)$$

$$g_2(x) = \frac{1}{1 + \left(\frac{x}{K}\right)^2}$$

其中 $x = |\nabla I|$，$K$ 是对比度参数，它的物理含义很直白：当梯度模长接近 K 时，扩散系数在 $0.3$ 到 $0.5$ 左右，介于“扩散”和“冻结”之间；当梯度模长远大于 K 时，扩散系数几乎为 0，像素被保留。说白了，K 决定了你把“多陡的梯度”算作边缘。

离散化后，迭代公式长这样：

$$I_{n+1} = I_n + \lambda \cdot \left( c_N \nabla_N I + c_S \nabla_S I + c_E \nabla_E I + c_W \nabla_W I \right)$$

每个方向上的 $c$ 都由该方向的梯度模长算出来，这就是为什么代码里要分别计算四个方向的梯度。三个核心参数的取值范围如下表：

2.3 原版模型的两个硬伤，也就是“修改代码”的动机

第一个硬伤是数学上的不适定性。在特定梯度区间内，$g$ 函数可能让扩散产生“反向”效果，边缘会被锐化过头，甚至把噪声尖峰放大成假边缘。这两个经典扩散系数在实际离散网格上表现敏感，稍微调大 λ 图像就花了。

第二个硬伤更实际：梯度估计本身对噪声敏感。含噪图像的噪声点也有很大的梯度模长，PM 会把噪声误判为边缘，结果就是噪声去不掉。原来的 $g_1$ 函数收敛慢，$g_2$ 函数去噪快但边缘也容易被磨掉。正因为这两个硬伤，大家拿到手里的“修改版”代码几乎都在同一个位置动手：改扩散系数函数、改梯度估计方式、改 K 的取值方式、加停止准则。后面三章就是这四个方向的完整拆解。

3. 修改版 PM 的四个改动方向：替换扩散系数、预平滑梯度、自适应 K 与自动停止

3.1 用 Tukey 鲁棒函数替换经典扩散系数，抑制逆扩散

最常见的修改是换掉扩散系数函数。两个经典函数虽然简单，但逆扩散区间让数值行为不可控。工程里我更偏好 Tukey 函数：

function c = diffusion_coeff(g, K, func) % g: 梯度模长 |∇I| % K: 对比度参数 % func: 扩散系数类型 switch func case 'classic1' % 原版 exp 形式 c = exp(-(g/K).^2); case 'classic2' % 原版倒数形式 c = 1 ./ (1 + (g/K).^2); case 'tukey' % 鲁棒统计中的 Tukey 权重 c = zeros(size(g), 'like', g); mask = g <= K; c(mask) = (1 - (g(mask)/K).^2).^2; case 'charbonnier' % 凸势函数对应的系数 c = 1 ./ sqrt(1 + (g/K).^2); otherwise error('Unknown func: %s', func); end end

Tukey 函数的关键行为是：梯度模长一旦超过 K，扩散系数直接置 0，完全不扩散。这从数学上掐掉了逆扩散区间，大梯度不可能被二次放大。用它之后，图像迭代稳定性明显提升，边缘也不容易出现“振铃”。

3.2 Catté 正则化：进扩散系数前先对梯度做一次高斯预平滑

第二处高频修改是梯度估计的正则化。原版直接用原始图像算梯度，噪声点的梯度会被当成边缘。常见做法是先对当前图像做一次高斯平滑，再拿平滑后的图算梯度，这就是 Catté 在 1992 年提出的修正思路：

% 在每次迭代计算梯度之前，先对 u 做高斯预滤波 if opts.sigma > 0 u_smooth = imgaussfilt(u, opts.sigma); else u_smooth = u; end % 四方向差分，全部基于 u_smooth 计算 dI_N = u_smooth([1, 1:end-1], :) - u_smooth; dI_S = u_smooth([2:end, end], :) - u_smooth; dI_E = u_smooth(:, [2:end, end]) - u_smooth; dI_W = u_smooth(:, [1, 1:end-1]) - u_smooth; % 扩散系数用带符号差分的绝对值来计算 cN = diffusion_coeff(abs(dI_N), K, 'tukey'); cS = diffusion_coeff(abs(dI_S), K, 'tukey'); cE = diffusion_coeff(abs(dI_E), K, 'tukey'); cW = diffusion_coeff(abs(dI_W), K, 'tukey');

这里只对梯度估计做平滑，平滑图本身不直接参与迭代混合，所以不会把整幅图模糊掉。σ 一般取 0.5～2 像素。取太小没用，取太大则会把真实边缘的梯度也压下去，造成最后边缘变钝。一个经验值是：噪声强的时候 σ 往 1.5 靠，噪声弱的时候 σ 设 0.5 左右。

3.3 用 MAD 噪声估计自动算 K 值，不再手填魔法数字

K 值是 PM 模型里最玄学的一个参数。同一段代码换个亮度分布，固定 K 的效果天差地别。修改版通常会用中值绝对偏差法自动估计噪声水平，再换算成 K：

function K = estimate_K(u) % 用梯度模长的 MAD 估计噪声标准差，再换算为 K [Gx, Gy] = imgradientxy(u); G = sqrt(Gx.^2 + Gy.^2); G_med = median(G(:)); noise_sigma = median(abs(G(:) - G_med)) / 0.6745; K = 1.5 * noise_sigma; % 乘 1.5 或 2 都可，视保留边缘强度而定 if K < 0.01 K = 0.01; % 下限保护，防止除零 end end

0.6745 是高斯分布下 MAD 到标准差的换算系数，这个做法来自小波去噪的 Donoho 阈值思想。如果你的环境没有imgradientxy，可以用内置gradient函数代替，得到的 K 值差别不大。K 设为噪声标准差的 1.5 倍左右，基本保证噪声梯度被平滑掉，而真实边缘的梯度因为更大会被保留。

3.4 用 RMS 变化量做停止准则，让迭代次数自己决定

固定迭代次数去跑 PM，容易出现“看着还行的时候已经过了最佳点”。修改代码里往往加一个停止准则：当相邻两次迭代的像素 RMS 变化小于阈值，认为扩散已经收敛，直接退循环。

for t = 1:iter_max u_prev = u; % ... 执行一次完整迭代 ... rms_change = sqrt(mean((u(:) - u_prev(:)).^2)); if rms_change < stop_tol fprintf('迭代 %d 步后收敛，RMS 变化 = %.2e\n', t, rms_change); break; end end

stop_tol 在图像归一化到 0～1 的情况下，取 1e-3 比较合适。想要更精细就取 1e-4，迭代次数会增加，视觉效果提升有限。这个机制比固定迭代次数更实用，因为不同图像收敛速度差很多：纹理多的图收敛慢，平坦图几步就停了。

4. 把修改版跑起来：pm_denoise 主函数、四方向差分与调用示例

4.1 pm_denoise 函数完整结构：参数入口与迭代主循环

这类压缩包解开后，常见结构是一个主函数加几个子函数。下面这套代码融合了前两章提到的全部四个修改方向，你可以直接当成修改版 PM 的主干来用：

function out = pm_denoise(img, opts) % pm_denoise Perona-Malik 扩散模型修改版 % 改动点: % 1. diffusion_coeff 支持 tukey / charbonnier % 2. Catté 正则化: 梯度估计前高斯预平滑 % 3. 自动 K 值: 基于梯度 MAD 的噪声估计 % 4. 停止准则: RMS 变化量低于阈值自动收敛 % % 输入: % img - 灰度图像，数值范围任意 % opts - 结构体参数，可选字段见下表 % 输出: % out - 去噪后的 double 图像，范围 0~1 % 默认参数 if nargin < 2 || isempty(opts) opts = struct(); end lambda = get_opt(opts, 'lambda', 0.2); % 时间步长，<=0.25 iter_max = get_opt(opts, 'iter_max', 200); % 最大迭代次数 K_mode = get_opt(opts, 'K_mode', 'auto'); % 'auto' 或 'manual' K_manual = get_opt(opts, 'K', 0.15); % manual 模式下的固定 K sigma = get_opt(opts, 'sigma', 1.0); % Catté 预平滑标准差 func = get_opt(opts, 'func', 'tukey'); % 扩散系数函数类型 stop_tol = get_opt(opts, 'stop_tol', 1e-3); % 停止阈值 u = im2double(img); if strcmp(K_mode, 'auto') K = estimate_K(u); fprintf('自动 K = %.4f\n', K); else K = K_manual; end for t = 1:iter_max u_prev = u; % Catté 正则化：只用于梯度估计 if sigma > 0 u_smooth = imgaussfilt(u, sigma); else u_smooth = u; end % 四方向带符号差分 dI_N = u_smooth([1, 1:end-1], :) - u_smooth; dI_S = u_smooth([2:end, end], :) - u_smooth; dI_E = u_smooth(:, [2:end, end]) - u_smooth; dI_W = u_smooth(:, [1, 1:end-1]) - u_smooth; % 各方向扩散系数 cN = diffusion_coeff(abs(dI_N), K, func); cS = diffusion_coeff(abs(dI_S), K, func); cE = diffusion_coeff(abs(dI_E), K, func); cW = diffusion_coeff(abs(dI_W), K, func); % 显式离散化主迭代 u = u + lambda * (cN .* dI_N + cS .* dI_S + cE .* dI_E + cW .* dI_W); % 停止准则 rms_change = sqrt(mean((u(:) - u_prev(:)).^2)); if rms_change < stop_tol fprintf('第 %d 步收敛, RMS=%.2e\n', t, rms_change); break; end end out = u; end

迭代主循环里的u = u + lambda * (cN .* dI_N + ...)就是离散化散度的展开式。每个方向的扩散系数 c 与该方向的差分相乘，再四方向求和，这一步等价于计算 $div(g(|\nabla I|)\nabla I)$。lambda 必须小于等于 0.25，否则显式格式不稳定，图像会发散成雪花点。

子函数get_opt和diffusion_coeff在前一章已经给过，estimate_K也一并放进同一个 m 文件即可。如果你机器上没有 Image Processing Toolbox，imgaussfilt可换成imfilter(u, fspecial('gaussian', [5 5], sigma))，结果几乎一致。

4.2 四方向差分与边界复制：这里最容易写错

四方向差分用的是“索引复制”而不是circshift，这是有讲究的。circshift是循环移位，图像左边缘会拿到右边缘的像素，对于自然图像来说会在边界产生虚假梯度。索引复制的写法是：

dI_N = u_smooth([1, 1:end-1], :) - u_smooth; dI_S = u_smooth([2:end, end], :) - u_smooth; dI_E = u_smooth(:, [2:end, end]) - u_smooth; dI_W = u_smooth(:, [1, 1:end-1]) - u_smooth;

以 N 方向为例，[1, 1:end-1]把第一行复制了一遍，形成“第一行外还有一个第一行”的假想邻居。此时dI_N在第一行为 0，相当于 Neumann 边界条件：边界外侧没有梯度通量。不需要额外写边界处理代码，也不用担心 NaN 或越界。

另外注意扩散系数必须用梯度模长abs(dI_N)，而迭代公式用的是带符号差分。两者不要混用，否则等号右边的扩散系数方向会乱掉。新手翻车往往是abs位置放错，导致扩散方向反转，图像越迭代越脏。

4.3 调用 Demo：sigma=0.1 高斯噪声下的参数配置

跑通最小示例很简单：

% 读入标准测试图并加高斯噪声 I = im2double(imread('cameraman.tif')); rng(2025); In = imnoise(I, 'gaussian', 0, 0.01); % 标准差约 0.1 % 修改版 PM 参数 opts = struct(); opts.lambda = 0.2; opts.iter_max = 200; opts.K_mode = 'auto'; opts.sigma = 1.0; opts.func = 'tukey'; opts.stop_tol = 1e-3; out = pm_denoise(In, opts); fprintf('PSNR: 噪声图 %.2f dB, 去噪后 %.2f dB\n', psnr(In, I), psnr(out, I)); fprintf('SSIM: 噪声图 %.3f, 去噪后 %.3f\n', ssim(In, I), ssim(out, I)); % 显示对比 subplot(1,3,1); imshow(I); title('原图'); subplot(1,3,2); imshow(In); title('噪声图'); subplot(1,3,3); imshow(out); title('PM 修改版去噪');

加噪时imnoise(I,'gaussian',0,0.01)的最后一个参数是方差，所以实际噪声标准差是 0.1，不是 0.01。这是 Matlab 里最容易看错的细节之一。如果噪声更强，比如标准差 0.2，建议把sigma提到 1.5，stop_tol保持 1e-3，迭代次数上限 300。

5. 常见问题与排查：五个反复出现的坑及解决办法

5.1 图像出现块状伪影/油彩感

现象：平滑区域不是细腻过渡，而是变成一块一块的色块或棋盘格。

原因：lambda 取值过大，显式迭代进入不稳定区；或者停止阈值放太宽，迭代在震荡区间提前退出，留下振荡残影。经典求稳让许多新手把 lambda 调到 0.25 以上，结果就是块状伪影。

解决：先把 lambda 降到 0.15～0.2 再跑。如果伪影消失，说明就是步长问题。若 lambda 已经小于 0.2 还有色块，查一下扩散系数，把func从'classic1'换成'tukey'，逆扩散区间消失后块状伪影通常也会跟着消失。

5.2 像素值发散，迭代过程爆炸

现象：迭代几步后，图像出现单个亮点或整片高亮，像素值超出 0～1 范围，甚至变成 NaN。

原因：最直接的是 lambda 大于 0.25，显式格式不稳定。另一类是输入图像没有经过im2double，还是 uint8，迭代时像素值溢出。扩散系数选 classic1/classic2 时逆扩散放大噪声尖峰，也会导致局部爆炸。

解决：检查输入是 double 且范围在 0～1；lambda 强制限制到 0.2 左右；扩散系数优先选 tukey，它在梯度超过 K 时直接置 0，没有逆扩散区间。改完这三处，爆炸问题基本不会出现。

5.3 边缘被磨平，去噪图看起来像高斯模糊

现象：噪声没了，但边缘也钝了，人物轮廓发虚。

原因：K 值太大，把大部分真实边缘梯度当成了噪声背景，扩散系数在边缘处依然很大；或者 Catté 预平滑的 sigma 太大，把边缘梯度在估计阶段就削平了。

解决：换成K_mode='auto'，让 MAD 估计接管 K 值。如果手动设 K，从噪声标准差的 1 倍起调。sigma 从 0.8 往下试，0.5 起步更保守。记住一个原则：Catté 预平滑只影响梯度估计，不代表要对整个图像做大尺度模糊，sigma 超过 2 多半就是边缘变糊的元凶。

5.4 运行极慢，小图都要跑几十秒

现象：256×256 的灰度图，200 次迭代跑起来像卡住一样。

原因：循环内反复调用imgradientxy或gradient，每次都会创建完整梯度矩阵；或者输出函数每步都打印调试信息；也有把扩散系数的计算写成多重 for 循环的版本，速度当然上不去。

解决：把梯度估计改成索引差分，就是 4.2 节那四行代码，不调用gradient函数。预分配变量、避免循环内fprintf。如果 512×512 的图还是要跑几十秒，把iter_max从 500 降到 200，加停止准则让它在达到稳态后提前退出。绝大多数情况下，停止准则能让实际迭代数比上限少一半。

5.5 整体亮度漂移，均值不守恒

现象：去噪前后图像平均灰度明显变化，直方图整体左移或右移。

原因：PM 理论上是均值保持的，但离散化时边界方向差分不对称、扩散系数不连续，会引入系统性偏差。边界复制方式不一致也会造成亮度泄漏。

解决：检查四个方向的边界索引是否都采用同一种复制策略，不要 N 向用复制、E 向用零填充。如果仍然漂移，在每次迭代结束后补一个均值校准：

u = u + (mean(u_prev(:)) - mean(u(:)));

这相当于把扩散过程里丢失的直流分量拉回来。轻微漂移时加这行就能解决，但当你发现需要频繁校准，还是回头检查边界的写法和梯度符号。

6. 验证修改是否有效：PSNR、SSIM 与 EPI 组合判断

6.1 三个指标各看什么

PSNR 衡量像素级误差，SSIM 衡量结构和亮度相似度，但它们对“边缘是否尖锐”并不敏感。这里补一个边缘保持指数 EPI：

gx_noisy = imgradientxy(In); gx_den = imgradientxy(out); epi = sum(sqrt(gx_den.x.^2 + gx_den.y.^2), 'all') / ... sum(sqrt(gx_noisy.x.^2 + gx_noisy.y.^2), 'all');

EPI 接近 1，代表去噪后梯度能量几乎没损失。若 PSNR 上升了但 EPI 掉到 0.6 以下，说明算法在靠糊边缘换噪声降低，不符合 PM 的本意。修改版一般希望 EPI 在 0.85 以上。

6.2 局部放大与残差图：看肉眼效果

数值指标之外一定要看残差图：imshowpair(I, out, 'diff')，残差里若有大片边缘轮廓，说明边缘区域被过度修改。常规做法是截取一个人物脸部或建筑边缘的局部区域，放大后对比原图和去噪图，观察边缘两侧的灰度过渡是否仍然陡峭。

6.3 一个脚本把原版和修改版放一起比对

% 原版 PM：固定 K，经典扩散系数，固定迭代次数 opts_orig = struct('lambda', 0.2, 'iter_max', 50, ... 'K_mode', 'manual', 'K', 0.1, ... 'sigma', 0, 'func', 'classic2'); out_orig = pm_denoise(In, opts_orig); % 修改版 PM：自动 K，Tukey 系数，Catté 预平滑 opts_mod = struct('lambda', 0.2, 'iter_max', 200, ... 'K_mode', 'auto', 'sigma', 1.0, ... 'func', 'tukey', 'stop_tol', 1e-3); out_mod = pm_denoise(In, opts_mod); fprintf('原版 : PSNR=%.2f dB SSIM=%.3f\n', psnr(out_orig, I), ssim(out_orig, I)); fprintf('修改版: PSNR=%.2f dB SSIM=%.3f\n', psnr(out_mod, I), ssim(out_mod, I));

我自己的习惯是每次都把修改版与原版、再和一幅高斯滤波结果做三方对比。高斯滤波保证平滑，原版 PM 代表经典方案，修改版看提升。曾经有一版代码 PSNR 高出原版 0.5 dB，但局部放大后边缘出现了轻微光晕，最后是 EPI 指标把问题暴露出来的。调 PM 这类偏微分方程算法，单看数字容易误判，边缘指标和目视检查必须和 PSNR 一起用。希望这套验证流程能帮你在自己的数据上少走弯路。

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