Activity 4: Enhancement by Histogram Manipulation

By looking at the histogram of an image, one can can determine wether the image is too bright or too dark, has poor or good contrast, i.e, if it has good dynamic range. Usually, we want an image to span all the gray levels resulting to a high contrast image. By manipulating the histogram of an image we can improve the quality of the image highlighting or enhancing certain features, or depending on the user, can mimick other imaging system[1].

The histogram of an image when normalized is equal to graylevel probability distribution function (PDF). We can alter the PDF of our image to the desired PDF by backprojecting the grayscale level values of the desired cumulative distribution function (CDF) to the grayscale values of the original CDF.

The steps in overview are as follows:
1. Per pixel in the image having grayscale r, find its CDF value T(r).
2. Find the value of T(r) in the y-axis of the desired CDF G(z).
3. Find the z value of this G(z).
4. Replace the image with grayscale value r to grayscale value z.
This process is called HISTOGRAM EQUALIZATION.

I have provided two methods in solving the problem, the difference is essentially in obtaining the PDF and CDF of the image. In the first solution, we use binning (middle of bin, i.e, 0.5, 1.5...254.5) in obtaining the histogram of the image. In the second solution, we use the function tabul of scilab which tabulates the frequency of distinct grayscale level values.

1st solution:

directory = 'C:\Documents and Settings\Orly\Desktop\AP186\Activity 4';
chdir(directory);
//load image in grayscale
img1=gray_imread('img2.png');
img2=img1; //holder
//obtain the histogram
bits=255;
bins = 0.5:bits; bins = bins./bits; //bins for [0,1] image.
for i=1:bits,
frequency(i) = sum(img1>=(bins(i)-1/(2*bits))&img1<(bins(i)+1/(2*bits))); end; frequency = frequency./max(cumsum(frequency)); f1=scf(1); plot(bins.*bits, frequency); //normalized histogram f2=scf(2); cdf = cumsum(frequency); //cumulative distribution function plot(bins.*bits,cdf); //desired cdf is a straight line, can be change descdf = (0:1.0/(length(cdf)-1):1); f3=scf(3); plot(bins.*bits, descdf); //mapping for i=1:bits-1, index = find(descdf>=cdf(i) & descdf=1 then //satisfied condtion above
//get first index
imgval = bins(index(1)); //desired pixel value
imgcur = bins(i); //current pixel value of image
img2(img2>=(imgcur-1/(2*bits))&img2<(imgcur+1/(2*bits))) = imgval; end, end; //for comparison of images f4=scf(4); imshow(img1); f5=scf(5); imshow(img2); //check the resulting CDF and PDF of the equalized image for i=1:bits, frequency2(i) = sum(img2>=(bins(i)-1/(2*bits))&img2<(bins(i)+1/(2*bits))); end; f6=scf(6); plot(bins.*bits, frequency2./max(cumsum(frequency2))); f7=scf(7); plot(bins.*bits, cumsum(frequency2)./max(cumsum(frequency2)))

The image below is the original image which we will use to demonstrate contrast enhancement.

Poor contrast imageHistogram of the original image, min = 0.5, max =181.5

CDF of the original image.
Desired CDF which is just the CDF of a uniform distribution

Using the above code, the resulting image and histograms is as follows.

Resulting image: Note that in our process, we span all the gray levels and redistribute them in the image in a uniform process. The midtones are enhanced but the problem is the image got saturated on parts with high gray level and too much shadows on parts with low gray level value. Also, the result is a bit noisy compared to the original.

Histogram of the resulting image, min = 0.5, max =253.5. The contrast of the resulting image base on the minimum and maximum gray values increased dramatically.

The CDF of the resulting image follows the shape of the desired CDF above.

Main problem encountered with the above solution is that for those images with limited dynamic range (i.e, 100-150), the resulting image would be too saturated for those pixels with higher pixel value (130-150) and too dark for those pixels with lower pixel value (100-120) because essentially, we are uniformly stretching the image into all the gray levels (0-255). The resulting image might not be too pleasing but it would have better dynamic range.

Another solution for this is to uniformly distribute the image (i.e, linear CDF) in the range of the original histogram. (i.e., from 100-150). In essence, we LOCALIZED the histogram equalization.

2nd Solution
directory = 'C:\Documents and Settings\UP DILIMAN\Desktop\ORLY\Acads\AP186\Activity 4';
chdir(directory);

//load image in grayscale
img1=gray_imread('img2.png');
info = imfinfo('img2.png');
img2=img1; //holder
bits=(2^info.Depth)-1;
//find frequency of distinct values
hist=tabul(img1,'i'); //arrange in increasing order
frequency = hist(:,2)./max(cumsum(hist(:,2))); //normalized frequency
bins=hist(:,1);

