contourenhance_main.m

An efficient way to implement image enhancement with contourlet transform. Welcome

function contourenhance_main
clc
clear all

in=imread('212.bmp');
im=rgb2gray(in);
figure,imshow(in);title('原始图像');axis on;
im = double(im) / 256;
pfilt='9-7';                % choose LP decomposition filter
dfilt='pkva';               % choose DFB decomposition filter
rho = 3;                    % noise level
nlevs = [3,4,4,5,5];       %   nlevs:  vector of numbers of directional filter bank decomposition levels 
                            %           at each pyramidal level (from coarse to fine scale).
%assure the values of J and l(j)
[nlrow,nlcol] = size(nlevs);
k1=nlcol; % k1 record the level of LP
for i=1:nlcol
    if nlevs(i) >= 2
       J = nlcol-i+1;
    break
    else i=i+1;
         J=0;
    end
end
if (J) 
   for j=J:-1:1
       l(j) = nlevs(i);
       i=i+1;
   end
end
% Generate noisy image. 
sig = std(im(:));
sigma = sig / rho;
sigmap=sigma^2;
nim = im + sigma * randn(size(im));
figure,imshow(nim);title('噪声图像');axis on;
%sigmaj = (4/3) * sigma *std(sqrt(nvar));
% Contourlet transform

     y = pdfbdec(nim, pfilt, dfilt, nlevs);
     wl = 0 ; %count the number of zero in nlevs,that is the level of wavelet decomposition
     for i=1:nlcol
         if nlevs(i)==0
         wl = wl+1;
         end
     end

        
    
% if wl=1,y{1}is the lowpass image ,and y{2}is the  horizontal,vertical and diagonally subbands                                     

%cope with the contourlet coefficients
cs = 5; % cs is a parameter ranging from 1 to 5
ap=0.8;
for j=1:J
     p=J+2-j+wl;
     for k=1:2^(l(j)-1)
          sigmaj=sigmap/std(y{p}{k}(:));       % essig is a function to estimate the noise variance of this subband
         [m,n]=size(y{p}{k}); % the size of directional subbands of scale 2^j,direction k
         for f=1:m 
            for g=1:n
                i=0;
                z=0;
                for nu=1:2^(l(j)-1)
                     i=i+1;
                     z(i)=y{p}{nu}(f,g);
                end
                   px=max(abs(z));   % the max coefficients of this subbands
                   pm=mean(z);
                  
                if  (px>=cs*sigmaj)
                           y{p}{k}(f,g) =max(((cs*sigmaj)/abs(y{p}{k}(f,g)))^ap,1)*y{p}{k}(f,g);
                elseif  (px>=2*&px<cs*sigmaj) 
                           y{p}{k}(f,g) =max(((cs*sigmaj)/abs(y{p}{k}(f,g)))^ap,1)*y{p}{k}(f,g);
                       elseif(px)
                   end
            end
         end
     end
     for k=2^(l(j)-1)+1:2^l(j)
         [m,n]=size(y{p}{k});
         sigmaj=sigmap/std(y{p}{k}(:));
          for f=1:m   
             for g=1:n
                 i=0;
                 for nu=2^(l(j)-1)+1:2^l(j)
                     i=i+1;
                     z(i)=y{p}{nu}(f,g);
                 end
                     px=max(abs(z));   % the max coefficients of this subbands
                     pm=mean(z);
                     if (px>=cs*sigmaj)
                            y{p}{k}(f,g)=max(((cs*sigmaj)/abs(y{p}{k}(f,g)))^ap,1)*y{p}{k}(f,g);
                     elseif (px<cs*sigmaj) 
                        y{p}{k}(f,g)=0 ;
                     end
             end
         end
    end
end
imrec = pdfbrec( y, pfilt, dfilt ) ;% use the modified cofficients to reconstruction
%psnr=PSNR(imrec, im);
g=1;
newim=adjgamma(imrec,g);
figure,imshow(newim);title(sprintf('增强后的图像(SNR = %.2f dB)', ...
              SNR(im, imrec)));axis on;

%figure,imshow(imdiff);title('差图像');axis on;