NSCT:Nonsubsampled Contourlet变换算法以及matlab实现

1、Nonsubsampled Contourlet变换算法介绍：
    对信号的稀疏表示是许多信号处理及应用的基础，2004年Minh N Do、Martin Vetterli提出了一种能够较好表示二维信号的数学工具--Contourlet变换。Contourlet是用金字塔方向滤波器组（PDFB）来将图像分解成不同尺度下的方向子带的。根据PDFB的结构，PDFB是一个拉普拉斯金字塔滤波器Laplacian Pyramid （LP）和一个方向滤波器组的叠加。实验证明，Contourlet变换在图像降噪，纹理，形状的特征提取方面的性能比2-D离散小波变换有了明显的提高。
为了获得平移不变性，本章所用的Nonsubsampled Contourlet变换（NSCT）是基于Nonsubsampled金字塔(NSP)和Nonsubsampled方向滤波器(NSDFB)的一种变换。首先由NSP对输入图像进行塔形分解，分解为高通和低通两个部分，然后由NSDFB将高频子带分解为多个方向子带，低频部分继续进行如上分解。NSCT是一种新型平移不变，多尺度，多方向性的快速变换。
1. Nonsubsampled Pyramid（NSP）：
Nonsubsampled Pyramid（NSP）和Contourlet的Laplacian Pyramid(LP)多尺度分析特性不同。图像通过Nonsubsampled Pyramid（NSP）进行多尺度分解，NSP去除了上采样和下采样，减少了采样在滤波器中的失真，获得了平移不变性。NSP为具有平移不变性滤波结构的NSCT多尺度分析，可以得到与LP分解一样的多尺度分析特性。图2.4(a)处分为3个尺度。
2. Nonsubsampled方向滤波器（NSDFB）
Nonsubsampled方向滤波器(NSDFB)是一个双通道的滤波器，将分布在同方向的奇异点合成NSCT的系数。方向滤波器（DFB）是Bamberger and Smith提出的。其通过一个l层的树状结构的分解，有效的将信号分成了 个子带，其频带分割成为锲形。Nonsubsampled DFB（NSDFB）为非采样，减少了采样在滤波器中的失真，获得了平移不变性。并且每个尺度下的方向子图的的大小都和原图同样大小，Contourlet变换为所有子带之和等于原图。NSCT有更多的细节得以保留，变换系数是冗余的。图2.4（b）为三个尺度下对图像频域的分割图，其中每个尺度的方向子带数目以2倍递增，以在1，2，3尺度下的方向子带数目分别为2，4，8个。

 

NSCT内部变量：

下面列出程序中主要的变量和结构体定义说明：

图像信息结构体:

typedef struct        

{  int high; //-----------------------------------------图像高度

   int width; //---------------------------------------图像宽度

  double max; //--------------------------------------最大像素值

  double min; //---------------------------------------最小像素值

  double *data; //--------------------------------------像素数据区域

}ImgData;  

 

NSCT变换空间结构体：（用于存放NSCT变换滤波器，尺度、自带空间等基本信息）

typedef struct

{ unsigned int  *nlevels;//子带方向数  必须为偶数  sublevels[0]为最高级方向数

sublevels[4]为最低级方向数   为负时为无效

  unsigned int  clevels; //------------------尺度层数

  ImgData ***high;//------------------------高通子带区域

  ImgData low;//-----------------------------低通子带区域

  ImgData H1;//------------------------------高通滤波器组

  ImgData H2; //------------------------------高通滤波器组

  ImgData G1; //------------------------------低通滤波器组

  ImgData G2; //------------------------------低通滤波器组

  ImgData ***filters_dec; //-----------------分解方向滤波器组

  ImgData ***filters_rec; //-----------------重构方向滤波器组

  BOOL sign;

}nsct_t;

尺度与分解层数滤波器：

const int  dlevels[3]={3,3,3};//--------------尺度为三，个尺度方向数为

const int   Qunx[4]={0,0,0,0};//--------------------------Q梅花采样矩阵

BMGImageStruct SourceBmp;//---------------------------位图操作结构体

int pfilter_type,dfilter_type; //滤波器类型选择,  两者都默认为,'maxflat'&'dmaxflat7'滤波器

