Lecture - Image Processing

Images are characterized by a certain intensity distribution, i.e. the intensity histogram (that often has an uneven shape). Astronomical images have in many cases an abundance of dark tones and few bright ones, representing the background with a few bright stars here and there, respectively. We have already seen that a manually applied transfer function can help increasing the contrast and, in turn, use more of the available range of brightness values that are shown on the screen or the printout. But, what if we want to go a step further and not only change the dynamic range and / or contrast of the image but the intensity distribution itself in order to suppress or promote different intensity levels? The process of histogram specification makes exactly this possible. Often also named histogram matching or histogram stretching, this method does exactly what it says: stretching the image intensity histogram via a certain transfer function so that the resulting histogram of the new image resembles the desired distribution. We will look at several examples of some distributions which are usually applied for deep-sky objects or lunar images.

The process works as follows. Below each step a python code sample is provided.

1) Compute the intensity histogram \(h(p)\) and the corresponding cumulative histogram \(h\_sum(p)\) of the original image, with \(p\) being the pixel value. On a programming level, these histograms are simply represented by arrays with the length of 256 for a standard 8-bit grayscale image. Each value indicates the histogram height at that brightness value (0...255).

# build histogram h
for y in range(0, ymax-1):
	for x in range(0, xmax-1):
		p = image(x,y)
		h[p] = h[p] + 1

# build cumul. histogram h_sum
for i in range(0, 255):
	h_sum[i] = h_sum[i] + h[i]

2) Generate a target histogram \(th(p)\), i.e. the desired shape of the intensity histogram of the transformed image. For this to happen we apply a transfer function \(a(p)\), which can be a for example a Gaussian distribution. The target histogram becomes 

\(th(p) = a(p/256)\)

for a function that runs from 0...1. This operation is executed within a loop over all pixel values.

# build target histogram th
# call a target function, use a Gaussian function here
for i in range(0, 255):
	th[i] = gaussian[i/256]

3) Compute the cumulative target histogram \(th\_sum(p)\). Note that the sums of both cumulative histograms are different. To make them equal, each element in \(th\_sum(p)\) must be multiplied with the ratio of the sums. Both cumulative histograms are now identical, except the distribution of pixel values.

# build target cumul. histogram th_sum
for i in range(0, 255):
	th_sum[i] = th_sum[i] + th[i]

# make sums of cumul. target histograms equal
ratio = h_sum[255] / th_sum[255]
for i in range(0, 255):
	th_sum[i] = th_sum[i] * ratio

4) Determine the look-up table (LUT) which lets \(h\_sum(p)\) track \(th\_sum(p)\) (the LUT is used in the next step as a transfer function for the pixel values of the original image). The process works as follows: a loop runs over every pixel value and calculates LUT(\(p\)), i.e. at every pixel value (0...255) the algorithm compares \(h\_sum(p)\) and \(th\_sum(p)\). If the value at the target histogram is larger that the original histogram LUT(\(p\)) we define a new pixel value \(p\_new\) (initially 0) that increases by 1. The value increases until the image histogram count matches the target histogram count. If the target histogram is smaller that the original image histogram \(p\_new\) stays the same. At the end of each iteration we set the LUT(\(v\)) equal to \(p\_new\).

# create LUT that makes the cumul. histogral track the target cumul. histogram
p_new = 0
for i in range(0, 255):
	while th_sum[p_new] < h_sum[i]
		p_new = p_new + 1
	LUT[i] = p_new	

5) Apply LUT to image. With the previous step done we have defined out transfer function that carries out the histogram specification. We apply the transfer function to the image at every pixel coordinate and obtain the new image.

# apply LUT
for y in range(0, ymax-1):
	for x in range(0, xmax-1):	
		image_new(x,y) = LUT(image(x,y))

The transformation of the original image simply applies the LUT to determine the target pixel value in the new image. The histograms of the newly created images are often characterized by odd-looking gaps between the histogram bars. This is due to the fact that a certain interval of intensities in the original have a limited number of gray levels. The algorithms that determines the LUT skips certain target pixel values until the target catches up with the original. 

Especially deep-sky images profit from histogram specification. The histogram of a typical deep-sky image, like in the example below, has a sharp peak at low intensity, corresponding to the dark background. The second component is a peak with exponential tail that describes the object itself. Some frequently used functions \(a\) in image processing that dictate the target histogram are for example:

- Flat-line / constant histogram, also called histogram equalization. The target histogram is a flat line, and the cumulative histogram consequently a wedge shape:

\(h(x) = c\),

where \(c\) is a constant. Because of its simplicity we list this function here, although it has limited use for deep-sky imagery. Many pixels that are originally dark or in the middle of the intensity scale are pushed to high intensities, letting the image look bright and washed-out. Images with a bi-modal histogram, like dark backgrounds and very bright objects (e.g. planets) are seemingly transformed into a over-exposed variant with gray backgrounds and extremely bright objects.

- Gaussian shape. This distribution assumes a variation of tones around a central mean. However, astronomical images often feature an exponential decline towards high intensities. Consequently, a Gaussian distribution that is shifted somewhat towards low intensities produces the most realistic and very effective images. The corresponding function is

\(h(x) = \frac{1}{\sigma \sqrt{2\pi}} e^{-\frac{1}{2}(\frac{x-\overline{x}}{\sigma})^2} \),

with the base of natural logarithms \(e\), the standard deviation \(\sigma\), and the mean value \(\bar{x}\) of \(x\).

- Exponential fall-off. Although many astronomical images appear like this histogram shape, they actually have low numbers of pixels at very low intensity values. Transforming the shape into a true exponential function lets many dark pixels appear even darker, if not black. However, images of bright extended objects display greater detail in the faint outer parts after applying this function, described by

\(h(x) = e^{-kx}\),

where \(k\) is a constant value.

Below we see the familiar image of the Horsehead Nebula. Apply the different versions of histogram specification and observe how the image and the corresponding histograms and cumulative histograms are changing.