一文搞懂NSCT(Nonsubsampled Contourlet)变换，matlab程序实现并讲解NSCT_TOOLBOX的使用方法

一、前言

在这篇文章中我将从图像分解领域的发展方面为大家讲解NSCT的出现背景和实现原理，最后附上NSCT的matalb实现代码。希望用最短的篇幅和最通俗的语言帮助大家搞懂NSCT（我在学习的时候花了好长时间，非常痛苦，希望痛苦到我为止）。上正文！

二、NSCT理论背景

传统的小波变换虽然具有良好的时频局部化特性，能将图像分解为水平、垂直、对角线三个方向上的高频子带，实现图像在不同频带、不同时段的分离，但由于其支撑基呈正方形，因而对二维图像中更高维的奇异点(如线、边缘、轮廓等特征)无法实现最优稀疏表示，也不能很好地表示图像的方向信息。同时小波变换的实现过程需要通过卷积完成，计算复杂，且不具备平移不变性，易导致重构图像出现明显的吉布斯(Gibbs)现象。为解决上述问题，1998年，Candes提出了Ridgelet变换，该理论有效解决了小波变换在处理二维图像时对方向信息表达不力的问题。脊波是通过对小波基函数添加一个表述方向的参数而得到，具有很强的方向选择和识别能力，并能非常有效地表示图像中诸如线性轮廓、直线信息等方向性特征。然而，对于图像中存在的大量不规则曲线，Ridgelet转换却表现的效果不佳，相比之下，Donoho等人在1999年推出的曲波(Curvelet)变换理论具有很好的方向敏感性。但由于Curvelet是在频域上定义和实现的，因而其在空域上的采样特性表现得并不明显，而且数字实现的计算复杂度通常很高，这促使了Do等人提出了轮廓波(Contourlet)变换理论。轮廓波变换是一种严格意义上的图像多尺度几何转换工具，它成功将小波的优点延伸到了高维空间，仅用少数非零系数就能有效地捕捉图像的光滑轮廓，但由于Contourlet在对图像滤波分解的过程中执行了下采样操作，导致其不具备平移不变性，因而在频域空间存在较严重的频率混叠，当应用于图像融合时同样会在融合图像中出现伪Gibbs效应。2006年，Cunha与Zhou等人将α ̀ trous算法应用到Contourlet的变换当中，发展提出了NSCT理论，该理论很好地解决了上述几种图像多尺度分解工具存在的缺陷，不但继承了轮廓波对图像多尺度、多方向分解的优良特性，还具备了平移不变性。

三、NSCT图像表述

3.1 图像变换综述

NSCT是一种优异的二维图像多尺度分解工具，它不仅能够处理小波变换所不能有效表示的更高维的奇异点，而且能够弥补Contourlet变换中所缺失的平移不变性(Shift-invariance)。NSCT变换中先后使用了非下采样的金字塔(Non-Subsampled Pyramid, NSP)分解滤波器和非下采样的方向滤波器组(Non-Subsampled Filter Banks, NSDFB)对二维图像进行剖分，其过程是首先利用NSP对源图像进行多尺度分解以有效“捕获”图像中的奇异性特征，然后采用NSDFB进一步对高频分量进行多方向分解，从而获得不同尺度、不同方向下的子带图像。NSP每次将图像分解为一个低通子图像和一个带通子图像，随后的分解都是在低通子带上迭代进行。NSDFB将NSP每级分解得到的带通子带分解成多个不同方向的子图像，以提取更加精确的方向性信息。NSCT没有对NSP以及NSDFB分解后的信号分量进分析滤波后的下采样以及综合滤波前的上采样操作，而是先对相应滤波器进行上采样，然后再对信号进行分析滤波和综合滤波。由于NSCT的分解过程中使用的是非下向采样的金字塔滤波器，因而具有平移不变性，这使得分解后的所有子图像大小均相等。