int pralevels; //------------图像金字塔层数（不能大于4，尺度在4以内足以满足一般图像处理需求）

nsct_t  nsct; //nsct-------数据计算区域

int shift[2]; //---------------------延迟分量
2.3.6 NSCT内部函数：

下面列出程序中主要的函数定义说明：

void ImgDataM_Var_3m(ImgData *x,float *mean,float *var,float *mom);

//----------------------------------计算nsct子带系数的均值和三阶中心矩

void ImgDataIFFT_2D(complex<double> * pCFData, complex<double> * pCTData, int nWidth, int nHeight);

//----------------------------------2D-FFT反变换

void ImgDataFFT_2D(complex<double> * pCTData, int nWidth, int nHeight, complex<double> * pCFData);

//----------------------------------2D-FFT反变换

void ImgDataFConv2(ImgData *out,ImgData *img,ImgData *mask);

//----------------------------------频域卷积函数

void Nsct_Extend2(ImgData *in,int ru,int rd,int cl,int cr,ImgData *out,const int flag=1);

//----------------------------------nsct镜像扩展函数

void Nsct_Zconv2(ImgData *in,ImgData *out,ImgData *filter,const int *mup);

//----------------------------------时域滤波器卷积函数（快速算法）

void Nsct_Efilter2(ImgData *in,ImgData *out,ImgData *filter);

//----------------------------------边缘扩展的2D滤波

void Nsct_Nssfbdec(ImgData * input,ImgData * out1,ImgData * out2,ImgData * filter1,ImgData * filter2,const int * Qx=Qunx);

//----------------------------------2通道非采样子带分解

void Nsct_Mirror2(ImgData *tomirror);

//----------------------------------镜像变换

void ImgDataConv2(ImgData *out,ImgData *img,ImgData *mask,ImgData *imgunext);

//----------------------------------相关卷积

void Nsct_Upsample2df(ImgData *h0,int lev,ImgData *ho);

//----------------------------------频域2D上采样

void Nsct_Symext(ImgData *imglow,ImgData *upedHi,ImgData *pp );

//----------------------------------2D平衡矩阵生成

void Nsct_Atrousc(ImgData * out,ImgData * signal,ImgData *filter,int (*matrix)[2]);

//----------------------------------包含抽取的滤波器卷积

void Nsct_Nsdfbdec(ImgData *imghigh,int currentlevel);

//----------------------------------非抽取方向子带分解

void Nsct_iTtransform();

//----------------------------------反变换主函数

void Nsct_New(HBITMAP hBitmap,const int *Nlevels=dlevels,const unsigned int Level=3,const int dfilter=1,const int pfilter=1);

//-----------------------------用来初始化原始目标图像，申请nsct计算空间和存储空间

void Nsct_Transform();

//-----------------------------变换主函数

void Nsct_NsfbDec(ImgData *imglow,ImgData *imghigh,int i);

//-----------------------------nsct图像金字塔分解

void NsctDelete();

//-----------------------------nsct析构

void ImgDataFreeMemery(ImgData  *image);

//-----------------------------内存空间释放

BOOL ImgDataAllocMemery(ImgData * image);

//-----------------------------内存空间申请

一个示例程序，和显示的样图：

