% imstretch_limit.m
% Explain the details of finding the limits for contrast stretch
% in functions "stretchlim" and "imadjust" used in imtool or imcontrast.
% >> edit stretchlimit

I = imread('Fig0316_low_contrast.tif');
[c, b] = imhist(I);
figure; plot(c);
cdf = cumsum(c)/sum(c);
figure; plot(cdf);

% By default, imadjust saturates the bottom 1% and the top 1%
% of all pixel values, thus the 2% outlier in imtool 
low = find (cdf > 0.01,1, 'first');
high = find (cdf >= 0.99, 1, 'first');

% Note the index starts from 1 in Matlab, correponding to 0 (pixel value) 
low_limit = low - 1
high_limit = high  - 1

% Using directly the Matlab function: same result
stretchlim(I)
low_scaled = (low - 1)/255
high_scaled = (high - 1)/255


% Mapping using the low_limit and high_limit

% Direct Mapping (slightly different result due to rounding)
J = uint8(double(I - low_limit)/double(high_limit - low_limit)*255);
% Saturation of outliers
J (I < low_limit) = 0;
J (I > high_limit) = 255;

% Mapping using the scaled low_scaled and high_scaled
I_scaled = double(I)/255;       % scale the pixel values to [0, 1]
J_scaled = (I_scaled - low_scaled)/(high_scaled - low_scaled);
% Saturation of outliers
J_scaled (I_scaled < low_scaled) = 0;
J_scaled (I_scaled > high_scaled) = 1;
J2 = im2uint8(J_scaled);

% Comparison with imadjust function: getting the same result
K = imadjust(I);
isequal(J2,K)

% Usage: with Gamma correction
I = imread('Fig0308(a)(fractured_spine).tif');
P = imadjust(I,[],[],0.4);