//displaying histogram
histogram=fhistogram(img1, info.Depth);
frequency2 = histogram(:,2)./max(cumsum(histogram(:,2)));
f1=scf(1);
plot(histogram(:,1).*bits, frequency2);

//find the minimum and maximum grayscale values
index=find(frequency~=0);
minB = bins(index(1))*255;
maxB = bins(index(length(index)))*255;

f2=scf(2);
cdf = cumsum(frequency);
plot(bins.*bits,cdf);
//desired cdf is a straight line, can be change
descdf = (min(cdf):(max(cdf)-min(cdf))/(length(cdf)-1):max(cdf));
desbins = (min(bins):(max(bins)-min(bins))/(length(bins)-1):max(bins));
f3=scf(3);
plot(desbins.*255, descdf);

//mapping
for i=1:length(bins)-1,
//returns the index in descdf where cdf(i) belongs
index = find(descdf>=cdf(i)); //returns a list of indices
if sum(index)>=1 then
//get first index
imgval = desbins(index(1)); //desired pixel value
imgcur = bins(i); //current pixel value
img2(img2==imgcur) = imgval;
else
x_message('cdf index is' + string(i));
end,
end;

f4=scf(4);
imshow(img1);
f5=scf(5);
imshow(img2);
f6=scf(6);
histogram2=fhistogram(img2, info.Depth);
frequency3 = histogram2(:,2)./max(cumsum(histogram2(:,2)));
plot(histogram2(:,1).*bits, frequency3);
f7=scf(7);
plot(histogram2(:,1).*bits, cumsum(frequency3));

code of filename: fhistogram.sce

function histogram = fhistogram(img, depth)
bits=(2^depth)-1;
histogram = zeros(bits,2);
bins = 0.5:bits; bins = bins./bits; //bins for [0,1] image.
for i=1:bits,
frequency(i) = sum(img>=(bins(i)-1/(2*bits))&img<(bins(i)+1/(2*bits))); end; histogram(:,1)=bins'; histogram(:,2)=frequency; endfunction
Using tabul, the bins are the distinct grayscale values and the frequency is just the occurrence of that distinct grayscale value.

Resulting image after mapping. Note that the image has better contrast compared to the reconstruction using the first method. The advantage of this reconstruction is that, the image would not be too saturated on some areas and too dark on other areas because we uniformly distributed the image in the range of the min and max grayscale values.

Resulting histogram, min=0, max=182

Resulting CDF, note that it is almost similar to the desired CDF.

By defining another set of bins ('desbins') that are also uniform in the min and max grayscale level values, we have essentially equally redistributed the histogram of the image.

Using other images:

Original: left, Contrast enhanced: right. The reconstructed image has better contrast compared to the original image. Image obtained from source [2].
Original: left, Reconstructed: right. Note that the effect of the histogram equalization is not quite visible. This is because the dynamic range of the image is limited as can be seen in the histogram and CDF of the reconstructed image. Image was obtained from source [3].
Note that in this histogram and CDF, we can see clearly how the method works. The histogram was equally distributed in the range of the histogram (min - max) hence the resulting CDF is linear in that range only.

USING OTHER CDF:
As a proof of principle, we use an exponential CDF to demonstrate what happens on the image when using nonlinear CDF. Depending on the user, we can opt to highlight the darker areas (low grayscale values), the whiter areas (high grayscale values) or the midtones (middle grayscale values). When highlighting the darker areas, we can use logarithmic CDF, when highlighting whiter areas, we can use exponential CDF and when highlighting the midtones we can use S-curve CDF.

The code above can be altered by defining the appropriate desired CDF. We just inserted the following code:
desbins = (min(bins):(max(bins)-min(bins))/(length(bins)-1):max(bins)); sigma =2; //strength of exponential descdf = exp(sigma*desbins);
descdf = (descdf - descdf(1))./max(descdf-descdf(1));

Resulting image using exponential CDF. Note that the resulting image highlighted the higher grayscale values (whiter areas) as expected from the shape of an exponential curve.

Left: Desired CDF, Right: CDF of the resulting image.

We have provided a solution for contrast enhancing an image by histogram manipulation. Our solution involves locally equalizing the histogram of an image using the available min and max pixel values of the image.

It must be noted that there is no single CDF that can enhance the contrast of every image. Each CDF has its own limitations and use. Care must be taken when contrast enhancing an image and knowledge of how each CDF alters the image must be understood.

Please email me (orly.tarun@gmail.com) if you plan to use the code above! Please leave your comments and suggestions. =)

For this activity, I give myself a grade of 10 for performing acceptable contrast enhancement algorithms.

References:
[1] Activity 4 Manual.
[2] http://myweb.lsbu.ac.uk/dirt/museum/margaret/08--452-1000120.jpg
[3] http://www.andrew.cmu.edu/user/timothyz/hw3/low_contrast.jpg