Ch-03 intensity transformations and spatial filtering

Chapter 3 …

What is this lesson about?

Intensity transformations operate on single pixels to effect contrast manipulation and image thresholding.
Spatial filtering works on a neighborhood centered around every pixel to perform operations such as image sharpening.
We study a number of intensity transformation and spatial filtering techniques in this lesson.

Spatial domain processing

Spatial domain techniques are computationally more efficient and require less processing resources to implement.
They can be denoted by the expression:
\[g(x,y) = T [f(x,y)]\]
For neighborhood operations, the border may need to be padded.
The neighborhood along with the operation is called a spatial filter, spatial mask, kernel, template, or window.

A $3 \times 3$ neighborhood about a point $(x,y)$

Intensity transformation

When the window is of size $1 \times 1$, $g$ depends on the value of $f$ at a single point $(x,y)$
$T$ becomes an intensity transformation function:
\[s = T(r)\]
where $s$ and $r$ represent the intensity of $g$ and $f$ at any point $(x,y)$.
Two transformations functions are shown below.
The first does contrast stretching.
The second is a limiting case of the first, and performs thresholding.
Intensity transformation functions

Intensity Transformation Functions

Intensity transformation functions are implemented as table lookups.
For 8-bit intensity images, size of lookup table is 256.
Some basic intensity transformation functions are shown figure below

Basic intensity transformation functions

Image negatives

The transformation function for photographic negative is:
\[s = (L - 1) - r\]
Shown below is original digital mammogram (showing a small lesion) and its negative.
Though the visual content is same in both figures, it is easier to analyze the breast tissue in the negative image.

Original mammogram image and its negative

Log Transformations

General form of a log transformation is:
\[s = c \log(1 + r)\]
where $c$ is a constant, and it is assumed that $r \ge 0$.
Log transformation maps a narrow range of low intensity values into a wider range of output levels.
The opposite is true of higher values of input levels.
An important characteristic of this transformation is that it compresses the dynamic range of images with large variations in pixel values.

Fourier spectrum and $\log$ transformation

Fourier spectrum values range from 0 to $10^6$ or higher.
On display systems, it is difficult to effectively display such a wide range of intensity values.
Significant degree of intensity detail can be lost.
For the left image shown below, intensity values range from 0 to $1.5 \times 10^6$.
Transforming this wide range of intensity values using the equation above will reduce the intensity values range to 0 to 6.2.
Scaling these intensity values linearly and displaying the spectrum on the same 8-bit display is shown in the right image below.

Fourier spectrum and its $\log$ transformation

Power-law ($\gamma$) transformations

Power-law ($\gamma$) transformations have the basic form
\[s = cr^\gamma\]
where $c$ and $\gamma$ are positive constants.
Sometimes the equation above is written as $s = c(r + \epsilon)^\gamma$ to account for measurable output when the input is zero.
Plots of $s$ versus $r$ for various values of $\gamma$ while keeping $c=1$ are shown below.
Power-law curves for fractional values of $\gamma$ behave like the $\log$ transformation (the opposite is true for higher values of $\gamma$).

Plots of the equation $s = cr^\gamma$ for various values of $\gamma$

Gamma $(\gamma)$ correction

Devices used for image capture, printing, and display respond according to a power law.
For example, CRT devices have an intensity-to-voltage response that is a power function, the exponent (i.e., $\gamma$) varies from 1.8 to 2.5.
As shown in above figure, such systems would produce images darker than intended.
Images shown below illustrates this effect.
$\gamma$ correction is effected by preprocessing the input image before inputting into the monitor by performing the following transformation:
\[s = r^{\frac{1}{2.5}} = r^{0.4}\]
The $\gamma$ value is device-dependent.

Intensity ramp image with $\gamma$ correction

Power-law transformations for contrast manipulation

In addition to $\gamma$ correction, power-law transformations are useful for contrast manipulation.
Shown in slide~\ref{sld:MRISpine} is an MRI image of a fractured human spine – results of applying $s = cr^\gamma$ with $c=1$ and $\gamma = 0.6, 0.4,$ and 0.3.
The second row, left image is the best and is produced with $\gamma = 0.6$. Further decrease in $\gamma$ value results in reduced contrast (washed-out appearance).

Contrast enhancement of an MRI image

Power-law transformations for contrast manipulation

The first row, left image shown below has a washed-out appearance.
Compression of intensity levels is desirable.
The remaining images in that slide show the results of applying equation with $\gamma = 3.0, 4.0, 5.0$ on the original image.
Results look good for $\gamma = 3.0 \text{ and } 4.0$; $\gamma = 5.0$ results in certain areas becoming quite dark.

An aerial image and $\gamma$ manipulation for improving contrast

Piecewise-Linear Transformation Functions

Contrast stretching

Low contrast images due to poor illumination, lack of dynamic range in the imaging sensor, or even the wrong setting of aperture.
Contrast stretching expands the range of intensity levels in an image so that the range spans the full intensity range of the display device or recording medium.
Shown in the top row, left image below is a typical transformation used for contrast stretching.
The original image, contrast-stretched image, and thresholded image are shown.

Contrast stretching and thresholding

Contrast stretching

Values of $(r_1,s_1)$ and $(r_2,s_2)$ control the shape of the transformation function.
$r_1 = s_1$ and $r_2 = s_2$ results in a linear function that produces no changes.
$r_1 = r_2, s_1 = 0$ and $s_2 = L - 1$ results in a thresholding function.
Intermediate values of $(r_1,s_1)$ and $(r_2,s_2)$ produces various degrees of spread in the intensity levels of output image.
In general, $r_1 \le r_2$ and $s_1 \le s_2$ so that the function is single valued and monotonically increasing.

Contrast stretching

The top row, right image is an scanning electron microscope image of pollen, magnified 700 times.
This is a low contrast image.
The bottom row, left image is contrast stretched with $(r_1,s_1) = (r_{min}, 0)$ and $(r_2,s_2) = (r_{max}, L - 1)$, where $r_{min}$ and $r_{max}$ are the minimum and maximum intensity values in the original input image.
The bottom row, right image is thresholded image with $(r_1, s_1) = (m, 0)$ and $(r_2, s_2) = (m, L - 1)$, where $m$ is the mean intensity level in the original input image.

Intensity-level slicing

Goal is to highlight specific range of intensities.
For example, enhancing features corresponding to masses of water in satellite imagery; enhancing flaws in X-ray images.
One approach is to set all values in the range of interest to just one value and all other values (i.e., those outside the range of interest) to another value (see left figure below).
In the second approach, brighten (or darken) the desired range of intensities and leave all other intensity values in tact (see right figure below).

Transformation functions for intensity-level slicing

Highlight major blood vessels

The left image below is an aortic angiogram near the kidney area.
The center image shows intensity-level slicing using a transformation of the form shown in left image above.
This type of enhancement produces a binary image and is useful for studying shape of flow of the contrast medium — detecting blockages, for example.
The right image below is the result of applying a transformation like the right image above.
A band of intensities in the mid-gray region around the mean intensity is set to black, and all others are left in tact.
Such an image might be useful in measuring the actual flow of the contrast medium as a function of time in a series of images.

