★ FEATURED ARTICLE 简介面向医学影像处理方向的初学者和研究人员围绕NSST与PCNN设计了一整套多模态医学图像融合MATLAB仿真流程重点解决CT、MRI等异源图像特征提取与融合决策问题。压缩包共40个文件、约2.31MB主体为27个.m源码文件含NSST工具箱、主程序及参数配置另有MEX动态库mexw32/mexw64/mexmac加速关键运算jpg、png、gif示例图与avi仿真录像用于对照运行结果、复现实验。已有960人学习下载。通过源码可理解非亚采样剪切变换的多尺度方向分解、PCNN脉冲耦合系数选择及逆变换重建的完整链路录像直观展示从图像输入到融合结果生成的全过程示例图像便于验证不同模态组合下的融合效果。整体体量精简适合在MATLAB 2021a中边看边练也是扩展医学图像融合实验的实用参考资料。1. 基于NSST和PCNN的多模态医学图像融合算法这套仿真在解决什么以及你大概要付出多少成本多模态医学图像融合这几年从论文选题走向了实际工程验证场景最常见的诉求就是把 CT 的骨骼细节和 MRI 的软组织信息拼到一张图里或者把 MRI 的解剖结构当成底图叠加 PET/SPECT 的功能代谢区域。标题里这条路线——NSST 非下采样剪切波变换负责多方向分解PCNN 脉冲耦合神经网络负责高频子带系数选取——是目前这类仿真里相当典型的组合。用 MATLAB 实现的好处是能快速验证想法并且配合仿真录像把迭代过程直观展示出来适合正在做课程设计、毕业课题以及刚接触多模态融合的从业者。这套方案解决的核心问题可以概括成一句话如何在保持图像边缘和纹理的同时从两个甚至多个源图像中挑选出“可信度更高”的细节。单纯做加权平均会让边缘模糊单纯取最大又会引入噪声NSSTPCNN 的组合则在方向上做了更精细的切分再通过神经元点火次数判断该信任谁。代价也很直接NSST 工具箱要自己找、参数要调、PCNN 有多个变体代码跑通不难跑出既清晰又自然的融合图需要一点排查功底。如果你只想要一张能交差的融合效果图这套组合可能有点重但如果你需要的是“能解释为什么这样融合”的完整仿真链路它确实值得投入。下面从原理、代码、录像和踩坑四个层面把这条线讲透。2. NSST分解与PCNN神经元先搞清楚为什么要这么组合2.1 从二维小波到非下采样剪切波平移不变性和多方向性从哪里来二维离散小波变换只有水平、垂直和对角三个方向这对自然图像可能够用但医学图像里的血管、骨折线、软组织边界方向是连续变化的。早期方案用非下采样轮廓波 NSCT 解决方向性问题但 NSCT 的扇形滤波器组实现复杂在多方向分辨率的灵活性上也不够好。NSST 的思路是把 Shearlet 变换改成非下采样结构通过非下采样金字塔获得多尺度再用剪切滤波器组获得多方向每个方向的信息分开处理。非下采样意味着去掉了下采样过程系数的尺寸跟原图保持一致平移变化时系数的幅值不会剧烈抖动医学图像配准如果存在亚像素残留这种平移不变性就能减少融合产生的伪影。另一个直观收益是NSST 对边缘的方向描述非常精确几个像素宽的毛细血管在传统小波里会被稀释放进对角方向而 NSST 可以把它分配到最接近实际角度的子带里后续 PCNN 选择系数时就更能保留真实结构。在具体参数上我一般会把分解层数 nlevels 设成 3 层方向数按 [2 3 3] 布置也就是第一层两个方向、第二层和第三层各三个方向。方向数越多图像细节切得越碎但计算量和存储也会线性增长。CT/MRI 这类灰度图[2 3 3] 已经能看出明显优势如果是血管造影之类的方向性极强图像可以改成 [3 3 4]。对 512×512 的输入三到四层分解在普通台式机上都能做到秒级不会成为仿真瓶颈。2.2 PCNN作为融合规则点火次数是怎么把“哪个源图系数更可靠”量化出来的PCNN 是从猫视觉皮层神经元同步放电特性抽象出来的网络模型每个像素位置对应一个神经元神经元接收外部刺激这里就是 NSST 子带系数同时通过链接权重接收周围神经元的点火反馈。当一个神经元的内部活动量超过动态阈值时它就点火点火后阈值被抬高过一阵阈值又慢慢下降于是出现周期性脉冲发放。用在图像融合上我们通常不关心脉冲的精确时间而是统计固定迭代步数内每个位置的点火总次数。某个位置的系数如果对比度突出、周围又跟着点火它的累计次数会明显更高。把这个逻辑套到高频子带A 图在某个位置点火次数高意味着该位置的局部结构更显著我们就优先选 A 图的系数。这个行为本质上是一种带邻域上下文信息的局部对比度投票比单纯“绝对值取大”更稳健因为单个像素的噪声不会轻易改变邻域整体的投票结果。PCNN 的简化模型通常关注五个参数链接强度 β、阈值衰减系数 αθ、阈值补偿幅度 Vθ、链接核大小以及最大迭代次数。β 控制邻域点火对当前神经元的影响程度太小就退化成单像素比较太大则会让亮斑向外扩散。Vθ 要跟归一化后的刺激量纲匹配否则全图点火或者全图不点火。这两个参数是仿真里最需要反复试的“玄学”区域后面我会给一组能直接起步的数值。2.3 低频和高频的差异化处理为什么不能一个规则打天下很多初学者把 NSST 分解到的所有系数都扔进 PCNN这是最容易翻车的做法。低频子带近似于原图的能量分布和亮度底色包含的是器官轮廓、大块组织层次这类信息用 PCNN 的点火投票来选很容易在亮度跳变区域产生不连续感融合图会显得很“脏”。低频更适合用区域能量加权因为低频本身就代表能量谁的局部能量大谁在那个区域的信息占主导。高频子带承载的是边缘、纹理、斑点这些细节直接比较系数绝对值很直观但会牺牲邻域一致性。PCNN 在这里有天然优势点火次数本身就考虑到了邻域耦合可以把方向子带里“既亮又连续”的结构挑出来而不是只挑最亮的孤点。所以常见做法是低频用固定窗口的区域能量计算权重高频用 PCNN 点火次数决定系数来源最后重建成像。这套分工既计算简单又能在视觉上保持整体亮度自然、边缘锐利。如果你的源图像有一张特别亮而另一张特别暗低频权重里还要再加一个全局亮度匹配项否则融合结果会一边倒地偏暗。3. 用MATLAB把NSSTPCNN融合跑通核心代码与参数设计3.1 融合框架与主流程整个融合流程可以分成六步读图与预处理、NSST 分解、低频融合、高频 PCNN 融合、逆 NSST 重建、指标评价。我习惯把它拆成几个独立函数而不是全部堆在一个 main 脚本里这样后面做参数扫描和录像时都能复用。% main_fusion.m clear; close all; clc; addpath(./nsst_toolbox); % 换成你本地NSST工具箱路径 addpath(./utils); % 自定义函数目录 % 1. 读图并转double imgA im2double(imread(CT_slice.tif)); imgB im2double(imread(MRI_slice.tif)); % 2. 统一尺寸NSST滤波器对尺寸有整除要求 if ~isequal(size(imgA), size(imgB)) imgB imresize(imgB, size(imgA)); end % 3. NSST分解 nlevels [2 3 3]; coeffsA nsst_dec1(imgA, nlevels); coeffsB nsst_dec1(imgB, nlevels); % 4. 低频融合 coeffsF{1} fuse_low(coeffsA{1}, coeffsB{1}); % 5. 高频PCNN融合逐层逐方向处理 params.beta 0.2; params.alpha_theta 0.9; params.v_theta 0.5; params.wsize 5; params.maxIter 30; for layer 2:numel(coeffsA) for dirIdx 1:numel(coeffsA{layer}) HA coeffsA{layer}{dirIdx}; HB coeffsB{layer}{dirIdx}; coeffsF{layer}{dirIdx} fuse_pcnn(HA, HB, params); end end % 6. 逆变换重建 imgF nsst_rec1(coeffsF); imshow(imgF);这段代码的逻辑很直白先把两幅图统一到同一尺寸然后分解低频调用一个函数高频进入 PCNN 融合函数。注意如果 nsst_dec1 返回两个输出比如第二个是剪切滤波器组逆变换时要把这个滤波器组也传给 nsst_rec1否则不同工具箱会报参数个数错误。不同来源的 NSST 工具箱封装略有差异但调用思路一致。3.2 NSST分解与逆变换的调用方式NSST 工具箱里的分解函数命名在不同版本里不完全相同最常见的是 nsst_dec1 和 nsst_rec1。分解结果一般用 cell 数组存放第一个 cell 是低频系数从第二个 cell 开始每层是一个 cell里面放该层所有方向的子带。方向子带本身还是 cell因为每层的方向数并不相同。% 查看分解结果结构 coeffsA nsst_dec1(imgA, [2 3 3]); disp(size(coeffsA{1})); % 低频与原图同尺寸 disp(numel(coeffsA{2})); % 第2层cell里的方向子带个数关键在于记住低频系数是全尺寸的高频子带也跟原图同尺寸这是 NSST 非下采样的特点。如果你以前用过小波工具箱很容易下意识认为高频子带会比原图小这会直接导致后面矩阵运算维度不匹配。逆变换时只要把融合后的 coeffsF 按相同结构传回去就行不需要手动指定方向数工具箱会从 cell 结构里推断。有一点需要特别提醒如果你的 MATLAB 版本较老某些 NSST 工具箱的 mex 文件需要重新编译。报错信息如果是“找不到已编译的函数”先确认工具箱里有对应的 .mex64 文件没有就把源文件里的 .c 文件用 mex 命令编译一遍。这个问题跟融合算法本身无关但卡住的人最多。3.3 低频融合区域能量与加权平均低频系数的融合规则我习惯用局部区域能量加权分别计算两个低频系数在 5×5 或 7×7 窗口内的平均能量能量高的那个在这个区域占更大权重。这样既能保留整体亮度又能让病灶区域或高密度区域自然过渡。function LF fuse_low(LA, LB, wsize) if nargin 3, wsize 5; end k ones(wsize) / (wsize^2); % 局部能量估计 EA conv2(LA.^2, k, same); EB conv2(LB.^2, k, same); % 权重归一化 wA EA ./ (EA EB eps); wB 1 - wA; LF wA .* LA wB .* LB; endconv2 的 same 参数能保证输出尺寸和输入一致边界处会隐式补零对 512×512 图像来说边界一小圈的影响可以忽略。eps 是为了防止两张图都是平坦黑区时除零。wsize 太小会让权重视觉上出现颗粒感太大又会让亮度过渡平滑、细节变得模糊CT/MRI 融合用 5 是比较折中的选择。关于低频还有一条经验如果两张源图整体亮度差异很大比如 MRI 比 CT 整体暗很多直接在权重里做归一化会导致暗图几乎不被采用。这种情况下可以先把两张低频系数分别归一化到相同均值和方差再做能量加权最后在重建前把融合结果调回合理灰度范围。这一步虽然简单但对最终视觉效果影响巨大。3.4 高频融合PCNN迭代统计点火次数高频子带可能包含正负系数所以送入 PCNN 之前我习惯先取绝对值然后统一归一化到 [0,1] 区间让不同层、不同方向子带的刺激量处于同一量纲。点火次数的计算用卷积实现链接输入避免写三重 for 循环速度能快一个数量级。function T pcnn_firing(S, params) [m, n] size(S); S double(S); beta params.beta; alpha_theta params.alpha_theta; v_theta params.v_theta; wsize params.wsize; maxIter params.maxIter; % 链接权重用归一化均值核 W ones(wsize) / (wsize^2); Y zeros(m, n); % 点火状态 L zeros(m, n); % 链接输入 U zeros(m, n); % 内部活动 Theta ones(m, n) * v_theta; % 初始阈值 T zeros(m, n); % 累计点火次数 for t 1:maxIter L conv2(Y, W, same); U S .* (1 beta * L); Y double(U Theta); T T Y; Theta alpha_theta * Theta v_theta * Y; end end这段代码里的阈值更新用的是离散恢复模型未点火时阈值按 αθ 比例缓慢衰减点火后通过 Vθ 把阈值重新抬升。αθ 取 0.9 时阈值衰减较慢神经元不会疯狂重复点火取 0.5 时阈值下降快点火次数普遍偏高适合细节特别稀疏的图像。β 的作用体现在链接输入上β0 时 PCNN 退化成纯阈值比较邻域信息完全没用所以 β 至少要保留 0.1 以上。融合时比较两幅图对应子带的点火次数图function HF fuse_pcnn(HA, HB, params) maxAbs max(abs(HA(:)), abs(HB(:))) eps; SA abs(HA) ./ maxAbs; SB abs(HB) ./ maxAbs; TA pcnn_firing(SA, params); TB pcnn_firing(SB, params); % 点火次数高的一方系数直接采用 HF HA .* (TA TB) HB .* (TA TB); end注意这里融合用的是原始系数 HA 和 HB而不是归一化后的 SA、SB。归一化只影响 PCNN 的点火决策决定权重归属后还是要把真实的子带系数填回重建通道。如果直接用绝对值参与重建图像会丢失正负交替的边缘结构视觉上出现明显的过冲。PCNN 参数建议按下表起步再根据你的图像特点微调参数含义常用范围初始值beta链接强度0.1~0.40.2alpha_theta阈值衰减系数0.6~0.990.9v_theta阈值基值/点火恢复量0.1~10.5wsize链接核尺寸3 或 5 或 75maxIter最大迭代次数20~8030调参有一个快速诊断手段跑完 pcnn_firing 后直接画 T 的直方图。如果 T 几乎全为 0说明阈值太高或迭代太少降低 v_theta 或增加 maxIter如果 T 全为 maxIter说明每次迭代都在点火需要增大 v_theta 或降低 β。一个健康的点火直方图应该是从 0 到 maxIter 都有分布集中在中间段。3.5 重建与客观评价指标融合完成后 nsst_rec1 会自动把低频和高频组合回灰度图。这时候不要忙着截图先用客观指标确认结果没有明显退化。医学图像融合没有标准参考图一般不直接算 PSNR而是看这几项指标衡量内容是否需要参考图空间频率 SF图像清晰度与边缘强度不需要信息熵 Entropy融合图保留的信息量不需要互信息 MI融合图与各源图的共享信息需要源图边缘保持度 QABF边缘方向与强度的保持程度需要源图SF 的计算可以简化为行、列方向差分后的均方根function sf spatial_frequency(img) RF sqrt(mean(mean(diff(img, 1, 1).^2))); CF sqrt(mean(mean(diff(img, 1, 2).^2))); sf sqrt(RF^2 CF^2); end如果仿真目标是证明算法在单一指标上的优势我建议把 MI 和 QABF 作为核心指标SF 作为辅助。SF 对噪声很敏感PCNN 融合图如果边缘提升明显SF 自然会涨但如果噪声也被保留了SF 可能虚高。这时候结合 QABF 看边缘方向是否被正确保留比单看 SF 更有说服力。4. 让仿真过程可录制从脚本到仿真录像的完整落地4.1 仿真录像要录什么不建议只录最终结果很多仿真录像只把最后的融合图摆在那里配上几行说明就完了这在实际交流中说服力很弱。真正有价值的是把“PCNN 迭代过程中点火次数如何变化”“融合图如何从模糊到清晰”这个过程录下来。评审或同事看到动态过程才能直观理解 NSST 分解出来的方向子带和 PCNN 的投票机制到底起了什么作用。我习惯在画面里放三个窗口左侧是 CT 源图或 MRI 源图中间是正在迭代的低频融合结果右侧实时显示某个代表性高方向子带的点火次数图。每迭代几次更新一次右侧画面点火次数图会从稀疏散点逐渐聚合成连续结构这个过程非常能说明算法在做“邻域投票”而不是单像素蛮力选择。类似地也可以把最终融合结果做成从低频重建逐步叠加高频的动画那更像是“剥洋葱”的反向过程。录制时长控制在 30 到 60 秒比较合适。帧率太低会很卡帧率太高文件会非常大。对于 MATLAB 的 VideoWriter5 到 10 帧每秒就够因为图像算法过程本身是平滑变化的不需要游戏级流畅度。4.2 用VideoWriter录制动态过程VideoWriter 是 MATLAB 里最直接的录像工具。核心思路是在循环里刷新图像显示然后用 getframe 抓取当前 figure 内容写入视频对象。注意必须先把图像句柄提前创建好循环内只更新 CData而不是反复调用 imshow否则窗口会闪烁录像也会出现花屏。% record_fusion.m fig figure(Position, [50 50 1200 400], Color, w); ax1 subplot(1, 3, 1); imshow(imgA); title(Source A: CT); ax2 subplot(1, 3, 2); imshow(imgB); title(Source B: MRI); ax3 subplot(1, 3, 3); hF imshow(zeros(size(imgA))); title(Fused Result); writer VideoWriter(fusion_process.avi); writer.FrameRate 5; % 5帧/秒 writer.Quality 90; % 视频质量 0~100 open(writer); for t 1:params.maxIter % 更新融合图这里用简化示意实际每次迭代需要重新计算融合结果 fusedNow update_fusion_once(coeffsA, coeffsB, params, t); set(hF, CData, fusedNow); title(ax3, sprintf(Fused Result (iter%d), t)); drawnow; % 强制刷新图形 frame getframe(fig); writeVideo(writer, frame); end close(writer);getframe 之前必须调用 drawnow否则 MATLAB 可能直接把上一次的图形缓存抓进去导致视频里大量重复帧。如果把 figure 创建时的 Position 设到屏幕外或者窗口被其他程序遮挡录像也可能黑屏。保险做法是把 Position 设置成屏幕可见区域并且把分辨率控制在 800×600 到 1024×768 之间。分辨率越高视频文件越大但融合细节会清楚一些如果只是演示800 宽完全足够。还有一个小习惯每迭代一步给标题加上当前迭代次数这样观看者能明确知道算法运行到哪个阶段。视频编码方面Motion JPEG 的 AVI 文件通用性最好MPEG-4 文件小但要确认系统安装了对应编码器。如果 writeVideo 报编码错误把 VideoWriter 构造改成 VideoWriter(fusion_process.avi)也就是默认 Motion JPEG AVI基本能解决兼容问题。5. 避坑与排查配准、取整、黑屏和三色归一化5.1 尺寸不匹配导致 NSST 分解直接报错注意“整除”和偶数尺寸现象nsst_dec1 执行到中途报错“数组维度不匹配”或者分解出来的低频系数尺寸跟原图不一致。原因很简单NSST 的非下采样金字塔需要输入图像尺寸能被特定因子整除某些工具箱还要求长宽是偶数。两张源图虽然后台能显示在同一坐标系里但像素尺寸往往不一样医学数据里 DICOM 导出后尤其常见。解决在分解前统一尺寸。我一般用 imresize 把两幅图都转成相同的正方形比如 512×512 或 256×256。如果不想缩放带来的细节损失可以先用 padarray 补零到公共尺寸记录补零区域重建后裁剪回来。注意补零前要把数据类型转成 doubleuint8 在求差和卷积时经常溢出这个细节能同时踩两个坑。5.2 融合结果灰蒙蒙、PCNN 不点火先查数据类型和归一化现象融合图边缘全没了像盖了一层磨砂玻璃或者 PCNN 点火次数 T 全为 0、全为 maxIter融合结果呈现出明显的人工“拼缝”。原因大多是读图后直接用了 uint8或者高频系数没有做归一化。uint8 在 0~255 范围内做 diff 会溢出卷积后局部能量偏高PCNN 阈值 Vθ 很难匹配到合适区间。另外高频系数有正有负负值绝对值大的位置恰恰是强边缘如果只取原始负值做比较PCNN 会被负值直接压住。解决读图立即 im2double高频输入 PCNN 前先取绝对值再除以整组最大值融合时用原始系数回填。诊断办法是查看 T 的直方图如果直方图呈现两端堆积的极端分布大概率是归一化没做到位。5.3 彩色功能图融合后颜色发黑多通道直接分解的翻车现象PET/SPECT 这类彩色功能图和 MRI 灰度图融合很多人直接对彩色图的三通道分别做 NSSTPCNN最后合成 RGB结果颜色发黑、色相漂移病灶区域看起来是一团暗斑。原因在于彩色图的 R、G、B 三个通道亮度相关性很强独立分解后三个通道的融合权重不一致重建出来的颜色自然混乱。解决先把彩色功能图转到 HSV 或 IHS 空间提取 V 通道亮度与解剖灰度图做融合融合结果再跟原始的 H 和 S 通道合成彩色图。V 通道融合可以用前面整套 NSSTPCNN 流程H 和 S 保持功能图原样。这样既能保留功能图的色彩定位信息又能用解剖图补充结构细节。如果 H/S 通道有明显噪声可以先做一次轻微高斯滤波但不要用 NSST 去分解它们。5.4 录像黑屏、文件巨大getframe 与图形刷新的配合问题现象录制完成的 AVI 打开后全黑或者视频文件几十 MB 只有几秒画面。原因通常有两个getframe 抓取时图形窗口没有刷新MATLAB 返回了后台缓存另一个是分辨率设得太大比如直接用默认 figure 尺寸往往接近 1280×960帧率还设置到 30几秒视频就上百 MB。解决循环里每更新一次 CData 就立刻执行 drawnowfigure 创建时要显式设置可见位置不要用 visible off。编码格式首选默认 Motion JPEG AVI如果文件仍然过大把分辨率从 1200 宽降到 800 宽帧率从 10 降到 5。另外如果录像内容只有灰度图最好在 figure 显示时用 colormap(gray) 统一否则视频里颜色会偏绿或偏橙。5.5 PCNN 参数调不动的“假死”状态现象β 改了 10 倍、迭代次数从 20 改成 80融合结果几乎没变。原因可能是链接核 W 没有归一化导致邻域影响被某个超大权重点主导或者阈值 Vθ 相对 S 的量纲相差太大神经元根本没有进入正常点火节奏。解决把 W 设为 ones(wsize)/wsize^2确保 L 的量级和 S 一致检查 S 是否已经归一化到 [0,1]用直方图确认 T 的分布。还有一个隐蔽坑PCNN 融合函数里如果对每个方向子带都重新做归一化方向子带的相对强度关系就被抹平了应该先算出整组高频系数的最大绝对值再统一归一化而不是每个子带各自归一。6. 用一盏“路灯”验证融合效果单次实验不足以证明算法有效如果你只跑通一张图就开始接下一个项目后面一定会被评审或导师问住为什么这个参数在这个图上有效换一组数据还行不行我自己的习惯是每次仿真至少准备三组不同类型的源图像对一组 CT-MRI、一组 MRI-PET 灰度模拟、一组彩色 PET-MRI全部记录 SF、熵、MI 三个指标存成一张表用均值±标准差说话。单张融合图“眼前一亮”经常是靠运气因为自然图像里的偶然结构恰好碰对了参数。验证时还要做参数敏感性扫描。固定其他参数只让 β 在 0.1、0.2、0.3 三个值里变化迭代次数在 20、30、50 里变化输出同一组图的结果for beta [0.1, 0.2, 0.3] for iter [20, 30, 50] params.beta beta; params.maxIter iter; fused run_fusion_pipeline(imgA, imgB, params); fprintf(beta%.1f iter%d SF%.4f MI%.4f\n, ... beta, iter, spatial_frequency(fused), mutual_info(imgA, imgB, fused)); end end如果某个指标在参数边界处剧烈跳变说明融合规则对参数过于敏感这时要回头检查是不是 PCNN 点火次数已经饱和而不是继续调参数。指标变好不等于视觉变好融合图里病灶边缘是否连续、灰度是否自然最终必须在显示器上用肉眼确认这一点比任何数值都重要。我现在做仿真时还会顺手把每版参数写进文件名比如 fused_b020_iter30.qabf.png这样回溯起来非常省力。图像融合这行没有银弹把每一次“看起来不错”都变成可复现的记录比追求所谓最优参数有用得多。希望以上的参数起步值和排查路径能帮到你让你少花几天在工具箱报错和灰蒙蒙的结果上。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
