（原创）极化SAR精致Lee滤波matlab…

      关于精致Lee滤波的原创matlab代码，只限仿真使用。是自己在初学的时候自己写的代码，后来才知道有官方的代码，只需运行PolSARpro软件即可。不过做为初学，还是有一定纪念意义，贴上来分享交流一下。
     
      仿真论文是：
         L. Jong-Sen, M. R. Grunes, and G. De Grandi, "Polarimetric SAR speckle filtering and its  implication for classification," Geoscience and Remote Sensing, IEEE Transactions on, vol. 37, pp. 2363-2373, 1999.
       
      关于数据 SanFrancisco_C.mat， SanFrancisco_T.mat的获取，可以先访问欧空局的网站 http://earth.eo.esa.int/polsarpro/datasets.html，下载到原始数据，然后解压成下面图片所示形式。之后用图片下面蓝色代码读出即可。



0.下面蓝色代码读出T矩阵和C矩阵

clc;
clear all;
addpath .\T3

file_idt11 = fopen('T3\T11.bin','rb');
file_idt22 = fopen('T3\T22.bin','rb');
file_idt33 = fopen('T3\T33.bin','rb');
file_idt12r = fopen('T3\T12_real.bin','rb');
file_idt12i = fopen('T3\T12_imag.bin','rb');
file_idt13r = fopen('T3\T13_real.bin','rb');
file_idt13i = fopen('T3\T13_imag.bin','rb');
file_idt23r = fopen('T3\T23_real.bin','rb');
file_idt23i = fopen('T3\T23_imag.bin','rb');
h=900;
w=1024;
col=w*h;
t11=zeros(h,w);
for i=1:h
    t11(i,:)=fread(file_idt11,w,'float32');%tmp=fread(file_idt11,12,'int8');
    t22(i,:)=fread(file_idt22,w,'float32');%tmp=fread(file_idt22,12,'int8');
    t33(i,:)=fread(file_idt33,w,'float32');%tmp=fread(file_idt33,12,'int8');
    t12i(i,:)=fread(file_idt12i,w,'float32');%tmp=fread(file_idt12i,12,'int8');
    t12r(i,:)=fread(file_idt12r,w,'float32');%tmp=fread(file_idt12r,12,'int8');
    t13i(i,:)=fread(file_idt13i,w,'float32');%tmp=fread(file_idt13i,12,'int8');
    t13r(i,:)=fread(file_idt13r,w,'float32');%tmp=fread(file_idt13r,12,'int8');
    t23i(i,:)=fread(file_idt23i,w,'float32');%tmp=fread(file_idt23i,12,'int8');
    t23r(i,:)=fread(file_idt23r,w,'float32');%tmp=fread(file_idt23r,12,'int8');
end
t12=complex(t12r,t12i);
t13=complex(t13r,t13i);
t23=complex(t23r,t23i);

save SanFrancisco_T t11 t22 t33 t12 t23 t13





1.主函数：
clear all,close all

load SanFrancisco_C.mat
load SanFrancisco_T.mat

span = c11 + c22 + c33;

m_orig = 900;
n_orig = 1024;

span = span(1:m_orig,1:n_orig);
t11 = t11(1:m_orig,1:n_orig);
t12 = t12(1:m_orig,1:n_orig);
t13 = t13(1:m_orig,1:n_orig);
t22 = t22(1:m_orig,1:n_orig);
t23 = t23(1:m_orig,1:n_orig);
t33 = t33(1:m_orig,1:n_orig);

span = double(span);
[msize,nsize]  = size(span);
A = flipdim(span,1);
A1 = [A;span;A];
A2 = flipdim(A1,2);
B = [A2,A1,A2];
expan_span = B;

deltau = 0.5;
blockSize = 7;
sub_blockNum = 9;
sub_blockWide = sub_blockNum/3;

 t11_filtered = zeros(msize,nsize);
 t12_filtered = zeros(msize,nsize);
 t13_filtered = zeros(msize,nsize);
 t22_filtered = zeros(msize,nsize);
 t23_filtered = zeros(msize,nsize);
 t33_filtered = zeros(msize,nsize);