个人见解：所谓实现非下采样的…，其实就是直接对滤波器进行上采样操作。传统的金字塔是对图像进行下采样，所有的是在图像上动(滤波器是不动的)，现在这个就是不动图像，通过动滤波器直接实现原来的“滤波+下采样”的工作。这也是NSCT区别于剪切波变换的地方。

图3-1是应用NSCT对图像分解的示意图，二维图像经NSCT分解后，频域空间被划分成楔形形状的方向子带，如下图所示:

图3-1 应用NSCT对图像分解示意图

3.2 非下采样的金字塔分解

NSCT采用双通道非下采样的二维滤波器组来实现NSP分解机制，分解过程如图3-2所示。

图3-2非下采样金字塔的双通道非下采样塔型滤波器组

其中，塔式分解滤波器$H_0(z),H_1(z)$和综合滤波器$G_0(z),G_1(z)$要求满足$Bezout$约束条件：
$$H_0(z)G_0(z)+H_1(z)G_1(z)=1$$
即满足上式所示的恒等式（此方程式中，$H_0(z)$表示低通分解滤波器，$H_1(z)$为高通分解滤波器，$G_0(z)$为低通重构滤波器，$G_1(z)$是高通重构滤波器）。这种约束方式保证了NSP满足完全重构(Perfect Reconstruction)的条件。上述变换过程与α ̀ trous算法的一维非下采样小波变换十分类似，并且当分解级数为$N$时，冗余度为$N+l$。为实现NSCT对图像的多尺度分解，在执行多级分解的过程中，NSP每一级均需要对前一级滤波器按照采样矩阵$D=2I$($I$为二阶单位矩阵)进行上采样，即对上一级尺度低频信号经上采样后的低通滤波器进行低通滤波，得到塔式分解后的低频信号(低频子带图像);同时还对上一粗糙尺度下的高通滤波器进行高通滤波，得到塔式分解后的高频信号(带通方向子带图像)。图像被NSP分解后，第$j$尺度下对应滤波器的理想传输频带的支撑区间是${\left[-\left(\pi /2^j\right),\left(\pi /2^j\right)\right]}^2$,而高通滤波器所对应的理想传输频带的支撑区间则为上一级低通滤波器的补集，即${\left[\left(-\pi /2^{j-1}\right),\left(\pi /2^{j-1}\right)\right]}^2/{\left[-\left(\pi /2^j\right),\left(\pi /2^j\right)\right]}^2$。随后的每一级滤波器只需对上一级滤波器进行上采样便可得到，而无需进行额外的滤波器设计。二维图像每经过一级NSP分解便可得到一幅低频子带图像和一幅带通子带图像，以后的每一级NSP分解都是在低通分量上通过迭代的方式进行以获取图像中的奇异点。因此，一幅二维图像通过k级NSP分解后，可产生k+1幅与源图像尺寸大小相同的子图像，其中包括1个低通子图像与k幅不同尺度下的带通子带图像。图3-3 (a)描述了一幅图像经三级金字塔分解的过程，其中，图中的$2I$表示对$H_0(z)$进行上采样操作。

图3-3 非下采样金字塔分解示意图

自行补充1：
如果利用一个通用公式来描述K级分解的NSCT变换中任意一级所使用的滤波器（H0和H1是什么），可以用如下公式表达：
K级 级联NSP的等效滤波器:

个人理解：第k级的级联滤波器的意思可以理解为此级下所使用的H0和H1。这个公式的意义在于：用一个公式就说明了金字塔分解的过程。

