简介NSCT程序实现的MATLAB工具包面向图像处理研究者与工程师提供非下采样剪切波变换的完整算法实现。该变换融合多尺度分析与方向敏感性可高效捕捉复杂图像的边缘与几何结构适用于去噪、压缩、边缘检测及特征提取等场景。RAR压缩包共48个文件以28个m函数脚本承载核心变换与重构算法6个tif图像作为标准测试样例配合c源码、mexw64/mexmac/dll等编译动态库实现跨平台加速另有txt与asv文件辅助说明整体仅293KB结构紧凑、易于扩展。已有1386人学习下载。借助该工具包读者可直接调用nsctdec执行正变换、isct2完成逆变换利用演示脚本加载示例图像观察系数分解与重建效果从阈值去噪到压缩编码各环节均有函数封装代码注释清晰适合科研入门、算法对比与工程二次开发。1. 从 NSCT 到 MATLAB 工具包为什么图像处理先跑通它做图像融合和去噪时我先后把手里的 MATLAB 工具包换过三轮最后还是回到 NSCT 程序实现这套非下采样轮廓波变换工具包。NSCT 在保留边缘细节上明显强于小波和 Contourlet它有平移不变性方向子带可分区调整最适合多尺度图像融合、阈值去噪和纹理特征提取。下面内容基于我实际拆过的工程写法讲清楚怎么调用核心函数、分解层数和方向数怎么设、滤波器选 9-7 还是 pkva以及最容易踩的五个坑。看完你就能把工具包跑通并把它嵌进自己的融合或去噪流程。2. 非下采样轮廓波变换分解结构与工具箱核心函数2.1 分解原理NSP 与 NSDFB 的分工NSCT 的分解链路分两步走。第一步是非下采样金字塔 NSP第二步是非下采样方向滤波器组 NSDFB。NSP 负责“尺度”NSDFB 负责“方向”。NSP 的结构类似拉普拉斯金字塔但去掉了下采样。每一级通过一对分解滤波器得到低频逼近和带通高频。下一级的输入是上一级低频但滤波前会先对滤波器系数进行插值上采样。这样一来每一级输出的子带尺寸都和原始图像一样大而且目标在图像中的位置只要平移各级子带系数也会平移相同像素不会像传统金字塔那样出现下采样造成的周期搬移。这个特性在图像配准、融合、去噪中非常宝贵。NSDFB 的作用是把每个高频带通子带按方向切分。它的基础模块是扇形滤波器组通过调制和重采样操作将二维频谱划分成楔形方向块。每做一次方向分解方向数翻倍分解级数为 L 时方向子带数是 2^L。所以工具包里nlevels给的是级数不是方向个数。在第三层给3代表 8 个方向子带。两层结构组合起来就是 NSCT。以三层为例NSP 先分出三层低频与三层高频对第一层高频最细尺度用较少方向级数对第二层、第三层高频逐渐增加方向级数最后得到类似楔形剖分的频率结构。因为每一步都保持了下采样不变NSCT 的冗余度大约在 4/3 到 8/3 之间具体取决于方向和层数。冗余带来的好处是稳定代价是内存和计算量成倍上升这是后面章节的踩坑核心。2.2 工具箱文件组成从主函数到滤波器组拿到工具包解压后先看目录结构。绝大多数版本的 NSCT 代码基于作者公开的测试框架文件命名规律如下文件名作用nsctdec.mNSCT 分解主函数封装 NSP 和 NSDFBnsctrec.mNSCT 重构主函数atrousdec.m非下采样金字塔分解atrousrec.m非下采样金字塔重构ddfbdec.m方向滤波器组分解ddfbrec.m方向滤波器组重构dfilters.m生成 9-7 与 pkva 的滤波器系数modulate2.m滤波器调制用于方向频谱搬移resampc.m重采样操作配合方向滤波zconv2.m / upfir.m二维卷积与上采样滤波基础函数这些文件之间存在明确的调用关系nsctdec先调atrousdec得到低频和高频金字塔再对每个高频层调ddfbdec完成方向分解。nsctrec反着来先调ddfbrec合成方向子带再调atrousrec重建低频。我建议你拿到工具包后第一件事就是打开nsctdec.m看前 20 行注释。因为不同版本扩展标记差很多有些函数名是nsctdec有些是nsctdec2还有的用 C 编译的 mex 文件替代zconv2。先确认主函数名和输入参数避免后面用错。如果目录里有mex_*.c说明需要先编译 mex 才能保证速度如果没有说明纯 m 代码也能跑只是慢一些。2.3 核心函数 nsctdec / nsctrec 的调用接口与参数含义核心调用代码[y, DFB] nsctdec(x, nlevels, dfilters, ext); xr nsctrec(y);四个输入参数的含义如下x二维 double 灰度矩阵尺寸建议偶数行偶数列。非 double 会先被强制转换但转换后可能丢失精度所以最好自己做转换。nlevels正整数行向量长度等于 NSCT 分解层数。每个元素表示该层方向滤波器组的分解级数实际方向数为2^nlevels(i)。例如[2 3]表示两层第一层 4 个方向第二层 8 个方向。dfilters滤波器类型。最常用的是9-7和pkva也可以给元胞数组按层指定比如{9-7, pkva}。ext边界扩展方式per是周期扩展sym是对称扩展。默认值在不同版本里不一样显式传入最稳。返回值y是嵌套元胞数组。y{1}保存低频子带是个和x同尺寸的矩阵y{2}是元胞数组y{2}{i}对应第 i 层y{2}{i}{j}是该层第 j 个方向子带。DFB是可选的第二个输出返回方向滤波器组设计信息初学者基本用不上但调试时看里面的h和g滤波器系数能知道每个方向是如何划分的。给出一个查看系数结构的代码x double(rand(128, 128)); y nsctdec(x, [2 3], 9-7, per); whos y fprintf(低频尺寸: %dx%d\n, size(y{1})); for i 1:2 fprintf(第%d层方向数: %d\n, i, length(y{2}{i})); end这段代码用随机小图测试分解whos会列出每个变量的尺寸。低频矩阵依然是 128x128每个方向子带也是 128x128这正是非下采样的直接证据。如果你看到某个子带尺寸不同说明工具包版本可能不是真正的 NSCT或者参数里有下采样选项。注意x需要是 double不要传 uint8。如果直接传imread的结果工具包内部虽然能跑但部分滤波操作会被隐式转换导致重构误差变大。所以统一用double(imread(...))。对 RGB 图像不能直接把三维矩阵传给nsctdec。我一般拆成三个通道分别分解融合function F nsct_fuse_rgb(A, B, nlevels, dfilters, ext) F zeros(size(A)); for c 1:3 yA nsctdec(double(A(:,:,c)), nlevels, dfilters, ext); yB nsctdec(double(B(:,:,c)), nlevels, dfilters, ext); % 融合规则在前面章节写这里需要补全 yF % F(:,:,c) nsctrec(yF); end end注意这里的yF需要按规则生成上面只是示意。逐通道处理虽然循环三次但换来了稳定的结果也比构造三维延展去钻工具包内部可靠得多。还有一个容易忽略的细节nsctrec(y)只接受一个参数它会从y的结构中反推分解层数和方向数。所以你在分解后如果手动删除了y{2}{i}中的某个方向子带重构通常会报错或产生错误尺寸因为滤波器组要求和分解时一致。如果不想保留全部方向子带建议用零矩阵填补而不是删除元素保证结构完整。3. 跑通第一个 NSCT 实验图像融合与去噪的完整流程3.1 准备工作安装工具包与设置路径下载解压后把文件夹放到一个无中文路径的位置因为 MATLAB 在某些中文路径下对 mex 调用会出问题。然后在命令窗口执行addpath(genpath(D:\work\matlab\nsct_toolbox)); savepath; % 保存路径避免每次重启都设置如果文件夹里包含.c源文件通常还需要编译 mexmex -setup mex mex_conv2.cmex -setup选择 C 编译器然后在工具包目录下编译所有.c文件。若使用的是纯 m 版本这一步可以跳过。经验上有 mex 的版本速度快 5~10 倍做融合实验时强烈建议编译。验证是否到位help nsctdec which nsctdec如果显示路径和帮助信息说明添加成功。如果提示未找到检查genpath是否用对了或者路径名里是否包含空格。注意不要用addpath(nsct_toolbox)加根目录而不带子目录工具包的内部子函数都在nsct_toolbox/下的多个子文件夹里根目录只有主函数找不到滤波器函数时会被误认成“函数未定义”。提示换电脑或换 MATLAB 版本后旧编译过的 mex 文件可能失效运行时会报Invalid MEX-file。此时删除mex_*.mex*重新编译即可。3.2 分解与重构验证往返误差在正式做融合前先跑往返测试确认工具包在你的环境下没问题。代码x double(imread(cameraman.tif)); % 换成自己的灰度图 y nsctdec(x, [2 3], 9-7, per); xr nsctrec(y); mse mean((x(:) - xr(:)).^2); fprintf(MSE %.6f\n, mse);这里用cameraman.tif是内建示例图实际上你可以用任意灰度图。[2 3]是两层分解配置第一层 4 方向第二层 8 方向速度和内存都适中。x转成 double 是因为滤波卷积需要浮点运算如果不转MATLAB 会路过整数运算重构误差会被量化放大。逻辑说明nsctdec把图像分解为低频和高频子带nsctrec原路重建。理想情况下这两个函数是完全自洽的MSE 应该接近 1e-10 甚至 0。如果你测出来 MSE 大于 1e-5优先怀疑输入是否 double、是否改了y的某个子带。这一步就像测试一台黑匣子只有往返误差通过了后面你才敢对系数做修改。参数说明mse的阈值我一般允许到 1e-8。如果到不了就检查扩展模式是否支持你的图像尺寸。例如[1 2 3]三层分解时高层滤波器上采样倍数大某些边界扩展在奇数尺寸下会产生不对称填充误差会明显偏大。最省事是先把图像裁剪成偶数尺寸。3.3 基于 NSCT 系数的图像融合实现融合是 NSCT 最经典的应用。假设有两幅配准好的灰度图像 A 和 B我们要得到一张包含两者清晰内容的图像。常见规则是低频取平均、高频取模极大A double(imread(A.tif)); B double(imread(B.tif)); nlevels [2 3]; d 9-7; ext per; [yA, DFB] nsctdec(A, nlevels, d, ext); [yB, ~] nsctdec(B, nlevels, d, ext); yF yA; yF{1} (yA{1} yB{1}) / 2; for i 1:length(nlevels) for j 1:length(yA{2}{i}) Acoef yA{2}{i}{j}; Bcoef yB{2}{i}{j}; Fcoef Acoef; mask abs(Bcoef) abs(Acoef); Fcoef(mask) Bcoef(mask); yF{2}{i}{j} Fcoef; end end F nsctrec(yF); F uint8(F);逻辑说明yA{1}是低频逼近反映了图像的整体亮度分布。两幅图的低频直接取平均可以防止融合结果偏向某一幅太亮或太暗。高频方向子带的系数绝对值代表边缘强度取绝对值大的系数可以保留更锐利的轮廓。mask是布尔矩阵选出 B 的系数比 A 强的地方把该位置的 A 系数替换成 B 系数。这个过程对每个方向每个位置做点对点比较效率高效果直观。参数说明如果图像是 RGB需要逐通道重复上述过程或者转为灰度处理但那样会丢失色彩信息。多聚焦融合用灰度即可医学图像或者遥感图像建议逐通道避免因为转换导致色偏。这里的d和ext选择会影响方向子带质量想获得更平滑的方向过渡可以用pkva代价是边界振铃比9-7略大。3.4 基于 NSCT 的阈值去噪实现去噪的思路是把噪声看成高频分量对高频系数做阈值收缩。实现x double(imread(noisy.tif)); y nsctdec(x, [2 3], 9-7, per); thr 3 * median(abs(y{2}{1}{1}(:))) / 0.6745; for i 1:length(nlevels) for j 1:length(y{2}{i}) coef y{2}{i}{j}; coef sign(coef) .* max(abs(coef) - thr, 0); y{2}{i}{j} coef; end end xr nsctrec(y); xr uint8(xr);逻辑说明阈值thr基于第一层第一个方向子带的系数中位数median(abs(...))计算这是一种经典的高斯噪声方差估计。常数 0.6745 是标准正态分布的中位绝对偏差与标准差的换算关系。系数绝对值大于阈值的部分保留小于等于阈值的被置零这样能压制小幅噪声波动。软阈值函数sign(coef) .* max(abs(coef)-thr, 0)避免在阈值处产生硬跳变。参数说明3是缩放系数数值越大去噪越强但也越容易抹平细节。想保守一点可以改成2.5想激进用4。也可以对每一层分别设阈值因为不同尺度的噪声标准差不完全一样。通常第 1 层噪声最多高层阈值可以递减。4. 参数怎么调分解层数、方向数、滤波器选择的边界与组合4.1 分解层数 nlevels 对频带划分的影响先说最容易搞混的地方nlevels不是一个方向数列表。比如[2 3]是两层分解第一层方向分解级数为 2第二层为 3实际方向数分别是 4 和 8。很多论文写“方向数为 [2 3]”严格来说不准确工具包内部读的是级数。分解层数由nlevels的元素个数决定。元素个数越多低频被逐级抽稀每层高频占据的尺度范围越细。举个例子[1 2 3]是三层分解第一层方向数 2第二层 4第三层 8[2 3 4]也是三层但方向数 4、8、16。后者对斜向纹理的刻画更细但系数总量从 124815 个尺度带变成 1481629 个子带内存和耗时差不多翻倍。在实际任务里我的默认选择是融合用[2 3]或[2 3 4]去噪用[2 3]特征提取用[1 2 3]。分解层数超过 4 层后最粗尺度的低频几乎是一个纯背景灰度图高频也大多是零散边缘继续加深只会增加重构累积误差不会带来肉眼可见的收益。如果图像本身只有 256x2564 层 NSCT 的最粗尺度已经只有 16x16 的区域信息方向分解再细也没有意义。4.2 方向数怎么给从 [2 3 4] 到 [1 2 3 4] 的取舍方向数选择的核心矛盾是方向越多边缘角度分辨率越高但每个子带的能量越稀疏后续处理越容易受到单点噪声影响。我用一张含斜纹理的测试图做过对比[1 2 3]去噪后45 度斜线边缘明显发虚[2 3 4]下同样参数斜线边缘保真度好很多但耗时是前者的 3 倍内存占用也翻倍。推荐的配置组合[1 2 3]三层方向数 2、4、8。计算量最小适合快速验证流程。[2 3]两层方向数 4、8。融合和去噪的均衡选择。[2 3 4]三层方向数 4、8、16。适合对边缘方向性要求高的图像比如遥感目标识别前的预处理。[1 2 3 4]四层方向数 2、4、8、16。方向层次最丰富但内存压力大512x512 下慎用。给方向数时还要注意尺度匹配。第一层是最细尺度方向数给多了容易把像素级噪声也分解成方向纹理高层是粗尺度方向数太少会丢失曲线轮廓。所以常见做法是“细尺度少方向、粗尺度多方向”也就是[1 2 3]这样的递增模式而不是[3 2 1]。如果发现某层方向子带里有明显的棋盘格伪影大概率是方向数给得过高该层频谱采样不足。4.3 滤波器类型与扩展模式别忽视边界效应dfilters参数控制两个关键点滤波器的长度和频率响应形状。9-7是 9/7 双正交滤波器线性相位边缘响应平缓重构误差小pkva是短滤波器计算快方向选择性更强但频率选择性略差。我的经验是默认用9-7如果图像纹理密集、需要快速调试换pkva先跑通流程最后再用9-7出正式结果。ext参数管的是图像边界怎么延拓。per是周期延拓把图像左右上下当循环连接对尺寸和频谱均匀的图像效果好sym是对称延拓把边界镜像外扩更符合自然图像边缘不连续的实际情况。但在 NSCT 工具包的老版本里sym的实现偶尔有 bug会导致分解的高频子带尺寸比输入大 2 像素。如果你发现重构结果边缘有一圈亮边就把ext改回per再测。另一个容易忽略的点是输入尺寸的奇偶性。NSP 中的上采样滤波器会让低频尺寸与输入保持一致但高层滤波器长度经过多次上采样后等效长度很长在奇数尺寸下边界延拓不对称重构误差会明显变大。所以我一般在预处理里强制把图像变成偶数尺寸x x(1:floor(size(x,1)/2)*2, 1:floor(size(x,2)/2)*2);4.4 参数组合速查表应用场景nlevelsdfiltersext备注通用图像融合[2 3]9-7per平衡速度与质量高质量融合[2 3 4]9-7sym方向数多注意内存快速去噪[2 3]pkvaper计算快纹理保留一般特征提取[1 2 3]9-7sym低层方向少粗尺度方向多小尺寸图像[1 2]pkvaper128x128 以下推荐这张表不是死规矩。你可以在自己数据集上先跑[2 3] 9-7 per然后用第 6 章的定量指标对比[2 3 4] 9-7 sym挑出最适合的。每次改参数都先跑一遍第 3.2 节的往返误差确认重构质量没崩再继续后续处理。5. 避开五个高频坑NSCT 工具箱使用中的常见问题与排查5.1 问题一分解后尺寸对不上重构结果发虚现象输入 256x256 图像分解后个别子带是 260x260 或 252x252重构后的图整体发虚尺寸也变了。原因ext选择不当导致卷积过程中尺寸偏移。NSCT 工具包的滤波器长度比较长周期延拓时如果图像尺寸是奇数或者延拓方式与滤波器相位不匹配卷积结果就会多出一两行。解决先用size(x)确认输入是偶数行偶数列把ext显式改成per再用第 3.2 节的 MSE 脚本验证。如果确认是滤波器长度问题就在分解前用padarray手动扩展边缘重构后裁剪掉对应像素。裁剪尺寸要匹配滤波器长度不要拍脑袋。换个思路先做往返测试再动参数。5.2 问题二方向数给错函数直接报错现象调用nsctdec(x, [4 8], 9-7, per)报错nlevels must be a vector of nonnegative integers或者断言失败。原因我们把方向数 4 和 8 直接给了nlevels但工具包要求的是方向分解级数方向数等于2^nlevels。[4 8]会被解释成 16 和 256 个方向明显超出合理范围。解决把[4 8]改成[2 3]。记住换算关系方向数 2、4、8、16 分别对应级数 1、2、3、4。我第一次用这个工具包的时候就在这里翻车了对着报错查了半天最后翻nsctdec.m源码注释才明白。如果你不确定先跑nsctdec(x, [1 1], ...)看能不能通过。5.3 问题三默认参数引起的边界失真现象重构图像边缘出现规则亮带像是被强加了一圈边框。原因nsctdec的默认扩展模式在不同版本里不一致有的默认per但你的图像边缘亮度是突变的周期延拓会把左右边缘强制连接产生虚假的横向纹理。解决每次调用都显式传sym或per不要依赖默认值。如果sym下仍有伪影就在分解前用edgetaper函数弱化边缘跳变或者干脆裁剪掉图像四周的暗角区域。这是血泪经验NSCT 工具包不是一个面向生产环境的黑匣子默认参数通常只是“能跑”不是“最优”。5.4 问题四内存不足 / 运行时间过长现象512x512 图像用[2 3 4]分解直接Out of memory或者跑十几分钟不出结果。原因NSCT 非下采样每个方向子带都和原始图像一样大。三层方向数 4、8、16 时高频子带总数是 28 个每个 512x512 的 double 矩阵约 2MB加上金字塔中间结果内存占用轻松超过 300MB。方向数再往上就是指数增长。解决把输入转成single能省一半内存降低分解层数和方向数比如从[2 3 4]降为[2 3]或者把图像分块处理融合后再拼接。分块时相邻块要留重叠区重叠宽度至少等于滤波器半径避免块边缘断层。另外用whos监控变量内存及时clear掉不再用的中间变量。5.5 问题五不同版本工具包函数名不一致现象从网上下了一份工具包主函数叫nsctdec另一份叫nsctdec_nsct或者nsctrec跑起来报“未定义函数”。原因NSCT 工具包有多版本迭代早期版本使用 mex 文件后续版本改纯 m 文件函数命名也有调整。不同开发者二次打包时还顺手改了函数名。解决解压后先读README.txt或Contents.m按里面给的示例调用。不要拿一个版本的脚本直接套另一个版本。如果怕多个版本混用在代码里封装一层兼容层if exist(nsctdec, file) y nsctdec(x, nlevels, dfilters, ext); elseif exist(nsctdec_nsct, file) y nsctdec_nsct(x, nlevels, dfilters, ext); else error(no NSCT function found); end这样换工具包时只改这一处后面的融合、去噪流程不受影响。6. 进阶技巧用系数统计做特征提取与融合规则定制6.1 系数分布与活动度量选 max 还是选 mean高频子带中系数绝对值大小代表局部边缘强度。最常用的融合规则是绝对值取最大但 max 对噪声很敏感一个孤立噪声点就能改变融合结果。更稳的是区域能量在 3x3 邻域内求平方能量localE conv2(coef.^2, ones(3)/9, same);这反映的是系数周围的活跃度而不是单点强度。6.2 自定义融合规则的 MATLAB 实现把第 3.3 节的 max 规则换成区域能量function Fcoef fuse_coef(Acoef, Bcoef) EA conv2(Acoef.^2, ones(3)/9, same); EB conv2(Bcoef.^2, ones(3)/9, same); Fcoef Acoef; Fcoef(EB EA) Bcoef(EB EA); end区域能量比 max 保留了更多柔和纹理但对比度稍低。想兼顾两者可以设定一个阈值当abs(EA-EB)/(EAEB)大于某个值时取能量大的否则取两者均值。这样既保边缘又保平缓区域。6.3 验证结果用熵、互信息、边缘保留度做定量评估融合结果不能只看眼睛还要量化。常见指标E entropy(uint8(F)); % 信息熵 MI mutualinfo(A(:), F(:)) mutualinfo(B(:), F(:)); % 互信息之和 Q qabf(A, B, F); % 边缘保留度entropy是 MATLAB 自带函数mutualinfo需要自己写直方图互信息脚本qabf是融合质量评估的常用函数网上能找到实现。我一般对同一组参数跑完记录这三个值再按应用偏好选参数。从那以后我每次换数据集都会强制跑一遍往返误差、熵和互信息三个脚本再决定是否继续调参避免白白浪费时间。希望帮到你。本文还有配套的精品资源点击获取
