简介这是面向MATLAB用户的NSCT非下采样剪切波变换实现工具包目标使用者是图像处理研究者、工程师及刚接触多尺度分析的学习者作为一种多分辨率分析工具它比传统小波更擅长捕捉图像边缘等几何结构可用于理解非下采样金字塔与方向滤波器组原理并直接支撑图像去噪、融合、边缘检测和特征提取等任务。资源共48个文件压缩包约293KB主体为28个m函数和脚本同时包含6张tif测试图MRI、PET、SPECT等、mexw64/mexmac/dll预编译接口、C源码与说明文档兼顾跨平台可运行性和源码级可读性。已有1386人浏览学习。核心代码覆盖NSCT正逆变换完整流程如nsctdec/nsctrec、nsdfbdec/nsdfbrec、atrousfilters、zconv2等模块并配有source images测试图和decdemo演示脚本可完成从图像加载、多尺度方向分解、系数操作到重构对比的完整实验对照特别适合对照论文和教材进行算法复现与验证C源码帮助理解并修改预编译模块mex接口省去自行编译的繁琐步骤。整体结构紧凑、入口清晰是一份能直接上手并深入剖析NSCT算法细节的MATLAB工具包。1. 找NSCT的MATLAB实现前先认清这个变换的代价与回报搜 NSCT 程序实现和 matlab 工具包的人绝大多数最后都卡在同一个地方分解和重构都能跑但重构回来的图和原图差了一圈亮边或者误差大得离谱。NSCT非下采样轮廓波变换不是 MATLAB 内置函数它通常以社区流传的 MATLAB 工具包形式出现这套工具包把小波变换那种横平竖直的分解升级成多个楔形方向分解而且全程不下采样。因此它在图像融合、去噪、边缘增强里表现好代价是系数冗余高、速度和内存开销大。这篇笔记适合要做融合或去噪实验却被 NSCT 参数和系数结构卡住的人我把工具包的调用方式、参数选择和我在实践中踩过的坑一次讲清。2. NSCT的MATLAB实现骨架非下采样金字塔配合方向滤波器组到底做了什么事2.1 金字塔不做下采样细节子带才会和原图一样大NSCT 的第一个模块叫非下采样金字塔NSP。它的任务是把图像按尺度分开低频分量一层层变平滑每一层减去上一层的低频就得到带通细节。关键点是整个过程不做抽取也就是不丢弃任何像素。滤波器靠插入零值来“膨胀”支撑范围每一层用到的等效滤波器越来越长覆盖的频率区域越来越窄但输出图像始终和原图尺寸一致。这个设计解决了离散小波变换的一个老问题小波每分解一层子带尺寸减半平移图像之后再分解系数分布会剧烈变化因为下采样点对不准。NSCT 不做下采样图像平移几个像素系数只是跟着平移不会产生混叠噪声。这就是所谓的平移不变性。代价也很直观原来小波 1024×1024 分解三层后最小子带只有 128×128NSCT 里每一层每种子带都是 1024×1024数据量直接膨胀几十倍。2.2 方向滤波器组用扇形滤波器把频域切成楔形而不是横平竖直NSCT 的第二个模块叫非下采样方向滤波器组NSDFB。它把每一层的高频细节按角度进一步切分二维频域被分成一个个楔形区域每个楔形对应一个方向子带。方向数按 2 的幂次增长比如一个尺度 8 个方向下一层细尺度可能 16 个方向。和小波只有水平、垂直、对角三个方向相比NSCT 能沿任意角度贴合纹理和曲线边缘。工具包内部用的是扇形滤波器组加梅花采样来切换方向坐标。理解到这一步就够了你不需要自己去设计那些扇形滤波器和梅花矩阵工具包已经把方向划分封装成黑匣子。需要做的是选择方向数然后记住一件事——方向子带和细节子带一样尺寸都等于原图尺寸。这个特性在做图像融合时特别方便系数之间不需要插值对齐直接逐像素比较就行。2.3 别走弯路在Contourlet上去掉下采样结果不是NSCT很多第一次接触 NSCT 的人会想Contourlet 和 NSCT 不就是差一个下采样吗把 Contourlet 的代码里抽取那几步删掉不就行了我试过结果重建出来全是条纹。原因在于 Contourlet 的下采样和它的滤波器组设计是一体的去掉抽取之后频域划分不再满足完全重构条件。NSCT 不是简单“不下采样的Contourlet”它的金字塔和方向滤波器都重新设计过才能在无下采样的情况下保持精确重构。这也是为什么我建议直接用社区流传的 MATLAB 工具包而不是自己从论文图纸开始写。NSCT 的滤波器组设计本身是论文级工作包含双通道滤波器、扇形滤波器、梅花采样矩阵的配合自己实现一遍的调试周期很长。工具包也不是十全十美但它的分解和重构代码已经调好拿来做实验比从零写可靠得多。2.4 NSP的atrous思路示意一段看懂“膨胀卷积”的代码理解 NSCT 的尺度分解看下面这段结构示意代码就够了。它展示的是 atrous 算法的核心即滤波器插零膨胀后与图像卷积function [lp, hp] nsp_step(lp_in, g0, g1, rate) % 单层非下采样金字塔atrous 思路示意 % lp_in : 当前低频图像 % g0/g1 : 低通/高通滤波器系数行向量 % rate : 膨胀间隔第 j 层传 2^(j-1) g0_up dyadup(g0, rate); % 在系数间插零撑大滤波器支撑 g1_up dyadup(g1, rate); % 先做行方向卷积再做列方向卷积保持图像尺寸 lp conv2(wextend(2D, sym, lp_in, length(g0_up)), g0_up, same); lp conv2(lp, g0_up, same); hp lp_in - lp; % 带通细节 输入低频 - 新低频 hp hp(1:size(lp_in, 1), 1:size(lp_in, 2)); end逻辑说明第一层 rate 为 1滤波器不膨胀第二层 rate 为 2滤波器每隔一个系数塞一个零等效支撑变长第三层 rate 为 4以此类推。这样每一层看到的频率范围更窄但输出尺寸不变。参数说明g0 和 g1 是一维滤波器系数这里用可分离卷积实现二维滤波wextend 的 sym 模式做对称延拓避免卷积后边缘收缩dyadup 来自小波工具箱。这段代码只是理解思路用的示意工具包里的真实实现还包含双通道重构条件不建议直接拿去替换 nsctdec。3. 把NSCT工具包跑起来nsctdec/nsctrec的调用顺序与参数选择3.1 安装验证把工具包放进路径后先确认nsctdec和nsctrec在不在下载得到的 NSCT 工具包一般是一个目录里面有一堆 .m 文件和子目录。第一步是把它加进 MATLAB 搜索路径。我习惯用 genpath 把整个目录树都加进去避免漏掉子目录里的滤波器初始化脚本dest fullfile(你的本地路径, NSCTToolbox); % 改成实际解压目录 addpath(genpath(dest)); if exist(nsctdec, file) 2 disp(nsctdec 已在路径上可以继续); else disp(nsctdec 未找到检查目录路径是否正确); end逻辑说明exist 返回 2 表示找到同名 .m 文件返回 0 表示不在路径上。这一步能过滤掉大量“明明下载了却报错”的情况。参数说明dest 路径不要带中文字符老工具包对非 ASCII 路径支持不好容易出现奇怪的加载失败。验证通过后我建议先运行工具包自带的 demo 脚本一般是以 demo 开头的 .m 文件跑通了再开始改自己的参数。3.2 最小分解与重构一组不会出错的初始参数初次上手不要一上来就追求 16 方向五层分解。先用一组保守参数跑通流程确认工具包在你机器上工作正常im im2double(imread(cameraman.tif)); % 任意灰度图都行 levels 4; % 金字塔层数 directions [8 4 2 1]; % 从细到粗各层方向数最后一个必须是 1 pfilter 9-7; % 金字塔滤波器 dfilter dmaxflat7; % 方向滤波器 dec nsctdec(im, levels, directions, pfilter, dfilter); rec nsctrec(dec, pfilter, dfilter); err max(abs(im(:) - rec(:))); fprintf(重构最大误差: %e\n, err);逻辑说明nsctdec 把图像分解成系数元胞nsctrec 从系数重构图像。重构误差在 1e-10 量级说明工具包安装正确分解和重构互为逆过程。参数说明im2double 把图像归一化到 [0,1]避免 uint8 的整型溢出问题levels4 是常见的折中既能看到多尺度效果计算量又不至于失控directions 数组长度必须等于 levels且从细尺度到粗尺度排列最后一个固定为 1代表最粗那层只有一个低频近似。3.3 四个必调参数levels、directions、pfilter、dfilter各自怎么改参数常用取值作用我一般怎么选levels2 到 5金字塔分解层数实验对比用 4快速验证用 2directions[4 2 1]、[8 4 2 1] 等每层方向子带数2 的幂次末位为 1纹理复杂用 8 起步平滑图像用 4pfilter9-7、maxflat金字塔滤波器默认 9-7速度和精度均衡dfilterdmaxflat7、vk方向滤波器论文复现优先 dmaxflat7逻辑说明directions 的每个值对应一个金字塔尺度的方向分解数不是随便填的。比如 directions[8 4 2 1]最细尺度 8 个方向第二尺度 4 个方向第三尺度 2 个方向最粗尺度 1 个。方向数从细到粗递减这是 NSCT 的标准设计。参数说明9-7 指的是类似 JPEG2000 的 9/7 双正交滤波器重建质量和速度都比较理想maxflat 支撑更长边缘保留更好但更慢dmaxflat7 方向选择性好纹理丰富时比 vk 稳。提示如果下载的工具包版本入口参数顺序不一样最快的确认方式是打开自带 demo 脚本看调用行我手头这份常用的是 nsctdec(x, levels, directions, pfilter, dfilter) 这个顺序。3.4 重构误差的量化检查为什么高频系数不能随便置零很多人拿到 NSCT 后第一件事就是把高频系数置零想看看低通效果结果发现重构图像出现条纹或亮度偏移。先做个基线实验认识高频系数的贡献decZero dec; % 复制上一小节的分解结果 for l 1:numel(dec) - 1 % 遍历所有方向尺度层保留最后的低频 for d 1:numel(dec{l}) decZero{l}{d} zeros(size(dec{l}{d})); % 高频全部清零 end end recZero nsctrec(decZero, pfilter, dfilter); errZero max(abs(im(:) - recZero(:))); fprintf(高频清零后重构误差: %e\n, errZero);逻辑说明高频清零后重结构相当于只保留低频近似原本的高频细节丢失误差自然巨大。这不是工具包有问题而是 NSCT 的系数不具备“随便改一个系数不影响其他区域”的性质。参数说明numel(dec)-1 是因为最后一个元胞元素是低频近似不能清零。这个实验还有个作用就是验证你对系数结构的理解是否正确——如果清零后重构出来的是一张模糊但正常的低通图说明结构索引用对了。4. 系数到底存在哪里NSCT输出的元胞结构、方向索引与平移不变性验证4.1 输出结构长什么样dec{1}到dec{end}的层级关系nsctdec 返回的是一个嵌套元胞数组这是新手最先懵的地方。结构并不复杂外层元胞按尺度分层dec{1} 是最细尺度dec{2} 是第二细尺度以此类推dec{end} 是最后一个低频近似。每个 dec{l} 又是一个元胞第 l 层的方向子带都装在这里面。内层元胞的元素个数就是 directions(l)。numel(dec) % 外层元素个数 levels 1加了最后的低频 numel(dec{1}) % 最细尺度方向子带数量应该等于 directions(1) size(dec{1}{1}) % 方向子带尺寸和原图完全一致 size(dec{end}) % 低频近似尺寸也和原图一致逻辑说明NSCT 由于全程不下采样从最细方向子带到最后的低频近似所有图像尺寸都等于原图尺寸。这是 NSCT 区别于小波最直观的一点。参数说明numel(dec) 等于 levels1是因为最后一层低频被单独放在 dec{end}而 dec{1} 到 dec{end-1} 才是真正的方向细节层。4.2 方向子带索引规则与尺寸恒定性取某个方向子带的写法是 dec{层号}{方向号}。比如 dec{1}{1} 是最细尺度第一个方向dec{2}{3} 是第二尺度第三个方向。和 directions 数组对应dec{l} 里的元素个数等于 directions(l)。低频近似单独放在 dec{end}不再套内层元胞。subFine dec{1}{1}; % 最细尺度第 1 方向 subMid dec{2}{2}; % 第二尺度第 2 方向 base dec{end}; % 低频近似 figure; subplot(1, 3, 1); imshow(mat2gray(abs(subFine))); title(细尺度方向); subplot(1, 3, 2); imshow(mat2gray(abs(subMid))); title(中尺度方向); subplot(1, 3, 3); imshow(mat2gray(base)); title(低频近似);逻辑说明方向子带系数有正有负直接 imshow 会把负值截断所以用 mat2gray(abs(...)) 把绝对值映射到 [0,1] 再显示。参数说明方向子带里绝对值大的位置代表该方向上的强边缘或纹理低频近似 base 是平滑后的图像保留了整体亮度分布。如果显示出来的方向子带全是噪声状的小颗粒先检查 directions 是不是设成了 [1 2 4] 这种从粗到细的顺序方向数排列反了会拿到错误的子带。4.3 一个验证平移不变性的最小实验NSCT 的平移不变性是它相对 Contourlet 的核心卖点但验证的时候有个常见错误图像平移后再分解直接拿两个分解系数逐像素相减发现差异巨大然后误以为 NSCT 不具备平移不变性。实际上图像平移了系数也跟着平移必须先平移回来再比。shiftRows 9; shiftCols 11; decShift nsctdec(circshift(im, [shiftRows, shiftCols]), levels, directions, pfilter, dfilter); sub1 circshift(dec{1}{1}, [shiftRows, shiftCols]); % 系数平移回来 sub2 decShift{1}{1}; % 平移后图像的分解 residual max(abs(sub1(:) - sub2(:))); fprintf(平移不变性验证残差: %e\n, residual);逻辑说明先把原图分解的系数平移回去再和平移后图像的分解系数比较。理想情况下残差应该非常小实际因为边界延拓差异会有一个不算大的残差但远小于子带本身的幅值。参数说明circshift 是循环移位边界处两边的差异会比内部大想看内部区域的平移不变性可以把比较范围裁掉边界几十个像素再算残差。这个方法比看论文里的示意图靠谱能直接量化你手上这个工具包的平移不变性水平。5. 避坑指南重构残差、方向数约束、内存爆炸与边界伪影的5条记录5.1 重构残差大、边缘一圈亮线先查边界延拓与图像尺寸现象用 nsctdec 分解再 nsctrec 重构中心区域没问题但图像四周出现一圈亮线或暗线重构误差在 1e-2 量级。原因NSCT 的滤波器支撑会随层数膨胀默认延拓宽度不够边缘像素在卷积时没有足够的对称数据可用另外图像尺寸太小的时候最粗那层的滤波器支撑已经比图像还宽边缘伪影会污染整个子带。解决分解前手动对称延拓重构后裁剪回来我一般加 20 像素左右ext 20; imPadded padarray(im, [ext ext], symmetric); decP nsctdec(imPadded, levels, directions, pfilter, dfilter); recP nsctrec(decP, pfilter, dfilter); rec recP(ext1:end-ext, ext1:end-ext); % 裁掉延拓区域逻辑说明延拓让边缘像素在卷积时有完整的邻域数据可用重构后再裁剪边缘伪影被隔离在延拓区域里。参数说明ext 的取值要大于最大滤波器支撑的一半20 像素对 4 层分解通常够用如果分解层数多或滤波器支撑长加大到 32。5.2 directions参数报错或子带对不上方向数数组要和金字塔层数严格配套现象设置 directions[8 4 2] 或 [16 8 4 1] 时报错或者跑完发现某个层方向子带数量不对。原因directions 数组长度必须等于 levels且每个值必须是 2 的幂次最后一个必须为 1。有些版本还会校验相邻尺度的方向数关系细尺度的方向数不能超过上一层方向数的合理倍数。解决先跑最小配置验证再逐步加方向数。levels 2; directions [2 1]; % 最小验证配置 decTest nsctdec(im, levels, directions, pfilter, dfilter); fprintf(方向子带数: %d, %d\n, numel(decTest{1}), numel(decTest{2}));逻辑说明方向数翻倍是有结构约束的不是想设多少就设多少。最粗尺度低频层只有 1 个方向往细的方向数按 2 的幂次增加而且相邻层之间要保持合理的倍数关系。参数说明从 [2 1] 验证通过后再跳到 [4 2 1]、[8 4 2 1]一次加一层能快速定位是哪个参数组合触发的报错。实际实验里 [8 4 2 1] 是绝大多数图像都能稳定跑通的配置。5.3 重建后偏暗或偏亮能量对不上检查图像类型与滤波器组合现象分解再重构图像整体变暗或变亮PSNR 掉到 20 以下但局部结构还在。原因最常见的是输入图像没有转成 double 或归一化uint8 图像直接送进 nsctdec卷积时溢出截断其次是 pfilter 和 dfilter 的组合来自不同工具包版本重构增益不为 1。解决统一走 im2double 再分解如果确认输入没问题换回 9-7 加 dmaxflat7 这个经过大量验证的组合。img imread(some_image.png); im im2double(img); % 先归一化再分解 if size(im, 3) 3 im rgb2gray(im); % NSCT 工具包默认按灰度图处理彩色要逐通道跑 end逻辑说明NSCT 工具包内部假设输入是 double 类型uint8 的整数卷积会产生截断。彩色图像如果直接送进去工具包会把第三维当作批量维度处理结果各通道混在一起。参数说明im2double 对所有格式的输入都安全彩色转灰度用 rgb2gray或者保持彩色但逐通道分解再合并后者适合需要颜色信息的融合实验。5.4 大图一跑就内存不足冗余变换的代价与降参数建议现象512×512 的彩色图levels5directions[32 16 8 4 1]运行到一半提示内存不足或者等了十几分钟还没结束。原因NSCT 是冗余变换系数总量约为原图像素数的 1sum(directions) 倍这个配置下接近 62 倍。512×512 灰度图就有约 1.6 亿个 double 值加上卷积过程的临时数组内存很容易爆。解决先降方向数再降层数或者逐通道处理。levels 4; directions [8 4 2 1]; % 比 32 方向温和得多 decR nsctdec(im(:,:,1), levels, directions, pfilter, dfilter); decG nsctdec(im(:,:,2), levels, directions, pfilter, dfilter); decB nsctdec(im(:,:,3), levels, directions, pfilter, dfilter);逻辑说明把彩色三通道分别分解内存峰值比一次性处理三维数据低很多。方向数是内存增长的主要因素每多一层方向数翻倍数据量就翻一倍。参数说明一般实验 8 个方向已经能体现 NSCT 的优势16 方向只在纹理非常密集时有必要4 层 8 方向是我测试多台机器都不会内存崩溃的配置。速度方面长支撑方向滤波器是大头vk 比 dmaxflat7 略快纹理要求不高时可以换。5.5 系数置零后重建出现条纹方向子带之间存在频域重叠现象把某一个方向子带全部置零重构图像出现均匀分布的方向性条纹而且条纹方向和置零子带的方向相关。原因NSDFB 的方向子带在频域有过渡带重叠方向子带并不是完全正交的。强行挖空一个楔形频带相当于在频域做了一个尖锐切口重构时产生 Gibbs 振荡。解决不要硬置零改用系数收缩或衰减如果只想抑制某个方向分量把系数乘以一个小系数而不是清零。decAtten dec; decAtten{1}{3} decAtten{1}{3} * 0.1; % 衰减而不是清零 recAtten nsctrec(decAtten, pfilter, dfilter);逻辑说明乘以 0.1 保留了该方向的大部分结构只是降低了能量频域的过渡带变化平缓不会产生明显条纹。参数说明衰减系数在 0 到 1 之间0.05 到 0.2 是常用区间如果确实需要彻底移除某个方向分量建议配合局部区域平滑而不是全局硬置零。这条经验在做去噪实验时特别有用很多去噪算法直接对高频系数做硬阈值结果出现条纹换成软阈值就好了。6. 进阶把NSCT系数用进图像融合实验的最小改动6.1 区域能量取大的融合规则以及它为什么比简单平均稳NSCT 最常见的落地场景是图像融合。两幅图像先各自做 NSCT 分解低频部分用平均高频方向子带用能量取大。低频决定整体亮度和结构平均融合能避免亮度突变高频子带里绝对值大的系数代表更强的边缘或纹理取大规则能保留更清晰的细节。decA nsctdec(A, levels, directions, pfilter, dfilter); decB nsctdec(B, levels, directions, pfilter, dfilter); decF decA; % 先按 A 的结构复制 for l 1:numel(decA) - 1 % 跳过最后的低频 for d 1:numel(decA{l}) eA sum(decA{l}{d}(:).^2); % 全子带能量 eB sum(decB{l}{d}(:).^2); if eA eB decF{l}{d} decA{l}{d}; else decF{l}{d} decB{l}{d}; end end end decF{end} (decA{end} decB{end}) / 2; % 低频平均保留整体亮度 F nsctrec(decF, pfilter, dfilter);逻辑说明全子带能量取大是最简版本它把整幅子带交给能量更高的一方边缘保留效果已经优于空间域平均。参数说明想要更精细把 sum(...)^2 改成局部滑窗求和比如 5×5 邻域内系数平方和这样融合边界会更自然。低频也可以用局部窗口平均替代全局平均对亮度差异大的图像更友好。注意融合前先把 A 和 B 的全局均值记录一下两幅图亮度差异过大时低频平均会把光照差带进结果里。我通常先做直方图匹配或者减去各自均值再分解融合完再加回来。这个预处理比调融合规则更影响最终效果。这块调整花掉的时间往往比调融合规则本身还多但做完之后你会对 NSCT 系数结构有更深的理解。另外一个经验是保存系数用 mat 文件不要用 imwrite 存子带图子带是 double 数组且有负值存成图片会丢信息。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/13 19:33:27 · 阅读 1.2k NSCT的MATLAB实现：工具包调用、参数选择与避坑指南
简介这是面向MATLAB用户的NSCT非下采样剪切波变换实现工具包目标使用者是图像处理研究者、工程师及刚接触多尺度分析的学习者作为一种多分辨率分析工具它比传统小波更擅长捕捉图像边缘等几何结构可用于理解非下采样金字塔与方向滤波器组原理并直接支撑图像去噪、融合、边缘检测和特征提取等任务。资源共48个文件压缩包约293KB主体为28个m函数和脚本同时包含6张tif测试图MRI、PET、SPECT等、mexw64/mexmac/dll预编译接口、C源码与说明文档兼顾跨平台可运行性和源码级可读性。已有1386人浏览学习。核心代码覆盖NSCT正逆变换完整流程如nsctdec/nsctrec、nsdfbdec/nsdfbrec、atrousfilters、zconv2等模块并配有source images测试图和decdemo演示脚本可完成从图像加载、多尺度方向分解、系数操作到重构对比的完整实验对照特别适合对照论文和教材进行算法复现与验证C源码帮助理解并修改预编译模块mex接口省去自行编译的繁琐步骤。整体结构紧凑、入口清晰是一份能直接上手并深入剖析NSCT算法细节的MATLAB工具包。1. 找NSCT的MATLAB实现前先认清这个变换的代价与回报搜 NSCT 程序实现和 matlab 工具包的人绝大多数最后都卡在同一个地方分解和重构都能跑但重构回来的图和原图差了一圈亮边或者误差大得离谱。NSCT非下采样轮廓波变换不是 MATLAB 内置函数它通常以社区流传的 MATLAB 工具包形式出现这套工具包把小波变换那种横平竖直的分解升级成多个楔形方向分解而且全程不下采样。因此它在图像融合、去噪、边缘增强里表现好代价是系数冗余高、速度和内存开销大。这篇笔记适合要做融合或去噪实验却被 NSCT 参数和系数结构卡住的人我把工具包的调用方式、参数选择和我在实践中踩过的坑一次讲清。2. NSCT的MATLAB实现骨架非下采样金字塔配合方向滤波器组到底做了什么事2.1 金字塔不做下采样细节子带才会和原图一样大NSCT 的第一个模块叫非下采样金字塔NSP。它的任务是把图像按尺度分开低频分量一层层变平滑每一层减去上一层的低频就得到带通细节。关键点是整个过程不做抽取也就是不丢弃任何像素。滤波器靠插入零值来“膨胀”支撑范围每一层用到的等效滤波器越来越长覆盖的频率区域越来越窄但输出图像始终和原图尺寸一致。这个设计解决了离散小波变换的一个老问题小波每分解一层子带尺寸减半平移图像之后再分解系数分布会剧烈变化因为下采样点对不准。NSCT 不做下采样图像平移几个像素系数只是跟着平移不会产生混叠噪声。这就是所谓的平移不变性。代价也很直观原来小波 1024×1024 分解三层后最小子带只有 128×128NSCT 里每一层每种子带都是 1024×1024数据量直接膨胀几十倍。2.2 方向滤波器组用扇形滤波器把频域切成楔形而不是横平竖直NSCT 的第二个模块叫非下采样方向滤波器组NSDFB。它把每一层的高频细节按角度进一步切分二维频域被分成一个个楔形区域每个楔形对应一个方向子带。方向数按 2 的幂次增长比如一个尺度 8 个方向下一层细尺度可能 16 个方向。和小波只有水平、垂直、对角三个方向相比NSCT 能沿任意角度贴合纹理和曲线边缘。工具包内部用的是扇形滤波器组加梅花采样来切换方向坐标。理解到这一步就够了你不需要自己去设计那些扇形滤波器和梅花矩阵工具包已经把方向划分封装成黑匣子。需要做的是选择方向数然后记住一件事——方向子带和细节子带一样尺寸都等于原图尺寸。这个特性在做图像融合时特别方便系数之间不需要插值对齐直接逐像素比较就行。2.3 别走弯路在Contourlet上去掉下采样结果不是NSCT很多第一次接触 NSCT 的人会想Contourlet 和 NSCT 不就是差一个下采样吗把 Contourlet 的代码里抽取那几步删掉不就行了我试过结果重建出来全是条纹。原因在于 Contourlet 的下采样和它的滤波器组设计是一体的去掉抽取之后频域划分不再满足完全重构条件。NSCT 不是简单“不下采样的Contourlet”它的金字塔和方向滤波器都重新设计过才能在无下采样的情况下保持精确重构。这也是为什么我建议直接用社区流传的 MATLAB 工具包而不是自己从论文图纸开始写。NSCT 的滤波器组设计本身是论文级工作包含双通道滤波器、扇形滤波器、梅花采样矩阵的配合自己实现一遍的调试周期很长。工具包也不是十全十美但它的分解和重构代码已经调好拿来做实验比从零写可靠得多。2.4 NSP的atrous思路示意一段看懂“膨胀卷积”的代码理解 NSCT 的尺度分解看下面这段结构示意代码就够了。它展示的是 atrous 算法的核心即滤波器插零膨胀后与图像卷积function [lp, hp] nsp_step(lp_in, g0, g1, rate) % 单层非下采样金字塔atrous 思路示意 % lp_in : 当前低频图像 % g0/g1 : 低通/高通滤波器系数行向量 % rate : 膨胀间隔第 j 层传 2^(j-1) g0_up dyadup(g0, rate); % 在系数间插零撑大滤波器支撑 g1_up dyadup(g1, rate); % 先做行方向卷积再做列方向卷积保持图像尺寸 lp conv2(wextend(2D, sym, lp_in, length(g0_up)), g0_up, same); lp conv2(lp, g0_up, same); hp lp_in - lp; % 带通细节 输入低频 - 新低频 hp hp(1:size(lp_in, 1), 1:size(lp_in, 2)); end逻辑说明第一层 rate 为 1滤波器不膨胀第二层 rate 为 2滤波器每隔一个系数塞一个零等效支撑变长第三层 rate 为 4以此类推。这样每一层看到的频率范围更窄但输出尺寸不变。参数说明g0 和 g1 是一维滤波器系数这里用可分离卷积实现二维滤波wextend 的 sym 模式做对称延拓避免卷积后边缘收缩dyadup 来自小波工具箱。这段代码只是理解思路用的示意工具包里的真实实现还包含双通道重构条件不建议直接拿去替换 nsctdec。3. 把NSCT工具包跑起来nsctdec/nsctrec的调用顺序与参数选择3.1 安装验证把工具包放进路径后先确认nsctdec和nsctrec在不在下载得到的 NSCT 工具包一般是一个目录里面有一堆 .m 文件和子目录。第一步是把它加进 MATLAB 搜索路径。我习惯用 genpath 把整个目录树都加进去避免漏掉子目录里的滤波器初始化脚本dest fullfile(你的本地路径, NSCTToolbox); % 改成实际解压目录 addpath(genpath(dest)); if exist(nsctdec, file) 2 disp(nsctdec 已在路径上可以继续); else disp(nsctdec 未找到检查目录路径是否正确); end逻辑说明exist 返回 2 表示找到同名 .m 文件返回 0 表示不在路径上。这一步能过滤掉大量“明明下载了却报错”的情况。参数说明dest 路径不要带中文字符老工具包对非 ASCII 路径支持不好容易出现奇怪的加载失败。验证通过后我建议先运行工具包自带的 demo 脚本一般是以 demo 开头的 .m 文件跑通了再开始改自己的参数。3.2 最小分解与重构一组不会出错的初始参数初次上手不要一上来就追求 16 方向五层分解。先用一组保守参数跑通流程确认工具包在你机器上工作正常im im2double(imread(cameraman.tif)); % 任意灰度图都行 levels 4; % 金字塔层数 directions [8 4 2 1]; % 从细到粗各层方向数最后一个必须是 1 pfilter 9-7; % 金字塔滤波器 dfilter dmaxflat7; % 方向滤波器 dec nsctdec(im, levels, directions, pfilter, dfilter); rec nsctrec(dec, pfilter, dfilter); err max(abs(im(:) - rec(:))); fprintf(重构最大误差: %e\n, err);逻辑说明nsctdec 把图像分解成系数元胞nsctrec 从系数重构图像。重构误差在 1e-10 量级说明工具包安装正确分解和重构互为逆过程。参数说明im2double 把图像归一化到 [0,1]避免 uint8 的整型溢出问题levels4 是常见的折中既能看到多尺度效果计算量又不至于失控directions 数组长度必须等于 levels且从细尺度到粗尺度排列最后一个固定为 1代表最粗那层只有一个低频近似。3.3 四个必调参数levels、directions、pfilter、dfilter各自怎么改参数常用取值作用我一般怎么选levels2 到 5金字塔分解层数实验对比用 4快速验证用 2directions[4 2 1]、[8 4 2 1] 等每层方向子带数2 的幂次末位为 1纹理复杂用 8 起步平滑图像用 4pfilter9-7、maxflat金字塔滤波器默认 9-7速度和精度均衡dfilterdmaxflat7、vk方向滤波器论文复现优先 dmaxflat7逻辑说明directions 的每个值对应一个金字塔尺度的方向分解数不是随便填的。比如 directions[8 4 2 1]最细尺度 8 个方向第二尺度 4 个方向第三尺度 2 个方向最粗尺度 1 个。方向数从细到粗递减这是 NSCT 的标准设计。参数说明9-7 指的是类似 JPEG2000 的 9/7 双正交滤波器重建质量和速度都比较理想maxflat 支撑更长边缘保留更好但更慢dmaxflat7 方向选择性好纹理丰富时比 vk 稳。提示如果下载的工具包版本入口参数顺序不一样最快的确认方式是打开自带 demo 脚本看调用行我手头这份常用的是 nsctdec(x, levels, directions, pfilter, dfilter) 这个顺序。3.4 重构误差的量化检查为什么高频系数不能随便置零很多人拿到 NSCT 后第一件事就是把高频系数置零想看看低通效果结果发现重构图像出现条纹或亮度偏移。先做个基线实验认识高频系数的贡献decZero dec; % 复制上一小节的分解结果 for l 1:numel(dec) - 1 % 遍历所有方向尺度层保留最后的低频 for d 1:numel(dec{l}) decZero{l}{d} zeros(size(dec{l}{d})); % 高频全部清零 end end recZero nsctrec(decZero, pfilter, dfilter); errZero max(abs(im(:) - recZero(:))); fprintf(高频清零后重构误差: %e\n, errZero);逻辑说明高频清零后重结构相当于只保留低频近似原本的高频细节丢失误差自然巨大。这不是工具包有问题而是 NSCT 的系数不具备“随便改一个系数不影响其他区域”的性质。参数说明numel(dec)-1 是因为最后一个元胞元素是低频近似不能清零。这个实验还有个作用就是验证你对系数结构的理解是否正确——如果清零后重构出来的是一张模糊但正常的低通图说明结构索引用对了。4. 系数到底存在哪里NSCT输出的元胞结构、方向索引与平移不变性验证4.1 输出结构长什么样dec{1}到dec{end}的层级关系nsctdec 返回的是一个嵌套元胞数组这是新手最先懵的地方。结构并不复杂外层元胞按尺度分层dec{1} 是最细尺度dec{2} 是第二细尺度以此类推dec{end} 是最后一个低频近似。每个 dec{l} 又是一个元胞第 l 层的方向子带都装在这里面。内层元胞的元素个数就是 directions(l)。numel(dec) % 外层元素个数 levels 1加了最后的低频 numel(dec{1}) % 最细尺度方向子带数量应该等于 directions(1) size(dec{1}{1}) % 方向子带尺寸和原图完全一致 size(dec{end}) % 低频近似尺寸也和原图一致逻辑说明NSCT 由于全程不下采样从最细方向子带到最后的低频近似所有图像尺寸都等于原图尺寸。这是 NSCT 区别于小波最直观的一点。参数说明numel(dec) 等于 levels1是因为最后一层低频被单独放在 dec{end}而 dec{1} 到 dec{end-1} 才是真正的方向细节层。4.2 方向子带索引规则与尺寸恒定性取某个方向子带的写法是 dec{层号}{方向号}。比如 dec{1}{1} 是最细尺度第一个方向dec{2}{3} 是第二尺度第三个方向。和 directions 数组对应dec{l} 里的元素个数等于 directions(l)。低频近似单独放在 dec{end}不再套内层元胞。subFine dec{1}{1}; % 最细尺度第 1 方向 subMid dec{2}{2}; % 第二尺度第 2 方向 base dec{end}; % 低频近似 figure; subplot(1, 3, 1); imshow(mat2gray(abs(subFine))); title(细尺度方向); subplot(1, 3, 2); imshow(mat2gray(abs(subMid))); title(中尺度方向); subplot(1, 3, 3); imshow(mat2gray(base)); title(低频近似);逻辑说明方向子带系数有正有负直接 imshow 会把负值截断所以用 mat2gray(abs(...)) 把绝对值映射到 [0,1] 再显示。参数说明方向子带里绝对值大的位置代表该方向上的强边缘或纹理低频近似 base 是平滑后的图像保留了整体亮度分布。如果显示出来的方向子带全是噪声状的小颗粒先检查 directions 是不是设成了 [1 2 4] 这种从粗到细的顺序方向数排列反了会拿到错误的子带。4.3 一个验证平移不变性的最小实验NSCT 的平移不变性是它相对 Contourlet 的核心卖点但验证的时候有个常见错误图像平移后再分解直接拿两个分解系数逐像素相减发现差异巨大然后误以为 NSCT 不具备平移不变性。实际上图像平移了系数也跟着平移必须先平移回来再比。shiftRows 9; shiftCols 11; decShift nsctdec(circshift(im, [shiftRows, shiftCols]), levels, directions, pfilter, dfilter); sub1 circshift(dec{1}{1}, [shiftRows, shiftCols]); % 系数平移回来 sub2 decShift{1}{1}; % 平移后图像的分解 residual max(abs(sub1(:) - sub2(:))); fprintf(平移不变性验证残差: %e\n, residual);逻辑说明先把原图分解的系数平移回去再和平移后图像的分解系数比较。理想情况下残差应该非常小实际因为边界延拓差异会有一个不算大的残差但远小于子带本身的幅值。参数说明circshift 是循环移位边界处两边的差异会比内部大想看内部区域的平移不变性可以把比较范围裁掉边界几十个像素再算残差。这个方法比看论文里的示意图靠谱能直接量化你手上这个工具包的平移不变性水平。5. 避坑指南重构残差、方向数约束、内存爆炸与边界伪影的5条记录5.1 重构残差大、边缘一圈亮线先查边界延拓与图像尺寸现象用 nsctdec 分解再 nsctrec 重构中心区域没问题但图像四周出现一圈亮线或暗线重构误差在 1e-2 量级。原因NSCT 的滤波器支撑会随层数膨胀默认延拓宽度不够边缘像素在卷积时没有足够的对称数据可用另外图像尺寸太小的时候最粗那层的滤波器支撑已经比图像还宽边缘伪影会污染整个子带。解决分解前手动对称延拓重构后裁剪回来我一般加 20 像素左右ext 20; imPadded padarray(im, [ext ext], symmetric); decP nsctdec(imPadded, levels, directions, pfilter, dfilter); recP nsctrec(decP, pfilter, dfilter); rec recP(ext1:end-ext, ext1:end-ext); % 裁掉延拓区域逻辑说明延拓让边缘像素在卷积时有完整的邻域数据可用重构后再裁剪边缘伪影被隔离在延拓区域里。参数说明ext 的取值要大于最大滤波器支撑的一半20 像素对 4 层分解通常够用如果分解层数多或滤波器支撑长加大到 32。5.2 directions参数报错或子带对不上方向数数组要和金字塔层数严格配套现象设置 directions[8 4 2] 或 [16 8 4 1] 时报错或者跑完发现某个层方向子带数量不对。原因directions 数组长度必须等于 levels且每个值必须是 2 的幂次最后一个必须为 1。有些版本还会校验相邻尺度的方向数关系细尺度的方向数不能超过上一层方向数的合理倍数。解决先跑最小配置验证再逐步加方向数。levels 2; directions [2 1]; % 最小验证配置 decTest nsctdec(im, levels, directions, pfilter, dfilter); fprintf(方向子带数: %d, %d\n, numel(decTest{1}), numel(decTest{2}));逻辑说明方向数翻倍是有结构约束的不是想设多少就设多少。最粗尺度低频层只有 1 个方向往细的方向数按 2 的幂次增加而且相邻层之间要保持合理的倍数关系。参数说明从 [2 1] 验证通过后再跳到 [4 2 1]、[8 4 2 1]一次加一层能快速定位是哪个参数组合触发的报错。实际实验里 [8 4 2 1] 是绝大多数图像都能稳定跑通的配置。5.3 重建后偏暗或偏亮能量对不上检查图像类型与滤波器组合现象分解再重构图像整体变暗或变亮PSNR 掉到 20 以下但局部结构还在。原因最常见的是输入图像没有转成 double 或归一化uint8 图像直接送进 nsctdec卷积时溢出截断其次是 pfilter 和 dfilter 的组合来自不同工具包版本重构增益不为 1。解决统一走 im2double 再分解如果确认输入没问题换回 9-7 加 dmaxflat7 这个经过大量验证的组合。img imread(some_image.png); im im2double(img); % 先归一化再分解 if size(im, 3) 3 im rgb2gray(im); % NSCT 工具包默认按灰度图处理彩色要逐通道跑 end逻辑说明NSCT 工具包内部假设输入是 double 类型uint8 的整数卷积会产生截断。彩色图像如果直接送进去工具包会把第三维当作批量维度处理结果各通道混在一起。参数说明im2double 对所有格式的输入都安全彩色转灰度用 rgb2gray或者保持彩色但逐通道分解再合并后者适合需要颜色信息的融合实验。5.4 大图一跑就内存不足冗余变换的代价与降参数建议现象512×512 的彩色图levels5directions[32 16 8 4 1]运行到一半提示内存不足或者等了十几分钟还没结束。原因NSCT 是冗余变换系数总量约为原图像素数的 1sum(directions) 倍这个配置下接近 62 倍。512×512 灰度图就有约 1.6 亿个 double 值加上卷积过程的临时数组内存很容易爆。解决先降方向数再降层数或者逐通道处理。levels 4; directions [8 4 2 1]; % 比 32 方向温和得多 decR nsctdec(im(:,:,1), levels, directions, pfilter, dfilter); decG nsctdec(im(:,:,2), levels, directions, pfilter, dfilter); decB nsctdec(im(:,:,3), levels, directions, pfilter, dfilter);逻辑说明把彩色三通道分别分解内存峰值比一次性处理三维数据低很多。方向数是内存增长的主要因素每多一层方向数翻倍数据量就翻一倍。参数说明一般实验 8 个方向已经能体现 NSCT 的优势16 方向只在纹理非常密集时有必要4 层 8 方向是我测试多台机器都不会内存崩溃的配置。速度方面长支撑方向滤波器是大头vk 比 dmaxflat7 略快纹理要求不高时可以换。5.5 系数置零后重建出现条纹方向子带之间存在频域重叠现象把某一个方向子带全部置零重构图像出现均匀分布的方向性条纹而且条纹方向和置零子带的方向相关。原因NSDFB 的方向子带在频域有过渡带重叠方向子带并不是完全正交的。强行挖空一个楔形频带相当于在频域做了一个尖锐切口重构时产生 Gibbs 振荡。解决不要硬置零改用系数收缩或衰减如果只想抑制某个方向分量把系数乘以一个小系数而不是清零。decAtten dec; decAtten{1}{3} decAtten{1}{3} * 0.1; % 衰减而不是清零 recAtten nsctrec(decAtten, pfilter, dfilter);逻辑说明乘以 0.1 保留了该方向的大部分结构只是降低了能量频域的过渡带变化平缓不会产生明显条纹。参数说明衰减系数在 0 到 1 之间0.05 到 0.2 是常用区间如果确实需要彻底移除某个方向分量建议配合局部区域平滑而不是全局硬置零。这条经验在做去噪实验时特别有用很多去噪算法直接对高频系数做硬阈值结果出现条纹换成软阈值就好了。6. 进阶把NSCT系数用进图像融合实验的最小改动6.1 区域能量取大的融合规则以及它为什么比简单平均稳NSCT 最常见的落地场景是图像融合。两幅图像先各自做 NSCT 分解低频部分用平均高频方向子带用能量取大。低频决定整体亮度和结构平均融合能避免亮度突变高频子带里绝对值大的系数代表更强的边缘或纹理取大规则能保留更清晰的细节。decA nsctdec(A, levels, directions, pfilter, dfilter); decB nsctdec(B, levels, directions, pfilter, dfilter); decF decA; % 先按 A 的结构复制 for l 1:numel(decA) - 1 % 跳过最后的低频 for d 1:numel(decA{l}) eA sum(decA{l}{d}(:).^2); % 全子带能量 eB sum(decB{l}{d}(:).^2); if eA eB decF{l}{d} decA{l}{d}; else decF{l}{d} decB{l}{d}; end end end decF{end} (decA{end} decB{end}) / 2; % 低频平均保留整体亮度 F nsctrec(decF, pfilter, dfilter);逻辑说明全子带能量取大是最简版本它把整幅子带交给能量更高的一方边缘保留效果已经优于空间域平均。参数说明想要更精细把 sum(...)^2 改成局部滑窗求和比如 5×5 邻域内系数平方和这样融合边界会更自然。低频也可以用局部窗口平均替代全局平均对亮度差异大的图像更友好。注意融合前先把 A 和 B 的全局均值记录一下两幅图亮度差异过大时低频平均会把光照差带进结果里。我通常先做直方图匹配或者减去各自均值再分解融合完再加回来。这个预处理比调融合规则更影响最终效果。这块调整花掉的时间往往比调融合规则本身还多但做完之后你会对 NSCT 系数结构有更深的理解。另外一个经验是保存系数用 mat 文件不要用 imwrite 存子带图子带是 double 数组且有负值存成图片会丢信息。希望帮到你。本文还有配套的精品资源点击获取