% ce script Matlab illustre la definition d'entropie, entropie conjointe et entropie conditionnelle
% Il utilise les pixels d'une image couleur, et calcule l'entropie des canaux rouge et vert, leur entropie conjointe
% et leur entropie conditionnelle
% Master SICOM, octobre 2006
% J. Rendas
%
close all
clear all

i1 = imread('Photo.jpg');
h=figure(1); set(h, 'Position', [661 535 560 420]);
image(i1);
title('complete image'); hold on; 

% select a small section of the image
disp('select an image window for analysis (not too big, otherwise program will go on forever...)')
h_t = text(size(i1,2)/2-300, size(i1,1)+150,'click on the left', 'Color','r'); 
[x,y] = ginput(1);
c_min = x;
line([x x], [1,size(i1,1)],'Color','r')
delete(h_t)
h_t = text(size(i1,2)/2-300, size(i1,1)+150, 'click on the bottom', 'Color','r');
[x,y] = ginput(1);
l_max = y;
h = line([1,size(i1,2)],[y y],'Color','r')
delete(h_t)
h_t = text(size(i1,2)/2-300, size(i1,1)+150,'click on the right', 'Color','r');
[x,y] = ginput(1);
c_max = x;
h = line([x x], [1,size(i1,1)],'Color','r')
delete(h_t)
h_t = text(size(i1,2)/2-300, size(i1,1)+150,'click on the top',  'Color','r');
[x,y] = ginput(1);
l_min = y;
h = line([1,size(i1,2)],[y y],'Color','r');
delete(h_t)

i11 = i1(ceil(l_min):fix(l_max),ceil(c_min):fix(c_max),:);
%
% image 
h=figure(2); set(h, 'Position',[18 527 560 420]);
clf
image(i11)
title('selected window')
axis equal
axis off

% couleurs
%
h=figure(3);set(h,'Position',[1046 57 219 420]);
clf
i11_r = i11(:,:,1);
subplot(3,1,1)
image(i11_r)
title('red')
colormap('gray')
axis equal
axis tight; axis off
i11_g = i11(:,:,2);
subplot(3,1,2)
image(i11_g)
title('green')
colormap('gray')
axis equal
axis tight;  axis off
i11_b = i11(:,:,3);
subplot(3,1,3)
image(i11_b)
title('blue')
colormap('gray')
axis equal
axis tight;  axis off
pause(0.1) % to allow display
% 
% compute the entropy of the red and green channels
%

i11r= double(i11_r(:));
i11g= double(i11_g(:));
m_r = min(i11r);
M_r = max(i11r);
m_g = min(i11g);
M_g = max(i11g);
for i=m_r:M_r
   h_r(i+1) = sum(i11r==i);
end
h_r = h_r/length(i11r); % normalize
for i=m_g:M_g
   h_g(i+1) = sum(i11g==i);
end
h_g = h_g/length(i11g); % normalize
%
H_r = entropy(h_r);
H_g = entropy(h_g);
%
disp('Entropy is bounded by the logarithm of the number of distinct values of the variable:')
disp(sprintf('H(r) = %.2f < log_2(%d) = %.2f', H_r, M_r-m_r+1, log2(M_r-m_r+1)));
disp(sprintf('H(g) = %.2f < log_2(%d) = %.2f', H_g, M_g-m_g+1, log2(M_g-m_g+1)));
disp('H(r): entropy of the red channel, H(g): entropy of the green channel')

%
close(1)
h=figure(4); set(h, 'Position', [661 535 560 420]);
clf
subplot(2,1,1)
bar(m_r:M_r, h_r(m_r+1:M_r+1))
axis tight
title('h_r: histogram of red channel')
subplot(2,1,2)
axis tight
bar(m_g:M_g, h_g(m_g+1:M_g+1))
title('h_g: histogram of green channel')
disp('press enter to continue')
pause
%
% joint distribution
%
for i=m_r:M_r
   for j=m_g:M_g
      h_rg(i+1,j+1) = sum(i11r==i & i11g==j);
      disp(i)
   end
end
h_rg =h_rg/length(i11r);
%
H_rg = entropy(h_rg(:))
disp('Joint entropy H(r,g) is bounded above by the sum of the entropies (H(r)+H(g)):')
disp(sprintf('%.2f = H(r,g)  < H(r)+H(g) = %.2f +%.2f = %.2f',...
   H_rg, H_r, H_g, H_r+H_g));
disp('and below by the maximum of the entropies')
disp(sprintf('%.2f = H(r,g)  > max(H(r),H(g)) = max(%.2f,%.2f)=%.2f',...
   H_rg, H_r, H_g, max(H_r, H_g) ) );
%
h = figure(5); set(h, 'Position', [661 535 560 420]);
imagesc(m_r+1:M_r+1,m_g:M_g+1,h_rg)
axis tight; axis equal; 
title('joint distribution')
disp('press enter to continue')
pause
disp('entering computation of conditional entropy...')
pause(0.1)
%
% conditional distributions
%
for j=m_g:M_g
   if (h_g(j+1) ~= 0)
      for i=m_r:M_r
         h_r_g(i+1,j+1) = h_rg(i+1,j+1)/h_g(j+1) ;
      end
      H_r_g(j+1) = entropy(h_r_g(:,j+1)); % conditional entropy of the red channel given i11_g == j
   else
      for i=m_r:M_r
         h_r_g(i+1,j+1) = NaN ;
      end
   end
   
end
%
h = figure(5); clf; set(h,'Position', [661 535 560 420]); hold on
clf
imagesc(m_r+1:M_r+1,m_g:M_g+1,h_r_g)
title('conditional distribution p(r|g)')
xlabel('red')
ylabel('green')
axis xy
text(40, 30, 'each line is a (color coded) probability law', 'Color','w')
text(10, 200,'you can see a particular conditional distribution, and its entropy','Color','y');
text(10, 190, 'by clicking in the row', 'Color','y')

[x,y] = ginput(1);
if ~isempty(x)
   figure(7);
   bar(h_r_g(:,fix(y)));
   text(10,0.9*max(h_r_g(:,fix(y))),sprintf('entropy H(r|g=%d)=%.2f', fix(y), H_r_g(fix(y))),...
      'Color', 'r');
end



%
% conditional entropy of the red channel given the green channel
%
CH_r_g=0;
for j=m_g:M_g
   if (h_g(j+1) ~= 0 )
      CH_r_g = CH_r_g + h_g(j+1)*H_r_g(j+1);
   end
end
disp('Conditional entropy H(r|g) is bounded above by the entropy H(r):')
disp(sprintf('H(r|g) = %f < %.2f =H(r)', CH_r_g, H_r))