% IEEE Image Processing, vol.9, No.6, 2000. 

%--------- 濾波器對圖像作增強提高視覺質量－－－－－－－

clear all;
% C=double(imread('t1.bmp'));%cameraman.tifSAR01.bmpBARB.BMP;5-51fft6649-060-1.tifmed01.jpg
C=double(imread('Car.JPG'));
h=[1,1,1;1,1,1;1,1,1]./9;
B=C(:,:,1);
%B=imfilter(B,h,'replicate');
S=size(B);H=S(1,1);W=S(1,2);

%%%%%%%%%%%%%%%%% CDF9/7
t=0.75;%0.730174;
h0 = (8*t*t*t-6*t*t+3*t)/(1+2*t)*(1/32);
h1 = (-16*t*t*t+20*t*t-12*t+3)/(1+2*t)*(1/32);
h2 = (2*t-3)/(1+2*t)*(1/8.);
h3 = (16*t*t*t-20*t*t+28*t+5)/(1+2*t)*(1/32.);
h4 = (-8*t*t*t+6*t*t+5*t+20)/(1+2*t)*(1/16.);
hL=[h0,h1,h2,h3,h4,h3,h2,h1,h0];

g0=-t/16;
g1=(1-2*t)/16;
g2=(t+4)/16;
g3=(3+2*t)/8;
hH=[-g0,g1,-g2,g3,-g2,g1,-g0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%  對淹沒在噪聲中的SAR圖像首先用Wiener濾波法自適應去噪
%Idenoise= medfilt2(B(:,:,1),[3 3]);%中值濾波
%Idenoise = wiener2(B(:,:,1),[5 5]);%自適應濾波
%Idenoise=double(Idenoise);
Idenoise=B;%imresize(B,[512,510]);

% 用多分辨分解法提取高頻邊緣BI

hL=[1/16,1/4,3/8,1/4,1/16];
for i=1:length(hL)
    for j=1:length(hL)
        hL2D(i,j)=hL(i)*hL(j);
    end
end
L0=imfilter(Idenoise,hL2D,'replicate');%低頻近似圖像L0
L0=double(L0);
L0max=max(max(L0));
c=0.3;s=2;
T1 = 0;
T2 = (1-c)*L0max;
k1=0.6; k2=s;

BI = Idenoise-L0;% 第一層高頻邊緣BI

 % 對BI作非線性插值后得到newBI
newBI=BI;

for i=1:H
    for j=1:W
        if abs(BI(i,j))<=T1
            newBI(i,j)=k1*BI(i,j);
        else if (abs(BI(i,j))>T1)&(abs(BI(i,j))<=T2)
                newBI(i,j)=sign(BI(i,j))*(k2*abs(BI(i,j))+T1*(k1-k2));
            end
        end
    end
end
%newBI=s*newBI;

% 對newBI作高通濾波得到增強后的高頻邊緣BI
hH=[-1/16,1/4,-3/8,1/4,-1/16];%[-1,2,-1]/2;
for i=1:length(hH)
    for j=1:length(hH)
        hH2D(i,j)=hH(i)*hH(j);
    end
end

BI=imfilter(newBI,hH2D,'replicate');%修正后的邊緣
%BI=newBI-double(imfilter(newBI,hL2D,'replicate'));

BI=double(BI);
Irecover=newBI+L0;

for i=1:H
    for j=1:W
        if Irecover(i,j)<0
            Irecover(i,j)=0;
        else if Irecover(i,j)>255
                Irecover(i,j)=255;
            end
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
EP1D=[-0.15,0.25,0.7,0.25,-0.15];
%EP1D=[-0.05,0.15,-0.25,0.3,0.9,0.3,-0.25,0.15,-0.05]./1.2
for i=1:length(EP1D)
    for j=1:length(EP1D)
        EP2D(i,j)=EP1D(i)*EP1D(j);
    end
end
Irepro = imfilter(Irecover,EP2D,'replicate');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%figure,imshow((100*BI)./255,[]);%imfilter(B,hH2D,'replicate')

figure,imshow(uint8(B),[]);
title('原圖');
%subplot(122);
subplot(121);imshow(uint8(Irecover),[]);
title(['高頻非線性增強后的圖像']);

subplot(122);imshow(uint8(Irepro),[]);
title(['非線性增強后再補償的圖像']);