% dip05_2.m             gr09mar15   gr10sep22
% Histogram manipulation
clear all; %clf(1); clf(2);

 f = double(imread('.\img\pentagon.tif')); % f = f(1:2:end,1:2:end);
 f = double(imread('.\img\cameraman.tif'));

h = hist(f(:),[0:255]) / length(f(:));      % one could use 'imhist'
figure(1);  imshow(f,[]);
figure(2);  subplot(221); bar(h); axis([0 255 0 max(h)]); 
pause

cdf = cumsum(h);        % cumulative distribution, to be used for equalization
 figure(2);  subplot(222); plot([0:255],cdf);
 pause

feq = f;
for lev = 0:255;                        %  equalization
    feq(f==lev) = cdf(lev+1)*255;
    end;
heq = hist(feq(:),[0:255]) / length(feq(:));     % hist. of the equalized image
figure(1);  imshow([f feq],[]);
figure(2);  subplot(222); bar(heq); axis([0 255 0 max(heq)]); 
        % note 'holes' in the histogram

% g = histeq(f/255,256); imshow(g,[]);      % one could just use 'histeq'

return

%% histogram specification 
% is not easily implemented, due to the necessity of inverting the second mapping, G.
%  Matlab operates by minimizing the error between the desired
% and actual cumulative histograms.

%  desired histogram (Matlab takes care of normalizing it):
 h_des = ones(1,256)/256;     % (simple equalization)
 h_des = [zeros(1,15)+.2   ones(1,45)   zeros(1,50)+.3  ones(1,70)   zeros(1,76)+.03];
% h_des = (1:256);
% h_des = (1:256).^0.3;

g = histeq(f/255,h_des)*256;         %  histogram specification
h_act = hist(g(:),[0:255]) / length(g(:));     %  actual histogram yielded
figure(1);  imshow([f feq g],[]);
figure(2);  subplot(223); bar(h_des); axis([0 255 0 max(h_des)]);
               subplot(224); bar(h_act); axis([0 255 0 max(h_act)]);
% pause;
% feqf = .4 * f + .6 * feq;             %  weighted average between original and equalized
% figure(1);  imshow([f feq; g feqf],[]);


%%  for color images too
% see first HSV color space:
% https://it.mathworks.com/help/images/understanding-color-spaces-and-color-space-conversion.html

f = double(imread('peppers.png'))/255;

 h_des = ones(1,256);     % (simple equalization)
 h_des = (1:256);
 h_des = (1:256).^0.3;
figure(30), plot(h_des);
 
RGB = f;
if size(RGB,3)==3
    HSV = rgb2hsv(RGB);
    f = HSV(:,:,3);
end

g = histeq(f,h_des);

if size(RGB,3)==3
    figure(31), montage({f,g});
    HSVeq = HSV;
    HSVeq(:,:,3) = g; %/255;
    g = hsv2rgb(HSVeq);
end

figure(32), montage({RGB,g});


%%  %%%%%  local histogram equalization
NT = 4;
g = adapthisteq(f/255,'NumTiles',[NT NT])*255;   % Try NumTiles from 2 to 32 - default 8
imshow([f g],[]);
return;


%%%%%%        see also dip05_2A_exact_histogram_specification     %%%%%%%%%%%%