%将T矩阵每个分量也相应扩展
[expan_t11,expan_t12,expan_t13,expan_t22,expan_t23,expan_t33] = expan_T(t11,t12,t13,t22,t23,t33);


for i = 1:msize
    for j = 1:nsize
        disp([ '  i=' num2str(i) '   j=' num2str(j)])
        original_T = [t11(i,j),       t12(i,j),       t13(i,j);
                      conj(t12(i,j)), t22(i,j),       t23(i,j);
                      conj(t13(i,j)), conj(t23(i,j)), t33(i,j)];
        %取出相应的小块
        [block,block_t11_resh,block_t12_resh,block_t13_resh,block_t22_resh,...
                block_t23_resh,block_t33_resh] ...
                = choose_block(expan_span,expan_t11,...
                expan_t12,expan_t13,expan_t22,expan_t23,expan_t33,i,j,msize,nsize,blockSize) ;        
                
       
        mean_sub_block = zeros(3,3);
       
        for p = 1:2:5%求均值矩阵
            for q = 1:2:5
                sub_block = block(p:p+sub_blockWide-1,q:q+sub_blockWide-1);
                x = (p+1)/2; y = (q+1)/2;
                mean_sub_block(x,y) = mean(sub_block(:));
            end
        end
        F1 = zeros(1,4);
        for n = 1:4
            grad_temp = choose_gradtemp(n);%用四个模板进行滤波
            F = imfilter(double(mean_sub_block),grad_temp,'symmetric');
            F1(1,n) = sum(abs(F(:)));
        end
       
        direct = choose_direction(F1,mean_sub_block);%每个7*7小块的最终方向
        [local_block,local_block_index] = choose_local(direct,block,blockSize);%选出的每个小块用于计算局部均值的点
        b = weightb(local_block,deltau);
        %计算每个小块中局域点的T矩阵的均值
        T_aver = compu_T_aver(local_block_index,block_t11_resh,block_t12_resh,block_t13_resh,...
                              block_t22_resh,block_t23_resh,block_t33_resh);
        T_filtered = T_aver + b*(original_T - T_aver);
        %滤波后的T矩阵
        t11_filtered(i,j)  = T_filtered(1,1);
        t12_filtered(i,j)  = T_filtered(1,2);
        t13_filtered(i,j)  = T_filtered(1,3);
        t22_filtered(i,j)  = T_filtered(2,2);
        t23_filtered(i,j)  = T_filtered(2,3);
        t33_filtered(i,j)  = T_filtered(3,3);
    end
end
save T_filtered t11_filtered t12_filtered t13_filtered t22_filtered t23_filtered t33_filtered
%画出RGB伪彩色图像
T_RGB = zeros(msize,nsize,3);
T_RGB(:,:,1) = sqrt(t22);
T_RGB(:,:,2) = sqrt(t33);
T_RGB(:,:,3) = sqrt(t11);
figure(1),imshow(T_RGB),title('original T RGB')
saveas(1,'original_T_RGB','fig')
filtered_T_RGB = zeros(msize,nsize,3);
filtered_T_RGB(:,:,1) = sqrt(t22_filtered);
filtered_T_RGB(:,:,2) = sqrt(t33_filtered);
filtered_T_RGB(:,:,3) = sqrt(t11_filtered);
figure(2),imshow(filtered_T_RGB),title('refined Lee filtered T RGB')
saveas(2,'refined Lee filtered_T_RGB','fig')


