% enh_pixel.m:  image enhancement based on pixel processing
% include image negatives, power law transform, log transform, averaging
% copyright (c) 1997-2005 by Yu Hen Hui
% created 10/13/97
% last modified: 9/2005
%
clear all, close all
echo on
% 1. load a gray scale image
[xi,ind]=imread('images\p64.bmp');
[nh,nv]=size(xi);  % size of the image
x=ind2gray(xi,ind);  % pixel value in [0 255]
figure(1),
imshow(x),title('original image')
posi=get(1,'position'); % a 1 x 4 vector gives lower left corner and size of image
set(1,'position',[0 540 posi(3:4)])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Image negatives:  T(r) = L-1-r   L=256
% Here the image is between 0 and 1 so we use 1-r
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pause
%
xneg=255-x;
figure(2),imshow(xneg),title('image negatives')
posi=get(2,'position'); % a 1 x 4 vector gives lower left corner and size of image
set(2,'position',[160 540 posi(3:4)])

pause
%
% Let us now consider some transformation function
%
r=[0:0.01:1]; 
s=[0.5*r(find(r<0.2)) 0.1+1.5*(r(find(0.2 <= r & r <= 0.7))-0.2)...
    1+0.5*(r(find(r>0.7))-1)];
figure(7),subplot(131),plot(r,s,'-'),title('mid-range strech')
set(7,'position',[275   265   712   194])
%
x1=double(x)/255; % x1 range to 0 and 1
y=(0.5*x1).*[x1<0.2]+(0.1+1.5*(x1-0.2)).*[0.2 <= x1 & x1 <= 0.7] ...
    +(1+0.5*(x1-1)).*[x1 > 0.7];
[xmid,indy]=cmunique(y);
figure(3),imshow(xmid,indy),title('mid-range stretched')
posi=get(3,'position'); % a 1 x 4 vector gives lower left corner and size of image
set(3,'position',[320 540 posi(3:4)])
pause

c = 255/log10(1+255); 
% c is chosen to ensure only 256 gray levels after transform
r=[0:255]; s=c*log10(1+abs(r));
figure(7),subplot(132),plot(r,s,'-'),title('logarithm compression')
axis([0 256 0 256])

xlog=uint8(round(c*log10(1+double(x)))); 
figure(4),imshow(xlog),title('log-transformed')
posi=get(4,'position'); % a 1 x 4 vector gives lower left corner and size of image
set(4,'position',[480 540 posi(3:4)]),

pause
%
% Next, we consider intensity slicing
%
r=[0:0.01:1];  s=r.*[r < 0.2 | r > .4] + .6*[0.2 <= r & r <=0.4];
figure(7),subplot(133),plot(r,s,'-'),title('intensity slicing')

y=x1.*[x1 < 0.2 | x1 > 0.4] + 0.6*[0.2 <= x1 & x1 <=0.4];
[xints,indy]=cmunique(y);
figure(5),imshow(xints,indy),title('intensity sliced')
posi=get(5,'position'); % a 1 x 4 vector gives lower left corner and size of image
set(5,'position',[640 540 posi(3:4)]),pause
%
% Now, let us examine the histogram of the original image:
% This can be conveniently done using the command "imhist" 

[N,grayscale]=imhist(x); 
figure(8),set(8,'position',[5    40   560   420])
subplot(211),
stem(grayscale,N),title('histogram of original image')
axis([0 256 0 max(N)])
pause

% Suppose we want to equalize the histogram. This can be done
% using the command "histeq"
xeq=histeq(x);  % equalized color map indy is a reshuffle

figure(6),imshow(xeq),title('histogram equalized')
posi=get(6,'position'); % a 1 x 4 vector gives lower left corner and size of image
set(6,'position',[800 540 posi(3:4)]),pause

[N,grayscale]=imhist(xeq);
figure(8),
subplot(212),stem(grayscale,N),title('equalized histogram')
axis([0 256 0 max(N)])
% YOu can see the histogram is not really quite flattened!