%
% Purpose: Histogram Equalization
% Written by: H.E.
% Date: 2/2009
%

close all;
clear;

% Fuzzify dark, gray, and whiter membership functions
L=256; z=0:L-1;
a=round((L-1)/4); b=round((L-1)/2);
figure;
subplot(2,1,1);
darkc=[ones(1,a+1) -(1:b-a)/(b-a)+1 zeros(1,L-1-b:L-1)];
grayc=[zeros(1,a+1) (1:b-a)/(b-a) -(1:b-a)/(b-a)+1 zeros(1,L-1-(a+b))];
brightc=[zeros(1,b+1) (1:b-a)/(b-a) ones(1,L-1-(a+b))];
plot(z,darkc,z,grayc,z,brightc); title('Input memberships');
ylabel('\mu(dark, gray, bright)'); xlabel('Intensity');

% Fuzzify darker, grayer, and whiter membership functions
vd=0; vg=127; vb=255;
subplot(2,1,2);
darker=[zeros(1,vd) 1 zeros(1,L-1-vd)];
grayer=[zeros(1,vg) 1 zeros(1,L-1-vg)];
brighter=[zeros(1,vb) 1 zeros(1,L-1-vb)];
plot(z,darker,z,grayer,z,brighter); title('Output memberships');
ylabel('\mu(darker, grayer, brighter)'); xlabel('Intensity');

% Read input
im=imread('Aerial_RGB.jpg');
figure; imagesc(im);
title('Low contrast aerial image');

[x,y,c]=size(im); im_HE=im; im_FL=im;
for k=1:3

    ['Processing color ...',int2str(k)]
    
    % Process color k
    imC=im(:,:,k);

    % Compute histogram equalization
    v=reshape(imC,1,x*y); 
    pdf=hist(v,z)/(x*y); cdf=pdf;
    figure; plot(pdf);
    for i=2:L
        cdf(i)=pdf(i)+cdf(i-1);
    end
    cdf=round((L-1)*cdf); imC_HE=imC; imC_FL=imC;
    for i=1:x
        for j=1:y
            
            % Histogram equalization step
            s=imC(i,j) == z;
            if(sum(s) ~= 0)
                imC_HE(i,j)=cdf(find(s));
            end
            
            % Fuzzy logic step
            m_dc=darkc(imC(i,j));
            m_gc=grayc(imC(i,j));
            m_bc=brightc(imC(i,j));
            imC_FL(i,j)=(m_dc*vd+m_gc*vg+m_bc*vb)/(m_dc+m_gc+m_bc);
            
        end
    end
    im_HE(:,:,k)=imC_HE; im_FL(:,:,k)=imC_FL;    
    
end    
figure; imagesc(im_HE);
title('Histogram equalized aerial image');

figure; imagesc(im_FL);
title('Intensity transformed (Fuzzy) aerial image');



%v=reshape(imR_n,1,x*y); pdf_n=hist(v,z)/(x*y);
%figure; subplot(2,2,1); plot(z,pdf); xlabel('rk'); ylabel('pr(rk)'); grid;
%subplot(2,2,2); plot(z,cdf); xlabel('rk'); ylabel('sr'); grid;
%subplot(2,2,3); plot(z,pdf_n); xlabel('sk)'); ylabel('ps(sk)'); grid;