###Date: 2018.4.30
===================================================
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; %計時結束 實驗結果圖：
參考：https://www.cnblogs.com/molakejin/p/5918976.html