Monday, January 4, 2010

[Review] A4 : Enhancement by Histogram Manipulation

A histogram of an image provides a quick glance of the tonal information that can be seen in the image. An image histogram can tell whether an image occupies a wide range of contrast or it is clumped in a relatively closer range [1].

Graphically, it plots the tonal (grayscale) value vs. the frequency of the said value. The extreme left side of the graph (small grayscale values) represent the darker (black) pixels while the extreme right represent the brighter pixels(white). The middle values represent the gray pixels. A histogram, when normalized, is called a Probability Distribution Function (PDF). It replaces the frequency of the grayscale values with its probability of appearing in the image. When the PDF is added cumulatively, the graph becomes a Cumulative Distribution Function (CDF). Examples of an image histogram, a PDF and CDF is shown below:

An Image Histrogram



A Probability Distribution Function (PDF) of a certain image [3]


A Cumulative Distribution Function (CDF) of a certain image [3]


Unfortunately, some images might be too dark (histogram is concentrated on the left side) or to bright (histogram concentrated on the right side) for some applications. An example is shown below:


Different histogram distribution [2]

To solve this problem, we could alter an image's histogram in order to redistribute its tonal information (grayscale values) to a wider range of values or to a specific range of interest. This process is called histogram manipulation. The steps for histrogram manipulation (also called CDF backprojection is enumerated below [3].
1. From the pixel value in the original image, find the corresponding CDF value.
2. Trace this CDF value in the desired CDF graph.
3. Find the corresponding new grayscale value with the given CDF value in the desired CDF graph (G(z)).
4. Replace all pixels in the original image with this new value.

To implement the activity, I have downloaded an image from a web showing an x-ray image with an uneven contrast contribution:


An X-ray Image

Below is a plot if its CDF and beside it is a plot of the desired histogram for it to have a more distributed contrast. We could see that many of the pixels are distributed in the darker shades as shown by the graph which is skewed to the left.


And below is the snippet of the code for the backprojection:

xdomain = linspace(double(min(image(:))),double(max(image(:))),256);
yrange = linspace(0,1,256);

for i=1:1:length(cdf)
% find the cdf value of the nonzero grayscale value in the image
location=find(image == i);
currentcdf = cdf(i);

% find the new grayscale value given the CDF value
% based on the linear CDF
new_grayscale_value = max(round(xdomain(find(yrange <= currentcdf))));
% replace old grayscale values with the new numbers
image2(location) = new_grayscale_value(1);
% end
end

Below is the resulting histogram equalized image from the backprojection process:
From the images above, we could see more details in the histogram equalized image as it "lights up" areas of the image that is too dim to see in the original image.

Shown below is the resulting CDF of the histogram equalized image. We could see that the new CDF now more similar to the desired CDF as it follows a more linear form.

The graphs below shows the PDF of the original image and the PDF of the modified image. The second PDF shows that the pixel values are distributed in a wider range of values.


Further, our eyes can mimic a nonlinear response. Examples of nonlinear response are exponential and logarithmic. Exponential and Logarithmic CDFs can actually redistribute pixel values either to the darker or brighter shades, whichever is necessary for a specific application. Shown below are the CDF graphs for a histogram equalized image, for an image that has more brighter pixels (exponentials) and one that has more darker pixels (logarithmic)


The following lines of code were executed to come up with the CDFs shown above:

degree = 4;
exponential = exp(degree*yrange)-1;
exponential = exponential./max(exponential);
logarithmic = 1 - exponential;
logarithmic = logarithmic(end:-1:1);

Exponential CDFs redistribute the pixels in a way that it subdues the darker pixels, and thus, in the process, highlights the brighter pixels. A logarithmic CDF do the exact opposite. This explanation is illustrated by using the histogram equalized image in the previous operations and let in undergo again back projections. But this time, the exponential and the logarithmic CDFs are now the desired CDFs.
As expected, the exponential CDF highlighted the brighter pixels and subdued the darker pixels while the opposite happens in the logarithmic CDF. The PDF and the CDF of the resulting images are shown below. The PDF shows that for the exponential PDF, the pixels are more distributed in the lighter side of the spectrum while for the logarithmic, the pixels are distributed in the other side.


Complete Code Listing:

function act4(image)
image = imread(image);
image2 = image;
[pdf cdf] = get_pdf_cdf(image);
% cdf_orig
% create x-axis and y-axis
xdomain = linspace(double(min(image(:))),double(max(image(:))),256);
yrange = linspace(0,1,256);

