基于小波包变换和高阶统计量的高斯判别准则的红外小目标的分割程序V2.0-MATLAB版

  基本思想：

1. 红外小目标的特点：处于复杂的环境中，目标信号能量很弱，信噪比很低，在一幅图像中主要处于高频段。其可以表示为：

   F(x,y)=f(x,y)+fg(x,y)+fn(x,y)。

   其中，f(x,y)为目标信号，主要是图像的高频部分；fg(x,y)为背景信号，主要是图像的低频部分；fn(x,y)为噪声，一般是在图像的高频部分，符合高斯分布的。

2. 检测的目标：是将背景和噪声去除，提取出小目标。

3. 小波包变换进行目标检测的主要思想:是在小波变换域中提取感兴趣的细节部分。将图像进行多尺度分解，每一次分解都分成四个大小相同的部分，左上角的低频部分，另三个的高频部分；递归的对这四个部分进行分解，直到所设定的分解层数，对分解后的高频的小波包系数和低频的小波包进行出理。

4. 小波包变换的特点：可以完全重构，可以提高高频的分辨率（小波变换只对低频进行递归分解，高频分辨率低，分解后的高频中有目标信息，高斯性噪声，还有部分背景）。

5. 分解后L层后，需要对分解后的小波包系数进行处理，那处理的原则是：将最低频系数置为0（因为它代表背景信息）；在高频域，将高斯性系数剔除，置为0，也就是将高斯性噪声去除；保留非高斯性系数，最后重构，可以得到只含目标的图像。

6. 在低频域处理很简单，只要分解后，将其置为0即可；主要是对高频域中系数处理，这里用到了基于高阶累计量的高斯判别准则，此准则主要用来判断变换后的小波包系数的高斯性。

7. 高阶统计量包括高阶矩和高阶累积量，高阶累积量有一个很重要的理论：高斯随机过程二阶以上的累计量恒等于0，也就是说含有加性高斯噪声的非高斯信号二阶以上的累积量理论上等于信号本身的累计量。

8. 峰度K是一个归一化的四阶统计量，理论上，当信号服从高斯分布时，其峰度K=0，当信号不服从高斯分布时，其峰度不等于0。

9. 理论上判断就是让信号的峰度值K和0相比。如果等于0，就是高斯性噪声，不等于0，就是非高斯性噪声。在实际中，由于误差的存在，信号的峰度K是与一个置信区间进行判断的，

   如果峰度K在这个区间了，则认为信号是服从高斯性分布的，高斯性噪声，否则认为非是高斯性的，

10. 具体的算法流程：

（1），递归地对图像进行小波包变换，

（2），将最低频带的小波包系数置为0，从最后一层自底向上对各高频带的小波包变换系数进行高斯性判断，如果是高斯性系数，就其置为0，如果不是，则保留。

（3），如果相邻的四个高频子带（同一父节点）小波包系数都为0，将这个四个子带合并，这样做的目的是要找到最优的小波包分解结构

（4），最后对于新的小波包系数进行重构，就可以重构目标图像

（5），得到的重构图像灰度值较低，采用自适应阈值的方式进行分割。

11，总结：打个比喻，就好像，用一个放大镜（小波包变换工具）在图像的每一个小区域（分解后的高频和低频）进行查找(查找依据高斯判别准则），将不是目标信息剔除，将是目标的信息保留。

此版本改动：增加了自下而上的将相邻四个（同一父节点）高斯性频带合并。

function infradDepartion
tic;
f=imread('G:\红外测试序列\112.bmp');
if ndims(f)>2
f=rgb2gray(f);
end; 
order=2;
depth=4;
f=double(f);
f=medfilt2(f);
T=wpdec2(f,order,'bior4.4'); //小波包分解
firstIndex=(order^depth-1)/3;          //depo2ind(order,[depth,0]);
lastIndex=((order^(depth+1)-1)/3)-1;   //depo2ind(order,[depth,order^depth])-1;
gaus=[];//用来存放节点是不是高斯性的值（是把1存进去，不是把0放进去）.
isGausIndex=[];//用来存放高斯性系数的节点值（也就是 indice value of nodes).
mgausIndex=[];//用来存放gaus数组中等于1的元素下标.
mindex=1;//记录gaus数组中等于1的元素下标
for i=firstIndex:lastIndex
 cfs=wpcoef(T,i);
 %cfs=cfs*1.5;
 igass=judgeGauss(cfs);
 gaus=[gaus,igass];
 if igass==1
     isGausIndex=[isGausIndex,i];
     mgausIndex=[mgausIndex,mindex];
 end;
 mindex=mindex+1;
 end;
 gaus;
 pargaus=[];//存放父节点的高斯性
 parIndex=[];//存放四个高斯性节点的父节点索引.
 //合并具有同一个父节点的四个高斯性节点.
 for j=1:order:lastIndex-firstIndex+1
     if j<(lastIndex-firstIndex-2)&gaus(j)==gaus(j+1)&gaus(j+1)==gaus(j+2)&gaus(j+2)==gaus(j+3)&gaus(j+3)==1
         parentIndex=(j+firstIndex-2)/4;
         parIndex=[parIndex,parentIndex];
         cfs=wpcoef(T,parentIndex);
         igass=judgeGauss(cfs);
         if igass==1
         pargaus=[pargaus,igass];
         isGausIndex=[isGausIndex,parentIndex];
         end;
         T=nodejoin(T,parentIndex);
     end;
 end;
  if numel(parIndex)~=0
  isGausIndex=deleteGausIndex(parIndex,isGausIndex); //这里要把isGausIndex中的高斯性系数的节点值给删除
  end;
nGausIndex=numel(isGausIndex);
 //抑制低频.
m=firstIndex;
cp=wpcoef(T,m);
T=write(T,'cfs',m,cp);
  //将高斯性系数置为0.    
 for j=1:nGausIndex
  m=isGausIndex(j);
  cp=wpcoef(T,m);
  [h,w]=size(cp);
  y=zeros(h,w);
  T=write(T,'cfs',m,y);
 end;

 f1=wprec2(T);
 means=mean2(f1);
 stds=std2(f1);
 v=means+stds*3;
 [l1,l2]=size(f1);
 for i=1:l1
     for j=1:l2
         if f1(i,j)<v
             f1(i,j)=0;
         else
             f1(i,j)=1;
         end;
     end;
 end;
figure(1);
 subplot(121);
 imagesc(f);
 title('原始红外图像');
 colormap('gray');
 subplot(122);
 imagesc(uint8(f1));
 title('分割后的结果');
toc;

//删除高斯性系数结点
function delGausInd=deleteGausIndex(parentIndex,isGausIndex)
for i=1:numel(parentIndex)
sonIndex=4*parentIndex(i)+1;
  for j=1:numel(isGausIndex)
    if isGausIndex(j)==sonIndex
        isGausIndex(j:j+3)=[];
        break;
    end;
  end;
  delGausInd=isGausIndex;
end;

//判断小波包系数是否是高斯性系数
function isGauss=judgeGauss(wpacketcoef)
L=numel(wpacketcoef);
[row,col]=size(wpacketcoef);
sum1=0;
sum2=0;
Confidence=0.9;//置信度
for i=1:row
    for j=1:col
        temp1=wpacketcoef(i,j)^4;
        temp2=wpacketcoef(i,j)^2;
        sum1=sum1+temp1;
        sum2=sum2+temp2;
    end;
end;
k=L*(sum1/(sum2^2))-3;
if abs(k)<sqrt(24/(L*(1-Confidence)))
    isGauss=1;
else
    isGauss=0;
end;

  实验结果如下：