Aortic angiogram and intensity-level transformations

Bit-plane representation of an 8-bit image

Bit-plane slicing

Shown below (top left) is an 8-bit grayscale image of size $500 \times 1192$ pixels.
Following figures are the bit-plane representations of bit planes 1 through 8 (plane 1 corresponds to the least significant bit).
Each bit plane is a binary image.
In terms of intensity transformation functions, $8^{th}$ bit plane of an 8-bit image is obtained by the following mappings:
- $[0..127] \rightarrow 0$ and $[128..255] \rightarrow 1$

An 8-bit image and its 8 bit planes (1 through 8)

Why bit-plane slicing?

Useful for determining relative importance of each bit.
This process aids in determining the adequacy of the number of bits used to quantize the image.
Bit-plane slicing is also useful in image compression.
For reconstruction, multiply the pixels in the $n^{th}$ plane by the constant $2^{n-1}$.
Add the selected plane values to construct the grayscale image.
As demonstrated below, storing the four highest-order bit planes allow us to reconstruct the original image in acceptable detail.

Reconstruction using bit planes

8, 7
8, 7, 6
8, 7, 6, 5

Histogram Processing

Intensity levels in the range $[0 \, .. \, L-1]$.
Histogram is a discrete function $h(r_k) = n_k$, where $r_k$ is the $k^{th}$ intensity value and $n_k$ is the number of pixels in the image with intensity $r_k$.
Normalized histogram is given by:
\[p(r_k) = \frac{n_k}{MN} \text{ for } k = 0, 1, 2, \ldots, L-1\]
$p(r_k)$ is the probability of occurrence of $r_k$ in an image.
The sum of all components of a normalized histogram is equal to 1.

Intensity distribution and histogram (a 3-bit, $64 \times 64$ image)

Image histograms

Image histograms are used for image enhancement, image compression, and segmentation.
Shown below are four figures and their histograms.
Histograms are plots of $h(r_k) = n_k \text{ versus } r_k$ or $p(r_k) = \frac{n_k}{MN} \text{ versus } r_k$.
An image whose pixels tend to occupy the entire range of possible intensity values, and tend to be distributed uniformly, will have an an appearance of high contrast and will exhibit a large variety of gray tones.

Dark, light, low-contrast, and high-contrast image types

Histogram Equalization

Histogram equalization

Histogram equalization automatically determines a transformation function which seeks to produce an output image that has a uniform histogram.
Sometimes you want to specify the exact shape of a histogram and like a transformation function automatically generated – histogram specification, histogram matching.
A function $T(r)$ is monotonically increasing if $T(r_2) \ge T(r_1) \text{ for } r_2 > r_1$.
A function $T(r)$ is strictly monotonically increasing if $T(r_2) > T(r_1) \text{ for } r_2 > r_1$.
Similar definitions apply to monotonically decreasing functions.
Note: Histogram equalization: ideally we want strictly monotonically increasing function transformation; due to discretization, we will get transformations given by the monotonically increasing function only.
Figure: Monotonically increasing functions

Assume that the intensity values are continuous, and consider transformations of the form: $s = T(r), \,\, 0 \le r \le L - 1$
Assume that $T(r)$ is:
- ( Condition 1) Monotonically increasing in the interval $0 \le r \le L - 1$; and
- ( Condition 2) $0 \le T(r) \le L - 1 \text{ for } 0 \le r \le L - 1$.
Condition 1 insures that for any two input values $r_1$ and $r_2$ such that $r2 \ge r_1$, the corresponding output values will have the property $s_2 \ge s_1$ — prevents artifacts from being created by reversals intensity.
Condition 2 insures that the range of output intensities is the same as input.
In some contexts, we use the inverse $r = T^{-1}(s), \quad 0 \le s \le L - 1$
In this case $T(r)$ is strictly monotonically increasing function in the interval $0 \le r \le L - 1$ (guarantees that the mapping from $s$ back to $r$ will be one-to-one).

Probability density function (PDF)

PDF is a fundamental descriptor of a random variable.
$p_r(r)$ and $p_s(s)$ denote PDFs of $r$ and $s$.
If $p_r(r)$ and $T(r)$ are known, and $T(r)$ is continuous and differentiable over a range of values of interest, then the PDF of $s$ (from probability distribution theory) is:
\[p_s(s) = p_r(r)\left|\frac{dr}{ds}\right|\]
A transformation function of interest in image processing has the form
\[s = T(r) = L - 1 \int_0^r p_r(w) dw \label{eqn:CDFX}\]
where $w$ is a dummy variable of integration.
The right side of equation above is the cumulative distribution function (CDF) of random variable $r$.
Note: $p_r(r)$ in the last line is the original histogram.

Cumulative distribution function (CDF)

Equation above is monotonically increasing (since the area under the function cannot decrease as $r$ increases) — Condition 1 is satisfied.
The maximum value of $s$ is $L-1$, so Condition 2 is also satisfied.
Use equation above to find $p_s(s)$.
\[\begin{aligned} \frac{d s}{d r} &=\frac{d T(r)}{d r} \\ &=(L-1) \frac{d}{d r}\left[\int_{0}^{r} p_{r}(w) d w\right] \\ &=(L-1) p_{r}(r) \end{aligned}\]

CDF

Substituting $\frac{ds}{dr}$ value in equation yields

\[\begin{aligned} p_{s}(s) &=p_{r}(r)\left|\frac{d r}{d s}\right| \\ &=p_{r}(r)\left|\frac{1}{(L-1) p_{r}(r)}\right| \\ &=\frac{1}{(L-1)}, \quad 0 \leq s \leq L-1 \end{aligned}\]

Recognize that $p_s(s)$ is an uniform probability density function.
Performing the intensity transformation in equation yields a random variable, $s$, characterized by a uniform PDF.
Though $T(r)$ depends on $p_r(r)$, the resulting $p_s(s)$ always is uniform, independently of the form of $p_r(r)$ as shown below.

Applying $s = T(r) = (L-1) \int_0^r p_r(w) \, dw$

Example to illustrate equations

Assume that the (continuous) intensity values in an image have the PDF
\[p_r(r) = \begin{cases} \frac{2r}{(L-1)^2} & \text{ for } 0 \le r \le L - 1 \\ 0 & \text{ otherwise} \end{cases}\]
Using equation
\[\begin{aligned} s=T(r) &=(L-1) \int_{0}^{r} p_{r}(w) d w \\ &=\frac{2}{L-1} \int_{0}^{r} w d w=\frac{r^{2}}{L-1} \end{aligned}\]
Consider an image in which $L=10$ and a pixel at an arbitrary location $(x,y)$ has intensity $r=3$.
This pixel value in the histogram equalized image is:
\[s = T(r) = \frac{r^2}{9} = 1\]
Compute $p_s(s)$ using equation
\[\begin{aligned} p_{s}(s)=p_{r}(r)\left|\frac{d r}{d s}\right| &=\frac{2 r}{(L-1)^{2}}\left|\left[\frac{d s}{d r}\right]^{-1}\right| \\ &=\frac{2 r}{(L-1)^{2}}\left|\left[\frac{d}{d r} \frac{r^{2}}{L-1}\right]^{-1}\right| \\ &=\frac{2 r}{(L-1)^{2}}\left|\frac{(L-1)}{2 r}\right|=\frac{1}{L-1} \end{aligned}\]