$$H_n^{eq}(z)=\{\begin{array}{ll}{H}_1\left(z^{2^{n-1}I}\right){\prod }_{k=0}^{n-2}{H}_0\left(z^{2^{k}I}\right),& 1\le n\le K\\ {\prod }_{k-0}^{n-2}{H}_0\left(z^{2^{k}I}\right),& n=K+1\end{array}$$
其中，$2I$仍表示上采样操作，所有的字母设定和前面的一致。$H_n^{eq}$就是第$n$级的级联滤波器，可以理解为第$n$级的$H0$和$H1$。公式中的$H0$和$H1$就是我们最初在第一级分解时使用的初始滤波器。

3.2 非下采样的方向滤波器组分解

非下采样的轮廓波转换中的NSDFB分解机制是在Bamberger和Smith所设计的扇形方向滤波器组的基础上，所构造的一组双通道非下采样滤波器组，如图3-4所示。

3-4 双通道非下采样扇形方向滤波器组

类似地，NSDFB的分析滤波器$U_0(z),U_1(z)$和综合滤波器$V_0(z),V_1(z)$也需满足Bezout恒等式:
$$U_0(z)V_0(z)+U_1(z)V_1(z)=1$$
NSDFB的工作流程如下:首先利用扇形滤波器组及象限滤波器组把图像分解为4个方向上的子带图像，然后再通过平行滤波器组将其分解为各个不同方向上的子带图像。具体地，采用理想的频域支撑区间为扇形的$U_0(z)$和$U_1(z)$滤波器实现双通道的方向性分解。在此基础上，对滤波器$U_0(z)$和$U_1(z)$使用不同的采样矩阵进行上采样，并对上一级粗糙层分解得到的方向子带进行高通滤波，完成频域中更为精确的方向分解。例如，可对滤波器$U_0(z)$和$U_1(z)$分别按采样矩阵$Q=[1,-1;1,1]$进行上采样得到滤波器$U_0(z^Q)$和$U_1(z^Q)$，再对前一级双通道方向分解后得到的子带图像进行滤波，实现四通道方向分解。如图3-5所示，NSDFB可将二维的频域平面分割为多个不同方向的楔形结构，每个楔形块表示图像对应在该方向上的细节特征，由此形成一个由多个双通道NSDFB组成的树形结构。同样地，也可对综合滤波器$V_0(z),V_1(z)$做滤波处理。若NSDFB将NSP变换得到的某一尺度下的带通子带再进行l级方向分解，则可得到$2^l$个与源图像大小相同的方向子带。因此，一幅图像经NSCT执行$k$级分解后，可得到一幅低通子带图像和${\sum }_j^k{2}^{l_j}$ 幅带通方向子带图像，其中，$l_j$为在第$j$尺度下的分解方向数。

3.5 非下采样方向滤波器组分解示意图

此图是理解NSCT方向滤波的重点！图(a)中的经过U1的上下扇形滤波后只剩下对应区域内的信息，在经过y0的二四象限滤波后，剩下的信息区域为(b)中的0区域（扇形和象限正方形取交集）。这就是所谓的支撑空间，前文在谈小波时说到的支撑基。

NSDFB滤波器进行上采样并未增加计算复杂度。假定上采样矩阵为S，二维滤波器为$H(z)$，使用滤波器$H(z^s)$对信号$x[n]$实施滤波，求得卷积输出$y[n]$:
$$y[n]=\sum _{k\in \mathrm{supp}(h)}h[n]x\left[n-S_k\right]$$
上述过程即为α ̀ trous滤波。因此，NSDFB中每个滤波器与方向扇形滤波器有相同的复杂度。同样，NSDFB中每个分解滤波器与第一级中采用的滤波器组的复杂度相同，即NSCT的计算复杂度同NSDFB的复杂度是一致的。NSCT对图像的分解和重构过程中没有进行上采样和下采样操作，这能消除Contourlet变换过程中易产生的频谱混叠现象。相对Contourlet变换，NSCT滤波器的设计方式更加灵活，其特有的设计方式使得NSCT具备Contourlet所没有的平移不变性，因此可以有效减少源图像误配准误差对融合性能的影响。同时图像经NSCT分解后得到的各子图像与源图像具有相同的尺寸大小，容易找到各分解子图像系数间的对应关系，从而有利于图像融合规则的制定。

四、NSCT的matlab程序实现

matlab中可以调用NSCT相关的工具包直接实现其分解功能，因此NSCT的函数便不再复现。关于matlab调用NSCT的方法如下：

我在学的时候发现网上好多的关于matlab实现NSCT都是付费教程，我想在我的博文里免费的教大家实现，基础知识不必为费用所累，希望大家站在这个免费的基础上去创造更多的更有价值的知识。做一名合格的引路人！

1. 下载NSCT的工具包

下载地址

2. 将压缩包解压并放置在你想要调用的程序项目文件序列中。以我的为例，我想创建一个NSCT_Learn项目：先创建一个NSCT_Learn文件夹，将这个压缩包解压放在文件夹里。
   
3. 打开matlab 利用 ”cd 你的项目文件夹的路径“的方法在matlab中打开这个项目文件。
   例如我的文件路径是”F:\Project_Code_any\Matlab_Project\NSCT_Learn“
   我在matlab底下的命令行窗口中输入：

cd F:\Project_Code_any\Matlab_Project\NSCT_Learn

然后创建一个脚本文件重命名为nsct_main，我们将在这个文件夹里实现调用nsct函数。

4. 将nsct包导入的matlab里，实现调用。
   在界面左上角选中NSCT_TOOLBOX_and_TEST IMAGES打开发现，它是由nsct和source images两个文件夹构成，此时选中NSCT_TOOLBOX_and_TEST IMAGES文件夹右键选择将此文件夹添加到路径，选择”选定的文件夹和字文件夹“。原因是该文件夹（工具包）里有我们要用的nsct的实现函数和相应的附属文件，将它导入matlab的过程就是这个添加到路径，因为这个文件夹下的所有子文件都可能是有用且相关的，所以要导入它文件名下的所有。这一过程就像python中的import。
   

5. 在nsct_main文件中写入以下代码：

clear
clc
close all;

image = imread('./lena.jpg');   %使用时替换成你所需要处理的文件的路径
image_gray = rgb2gray(image);
f=im2double(image_gray);

y=nsctdec(f,1,'pkva','9-7');

n = f - y{1};
figure;imshow(f);
figure;imshow(y{1});
figure;imshow(n,[]);
     

阅读上面代码不难发现，nsct的实现主要是靠y=nsctdec(f,1,‘pkva’,‘9-7’);这句中的nsctdec函数实现的。为了更好的理解这个函数怎么用，可以将鼠标光标单击这个函数，然后右键选择打开。

这是我们会发现这个nsctdec函数/文件存在了我们刚才导入的工具包中的nsct文件夹下（该文件夹里集成了很多其它的功能文件，有兴趣的话可以一一探索）。能打开这个说明我们已经把nsct包成功导入了（其实就是文件调用）。这个文件里是nsct的代码实现，有更深入学习原理的小伙伴可以逐行解读代码从而彻底弄清NSCT的工程实现，有助于后续对NSCT进行创新改进。我这里就只讲它怎么用了。
从注释中可以看出nsctdec函数有四个参数：

参数	含义
x	表示待处理的图像
levels	表示分解等级，浅显的理解就是充分分解多少次，前面原理部分提到过。特别注意：因为每次分解会产生一个细节图像，就是被分解/过滤掉的部分，所以此函数的返回值y是一个集合，它包含了前k-1次过滤时产生的图像和过滤后的图像（第k次）（以level=k为例）。因此在上面的程序中表示分解后的图像时采用的时y{1}。
dfilt	选取方向滤波器
pfilt	选取塔式滤波器

ps:滤波器具体都有什么l选择在注释里都有明确的解释，这里就不赘述了。本文示例程序选的是：‘pkva’和’9-7’

6. 程序运行结果展示：
   

五、总结

5.1 NSCT现存不足

由于NSCT中的固定方向滤波器难以处理复杂的方向表示，NSCT也可能产生晕伪影。

5.2 致谢

在本文关于NSCT理论介绍中参考了来自上电子科技大学博士论文：

《基于多分辨率分析的医学图像融合关键技术研究》

具体援引内容就不一一标注了，特此表示感谢！