% plot the original cdf and the desired cdf
% figure
% subplot(121), plot(cdf)
% xlabel('Possible Grayscale Values')
% ylabel('CDF')
% title('Original CDF')
% subplot(122), plot(xdomain, yrange)
% xlabel('Possible Grayscale Values')
% ylabel('CDF')
% title('Desired CDF')

for i=1:1:length(cdf)
% find the cdf value of the nonzero grayscale value in the image
location=find(image == i);
currentcdf = cdf(i);

% find the new grayscale value given the CDF value
% based on the linear CDF
new_grayscale_value = max(round(xdomain(find(yrange <= currentcdf))));
% replace old grayscale values with the new numbers
image2(location) = new_grayscale_value(1);
% end
end

% figure
% subplot(121), imshow(image)
% title('Original Image')
% subplot(122), imshow(image2)
% title('Histogram Equalized')

[pdf2 cdf2] = get_pdf_cdf(image2);

% figure
% subplot(121), plot(cdf)
% title('Original Image')
% subplot(122), plot(cdf2)
% title('Histogram Equalized')

% % for the nonlinear desired CDF
image3 = image2;
image4 = image2;
degree = 4;
exponential = exp(degree*yrange)-1;
exponential = exponential./max(exponential);
logarithmic = 1 - exponential;
logarithmic = logarithmic(end:-1:1);

% % plot the original cdf and the desired cdf
% figure
% subplot(131), plot(cdf2)
% xlabel('Possible Grayscale Values')
% ylabel('CDF')
% title('Original CDF')
% subplot(132), plot(xdomain, exponential)
% xlabel('Possible Grayscale Values')
% ylabel('CDF')
% title('Desired CDF')
% subplot(133), plot(xdomain, logarithmic)
% xlabel('Possible Grayscale Values')
% ylabel('CDF')
% title('Desired CDF')
%
%
for i=1:1:length(cdf2)
% find the cdf value of the nonzero grayscale value in the image
location=find(image == i);
currentcdf = cdf2(i);

% find the new grayscale value given the CDF value
% based on the linear CDF
new_grayscale_value3 = max(round(xdomain(find(exponential <= currentcdf))));
new_grayscale_value4 = max(round(xdomain(find(logarithmic <= currentcdf))));
% replace old grayscale values with the new numbers
image3(location) = new_grayscale_value3(1);
image4(location) = new_grayscale_value4(1);
% end
end
%
% figure
% subplot(131), imshow(image)
% title('Original Image')
% subplot(132), imshow(image3)
% title('Exponential Histogram')
% subplot(133), imshow(image4)
% title('Logarithmic Histogram')
%
%
% figure
% subplot(131), imshow(image2)
% title('Histogram Equalized Image')
% subplot(132), imshow(image3)
% title('Exponential Histogram')
% subplot(133), imshow(image4)
% title('Logarithmic Histogram')
%
[pdf3 cdf3] = get_pdf_cdf(image3);
[pdf4 cdf4] = get_pdf_cdf(image4);
%
% % % plot the histogram equalized cdf and the exponential cdfs
% figure
% subplot(131), plot(cdf2)
% xlabel('Possible Grayscale Values')
% ylabel('CDF')
% title('Original CDF')
% subplot(132), plot(cdf3)
% xlabel('Possible Grayscale Values')
% ylabel('CDF')
% title('Resulting CDF')
% subplot(133), plot(cdf4)
% xlabel('Possible Grayscale Values')
% ylabel('CDF')
% title('Resulting CDF')

figure
subplot(121), plot(pdf)
xlabel('Possible Grayscale Values')
ylabel('Probability')
title('Original PDF')
subplot(122), plot(pdf2)
xlabel('Possible Grayscale Values')
ylabel('Probability')
title('Resulting PDF')

figure
subplot(121), plot(pdf3)
xlabel('Possible Grayscale Values')
ylabel('Probability')
title('Exponential PDF')
subplot(122), plot(pdf4)
xlabel('Possible Grayscale Values')
ylabel('Probability')
title('Logarithmic PDF')

function [pdf cdf] = get_pdf_cdf(image)
[count, values] = imhist(image);
numpoints = sum(count);
pdf = count/numpoints;
cdf_orig = zeros(256,1);
cdf_orig(1) = pdf(1);
for i=2:1:256
cdf_orig(i) = pdf(i) + cdf_orig(i-1);
end
cdf = cdf_orig;

I give myself a grade of 9 for completing the requirements. Though I did have a hard time looking for a solution for back projecting the new values since the points in the graph have only discrete values. *The '<=' sign did the trick. :)



[1] http://en.wikipedia.org/wiki/Image_histogram
[3] Maam Maricor's Activity Handout