Discrete form of intensity transformations

For discrete values, we use probabilities (histogram values) and summations instead of PDFs and integrals.
As discussed earlier, the probability of occurrence of intensity level $r_k$ is approximated by:
\[p_r(r_k) = \frac{n_k}{MN} \quad k = 0, 1, 2, \ldots, L-1\]

Discrete form of intensity transformations

The discrete form of equation is:
\[\begin{aligned} s_{k}=T\left(r_{k}\right) &=(L-1) \sum_{j=0}^{k} p_{r}\left(r_{j}\right) \\ &=\frac{(L-1)}{M N} \sum_{j=0}^{k} n_{j} \quad k=0,1,2, \ldots, L-1 \end{aligned}\]
The transformation specified by equation above is called histogram equalization or histogram linearization.

Example to illustrate histogram equalization

Consider a hypothetical 3-bit image $(L=8)$ of size $64 \times 64$ pixels $(MN = 4096)$.
Range of intensity levels is $[0, L-1] = [0,7]$.
Shown below are intensity distribution, histogram, histogram equalization function, and equalized histogram.

Intensity distribution and histogram values

Example to illustrate histogram equalization

Calculate $s_1, s_2, \cdots, s_7$ using equation.
For example, $s_0 = T(r_0) = 7 \displaystyle\sum_{j=0}^{0} p_r(r_j) = 7 p_r(r_0) = 7 \times 0.19 = 1.33$.
Likewise, $s_1 = 3.08, s_2 = 4.55, s_3 = 5.67, s_4 = 6.23, s_5 = 6.65, s_6 = 6.86, \text{ and } s_7 = 7.00$

Example to illustrate histogram equalization

Round $s$ values to the nearest integer.
- $s_0 = 1.33 \rightarrow 1 \quad s_1 = 3.08 \rightarrow 3$
- $s_2 = 4.55 \rightarrow 5 \quad s_3 = 5.67 \rightarrow 6$
- $s_4 = 6.23 \rightarrow 6 \quad s_5 = 6.65 \rightarrow 7$
- $s_6 = 6.86 \rightarrow 7 \quad s_7 = 7.00 \rightarrow 7$
Equalized histogram is shown in the right side figure above.
Unlike the continuous case, it cannot be proved (in general) that discrete histogram equalization results in a uniform histogram.

Histogram equalized images

Histogram Matching (Specification)

Histogram equalization based on a specified histogram

Original histogram, specified histogram, transformation function, and resulting histogram.

Histogram matching (specification)

To generate a processed image that has a specified histogram.
Equalize the input image $r$: $s = T(r) = \int_0^r p_r(w) \, dw$
If the desired image $z$ were available, it could also be equalized: $s^\prime = G(z) = \int_0^r p_z(w) \, dw$, where $p_z(w)$ is the desired histogram.
The inverse of the above transform is: $z = G^{-1}(s^\prime) = G^{-1}(s)$
Since $s$ and $z$ have the same histogram, $z = G^{-1}(s) = G^{-1}(T(r))$
Both $T(r)$ and $G(z)$ can be found from the histogram of the input image and the specified histogram, respectively.
For continuous intensities $r$ (input) and $z$ (output), let $p_r(r)$ and $p_z(z)$ denote their PDFs.
You can estimate $p_r(r)$ from the input image, while $p_z(z)$ is the specified PDF that we wish the output image to have.
Let $s$ be a random variable with the property
\[s = T(r) = (L - 1) \int_0^r p_r(w) \, dw\]
Let $z$ be a random variable with the property
\[G(z) = (L - 1) \int_0^z p_z(t) \, dt = s\]
It follows from equations that $G(z) = T(r)$.
Then $z$ must satisfy the condition
\[z = G^{-1}[T(r)] = G^{-1}(s)\]

Histogram matching procedure

Obtain $p_r(r)$ from the input image and use equation~\eqref{eqn:CDFR} to obtain values of $s$.
Use the specified PDF in equation~\eqref{eqn:CDFR2} to obtain $G(z)$.
Obtain $z = G^{-1}(s)$.
Obtain the output image by first equalizing the input image using equation~\eqref{eqn:CDFR}; pixel values in this image are the $s$ values. For each $s$ compute $z$ using the inverse mapping $z = G^{-1}(s)$.

Example to illustrate histogram specification

An image has intensity PDF
\[p_r(r) = \frac{2r}{(L-1)^2} \text{ for } 0 \le r \le (L-1)\]
and $p_r(r) = 0$ for other values of $r$.
Find a transformation function that will produce an image whose intensity PDF is
\[p_z(z) = \frac{2z^2}{(L-1)^3} \text{ for } 0 \le z \le (L-1)\]
and $p_z(z) = 0$ for other values of $z$.
First find histogram equalization transformation for the interval $[0..L-1]$.
\[s = T(r) = (L-1)\int_{0}^{r} p_r(w) \, dw = \frac{2}{(L-1)} \int_0^r w \, dw = \frac{r^2}{(L-1)}\]
Next, find $G(z)$
\[G(z) = (L-1)\int_{0}^{z} p_z(w) \, dw = \frac{3}{(L-1)^2} \int_0^z w^2 \, dw = \frac{z^3}{(L-1)^2}\]
over the interval $[0..L-1]$; this function is 0 elsewhere by definition.
Finally, we require that $G(z) = s$. Therefore, $z = [(L-1)^2s]^\frac{1}{3}$.
Substituting for $s$ yields:
\[z = \left[(L-1)^2 s\right]^\frac{1}{3} = \left[(L-1)^2 \frac{r^2}{(L-1)}\right]^\frac{1}{3} = [(L-1)r^2]^\frac{1}{3}\]
The two steps are combined into a single transformation from $r$ to $z$.

Histogram specification

In practice, finding meaningful analytical expressions for $T(r)$ and $G^{-1}$ are difficult for the continuous case.
However, the problem is simplified for the discrete case.
The price paid is that only approximation to the desired histogram is achieved.

Histogram equalization: discrete formulation