2、NSCT用于图像融合的应用（matlab实现）

    %% 实现红外与可见光图像融合NSTC
    %计时开始%
    tic;%计时  
    path(path,'nsct_toolbox');%加载matlab自带的nsct工具箱
    path(path,'transform');%加载RGB和IHS转换用的函数的文件夹
     
    %% 1、输入图像和初始化参数
    kk=1.5;%参数
    [input_image_TV_int_GIF, input_image_TV_int_MAP] = imread('test_kejianguang.gif' );%输入gif格式的可见光图像%
    input_image_TV_int_RGB = uint8(256*ind2rgb(input_image_TV_int_GIF, input_image_TV_int_MAP));%将gif格式数据转换为RGB格式数据
     
    [input_image_IR_int_GIF, input_image_IR_int_MAP] = imread('test_hongwai.gif' ); %输入gif格式的红外图像%
    input_image_IR_int_RGB = uint8(256*ind2rgb(input_image_IR_int_GIF, input_image_IR_int_MAP));%将gif格式数据转换为RGB格式数据%
     
    [Ny,Nx] = size(input_image_TV_int_RGB);%求取输入图像的宽和高
     
    %将可见光图像的格式转换为IHS格式
    input_image_TV_int_IHS = rgb2ihs( input_image_TV_int_RGB );%IHS
    input_image_TV_IHS = double(input_image_TV_int_IHS);  %将输入图像的数据类型转换为双精度数据类型，便于处理
    input_image_TV_I =  256*input_image_TV_IHS(:,:,3);  %将输入图像的I分量提取出来以进行NSCT分解
     
    %将红外图像的格式转换为gray格式，灰度格式
    input_image_IR_int = rgb2gray( input_image_IR_int_RGB );%灰度化
    input_image_IR= double(input_image_IR_int); %将输入图像的数据类型转换为双精度数据类型
     
    input_image_TV_S =  input_image_TV_IHS(:,:,2);
    input_image_TV_H =  input_image_TV_IHS(:,:,1);
     
    figure;
    subplot(1,2,1);
    imshow(uint8(input_image_TV_I));title('可见光图像亮度分量I') %显示输入可见光图像的I分量
    subplot(1,2,2);
    imshow(uint8(input_image_IR));title('红外图像灰度化') %显示输入红外图像的灰度图像
     
    %金字塔表示参数%
    Nsc = ceil(log2(min(Ny,Nx)) - 7);    %分解尺度的数量 (自适应于图像尺寸)%
    Nor = 8;                             %每级分解的方向数%
     
    %% 2、NSCT子带分解
    %初始化NSCT子带分解参数%
    pfilter = 'maxflat' ;  %金字塔滤波器
    dfilter = 'dmaxflat7' ; %方向滤波器
     
    nlevels=zeros(1,Nsc+1);
    for i=1:Nsc+1
        nlevels(i)=log2(Nor); %初始化分解尺度%
    end
     
    %NSCT分解
    coeffs_TV_int = nsctdec( input_image_TV_I, nlevels, dfilter, pfilter ); %分解可见光图像
    coeffs_IR_int = nsctdec( input_image_IR, nlevels, dfilter, pfilter ); %分解红外图像
     
    %% 3、图像融合
    %3.1 低频分量采用自适应模糊逻辑算法
    u=mean2(coeffs_IR_int{1});%均值
    s=std2(coeffs_IR_int{1});%方差
    [lx,ly]=size(coeffs_IR_int{1});
    for x=1:lx
        for y=1:ly
            Nb=exp(-(coeffs_IR_int{1}(x,y)-u)^2/(kk*(2*s)^2));
            Nt=1-Nb;
            coeffs_rec{1}(x,y)=Nt*coeffs_IR_int{1}(x,y)+Nb*coeffs_TV_int{1}(x,y);
        end
    end
    %3.2 高频分量采用模值取大的方法
    for i=2:Nsc+2
        for j=1:Nor
            band=abs(coeffs_TV_int{i}{j})-abs(coeffs_IR_int{i}{j});
            band_coffes1=double(im2bw(band,0));
            band_coffes2=ones(size(band))-band_coffes1;
            coeffs_rec{i}{j} = coeffs_TV_int{i}{j}.*band_coffes1+coeffs_IR_int{i}{j}.*band_coffes2;  
        end
    end
     
    %% 4、NSCT子带重构
    out_image_end_I = nsctrec( coeffs_rec, dfilter, pfilter ) ;%重构灰度图像%  
    figure;
    imshow(uint8(out_image_end_I));title('融合图像亮度分量');
    out_image_end=zeros(Ny,Nx,3);   %初始化彩色图像输出矩阵
    out_image_end=cat(3,input_image_TV_H,input_image_TV_S,out_image_end_I/256);                      
    out_image_end_RGB=256*ihs2rgb(double(out_image_end)); %将输出彩色图像转换为RGB格式              
     
    figure;
    imshow(uint8(out_image_end_RGB));%显示融合后彩色图像%
    title('融合后的彩色图像')
    imwrite(uint8(out_image_end_RGB),'fusion_image.bmp','bmp');
    toc; %计时结束
---------------------
作者：SoaringLee_fighting
来源：CSDN
原文：https://blog.csdn.net/soaringlee_fighting/article/details/80150711
版权声明：本文为博主原创文章，转载请附上博文链接！