
% Contrast enhancement based on local statistics.
%
% This script implements the technique illustred in example 3.12, 
% at page 164 of the book:
%    R.C. Gonzalez & R.E. Woods, 
%    Digital image processing, 
%    Addison-Wesley (1992, 2002 2nd ed., 2008 3rd ed.). ISBN 9780131687288. 
% 


FONTSIZE = 14;
%I = imread('pout.tif');
I = imread('filament.png');
I = rgb2gray(I);
X = double(I);

radius = 1;
num_rows = size(X,1);
num_columns = size(X,2);

% apply a linear transformation for mapping the intensity values
% in the range [0..255] 
minX = min(X(:));
maxX = max(X(:));
X = (X-minX)/(maxX-minX)*255;

m = mean(X(:));
s = std(X(:));

tic;
for i = 1:num_rows,  % row
 for j = 1:num_columns, % column
 
     % select the neighborhood
     i_r = i+(-radius:radius);
     i_c = j+(-radius:radius);
     % trim the elements out of the legal range
     i_r = i_r(find(i_r>=1 & i_r<=num_rows)); 
     i_c = i_c(find(i_c>=1 & i_c<=num_columns));

     xn = X(i_r, i_c);
     
     if (isempty(xn)),
         disp('empty xn; hit a key'); pause;
     end;
     % compute statistics for pixel (i,j)
     l_m(i,j) = mean(xn(:));
     l_s(i,j) = std(xn(:));
     
 end;
end;
disp('local statistics computed with neighborhood selection');
toc;

%disp('hit a key'); pause;

% alternative implementation of the computation of the local
% statistics
tic;
W = X;
k = m; % instead of zero, the gloabal mean can be used for padding
for i = -radius:radius,  % row shift
    if (i<0),
        Y = [k*ones(abs(i), num_columns); X(1:end+i, :)];
        % instead of padding, mirroring can also be used, e.g.:
        % Y = [double(X(1:, :); X(1:end+i, :)]; 
    elseif (i>0),
        Y = [X(i+1:end, :); k*ones(abs(i), num_columns)];      
    else
        Y = X;
    end;    
    
    for j = -radius:radius, % column shift
        % generate the neighborhood image
        if (j<0),
            Y2 = [k*ones(num_rows, abs(j)) Y(:, 1:end+j)];
        elseif (j>0),
            Y2 = [Y(:, j+1:end) k*ones(num_rows, abs(j))];      
        else
            Y2 = Y;
        end;    
        if (i~=0 | j~=0),
            W(:,:, end+1) = Y2;
        end;
 end;
disp('local statistics computed with matrix shift');
end;
toc;

% compute statistics on W
l_m2 = mean(W,3);
l_s2 = std(W, [], 3);

m2 = mean(X(:));
s2 = std(X(:));

figure;
subplot(1,2,1);
imshow(uint8(l_m));
title('local mean with neighborhood selection', 'FontSize', FONTSIZE);
subplot(1,2,2);
imshow(uint8(l_m2));
title('local mean with matrix shift', 'FontSize', FONTSIZE);

% Computation of customized contrast enhancement
% based on local statistics
%
%   X:   original image
%   Y:   transformed image
%   m:   global mean
%   s:   global standard deviation
%   l_m: local mean
%   l_s: local standard deviation
%   k0:  darkness factor
%   k1:  low contrast factor
%   k2:  high contrast factor
%   E:   enhancement factor
%
%             /
%             | E*X(i,j),  if  l_m(i,j) < k0*m AND 
%   Y(i,j) = <                    k1*s <= l_s(i,j) <= k2*s
%             | X(i,j),    otherwise 
%             \
%
% In the book, the following constants has been provided:
%   E = 4; k0 = 0.4; k1 = 0.02; k2 = 0.4;
% but they have to be customized for each image:
E = 2; k0 = 0.4; k1 = 0.05; k2 = 0.4;

Y=X;
ind = find(l_m<k0*m & l_s>=k1*s & l_s<=k2*s);
Y(ind) = E*X(ind);

figure;
subplot(1,3,1);
imshow(uint8(X));
title('original', 'FontSize', FONTSIZE);
subplot(1,3,2);
imshow(histeq(uint8(X)));
title('equalized', 'FontSize', FONTSIZE);
subplot(1,3,3);
imshow(uint8(Y));
title('local stats', 'FontSize', FONTSIZE);

figure;
subplot(1,3,1);
imshow(l_m<k0*m);
title('dark region (l\_m < k0*m)', 'FontSize', FONTSIZE);
subplot(1,3,2);
imshow(l_s >= k1*s);
title('detail region (l\_s >= k1*s)', 'FontSize', FONTSIZE);
subplot(1,3,3);
imshow(l_s<=k2*s);
title('non-edge region (l\_s <= k2*s)', 'FontSize', FONTSIZE);