Recall the discrete formulation of histogram equalization discussed earlier:
\[\begin{aligned} s_{k}=T\left(r_{k}\right) &=(L-1) \sum_{j=0}^{k} p_{r}\left(r_{j}\right) \\ &=\frac{(L-1)}{M N} \sum_{j=0}^{k} n_{j} \quad k=0,1,2, \cdots, L-1 \end{aligned}\]
Given a specific value of $s_k$, the discrete formulation of equation~\eqref{eqn:CDFR2} involves computing the transformation function
\[G(z_q) = (L-1)\sum_{i=0}^{q}p_z(z_i)\]
for a value of $q$ such that $G(z_q) = s_k$, where $p_z(z_i)$ is the $i^{th}$ value of the specified histogram.
Find $z_q$ by obtaining the inverse transformation, which is a mapping from $s$ to $z$.
\[z_q = G^{-1}(s_k)\]

Summary of histogram specification procedure

Compute the histogram of $p_r(r)$ of the given image, and use it to find histogram equalization transformation in equation~\eqref{eqn:HistEq2}.

Round the resulting values ($s_k$) to the integer range $[0 .. L-1]$.
Compute all values of $G$ using the equation~\eqref{eqn:Gz} for $q = 0, 1, 2, \cdots, L-1$, where $p_z(z)$ are the values of the specified histogram.

Round $G$ values to integers in the range $[0 .. L-1]$ and store them in a table.
For each $s_k, k = 0, 1, 2, \cdots, L-1$, use the stored value of $G$ from step~\ref{stepX} to find the corresponding value of $z_q$ so that $G(z_q)$ is closest to $s_k$ and store these mappings from $s$ to $z$.

When there is more than one $z_q$ satisfies the given $s_k$, choose the smallest $z_q$ (convention).
Form the histogram-specified image by mapping the pixels of histogram equalized input image to $z_q$ pixel values in the histogram-specified image using the mappings of step~\ref{stepY}.

Example to illustrate histogram specification

Consider the $64 \times 64$ hypothetical image introduced in slide~\ref{sld:IntDist}.
Its histogram is shown in top left in slide~\ref{sld:HEExample2}.
You like to transform this histogram so that it will have the values as specified in the center column of the figure in slide~\ref{sld:HEExample3}.
The top right area of slide~\ref{sld:HEExample2} shows this histogram.

Histogram equalization based on a specified histogram

Original histogram, specified histogram, transformation function, and resulting histogram.

Specified and actual histograms

Example to illustrate histogram specification

First obtain the $s_k$ values as illustrated in the example beginning on slide~\ref{sld:HEExample}.
\[s_0 = 1, s_1 = 3, s_2 = 5, s_3 = 6, s_4 = 6, s_5 = 7, s_6 = 7, \text{ and } s_7 = 7\]
Compute all values of $G$ using equation~\eqref{eqn:Gz}.
\[\begin{array}{l} G\left(z_{0}\right)=7 \sum_{j=0}^{0} p_{z}\left(z_{j}\right)=0.00 \\ G\left(z_{1}\right)=7 \sum_{j=0}^{1} p_{z}\left(z_{j}\right)=7\left[p_{z}\left(z_{0}\right)+p_{z}\left(z_{1}\right)\right]=0.00 \\ G\left(z_{2}\right)=7 \sum_{j=0}^{2} p_{z}\left(z_{j}\right)=7\left[p_{z}\left(z_{0}\right)+p_{z}\left(z_{1}\right)+p_{z}\left(z_{2}\right)\right]=0.00 \end{array}\]
- $G(z_0) = 0.00 \rightarrow 0 \quad G(z_1) = 0.00 \rightarrow 0$
- $G(z_2) = 0.00 \rightarrow 0 \quad G(z_3) = 1.05 \rightarrow 1$
- $G(z_4) = 2.45 \rightarrow 2 \quad G(z_5) = 4.55 \rightarrow 5$
- $G(z_6) = 5.95 \rightarrow 6 \quad G(z_7) = 7.00 \rightarrow 7$
These values are summarized in slide~\ref{sld:tab3} and are sketched in bottom left area of slide~\ref{sld:HEExample2}.
$G$ is not strictly monotonic, use step~\ref{stepY} of the algorithm to handle this situation.

Example to illustrate histogram specification: $G$ values

Example to illustrate histogram specification

For every value of $s_k$, find smallest $z_q$ so that the value $G(z_q)$ is the closest to $s_k$.
For $s_0 = 1$, $G(z_3) = 1$ is perfect match. Therefore, we have the correspondence $s_0 \rightarrow z_3$
These mappings are shown in slide~\ref{sld:tab4}.
The third column of the table in slide~\ref{sld:HEExample3} lists the values of the resulting histogram.
The resulting histogram is sketched in the bottom right region of slide~\ref{sld:HEExample2}.

Example to illustrate histogram specification $s_k \rightarrow z_q$

Mars moon: Phobos

Shown on left in slide~\ref{sld:MarsMoon} is an image of Mars moon (Phobos) taken by NASA’s Mars Global surveyor.
Shown on the right is its histogram. The image is dominated by large dark areas.
The result of applying the transformation function of equation~\eqref{eqn:CDF} to the Phobos image is shown in slide~\ref{sld:MarsMoon2}.
Notice the shape of the transformation function and resulting washed-out image.
As shown by the histogram, all the intensity levels are biased toward the upper one-half of the gray scale.

Mars moon and its histogram

Histogram equalized Mars moon image

Mars moon: Phobos

Shown in the top left of slide~\ref{sld:MarsMoon3} is the manually specified function.
Sampling this function into 256 equally spaced discrete values produces the desired specified histogram.
Obtain $G(z)$ using equation~\eqref{eqn:Gz}, which is shown by the transformation labeled (1) in the center image in the left column of slide~\ref{sld:MarsMoon3}.
Obtain $G^{-1}(s)$ using equation~\eqref{eqn:Gq}, which is shown by the transformation labeled (2) in the center image in the left column of slide~\ref{sld:MarsMoon3}.

Mars moon: Phobos (enhanced image using mappings from curve (2))

Local Histogram Processing

Define a neighborhood and its center from pixel to pixel.
At each location, compute the histogram of the pixels in the neighborhood.
Obtain either histogram equalization or histogram specification transformation function.
Use this function to map the intensity of the pixel centered in the neighborhood.
Repeat this process.

Local histogram processing example

Shown in the left in slide~\ref{sld:Global} is an 8-bit, $512 \times 512$ image.
The center image shows global histogram equalization. Gained significant enhancement of noise.
Shown right in slide~\ref{sld:Global} is the image obtained using local histogram equalization with a neighborhood size of $3 \times 3$.
Notice the significant detail within the dark squares.

Global and local histogram equalizations

Image histogram statistics

Intensity levels range: $[0 .. L-1]$
$p(r_i)$ denotes normalized histogram component corresponding to $r_i$.
View $p(r_i)$ as an estimate of the probability of $r_i$ occurring in the image.
The $n^{th}$ moment of $r$ about its mean $m$

\[\begin{aligned} \mu_{n}(r) &=\sum_{i=0}^{L-1}\left(r_{i}-m\right)^{n} p\left(r_{i}\right) \\ m &=\sum_{i=0}^{L-1} r_{i} p\left(r_{i}\right) \end{aligned}\]

