简介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_{n1} 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$ 都由该方向的梯度模长算出来这就是为什么代码里要分别计算四个方向的梯度。三个核心参数的取值范围如下表参数常见取值范围作用K0.050.3图像归一化到 01边缘阈值越大越平滑边缘保留越弱λ0.10.25单步迭代步长超过 0.25 数值不稳迭代次数50500总平滑量越多越模糊但边缘仍会保留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 endTukey 函数的关键行为是梯度模长一旦超过 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.52 像素。取太小没用取太大则会把真实边缘的梯度也压下去造成最后边缘变钝。一个经验值是噪声强的时候 σ 往 1.5 靠噪声弱的时候 σ 设 0.5 左右。3.3 用 MAD 噪声估计自动算 K 值不再手填魔法数字K 值是 PM 模型里最玄学的一个参数。同一段代码换个亮度分布固定 K 的效果天差地别。修改版通常会用中值绝对偏差法自动估计噪声水平再换算成 Kfunction 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 end0.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 endstop_tol 在图像归一化到 01 的情况下取 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 Toolboximgaussfilt可换成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 调用 Demosigma0.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.5stop_tol保持 1e-3迭代次数上限 300。5. 常见问题与排查五个反复出现的坑及解决办法5.1 图像出现块状伪影/油彩感现象平滑区域不是细腻过渡而是变成一块一块的色块或棋盘格。原因lambda 取值过大显式迭代进入不稳定区或者停止阈值放太宽迭代在震荡区间提前退出留下振荡残影。经典求稳让许多新手把 lambda 调到 0.25 以上结果就是块状伪影。解决先把 lambda 降到 0.150.2 再跑。如果伪影消失说明就是步长问题。若 lambda 已经小于 0.2 还有色块查一下扩散系数把func从classic1换成tukey逆扩散区间消失后块状伪影通常也会跟着消失。5.2 像素值发散迭代过程爆炸现象迭代几步后图像出现单个亮点或整片高亮像素值超出 01 范围甚至变成 NaN。原因最直接的是 lambda 大于 0.25显式格式不稳定。另一类是输入图像没有经过im2double还是 uint8迭代时像素值溢出。扩散系数选 classic1/classic2 时逆扩散放大噪声尖峰也会导致局部爆炸。解决检查输入是 double 且范围在 01lambda 强制限制到 0.2 左右扩散系数优先选 tukey它在梯度超过 K 时直接置 0没有逆扩散区间。改完这三处爆炸问题基本不会出现。5.3 边缘被磨平去噪图看起来像高斯模糊现象噪声没了但边缘也钝了人物轮廓发虚。原因K 值太大把大部分真实边缘梯度当成了噪声背景扩散系数在边缘处依然很大或者 Catté 预平滑的 sigma 太大把边缘梯度在估计阶段就削平了。解决换成K_modeauto让 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 衡量结构和亮度相似度但它们对“边缘是否尖锐”并不敏感。这里补一个边缘保持指数 EPIgx_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自动 KTukey 系数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 一起用。希望这套验证流程能帮你在自己的数据上少走弯路。本文还有配套的精品资源点击获取 Perona-Malik扩散模型Matlab修改版：保边去噪算法与实现
简介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_{n1} 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$ 都由该方向的梯度模长算出来这就是为什么代码里要分别计算四个方向的梯度。三个核心参数的取值范围如下表参数常见取值范围作用K0.050.3图像归一化到 01边缘阈值越大越平滑边缘保留越弱λ0.10.25单步迭代步长超过 0.25 数值不稳迭代次数50500总平滑量越多越模糊但边缘仍会保留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 endTukey 函数的关键行为是梯度模长一旦超过 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.52 像素。取太小没用取太大则会把真实边缘的梯度也压下去造成最后边缘变钝。一个经验值是噪声强的时候 σ 往 1.5 靠噪声弱的时候 σ 设 0.5 左右。3.3 用 MAD 噪声估计自动算 K 值不再手填魔法数字K 值是 PM 模型里最玄学的一个参数。同一段代码换个亮度分布固定 K 的效果天差地别。修改版通常会用中值绝对偏差法自动估计噪声水平再换算成 Kfunction 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 end0.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 endstop_tol 在图像归一化到 01 的情况下取 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 Toolboximgaussfilt可换成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 调用 Demosigma0.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.5stop_tol保持 1e-3迭代次数上限 300。5. 常见问题与排查五个反复出现的坑及解决办法5.1 图像出现块状伪影/油彩感现象平滑区域不是细腻过渡而是变成一块一块的色块或棋盘格。原因lambda 取值过大显式迭代进入不稳定区或者停止阈值放太宽迭代在震荡区间提前退出留下振荡残影。经典求稳让许多新手把 lambda 调到 0.25 以上结果就是块状伪影。解决先把 lambda 降到 0.150.2 再跑。如果伪影消失说明就是步长问题。若 lambda 已经小于 0.2 还有色块查一下扩散系数把func从classic1换成tukey逆扩散区间消失后块状伪影通常也会跟着消失。5.2 像素值发散迭代过程爆炸现象迭代几步后图像出现单个亮点或整片高亮像素值超出 01 范围甚至变成 NaN。原因最直接的是 lambda 大于 0.25显式格式不稳定。另一类是输入图像没有经过im2double还是 uint8迭代时像素值溢出。扩散系数选 classic1/classic2 时逆扩散放大噪声尖峰也会导致局部爆炸。解决检查输入是 double 且范围在 01lambda 强制限制到 0.2 左右扩散系数优先选 tukey它在梯度超过 K 时直接置 0没有逆扩散区间。改完这三处爆炸问题基本不会出现。5.3 边缘被磨平去噪图看起来像高斯模糊现象噪声没了但边缘也钝了人物轮廓发虚。原因K 值太大把大部分真实边缘梯度当成了噪声背景扩散系数在边缘处依然很大或者 Catté 预平滑的 sigma 太大把边缘梯度在估计阶段就削平了。解决换成K_modeauto让 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 衡量结构和亮度相似度但它们对“边缘是否尖锐”并不敏感。这里补一个边缘保持指数 EPIgx_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自动 KTukey 系数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 一起用。希望这套验证流程能帮你在自己的数据上少走弯路。本文还有配套的精品资源点击获取