% Intensity transformations

L = 256; % uint8 images

X = imread('tire.tif');
%X = imread('pout.tif');
T = imread('teethwhitening_Woman_Smile_small.jpg');

%% negative transformation
Y = 255-X;
Y2 = imcomplement(X);

figure;
subplot(1,3,1);
imshow(X); title('original', 'FontSize', 14);
subplot(1,3,2);
imshow(Y); title('negative', 'FontSize', 14);
subplot(1,3,3);
imshow(Y2); title('imcomplement', 'FontSize', 14);

% Note: it works also for color images: it operates on each color channel
S = 255-T;
S2 = imcomplement(T);

figure;
subplot(3,4,1);
imshow(T); title('original', 'FontSize', 14);
subplot(3,4,5);
imshow(S); title('negative', 'FontSize', 14);
subplot(3,4,9);
imshow(S2); title('imcomplement', 'FontSize', 14);
for i = 1:3,
    subplot(3,4,1+i);
    imshow(T(:,:,i)); 
    subplot(3,4,5+i);
    imshow(S(:,:,i)); 
    subplot(3,4,9+i);
    imshow(S2(:,:,i)); 
end;

%% Logarithmic transformation

% note the castings: log doesn't work with uint8
c = (L-1)/log(L); % it maps L-1 into L-1
Y = uint8(c*log(double(X)+1)); 

figure;
subplot(1,2,1);
imshow(X); title('original', 'FontSize', 14);
subplot(1,2,2);
imshow(Y); title('logarithmic', 'FontSize', 14);


%% Gamma transformation
gamma = .3; % lightening
c = (L-1)/(L-1)^gamma; % it maps L-1 into L-1
% note the castings
Y = uint8(c*double(X).^gamma); 

figure;
subplot(1,2,1);
imshow(X); title('original', 'FontSize', 14);
subplot(1,2,2);
imshow(Y); title(['gamma = ' num2str(gamma)], 'FontSize', 14);

gamma = 3; % darkening
c = (L-1)/(L-1)^gamma; % it maps L-1 into L-1
% note the castings
Y = uint8(c*double(X).^gamma); 

figure;
subplot(1,2,1);
imshow(X); title('original', 'FontSize', 14);
subplot(1,2,2);
imshow(Y); title(['gamma = ' num2str(gamma)], 'FontSize', 14);

%% Contrast stretching transformation

% general shape: three piecewise linear transformations
%                defined by the arrays r and s
%                 - r: input intensities
%                 - s: corresponding output intensities
r = [0 50 200 L-1]; % note: r(1) always 0, r(4) always L-1 
s = [0 100 150 L-1];

% Even for simple transformations it is safer to use casting

R = uint8(0:L-1);
S = R; % will be mapped without casting
S2 = R; % will be mapped without casting, but with more parenthesis 
S3 = R; % will be mapped with casting

% intensities in [r(1)..r(2))
ind = find(R>=r(1) & R<r(2));
S(ind) = (R(ind)-r(1))*(s(2)-s(1))/(r(2)-r(1))+s(1);
S2(ind) = (R(ind)-r(1))*((s(2)-s(1))/(r(2)-r(1)))+s(1);
S3(ind) = uint8((double(R(ind))-r(1))*(s(2)-s(1))/(r(2)-r(1))+s(1));
% intensities in [r(2)..r(3))
ind = find(R>=r(2) & R<r(3));
S(ind) = (R(ind)-r(2))*(s(3)-s(2))/(r(3)-r(2))+s(2);
S2(ind) = (R(ind)-r(2))*((s(3)-s(2))/(r(3)-r(2)))+s(2);
S3(ind) = uint8((double(R(ind))-r(2))*(s(3)-s(2))/(r(3)-r(2))+s(2));
% intensities in [r(3)..r(4)]
ind = find(R>=r(3) & R<=r(4));
S(ind) = (R(ind)-r(3))*(s(4)-s(3))/(r(4)-r(3))+s(3);
S2(ind) = (R(ind)-r(3))*((s(4)-s(3))/(r(4)-r(3)))+s(3);
S3(ind) = uint8((double(R(ind))-r(3))*(s(4)-s(3))/(r(4)-r(3))+s(3));

figure;
subplot(2,2,1);
plot(r, s, 'LineWidth', 2); title('transformation', 'FontSize', 14);
axis equal; axis([0 L-1 0 L-1]);
subplot(2,2,2);
plot(R, S, 'LineWidth', 2); title('no casting 1', 'FontSize', 14);
axis equal; axis([0 L-1 0 L-1]);
subplot(2,2,3);
plot(R, S2, 'LineWidth', 2); title('no casting 2', 'FontSize', 14);
axis equal; axis([0 L-1 0 L-1]);
subplot(2,2,4);
plot(R, S3, 'LineWidth', 2); title('casting', 'FontSize', 14);
axis equal; axis([0 L-1 0 L-1]);