2.将蓝色部分保存为choose_block.m文件
 function [block,block_t11_resh,block_t12_resh,block_t13_resh,block_t22_resh,...
                            block_t23_resh,block_t33_resh]...
                            = choose_block(expan_span,expan_t11,expan_t12,expan_t13,...
                             expan_t22,expan_t23,expan_t33,i,j,msize,nsize,blockSize)        
       
        block = expan_span( i+msize - (blockSize-1)/2  : i+msize + (blockSize-1)/2 ,...
                            j+nsize - (blockSize-1)/2  : j+nsize + (blockSize-1)/2);
        block_t11 = expan_t11( i+msize - (blockSize-1)/2  : i+msize + (blockSize-1)/2 ,...
                            j+nsize - (blockSize-1)/2  : j+nsize + (blockSize-1)/2);
        block_t12 = expan_t12( i+msize - (blockSize-1)/2  : i+msize + (blockSize-1)/2 ,...
                            j+nsize - (blockSize-1)/2  : j+nsize + (blockSize-1)/2);
        block_t13 = expan_t13( i+msize - (blockSize-1)/2  : i+msize + (blockSize-1)/2 ,...
                            j+nsize - (blockSize-1)/2  : j+nsize + (blockSize-1)/2);
        block_t22 = expan_t22( i+msize - (blockSize-1)/2  : i+msize + (blockSize-1)/2 ,...
                            j+nsize - (blockSize-1)/2  : j+nsize + (blockSize-1)/2);
        block_t23 = expan_t23( i+msize - (blockSize-1)/2  : i+msize + (blockSize-1)/2 ,...
                            j+nsize - (blockSize-1)/2  : j+nsize + (blockSize-1)/2);               
        block_t33 = expan_t33( i+msize - (blockSize-1)/2  : i+msize + (blockSize-1)/2 ,...
                            j+nsize - (blockSize-1)/2  : j+nsize + (blockSize-1)/2);
       
        block_t11_resh = reshape(block_t11,1,blockSize*blockSize);
        block_t12_resh = reshape(block_t12,1,blockSize*blockSize);
        block_t13_resh = reshape(block_t13,1,blockSize*blockSize);
        block_t22_resh = reshape(block_t22,1,blockSize*blockSize);
        block_t23_resh = reshape(block_t23,1,blockSize*blockSize);
        block_t33_resh = reshape(block_t33,1,blockSize*blockSize);


3.将棕色部分保存为choose_direction.m文件
function direct = choose_direction(F1,mean_sub_block)
        [Y,ind] = max(F1);
        if ind==1
           diff1 = abs(mean_sub_block(2,3)-mean_sub_block(2,2));
           diff2 = abs(mean_sub_block(2,1)-mean_sub_block(2,3));
           if diff1>=diff2
               direct = 0;
           elseif diff1               direct = 4;
           end
        end
        if ind==2
           diff1 = abs(mean_sub_block(1,3)-mean_sub_block(2,2));
           diff2 = abs(mean_sub_block(3,1)-mean_sub_block(2,2));
           if diff1>=diff2
               direct = 1;
           elseif diff1               direct = 5;
           end
        end
        if ind==3
           diff1 = abs(mean_sub_block(1,2)-mean_sub_block(2,2));
           diff2 = abs(mean_sub_block(3,2)-mean_sub_block(2,2));
           if diff1>=diff2
               direct = 2;
           elseif diff1               direct = 6;
           end
        end
        if ind==4
           diff1 = abs(mean_sub_block(1,1)-mean_sub_block(2,2));
           diff2 = abs(mean_sub_block(3,3)-mean_sub_block(2,2));
           if diff1>=diff2
               direct = 3;
           elseif diff1               direct = 7;
           end
        end

4.将下面粉紫色部分保存为 choose_gradtemp.m文件
function grad_temp = choose_gradtemp(n)
if n==1
 grad_temp = [-1  0  1;
              -1  0  1;
              -1  0  1];
 elseif n==2
    grad_temp = [ 0  1  1;
                 -1  0  1;
                 -1 -1  0];

elseif n==3
    grad_temp = [ 1   1  1;
                  0   0  0;
                 -1 -1 -1];
 
elseif n==4
    grad_temp = [1  1  0;
                 1  0  -1;
                 0 -1  -1];
 
end


