amer.m

amer s noise level estimation implementation on MATLAB

%===========================================================================
%  Amer NLE fuction (by europium 2008)
% --------------------------------------------------------------------------
%  Output 
%     est_std  : estimated noise standard deviation
%     blknum  : total number of blocks
%     homobk :  # of homogeneous blocks
% 
%  Input
%     id                  : image id
%     std                : given noise standard deviation
%     w_size          : operation block size
%     th                  : threshold for measuring homogeneous blocks
%     Estitmation    : selection of noise estimation methods 
%     Adaptive        :  selection of adaptive NLE methods
%     showimg        :  display homogeneous blocks in the given noisy image
%===========================================================================

function [est_std, blknum, homoblk] = amer(id, std, w_size, th, Estimation, Adaptive, showimg)

%===========================================================================
% image read
%===========================================================================
s1 = 'D:\KODAK_PNG\noise\';
s2 = sprintf('img%02d', id);
s3 = '\';
s4 = sprintf('noise_std%02d',std);
s5 = '.png';    

originalfile = strcat(s1,s2,s3, s4,s5);
a = imread(originalfile,'png');
noise = (double(uint8(a)));
  
%===========================================================================
% Masking operation for finding homogeneous blocks
%===========================================================================
      
Cnt = 0;
blknum = 0;    

img_width = size(noise,1);
img_height = size(noise,2);
center = median(1:w_size);
center2 = median(1:w_size^2);
acc_dis = 0;

for i=1:w_size:img_width-w_size %x
    for j=1:w_size:img_height-w_size %y

        blknum = blknum + 1;
        current_pixel = w_size*noise(i+center-1,j+center-1,1);
    
        % mask1  (vertical)
        x_loc = i:i+w_size-1;
        y_loc = j+center-1;
        m1 = abs(current_pixel - sum(noise(x_loc,y_loc,1)));
        % mask2  (horizontal)
        x_loc = i+center-1;
        y_loc = j:j+w_size-1;
        m2 = abs(current_pixel - sum(noise(x_loc,y_loc,1)));    
   
        % mask3 (diagonal)
        x_loc = i:i+w_size-1;
        y_loc = j:j+w_size-1;
        m3 = abs(current_pixel - sum(diag(noise(x_loc,y_loc,1))));        
    
        % mask4 (revers diagonal)
        x_loc = i:i+w_size-1;
        y_loc = j:j+w_size-1;
        m4 = abs(current_pixel - sum(diag(fliplr((noise(x_loc,y_loc,1))))));      
  
        homo = m1+m2+m3+m4; % homogeneity measurement
    
       if homo <= th
        
            if showimg == 1
                noise(x_loc,y_loc,2) = 0; % To show homogeneous blocks in the noisy image
            end           
           
            Cnt = Cnt+1;
            homo_index(Cnt) = homo;
            x_loc = i:i+w_size-1;
            y_loc = j:j+w_size-1;
         
            if Estimation == 3 % Estimation = 2 : Histogram Approximation (68% -> sigma)
                distance = abs(noise(x_loc,y_loc,1)-noise(i+center-1,j+center-1,1));
                dis = reshape(distance, 1, w_size^2);
                adis = [dis(1:center2-1),dis(center2+1:w_size^2)];
                acc_dis = [acc_dis, adis];      
            else % Estimation = 1 & 2 : Averaging 
                h_std(Cnt) = std2(noise(x_loc,y_loc,1));
                
            end
        end
    end
end

homoblk = Cnt;

%===========================================================================
% Noise level estimation
%---------------------------------------------------------------------------
% Estimation = 1 : Averaging local standard deviations
% Estimation = 2 : Adaptive averaging local standard deviations(mean of 'local std - ref_std')
% Estimation = 3 : Histogram Approximation (68% -> sigma)
%===========================================================================

if Estimation == 1 % Estimation = 1 : Averaging local standard deviations
    Cnt
    est_std1 = mean(h_std);
elseif Estimation == 2 % Estimation = 2 : Adaptive averaging local standard deviations(mean of 'local std - ref_std')       
    [sort_homo, index] = sort(homo_index);
    if index >= 3
        sigma_ref = mean(h_std(index(1:3)));
    else
        sigma_ref = h_std(index(1));
    end
        sort_index = find(abs(h_std - sigma_ref) < 2);
        
    if length(sort_index) < 1
        sort_index = find(abs(h_std - sigma_ref) < 3);
    end
        est_std1 = mean(h_std(sort_index));
    
else % Estimation = 3 : Histogram Approximation (68% -> sigma)
    sig_max = th;
    [n,xout] = hist(acc_dis,0:sig_max);
    sigma68 = sum(n)*0.68;
    sum_n = 0;
    for k=1:sig_max+1
        sum_n = sum_n + n(k);
        if sum_n > sigma68
            sigma = xout(k)-1;
            break;
        end
    end

    if sigma < 0
        est_std1 = 0;
    else
        est_std1 = sigma;
    end    
end    

%===========================================================================
% Adaptive NLE (Amer -> filtering NLE)
%---------------------------------------------------------------------------
% Adaptive = 0 : No filtering NLE
% Adaptive = 1 : Amer(est_std) -> Gaussian Filtering (sigma = est_std, Homo = loc_std < th)
% Adaptive = 2 : Amer(est_std) -> Gaussian Filteing (sigma = est_std, Homo = abs(est_std - loc_std) < th)
%===========================================================================
% Adaptive = 0 : No filtering NLE
if Adaptive == 0
    est_std = est_std1;

else
    img = double(rgb2gray(imread(originalfile,'png')));
    filtering_img =GaussianFilter(img,est_std1);
    test_img = (img - filtering_img);

    sig_max = est_std1+5;

    clear dis;
    clear adis;
    clear acc_dis;
    acc_dis = 0;

    for i=1:w_size:img_width-w_size %x
        for j=1:w_size:img_height-w_size %y

            Cnt = Cnt+1;
            x_loc = i:i+w_size-1; 
            y_loc = j:j+w_size-1;
 
            s_y = std2(img(x_loc,y_loc));
    
            if Adaptive == 1
                if s_y < sig_max
                    IsHomoblk = 1;
                else
                    IsHomoblk = 0;
                end
            else
                if abs(s_y - est_std1) < 0.5
                    IsHomoblk = 1;
                else
                    IsHomoblk = 0;
                end                
            end
            
            if IsHomoblk == 1
                distance = abs(test_img(x_loc,y_loc,1)-test_img(i+center-1,j+center-1,1));        
                dis = reshape(distance, 1, w_size^2);
                adis = [dis(1:center2-1),dis(center2+1:w_size^2)];
                acc_dis = [acc_dis, adis];
            end
        end
    end

    [n,xout] = hist(acc_dis,0:sig_max);
    sigma68 = sum(n(2:length(xout)))*0.68;
    sum_n = 0;

    for k=2:sig_max+1
        sum_n = sum_n + n(k);
        
        if sum_n > sigma68
            sigma = xout(k)-1;
            break;
        end
    end

    if sigma < 0
        est_std = 0;
    else
        est_std = sigma;
    end        
    
    
end

if showimg == 1
    figure;imshow(uint8(noise));
    title(['given std='  ,num2str(std), ' , est std =',  num2str(est_std), ' ,blk number =', num2str(homoblk)]);
end