% contrast stretching with casting
Y = X;
% intensities in [r(1)..r(2)]
ind = find(X>=r(1) & X<r(2));
if ~isempty(ind), % if empty, it will produce an error
    Y(ind) = uint8((double(X(ind))-r(1))*(s(2)-s(1))/(r(2)-r(1))+s(1));
end;

% intensities in [r(2)..r(3)]
ind = find(X>=r(2) & X<r(3));
if ~isempty(ind), % if empty, it will produce an error
    Y(ind) = uint8((double(X(ind))-r(2))*(s(3)-s(2))/(r(3)-r(2))+s(2));
end;

% intensities in [r(3)..r(4)]
ind = find(X>=r(3) & X<=r(4));
if ~isempty(ind), % if empty, it will produce an error
    Y(ind) = uint8((double(X(ind))-r(3))*(s(4)-s(3))/(r(4)-r(3))+s(3));
end;

figure;
subplot(1,3,1);
plot(r, s, 'LineWidth', 2); title('transformation', 'FontSize', 14);
axis equal; axis([0 L-1 0 L-1]);
subplot(1,3,2);
imshow(X); title('original', 'FontSize', 14);
subplot(1,3,3);
imshow(Y); title('transformed', 'FontSize', 14);

% a linear transformation can be used to change the intensity
% resolution, i.e. the number of bits required for each pixel

r = [0 L-1]; 
k = 3;
s = [0 2^k-1]; % the transformed image will use k bits/pixel

Y = uint8((double(X)-r(1))*(s(2)-s(1))/(r(2)-r(1))+s(1));

figure;
subplot(1,3,1);
plot(r, s, 'LineWidth', 2); title('transformation', 'FontSize', 14);
axis equal; axis([0 L-1 0 L-1]);
subplot(1,3,2);
imshow(X); title('original', 'FontSize', 14);
subplot(1,3,3);
imshow(Y); title('transformed', 'FontSize', 14);


%% Intensity level slicing

A = 50; % [A..B] is the interval of intensity levels to be enhanced
B = 100;

R = uint8(0:L-1);
S = R; % for binary slicing
S2 = R; % for interval enhancing
S(find(R>=A & R<=B)) = 200;
S(find(R<A | R>B)) = 20;
S2(find(R>=A & R<=B)) = 200;

figure;
subplot(2,3,1);
imshow(X); title('original', 'FontSize', 14);
subplot(2,3,2);
plot(R, S, 'LineWidth', 2); title('binary slicing', 'FontSize', 14);
axis equal; axis([0 L-1 0 L-1]);
subplot(2,3,3);
plot(R, S2, 'LineWidth', 2); title('interval highlighting', 'FontSize', 14);
axis equal; axis([0 L-1 0 L-1]);

Y =X;
Y2 =X;
Y(find(X>=A & X<=B)) = 200;
Y(find(X<A | X>B)) = 20;
Y2(find(X>=A & X<=B)) = 200;

subplot(2,3,5);
imshow(Y); title('binary slicing', 'FontSize', 14);
subplot(2,3,6);
imshow(Y2); title('interval highlighting', 'FontSize', 14);

% color can be more effective for highlighting
YR = X;
YG = X;
YB = X;

YR(find(X>=A & X<=B)) = L-1;
YG(find(X>=A & X<=B)) = 0;
YB(find(X>=A & X<=B)) = 0;

Y3 = YR;
Y3(:,:,2) = YG;
Y3(:,:,3) = YB;

subplot(2,3,4);
imshow(Y3); title('color highlighting', 'FontSize', 14);


%% Bit-plane transformation

% Note: the division of an instance of the uint8 class involve a
% round-up of the result; this means that uint8(1)/2 = 1.
% Hence, the '/' operator cannot be used to emulate the
% bit-shifting effect of the division by two.
% Instead, the function idivide can be used for this goal.
% Besides, Matlab provides the function bitshift for explicitely
% realize this effect.
% Hence, two implementations are possible to realize the bit-plane
% transformation.

figure;
subplot(3,3,1);
imshow(X); title('original', 'FontSize', 14);
for i = 1:8,
    % Selection of the i-th bitplane 
    %    i = 1: least significant bit
    %    i = 8: most significant bit
    
    T = bitshift(X, 8-i); % discard the higher bitplanes
    T = bitshift(T, -7); % discard the lower bitplanes
    % now T is composed of 0- or 1- pixels
    T = T*255; % scale the 1-pixels to white (255)

    subplot(3,3,i+1);
    imshow(T);
    title(['bit-plane ' int2str(i)], 'FontSize', 14);
end;

figure;
subplot(3,3,1);
imshow(X); title('original', 'FontSize', 14);
for i = 1:8,
    % Selection of the i-th bitplane 
    %    i = 1: least significant bit
    %    i = 8: most significant bit
    
    T = idivide(X, 2^(i-1)); % discard the lower bitplanes
    % now T is odd if the i-th bitplane of X is 1 
    T = rem(T,2)*255; % selection of the first bit and scaling to white (255)
    
    subplot(3,3,i+1);
    imshow(T);
    title(['bit-plane ' int2str(i)], 'FontSize', 14);
end;