Written September 8th, 2021

Images are, of course, sets of visual data. With any sort of processing, it is important to differentiate between analog and digital types. Analog images are real life things such as photographs, paintings, medical images, et cetera. So then, digital images are things like .pngs, .jpgs, and so on (while these are usually called "photos" it is important to state their explicit name so not to confuse them). We can explicitly represent images mathematically by a function I(x, y).

I(x, y) is a function in which *x* represents the horizontal spatial dimension and *y* represents the vertical spatial dimension. Thus, the units are some form of length. Videos would add a time component to this. The function I(x, y) is a continuous function that, for our purposes, is sampled at regular intervals to produce a digital image. This sampling is best represented as pixels per inch (PPI). A PPI of 2000 means that an analog image is sampled 2000 times in one inch. For a 1" by 1" square, this would mean there are 1000^{2} samples.

The equation above is the full representation of a digital image. The image consists of *M* rows and *N* columns where any arbitrary location can be found between 1 ≤ m ≤ M and 1 ≤ n ≤ N. The values *T _{1}* is the sampling distance (inverse of the sampling rate such as PPI) in the horizontal direction and

Monochromatic images are typically displayed in grayscale which ranges from white to black in intensity. A sampled grayscale image would be represented as a two-dimensional array of numbers *I[m, n]*. For most digital images, the value of *I[m, n]* at any given point is always positive and finite. E.g:

This bounding of the function is from the physical nature of the image; it doesn't make sense to have an infinitely intense physical value of light. One's monitor wouldn't be able to reproduce it nor would a computer be able to store it. The bound at *I _{max}* comes from the storage of the data on the computer. This allows the intensity of the image to be relative to a certain level. This is very useful for some digital signal processes that expand the range of an image's intensity, for example. Most digital images use an 8-bit intensity value. This means that

Take, for example, the image of the moon displayed below. Next to it is a histogram of how many pixels are at what intensity.

Notice all how the image has an even selection of intensity values - it doesn't favor being too bright or too dark so the image is balanced or equalized adequately well. If an image isn't balanced well, it's said to have a low contrast. Low contrast is usually hard for us with human eyes to see, and equally as difficult for digital image processors to make use of. While this image has rather well contrast width, it could be a bit better equalized due to the value of intensities at r = 0. Now we need to start thinking in terms of transformations to improve the contrast equalization.

The most generic and arbitrary transformation is a mapping function *z = T(r)* where *r* is the input intensity and *z* is the output intensity at any given pixel at location *[m, n]*. There is generally no restriction to the transformation *T* other than it must be defined within the working range of the intensities *r* and *z*. One very common type of transformation is a thresholding function where the graph of *z(r)* is split into two distinct regions of *r*. Two points are defined at (r1, z1) and (r2, z2). If z1 = 0, z2 = L - 1 (where L is 2^{# of bits}), and r1 = r2. This results in the shape of the graph given here (desmos link). What this function does is it ignores any data below a certain input intensity (T(r) = 0) and intensifies any data above a certain intensity (T(r) = L - 1).

Another type of contrast stretching is histogram equalization. I know MATLAB has this function already available, but it's worth going into what exactly histogram equalization does. The first new function I must introduce is the Cumulative Distribution Function (CDF). This function is a summation of the number of pixels at each intensity. The explicit definition of the CDF is displayed below. Here is also a link to the wikipedia article.

The CDF of an image always sums up to equal the number of pixels in the image. For a 3x3, 2-bit intensity image with 6 black (intensity is zero) and 3 white (intensity is 3) pixels, the CDF would be [6, 0, 0, 3]. Notice how the CDF is a length of four, which corresponds to the 2-bit resolution (2^{2} levels of intensity). With this function in mind, the CDF of the moon image is displayed below. One important thing to take notice of is the vertical axis of the graph - it doesn't start at r_{0} = 0 - but rather r_{0} ~= 1.2E5. This means that the lower intensity values are heavily favored in the image. We can clearly see this with the black background of the original moon.png image.

This CDF is our transformation *T(r)*. This can be explained conceptually. We want to even out the full range of an image so that each intensity has an equal part in the image intensity histogram. So, if there is a spike in low intensities (like our case), then we should prioritize the high intensities so that the high intensities have equal footing. We know our image has an equalized histogram when the image post-transformation results in a CDF that is representative of a linear function as close as possible. (An example of this is under "full-sized image" on the wikipedia page). After normalizing the CDF for the moon image, below is what the transformation did to the image along with the intensity histogram and the CDF itself.

Sadly, the image loses some contrast width when we put it through this function, but the function did *exactly* what we wanted it to do. It equalized the intensity histogram and made the CDF as linear as possible. (Also, the image did decide to die a bit as I put it through my processor, I wonder where this bug lies within MATLAB).