5.将下面绿色部分保存为 choose_local.m文件
function [local_block,local_block_index] = choose_local(direct,w,blockSize)
       
        w_index = [1:blockSize*blockSize];
        w_index1 = reshape(w_index,blockSize,blockSize);
        local_block = [];%每个7*7小块所取出的值
        local_block_index = []; %每个7*7小块所取出值的下标
        if direct==0
            local_block = w(:,4:7);
            local_block_index = w_index1(:,4:7);
        elseif direct==4
            local_block = w(:,1:4);
            local_block_index = w_index1(:,1:4);
        elseif direct==1;
            for i = 1:blockSize
              local_block = [ local_block , w(i,i:end) ];
              local_block_index = [ local_block_index , w_index1(i,i:end) ];
            end
        elseif direct==5
            for i = 1:blockSize
              local_block = [ local_block , w(i,1:i) ];
              local_block_index = [ local_block_index , w_index1(i,1:i) ];
            end
        elseif direct==2
            local_block = w(1:4,:);
            local_block_index = w_index1(1:4,:);
        elseif direct==6
            local_block = w(4:7,:);
            local_block_index = w_index1(4:7,:);
        elseif direct==3
            for i = 1:blockSize
                local_block = [ local_block , w( i,1:end-(i-1) ) ];
                local_block_index = [ local_block_index , w_index1( i,1:end-(i-1) ) ];
            end
        elseif direct==7
            for i = 1:blockSize
                local_block = [ local_block, w( i, (blockSize+1-i:end) ) ];
                local_block_index = [ local_block_index, w_index1( i, (blockSize+1-i:end) ) ];
            end
        end


6.将下面（come on， 貌似蓝色） 部分保存为compu_T_aver .m文件
function T_aver = compu_T_aver(local_block_index,block_t11_resh,block_t12_resh,...
                               block_t13_resh,block_t22_resh,block_t23_resh,block_t33_resh)
        perNum = length(local_block_index(:));
        block_t11_sum = sum( block_t11_resh(1,local_block_index(:)));
        block_t12_sum = sum( block_t12_resh(1,local_block_index(:)));
        block_t13_sum = sum( block_t13_resh(1,local_block_index(:)));
        block_t22_sum = sum( block_t22_resh(1,local_block_index(:)));
        block_t23_sum = sum( block_t23_resh(1,local_block_index(:)));
        block_t33_sum = sum( block_t33_resh(1,local_block_index(:)));
       
        T_aver =  [block_t11_sum        block_t12_sum         block_t13_sum;
                   conj(block_t12_sum)  block_t22_sum         block_t23_sum;
                   conj(block_t13_sum)  conj(block_t23_sum)   block_t33_sum];
        T_aver = T_aver./perNum;



7.将下面菊黄色部分保存为expan_T.m文件
function [expan_t11,expan_t12,expan_t13,expan_t22,expan_t23,expan_t33] = expan_T(t11,t12,t13,t22,t23,t33)
A = flipdim(t11,1);
A1 = [A;t11;A];
A2 = flipdim(A1,2);
B = [A2,A1,A2];
expan_t11 = B;

A = flipdim(t12,1);
A1 = [A;t12;A];
A2 = flipdim(A1,2);
B = [A2,A1,A2];
expan_t12 = B;

A = flipdim(t13,1);
A1 = [A;t13;A];
A2 = flipdim(A1,2);
B = [A2,A1,A2];
expan_t13 = B;

A = flipdim(t22,1);
A1 = [A;t22;A];
A2 = flipdim(A1,2);
B = [A2,A1,A2];
expan_t22 = B;

A = flipdim(t23,1);
A1 = [A;t23;A];
A2 = flipdim(A1,2);
B = [A2,A1,A2];
expan_t23 = B;

A = flipdim(t33,1);
A1 = [A;t33;A];
A2 = flipdim(A1,2);
B = [A2,A1,A2];
expan_t33 = B;


8.将下面紫色部分保存为weightb.m文件
function b = weightb(local_block,deltau) 
            tt = local_block;
            mean_tt = mean(tt(:));
            var_x = (var(tt(:)) - (mean_tt.^2)*(deltau.^2))./(1+deltau.^2);
            b = var_x/var(tt(:));


9.结果