The $2^{nd}$ moment (also denoted by $\sigma^2$) of $r$ about its mean $m$

\[\mu_{2}(r)=\sum_{i=0}^{L-1}\left(r_{i}-m\right)^{2} p\left(r_{i}\right)\]

Mean ($m$) is a measure of average intensity, where as the variance (or $\sigma$) is a measure of contrast in the image.
If only mean and variance are desired, compute them directly without computing the histogram.

\[\begin{aligned} m &=\frac{1}{M N} \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x, y) \\ \sigma^{2} &=\frac{1}{M N} \sum_{x=0}^{M-1} \sum_{y=0}^{N-1}\left[(f(x, y)-m]^{2}\right. \end{aligned}\]

Image histogram statistics example

Consider the following 2-bit image of size $5 \times 5$

\[\left(\begin{array}{lllll} 0 & 0 & 1 & 1 & 2 \\ 1 & 2 & 3 & 0 & 1 \\ 3 & 3 & 2 & 2 & 0 \\ 2 & 3 & 1 & 0 & 0 \\ 1 & 1 & 3 & 2 & 2 \end{array}\right)\]

L = $2^2 = 4$, and intensity range is $[0..L-1] = [0..3]$
Compute $p(r_i): p(r_0) = \frac{6}{25} = 0.24, p(r_1) = \frac{7}{25} = 0.28, p(r_2) = \frac{7}{25} = 0.28, p(r_3) = \frac{5}{25} = 0.20$
Compute $m$ using equation~\eqref{eqn:CM} $\begin{aligned} m &=\sum_{i=0}^{3} r_{i} p\left(r_{i}\right) \\ &=(0)(0.24)+(1)(0.28)+(2)(0.28)+(3)(0.20) \\ &=1.44 \end{aligned}$
Compute $m$ using equation~\eqref{eqn:meanVal}

\[\begin{aligned} m &=\frac{1}{25} \sum_{x=0}^{4} \sum_{y=0}^{4} f(x, y) \\ &=1.44 \end{aligned}\]

Global (and local) mean and variance

Global mean and variance used for gross adjustments in overall intensity and contrast.
More powerful use is in local enhancement using local mean and variance.
Let $S_{xy}$ denote a neighborhood centered on $(x,y)$.
The mean value of pixels in the neighborhood is:

\[m_{S_{x y}}=\sum_{i=0}^{L-1} r_{i} p_{S_{x y}}\left(r_{i}\right)\]

where $p_{S_{xy}}$ is the histogram of the pixel in region $S_{xy}$.
Variance of the pixels in the neighborhood is

\[\sigma_{S_{x y}}^{2}=\sum_{i=0}^{L-1}\left(r_{i}-m_{S_{x y}}\right)^{2} p_{S_{x y}}\left(r_{i}\right)\]

Local enhancement using histogram statistics

Shown in slide~\ref{sld:SEMTungsten} is an SEM of tungsten filament wrapped around a support.
The filament in the center of the image is quite clear, but there is another one on the right side which is barely perceptible.
Enhance the dark area and leave the light areas unchanged — ideal situation for local enhancement.

Global and local histogram equalization

SEM image of a magnified tungsten filament, and results of global and local histogram equalization.

Local enhancement using histogram statistics

Let $m_G$ and $m_{S_{xy}}$ denote global and local mean intensities.
Compare $m_{S_{xy}}$ to $m_G$ to determine whether an area is relatively light or dark at point $(x,y)$.
Process pixel at $(x,y)$ if $m_{S_{xy}} \le k_0 m_G$ where $k_0$ is positive constant with value less than 1.0
A pixel at $(x,y)$ is a candidate for enhancement if $\sigma_{S_{xy}} \le k_2 \, \sigma_G$ where $\sigma_G$ is global standard deviation and $k_2$ is a positive constant.
$k > 1.0$ for enhancing light areas and $k < 1.0$ for dark areas.
We need to restrict the lowest values of contrast we are willing to accept.
Set a lower limit on the local standard deviation: $k_1\sigma_G \le \sigma_{S_{xy}}$ with $k_1 < k_2$.
Pixels that qualify for local enhancement are processed by multiplying the pixel value by a specified constant $E$.

\[g(x, y)=\left\{\begin{array}{ll} E \cdot f(x, y) & \text { if } m_{S_{x y}} \leq k_{0} m_{G} \text { and } k_{1} \sigma_{G} \leq \sigma_{S_{x y}} \leq k_{2} \sigma_{G} \\ f(x, y) & \text { otherwise } \end{array}\right.\]

Choosing the parameter values in equation~\eqref{eqn:Sum} requires a bit of experimentation.
For the figures in slide~\ref{sld:SEMTungsten}, $E = 4.0, k_0 = 0.4, k_1 = 0.02, \text{ and }, k_2 = 0.4$
$E = 4.0$ was chosen to preserve the general visual balance of the image (dark areas should be left relatively dark).
$S_{xy}$ should be as small as possible, so we chose $3 \times 3$.

Spatial Filtering

Linear Spatial filtering

A spatial filter consists of
- A neighborhood (typically a small rectangle)
- A predefined operation (performed on the pixels encompassed by the neighborhood)
If the operation performed on the pixels is linear, the filter is called a linear spatial filter; otherwise, it is a nonlinear spatial filter.

\[\begin{aligned} g(x, y)=& w(-1,-1) f(x-1, y-1)+w(-1,0) f(x-1, y)+\cdots \\ &+w(0,0) f(x, y)+\cdots+w(1,1) f(x+1, y+1) \end{aligned}\]

The center coefficient of the filter, $w(0,0)$, aligns with the pixel at location $(x,y)$.
Mask of size $m \times n$, assume that $m=2a+1$ and $n=2b+1$, where $a$ and $b$ are positive integers.

\[g(x, y)=\sum_{s=-a}^{a} \sum_{t=-b}^{b} w(s, t) f(x+s, y+t)\]

Spatial Correlation and Convolution

Correlation: process of moving a filter mask over the image and computing the sum of products at each location.
Convolution: same as correlation except that the filter is first rotated by $180^\circ$.

Correlation is a function of the displacement of the filter.
The first value of correlation corresponds to zero displacement of the filter, the second value of correlation corresponds to one unit displacement of the filter, and so on.
Correlating a filter $w$ with a function that contains all zeros and a single one ( discrete unit impulse function) yields a result that is a $180^\circ$ rotated copy of $w$.a
A fundamental property of convolution is that convolving a function with a unit impulse yields a copy of the function at the location of the impulse.
As shown in slide~\ref{sld:ImCorr}, correlation and convolution concepts apply to images as well.
If the filter mask is symmetric, both correlation and convolution yield the same result.

Correlation of a filter $w(x,y)$ of size $m \times n$ with an image $f(x,y)$, denoted $w(x,y) **\hollowstar** f(x,y)$:

\[w(x,y) **\hollowstar** f(x,y) = \sum_{s=-a}^{a} \sum_{t=-b}^{b} w(s,t) f(x+s,y+t)\]

where $a=(m-1)/2$, $b=(n-1)/2$, and $m$ and $n$ are odd integers.
Similarly, convolution of $w(x,y)$ and $f(x,y)$, denoted $w(x,y) **\star** f(x,y)$

\[w(x, y) \star f(x, y)=\sum_{s=-a}^{a} \sum_{t=-b}^{b} w(s, t) f(x-s, y-t)\]

Convolution is commutative, whereas correlation is not.

Vector representation of linear filtering

Convenient to write the sum of products as $R &= w_1z_1 + w_2z_2 + \ldots + w_{mn}z_{mn} \\ &= \sum_{k=1}^{mn} w_kz_k \\ &= \mathbf{w}^T\mathbf{z}$
For the filter mask shown in slide~\ref{sld:VectRep}

\[R &= w_1z_1 + w_2z_2 + \cdots + w_{9}z_{9} \\ &= \sum_{k=1}^{9} w_kz_k \\ &= \mathbf{w}^T\mathbf{z}\]

Generating spatial filter masks

We want to replace the pixels in an image by the average intensity of a $3 \times 3$ neighborhood centered on those pixels.

\[R = \frac{1}{9}\sum_{i=1}^{9} z_i\]

Note that $w_i$ are equal to $\frac{1}{9}$ and $z$ values are intensities.
Gaussian function of two variables has the basic form:

\[h(x,y) = e^{-\frac{x^2+y^2}{2\sigma^2}}\]

To generate a $3 \times 3$ mask from the Gaussian function, sample it about its center.
- \[w_1 = h(-1,-1), w_2 = h(-1,0), \cdots, w_9 = h(1,1)\]

Smoothing spatial filters

Used for blurring and noise reduction.
Response of a smoothing, linear spatial filter is simply the average of pixels contained in the neighborhood of the filter mask.
They are also called lowpass filters and averaging filters.
Effects reduced sharp transitions in intensities.
Random noise typically consists of sharp transitions, smoothing is used as a noise reduction technique.
False contours due to insufficient number of intensity levels can be smoothed as well.
However, edges are also characterized by sharp intensity transitions. Smoothing blurs edges — an undesirable effect.
Shown in slide~\ref{sld:SFilters} are two $3 \times 3$ smoothing filters.

\[R = \frac{1}{9}\sum_{i=1}^{9} z_i\]

Coefficients of the first filter are all 1s — computationally more efficient.
A spatial averaging filter in which all the coefficients are equal is called a box filter.

The second filter is called a weighted average filter.
The center pixel is weighted more, and other pixels are inversely weighted as a function of their distance from the center pixel.
The idea is to reduce blurring in the smoothing process.
Filtering an image of size $M \times N$ with a weighted averaging filter of size $m \times n$ ($m$ and $n$ are odd) is given by

\[g(x,y) = \frac{\displaystyle\sum_{s=-a}^{a} \displaystyle\sum_{t=-b}^{b} w(s,t) f(x+s,y+t)}{\displaystyle\sum_{s=-a}^{a} \displaystyle\sum_{t=-b}^{b} w(s,t)}\]

The denominator needs to be computed only once.

Smoothing with different filter sizes: example

Shown in slide~\ref{sld:DSFM} are the effects of smoothing an image with different filter sizes.
The original image (top, first column) is of size $500 \times 500$ pixels.
The figures resulting from smoothing with square filters of size $m=$ 3, 5, 9, 15, and 35 are shown (left to right, top to bottom order).
Blurring with filter sizes 15 and 35 results in aggressive blurring and is used to eliminate small objects from an image.
Notice the gradual but pronounced appearance of the black border due to padding of the image border with 0s.

Smoothing followed by thresholding: example

Shown on the left in slide~\ref{sld:Hubble1} is an image from Hubble telescope in orbit around the Earth.
The center image shows the result of applying a $15 \times 15$ averaging mask.
This smoothing has effected a number of objects to blend with the background or diminished their intensity significantly.
Thresholding the center image with a threshold value equal to 25\% of the highest intensity (in the blurred image) is shown on the right.
The result is a reasonable representation of the largest and brightest objects in the original image.

Order-statistic (nonlinear) filters

Response of nonlinear filters involve ordering (or ranking) the pixels in the neighborhood encompassed by the filter.
Best known filter in this category is the median filter.
For certain types of random noise, median filters provide excellent noise reduction with considerably less blurring than linear smoothing filters of similar size.
Median filters are particularly effective in the presence of impulse noise (aka salt-and-pepper noise).

Median filter**

The median, $\xi$, of a set of values is such that half the values in the set are less than or equal to $\xi$, and half are greater than or equal to $\xi$.
First sort the pixel values in the neighborhood encompassed by the filter window.
Once sorted, for a $3 \times 3$ window the median is the $5^{th}$ largest value, and for a $5 \times 5$ window the median is the $13^{th}$ largest value.
Median filters force pixels with distinct intensity values to be more like their neighbors.
Isolated clusters of pixels which are light or dark (wrt their neighbors) and whose area is less than $m^2/2$ are eliminated by $m \times m$ median filter.
Median filter is by far the most useful order-statistic filter.
The median represents the $50^{th}$ percentile of a ranked set of numbers.
Using the $100^{th}$ percentile results in the so-called max filter — useful for finding the brightest spots in the image.
$3 \times 3$ filter’s response is

\[R = \max\{z_k \mid k = 1, 2, \cdots, 9\}\]

The $0^{th}$ percentile filter is the min filter, which has the opposite effect of \emph{max filter}.

Median filter example

Shown in left in slide~\ref{sld:MedFilter} is an X-ray image of a circuit board corrupted by salt-and-pepper noise.
Shown in the center is the result of processing the corrupted image with a $3 \times 3$ neighborhood averaging mask.
The averaging filter blurred the image and noise reduction is poor.
Shown in right is the result of processing the corrupted image with a $3 \times 3$ median filter.
Compare and contrast the effectiveness of these two filters in removing salt-and-pepper noise.

Sharpening Spatial Filters

Sharpening is to highlight transitions in intensity.
Applications range from electronic printing and medical imaging to industrial inspection and autonomous guidance in military systems.
Averaging is analogous to integration, sharpening is spatial differentiation.
Strength of the response of a derivative operator is proportional to the degree of intensity discontinuity.
Image differentiation enhances edges and other discontinuities such as noise, and deemphasizes areas of slowly varying intensities.

Derivatives of a digital function

Defined in terms of differences.
Properties of first derivative:
- Must be zero in areas of constant intensity,
- Must be non-zero at the onset an intensity step or ramp, and
- Must be non-zero along ramps.
Properties of second derivative:
- Must be zero in areas of constant intensity,
- Must be non-zero at the onset and end of an intensity step or ramp, and
- Must be zero along ramps of constant slope.
The first-order derivative of a 1-D function $f(x)$ is the difference:

\[\frac{\partial f}{\partial x} = f(x+1) - f(x)\]

The second-order derivative of a 1-D function $f(x)$ is the difference:

\[\frac{\partial^2 f}{\partial^2 x} = f(x+1) + f(x-1) - 2 f(x)\]

These two definitions satisfy the conditions stated in slide~\ref{sld:DerCond}.
Similarities and difference between the first- and second-order derivatives are illustrated in slide~\ref{sld:FSDer}.

First- and second-order derivatives

Observe that the sign of a second derivative changes sign at the onset and end of a step or ramp.
In a step transition, the line joining these two values crosses the horizontal axis midway between the two extremes — zero crossing property.
First derivative gives thick edges since the derivative is non-zero along a ramp.
On the other hand, the second derivative would produce a double edge of one pixel thick, separated by zeros.
Second derivatives are much easier to implement than the first derivatives.

The Laplacian: a second derivativ**

Isotropic filters are rotation invariant.
Rotating the image and then applying the filter will give the same result as applying the filter first and then rotating the result.
The simplest isotropic derivative operator is the Laplacian.

\[\nabla^2f = \frac{\partial^2f}{\partial x^2} + \frac{\partial^2f}{\partial y^2}\]

As derivatives of any order are linear, the Laplacian is a linear operator.
Discrete form of equation~\eqref{eqn:TheLap} in $x-$ and $y-$directions

\[\frac{\partial^2f}{\partial x^2} &= f(x+1, y) + f(x-1,y) - 2f(x,y) \label{eqn:nablaOne} \\ \frac{\partial^2f}{\partial y^2} &= f(x, y+1) + f(x,y-1) - 2f(x,y) \label{eqn:nablaTwo}\]

Discrete Laplacian of two variables is

\[\nabla^2f(x,y) = f(x+1, y) + f(x-1,y) + f(x, y+1) + f(x,y-1) - 4f(x,y)\]

Implementing the discrete Laplacian

Equation~\eqref{eqn:TheLapDis} can be implemented using the top, left mask shown in slide~\ref{sld:LapFilters}. Yields isotropic result for rotations in increments of $90^\circ$.
Diagonal directions can be incorporated by adding two more terms to Equation~\eqref{eqn:TheLapDis} (note that each diagonal term also contains a $-2f(x,y)$).
The total subtracted from the difference terms now would be $-8f(x,y)$.
The top, right mask shown in slide~\ref{sld:LapFilters} incorporates diagonal directions as well. Yields isotropic result for rotations in increments of $45^\circ$.

In practice, you are likely to see the masks shown in the bottom row in slide~\ref{sld:LapFilters}.
They are the negatives of the discrete versions of the ones specified in equations~\eqref{eqn:nablaOne} and~\eqref{eqn:nablaTwo}.
The Laplacian highlights intensity discontinuities and deemphasizes regions with slowly varying intensity levels.
Produces figures that have grayish edge lines and other discontinuities superimposed on a dark, featureless background.
Recover background features by adding the Laplacian image to the original.

\[g(x,y) = f(x,y) + c [\nabla^2f(x,y)]\]

where $c = -1$ for the top row of the Lapacian filters in slide~\ref{sld:LapFilters} and $c = 1$ for the filters in the bottom row.

Image sharpening using the Laplacian: example

The image in the top row in slide~\ref{sld:LapSharpening} is a slightly blurred image of the north pole of the moon.
Shown in the left in the middle row is the result of filtering the image with the mask shown in the top, left in slide~\ref{sld:LapFilters}.
Since the Laplacian contains both positive and negative values, all negative values are clipped to 0 by the display — large sections of this image are black.
Bring the minimum value in this image to zero and scale the result to the full intensity range — right image in the center row in slide~\ref{sld:LapSharpening}.
The left image in the bottom row in slide~\ref{sld:LapSharpening} is the result of applying equation~\eqref{eqn:TheLapDis2} with $c=-1$.
The right image in the bottom row is the result of repeating the above procedure with the directional filter (top row, right filter of slide~\ref{sld:LapFilters}).

Unsharp masking and highboost filtering

Idea: sharpen an image by subtracting an unsharp (smoothed) version of an image from the original.
Let $f(x,y)$ and $\overline{f}(x,y)$ be the original and blurred figures.
Obtain the mask (unsharp or smoothed) image.

\[g_{\text{mask}}(x,y) = f(x,y) - \overline{f}(x,y)\]

Add a weighted portion of the mask back to the original image

\[g(x,y) = f(x,y) + k g_{\text{mask}}(x,y)\]

where $k \ge 0$, is a weight.
When $k=1$, we get unsharp masking.
When $k > 1$, we get highboost filtering.
$k < 1$ deemphasizes the contribution of the unsharp mask.

Image sharpening using unsharp masking

The top row image in slide~\ref{sld:USMasking} is a slightly blurred image.
The image in the second row is the result of applying a Gaussian smoothing filter of size $5 \times 5$ with $\sigma=3$.
The image in the third row is the unsharp mask using equation~\eqref{eqn:Gmask}.
The next row image is obtained using unsharp masking (equation~\ref{eqn:Gmask2} with $c=1$).
The last row image is obtained using unsharp masking (equation~\ref{eqn:Gmask2} with $c=4.5$).

The Gradient

First derivatives are implemented using the magnitude of the gradient.
$f(x,y)$ is the image, and the gradient of $f$ at coordinates $(x,y)$ is the column vector

\[\nabla f \equiv \text{grad}(f) \equiv \begin{bmatrix} g_x \\ g_y \end{bmatrix} = \begin{bmatrix} \frac{\partial f}{\partial x} \\ \\ \frac{\partial f}{\partial y} \end{bmatrix}\]

This vector points in the direction of the greatest rate of change of $f$ at location $f(x,y)$.
The magnitude of the vector $\nabla f$, denoted as $M(x,y)$

\[M(x,y) = \text{mag}(\nabla f) = \sqrt{(g_x^2 + g_y^2)}\]

$M(x,y)$ is called the gradient image.
The partial derivatives in equation~\ref{eqn:TheGradient} are not rotation invariant (isotropic), but the magnitude of the gradient vector is.
Computationally more suitable to approximate the squares and square root operations by absolute values.

\[M(x,y) \approx \,\, \mid \! g_x \! \mid + \mid \! g_y \! \mid\]

Isotropic property is lost due to the above approximation.
However, this is not significant as most popular masks used to approximate the gradient are isotropic at multiples of $90^\circ$.

Discrete approximations of the gradient

We use the notation depicted in top row in slide~\ref{sld:Roberts}.
$z_5$ corresponds to $f(x,y)$, $z_1$ corresponds to $f(x-1, y-1)$, and so on.
Simplest approximations to first-order derivatives:

\[g_x = (z_8 - z_5) \text{ and } g_y = (z_6 - z_5)\]

Approximations proposed by Roberts uses cross differences:

\[g_x = (z_9 - z_5) \text{ and } g_y = (z_8 - z_6)\]

Using equations~\eqref{eqn:TheGradientM} and~\eqref{eqn:Roberts}, gradient is computed as

\[M(x,y) = [(z_9 - z_5)^2 + (z_8 - z_6)^2]^\frac{1}{2}\]

If we use equations~\eqref{eqn:TheGradientApprox} and~\eqref{eqn:Roberts}, gradient is computed as

\[M(x,y) \approx \,\, \mid \!z_9 - z_5\!\mid + \mid\!z_8 - z_6 \!\mid\]

Roberts approximations to the gradient (equation~\eqref{eqn:Roberts}) can be implemented using the masks shown in the center row of slide~\ref{sld:Roberts} — Roberts cross-gradient operators.

Approximations to $g_x$ and $g_y$

Approximations to $g_x$ and $g_y$ using a $3 \times 3$ neighborhood centered on $z_5$ are:

\[g_x &= \frac{\partial f}{\partial x} = (z_7 + 2z_8 + z_9) - (z_1 + 2z_2 + z_3) \\ g_y &= \frac{\partial f}{\partial y} = (z_3 + 2z_6 + z_9) - (z_1 + 2z_4 + z_7)\]

$g_x$ and $g_y$ can be implemented using the masks shown in the bottom row of slide~\ref{sld:Roberts} — Sobel operators.
Substituting the values of $g_x$ and $g_y$ in equation~\eqref{eqn:TheGradientApprox} gives:

\[M(x,y) \approx \,\,\mid \! \! (z_7 &+ 2z_8 + z_9) - (z_1 + 2z_2 + z_3) \!\! \mid \\ &+ \mid \! \!(z_3 + 2z_6 + z_9) - (z_1 + 2z_4 + z_7) \! \! \mid\]

Gradient operators

The coefficients in all the masks shown in slide~\ref{sld:Roberts} sum to zero — filter response is zero in an area of constant intensity.
$g_x$ and $g_y$ are linear operators because they involve derivatives — implemented as sum of products.
The nonlinear aspect of sharpening with the gradient comes from the computation of $M(x,y)$, which involves squaring, square roots, or finding absolute values — nonlinear operations.

Use of the gradient for edge enhancement

Gradient is used frequently in industrial inspection to aid humans in the detection of defects, or as a preprocessing step in automated inspection.
Gradient can be used to enhance defects and to eliminate slowly changing background features.
Shown in left in slide~\ref{sld:ContactLens} is an optical image of a contact lens.
Shown in right is the gradient image of equation~\eqref{eqn:TheGradientApprox} implemented using the Sobel masks shown in slide~\ref{sld:Roberts}.

Combining spatial enhancement methods: example

This example illustrates how to combine several of the approaches developed thus far to address a difficult image enhancement task.
Shown in top, left in slide~\ref{sld:Combination1} is a nuclear whole body bone scan — to detect bone infection and tumors.
Narrow dynamic range of intensity values and high noise content makes enhancement by sharpening and bringing out more skeletal detail difficult.
Strategy: utilize the Laplacian to highlight fine detail, and the gradient to enhance prominent edges.

Top, right image is the Laplacian (using the bottom, right filter in slide~\ref{sld:LapFilters}).
Add the first row figures to obtain the bottom, left image — too noisy.
We can remove noise by using a median filter. However, median filtering is a nonlinear process capable of removing image features — unacceptable for medical applications.
The result of applying the Sobel gradient on the original image is shown in bottom, right. Edges are much more dominant now.
Shown in top, left in slide~\ref{sld:Combination2} smoothed gradient image by using an averaging filter of size $5 \times 5$.
The product of the Laplacian and smoothed gradient image is shown in top, right — mask image.
Adding the mask to the original image resulted in the bottom, left image. Notice the significant increase in sharpness of detail.
Next, increase the dynamic range of the above image. Apply power-law transformation with $\gamma = 0.5$ and $c=1$. The resulting image is shown in bottom, right. Notice the significant new detail.

Fuzzy Techniques

Summary

Intensity transformation methods.
Spatial filtering methods.
Though intensity transformation and spatial filtering methods are demonstrated for image enhancement task, they have other applications in DIP.

Back to course home

Ch-03 intensity transformations and spatial filtering

What is this lesson about?

Spatial domain processing

A \(3 \times 3\) neighborhood about a point \((x,y)\)

Intensity transformation

Intensity Transformation Functions

Basic intensity transformation functions

Image negatives

Original mammogram image and its negative

Log Transformations

Fourier spectrum and \(\log\) transformation

Fourier spectrum and its \(\log\) transformation

Power-law (\(\gamma\)) transformations

Plots of the equation \(s = cr^\gamma\) for various values of \(\gamma\)

Gamma \((\gamma)\) correction

Intensity ramp image with \(\gamma\) correction

Power-law transformations for contrast manipulation

Contrast enhancement of an MRI image

Power-law transformations for contrast manipulation

An aerial image and \(\gamma\) manipulation for improving contrast

Piecewise-Linear Transformation Functions

Contrast stretching

Contrast stretching and thresholding

Contrast stretching

Contrast stretching

Intensity-level slicing

Transformation functions for intensity-level slicing

Highlight major blood vessels

Aortic angiogram and intensity-level transformations

Bit-plane representation of an 8-bit image

Bit-plane slicing

An 8-bit image and its 8 bit planes (1 through 8)

Why bit-plane slicing?

Reconstruction using bit planes

Histogram Processing

Intensity distribution and histogram (a 3-bit, \(64 \times 64\) image)

Image histograms

Dark, light, low-contrast, and high-contrast image types

Histogram Equalization

Histogram equalization

Probability density function (PDF)

Cumulative distribution function (CDF)

CDF

Applying \(s = T(r) = (L-1) \int_0^r p_r(w) \, dw\)

Example to illustrate equations

Discrete form of intensity transformations

Discrete form of intensity transformations

Example to illustrate histogram equalization

Intensity distribution and histogram values

Example to illustrate histogram equalization

Example to illustrate histogram equalization

Histogram equalized images

Histogram Matching (Specification)

Histogram equalization based on a specified histogram

Histogram matching (specification)

Histogram matching procedure

Example to illustrate histogram specification

Histogram specification

Histogram equalization: discrete formulation

Summary of histogram specification procedure

Example to illustrate histogram specification

Histogram equalization based on a specified histogram

Specified and actual histograms

Example to illustrate histogram specification

Example to illustrate histogram specification: \(G\) values

Example to illustrate histogram specification

Example to illustrate histogram specification \(s_k \rightarrow z_q\)

Mars moon: Phobos

Mars moon and its histogram

Histogram equalized Mars moon image

Mars moon: Phobos

Mars moon: Phobos (enhanced image using mappings from curve (2))

Local Histogram Processing

Local histogram processing example

Global and local histogram equalizations

Image histogram statistics

Image histogram statistics example

Global (and local) mean and variance

Local enhancement using histogram statistics

Global and local histogram equalization