\documentclass[12pt]{article}
\usepackage{graphicx}
\usepackage{float}
\usepackage{hyperref}
\usepackage{multirow}
\usepackage{amsmath}
\usepackage{enumitem}
\usepackage{subcaption}
\usepackage{url}
\usepackage{array}
\usepackage[table]{xcolor}

\usepackage{listings}
\lstset{%
%backgroundcolor=\color{cyan!10},
basicstyle=\ttfamily,
numbers=left,numberstyle=\scriptsize
}

% adjust margins
\addtolength{\oddsidemargin}{-.5in}
\addtolength{\evensidemargin}{-.5in}
\addtolength{\textwidth}{1in}

\addtolength{\topmargin}{-.5in}
\addtolength{\textheight}{1in}

\renewcommand{\arraystretch}{2}


% define commands for tables
\makeatletter
\newcommand{\thickhline}{%
    \noalign {\ifnum 0=`}\fi \hrule height 1pt
    \futurelet \reserved@a \@xhline
}
\newcolumntype{"}{@{\hskip\tabcolsep\vrule width 1pt\hskip\tabcolsep}}
\makeatother

\DeclareMathOperator{\atantwo}{atan2}



\title{Filtering Digital Images}
\author{EE 224: Signals and Systems I}
\date{}

\begin{document}

\maketitle

\section{Overview}

In this lab, we introduce linear sliding window filtering of images. The filters work in two dimensions to enhance the image in ways such as averaging out noise or emphasizing edges.

\section{Learning Objectives}

By the end of this lab, students will be able to:
\begin{enumerate}
\item Describe spatial frequency in 2D digital images.
\item Design and implement linear, spatially invariant (LSI) filters using two-dimenional convolutional windows.
\item Use the properties of LSI filters to implement the Canny Edge Detection Protocol for images.
\end{enumerate}

\section{Pre-Lab Reading and Exercises}
\label{sec:prelab}

\subsection{Image Filtering and Spatial Frequency}

Image analysis deals with techniques for extracting information from images. The first step is generally to segment an image. Segmentation divides an image into its constituent parts or objects. For example, in a military air-to-ground target acquisition application one may want to identify tanks on a road. The first step is to segment the road from the rest of the image and then to segment objects on the road for possible identification. Segmentation algorithms are usually based on one of two properties of gray-level values: discontinuity and similarity. For discontinuity, the approach is to partition an image based on abrupt changes in gray level. Objects of interest are isolated points, lines, and edges.

\begin{figure}[h]
\centering
\includegraphics[width=0.6\textwidth]{edges.png}
\caption{Synthetic images with horizontal pixel intensity, first derivative of intensity, and second derivative of intensity.}
\label{fig:edge}
\end{figure}

An edge is the boundary between two regions with relatively distinct gray-level properties. The idea underlying most edge-detection techniques is the computation of a relative derivative operator. Figure~\ref{fig:edge} illustrates this concept. The picture on the left shows an image of a light stripe on a dark background, the gray-level profile of a horizontal scan line, and the first and second derivatives of the profile. The first derivative is positive when the image changes from dark to light and zero where the image is constant. The second derivative is positive for the part of the transition associated with the dark side of the edge and negative for the transition associated with the light side of the edge. Thus, the magnitude of the first derivative can be used to detect the presence of an edge in the image and the sign of the second derivative can be used to detect whether a pixel lies on the light or dark side of the edge.


\subsection{Sliding Window Filters}

One way of filtering is to first filter the rows with a one-dimensional filter and then filter the columns with a one-dimensional filter. Alternatively, two dimensional filters can be applied to the neighborhood of each pixel to filter a signal. Filtering can help remove noise from an image by averaging out the effects of random noise fluctuations. High pass filters help detect edges and sudden spatial changes in an image.

\begin{figure}
     \centering
     \begin{subfigure}[b]{0.45\textwidth}
         \centering
         \includegraphics[width=\textwidth]{spatialFilter1.png}
         %\caption{$y=3sinx$}
         %\label{fig:three sin x}
     \end{subfigure}
     \hfill
     \begin{subfigure}[b]{0.45\textwidth}
         \centering
         \includegraphics[width=0.8\textwidth]{spatialFilter2.png}
         %\caption{$y=5/x$}
         %\label{fig:five over x}
     \end{subfigure}
        \caption{Spatial filtering with a sliding window; (left) three by three neighborhood of pixel $(i,j)$; (right) three by three spatial mask.}
        \label{fig:spatial}
\end{figure}

The operations take place in the \emph{spatial domain} and operate directly on pixel values. The general form for operations is: $y(i, j) = T[x (i, j)]$ where $x(i,j)$ is the input image, $y(i,j)$ is the filtered output image, $T$ is an operator on $x$ defined over a neighborhood of point $(i,j)$. The operator, $T$, defines the size of the neighborhood around each pixel. The neighborhood is usually rectangular, centered on $(i,j)$ and is much smaller than the image. This neighborhood is moved around to each pixel, hence the name sliding window or spatial mask. Use of spatial masks for filtering is called spatial filtering and the masks can be linear or nonlinear. As a note, the linear filter is equivalent to a two-dimensional convolution of the image with the filter.

A three by three linear filter can be written in matrix form as:
$$\begin{bmatrix}
w_{1,1} & w_{1,0} & w_{1,-1} \\
w_{0,1} & w_{0,0} & w_{0,-1} \\
w_{-1,1} & w_{-1,0} & w_{-1,-1}
\end{bmatrix}$$
The output of applying a three by three mask to image $x$ is given by the two-dimensional convolution:
\begin{align*}
y(i,j) &= \sum_{k=-1}^{1} \sum_{l=-1}^{1} w_{k,l} x(i-k,j-l) \\
&= w_{1,1} x(i-1,j-1) + w_{1,0} x(i-1,j) + w_{1,-1} x(i-1,j+1) \\
& \quad + w_{0,1} x(i,j-1) + w_{0,0} x(i,j) + w_{0,-1} x(i,j+1) \\
& \quad + w_{-1,1} x(i+1,j-1) + w_{-1,0} x(i+1,j) + w_{-1,-1} x(i+1,j+1) \\
\end{align*}

This masking process continues until all pixels in the image are covered. The standard practice is to create a new image with the new values. For pixels near the boundary of the image, the output may be computed using partial neighborhoods or by padding the input appropriately by adding in repeated edge rows/columns or by a frame of black or white pixels.
 
These types of spatial filters are linear and spatially invariant (LSI). They can also be interpreted in the frequency domain. Lowpass filters attenuate (or eliminate) high frequency components characterized by edges and sharp details in an image with the net effect of blurring the image. Highpass filters attenuate (or eliminate) low frequency components such as slowly varying characteristics. The overall effect is a sharpening of edges. Bandpass filters attenuate (or eliminate) a given frequency range of frequencies. This is primarily used for image restoration.


\subsection{Technical Notes for Implementing Filter Masks}

In this lab, the MATLAB routine \texttt{imfilter} will be used: \texttt{Y = imfilter(X,h);} where $X$ is the image to be filtered, $h$ is the filter in matrix form, and $Y$ is the filtered image.
A simple example for a three by three averaging filter is given below:

\begin{lstlisting}
>>X=imread(‘moon.tif’);
>>h = ones(3)/9;% or [1,1,1;1 1 1;1 1 1]/9;
>>Y = imfilter(X,h);
\end{lstlisting}

There are two things to keep in mind with \texttt{imfilter}. First, it returns an image in the same format as the input. This means that if, for a uint8 image, the function gets a non-integer value, it rounds it to the closest integer. Second, it clips pixel values outside of the range [0, 255] to 0 or 255. These are non-linear operations which effect whether the LSI operations can be done in any order. For pixels near the boundary of the image, the output may be computed using partial neighborhoods or by padding the input appropriately by adding in repeated edge rows/columns or by a frame of black or white pixels. The default in \texttt{imfilter} is to assume a black border. Details can be found by reading about padding options on the help page.

\subsection{Pre-Lab Questions}

\begin{enumerate}
\item{Which part of the image in Figure~\ref{fig:q1} has the lowest spatial frequency: specify top, bottom, left, right, middle, or a corner.
\begin{figure}[H]
\centering
\includegraphics[width=0.35\textwidth]{rover.png}
\caption{Image for question 1.}
\label{fig:q1}
\end{figure}
}
\item{Compute the spatial period of the image below in terms of pixels/cycle. Give a numerical answer.
\begin{figure}[H]
\centering
\includegraphics[width=0.35\textwidth]{horizontalSinusoid.png}
\caption{Image for question 2.}
\label{fig:q2}
\end{figure}
}
\item{Given the 2D filter mask below and the intensity values from an image, compute the output for pixel at location $(3,3)$ in the image.
\begin{center}
%\begin{table}[]
\begin{tabular}{lll}
\hline
\multicolumn{1}{|c|}{-1}                & \multicolumn{1}{c|}{0}                       & \multicolumn{1}{c|}{1}                       \\ \hline
\multicolumn{1}{|c|}{-2}                & \multicolumn{1}{c|}{0}                       & \multicolumn{1}{c|}{2}                       \\ \hline
\multicolumn{1}{|c|}{-1}                & \multicolumn{1}{c|}{0}                       & \multicolumn{1}{c|}{1}                       \\ \hline\end{tabular}
%\end{table}
\end{center}

\begin{center}
%\begin{table}[]
\begin{tabular}{lllllll}
\hline
\multicolumn{1}{|l|}{\textbf{Pixel \#}} & \multicolumn{1}{l|}{\textbf{1}} & \multicolumn{1}{l|}{\textbf{2}} & \multicolumn{1}{l|}{\textbf{3}} & \multicolumn{1}{l|}{\textbf{4}} & \multicolumn{1}{l|}{\textbf{5}} & \multicolumn{1}{l|}{\textbf{6}} \\ \hline
\multicolumn{1}{|c|}{\textbf{1}}  & \multicolumn{1}{c|}{5}  & \multicolumn{1}{l|}{25}  & \multicolumn{1}{c|}{25}  & \multicolumn{1}{l|}{30}  & \multicolumn{1}{c|}{20}  & \multicolumn{1}{l|}{10}                       \\ \hline
\multicolumn{1}{|c|}{\textbf{2}}  & \multicolumn{1}{c|}{20}  & \multicolumn{1}{l|}{20}  & \multicolumn{1}{c|}{30}  & \multicolumn{1}{l|}{60}  & \multicolumn{1}{c|}{25}  & \multicolumn{1}{l|}{30}                       \\ \hline
\multicolumn{1}{|c|}{\textbf{3}}  & \multicolumn{1}{c|}{35}  & \multicolumn{1}{l|}{30}  & \multicolumn{1}{c|}{\cellcolor{lightgray}40}  & \multicolumn{1}{l|}{70}  & \multicolumn{1}{c|}{40}  & \multicolumn{1}{l|}{40}                       \\ \hline
\multicolumn{1}{|c|}{\textbf{4}}  & \multicolumn{1}{c|}{80}  & \multicolumn{1}{l|}{80}  & \multicolumn{1}{c|}{70}  & \multicolumn{1}{l|}{70}  & \multicolumn{1}{c|}{60}  & \multicolumn{1}{l|}{50}                       \\ \hline
\end{tabular}
%\end{table}
\end{center}

}
\end{enumerate}


\section{Lab Exercises}
\label{sec:exercises}

%The images will change across semesters. 

This part of the lab will be done using the three images that were provided with the lab. They are called $I$, $I_g$, and $I_{sp}$. For each piece of the lab, you will need to display the entire image before and after filtering. The command for this is:
\begin{lstlisting} 
>> imshowpair(I,filtI,`montage')
\end{lstlisting}
You will also need to display a close up of parts of the image that illustrate the filter effect well. To do this, it helps to display the image with the pixels numbered to help select the proper subimage. If the routine \texttt{imshow} does not do this automatically, then use the following code to display image $I$, for example:
\begin{lstlisting}
>> RI = imref2d(size(I)); figure; imshow(I,RI);
\end{lstlisting}
Load the Matlab file \emph{matlab.mat} with the images.

\subsection{Low-Pass Spatial Filtering with a Moving Average Filter}
\label{sec:lpf}

A simple three by three or $N=3$ averaging filter has the form:
\begin{center}
%\begin{table}[]
\begin{tabular}{lll}
\hline
\multicolumn{1}{|c|}{$\frac{1}{9}$}                & \multicolumn{1}{c|}{$\frac{1}{9}$}                       & \multicolumn{1}{c|}{$\frac{1}{9}$}                       \\ \hline
\multicolumn{1}{|c|}{$\frac{1}{9}$}                & \multicolumn{1}{c|}{$\frac{1}{9}$}                       & \multicolumn{1}{c|}{$\frac{1}{9}$}                       \\ \hline
\multicolumn{1}{|c|}{$\frac{1}{9}$}                & \multicolumn{1}{c|}{$\frac{1}{9}$}                       & \multicolumn{1}{c|}{$\frac{1}{9}$}                       \\ \hline\end{tabular}
%\end{table}
\end{center}

Write out the 2D filter for the $N = 5$ case. For averaging filters, the sum of the weights over the filter should be equal to 1. Use filters with dimensions $N=5$, $9$, and~$13$ to filter the image $I$. What is the effect of increasing the filter size in parts of the image with small details? Include your output images in your report.


\subsection{High-Pass Spatial Filtering}
\label{sec:hpf}

The image filters shown below respond more strongly to lines of different directions in an image. The first mask responds very strongly to horizontal lines with a thickness of one pixel. The maximum response occurs when a line passed through the middle row of the mask. The direction of a mask can be established by noting that the preferred direction is weighted with a larger coefficient than the other possible directions.

\begin{center}
%\begin{table}[]
\begin{tabular}{llllllllllll}
\hline
\multicolumn{3}{|c|}{(a)} & \multicolumn{3}{|c|}{(b)} & \multicolumn{3}{c|}{(c)} & \multicolumn{3}{c|}{(d)}     \\ \thickhline
\multicolumn{1}{|c|}{-1} & \multicolumn{1}{c|}{-1} & \multicolumn{1}{c|}{-1} & \multicolumn{1}{|c|}{-1} & \multicolumn{1}{c|}{-1} & \multicolumn{1}{c|}{2} &
\multicolumn{1}{|c|}{-1} & \multicolumn{1}{c|}{2} & \multicolumn{1}{c|}{-1} & \multicolumn{1}{|c|}{2} & \multicolumn{1}{c|}{-1} & \multicolumn{1}{c|}{-1}     \\ \hline
\multicolumn{1}{|c|}{2} & \multicolumn{1}{c|}{2} & \multicolumn{1}{c|}{2} & \multicolumn{1}{|c|}{-1} & \multicolumn{1}{c|}{2} & \multicolumn{1}{c|}{-1} &
\multicolumn{1}{|c|}{-1} & \multicolumn{1}{c|}{2} & \multicolumn{1}{c|}{-1} & \multicolumn{1}{|c|}{-1} & \multicolumn{1}{c|}{2} & \multicolumn{1}{c|}{-1}     \\ \hline
\multicolumn{1}{|c|}{-1} & \multicolumn{1}{c|}{-1} & \multicolumn{1}{c|}{-1} & \multicolumn{1}{|c|}{2} & \multicolumn{1}{c|}{-1} & \multicolumn{1}{c|}{-1} &
\multicolumn{1}{|c|}{-1} & \multicolumn{1}{c|}{2} & \multicolumn{1}{c|}{-1} & \multicolumn{1}{|c|}{-1} & \multicolumn{1}{c|}{-1} & \multicolumn{1}{c|}{2}     \\ \hline
\end{tabular}
%\end{table}
\end{center}

Use the image $I$. Filter the image using the four filters above. Describe your results for each case. Which edges were highlighted in each case? The MATLAB commands for the first case are given below.
\begin{lstlisting}
>> wa=[-1,-1,-1;2,2,2;-1,-1,-1];
>> I_a=imfilter(I,wa);
\end{lstlisting}
Display your results by plotting the four output images together. Comment on your results and what edges were accented in each case.


\subsection{Filtering Noise}
\label{sec:nfilt}

Noise often has high spatial frequency so that it looks like salt and pepper in the picture. Most objects in the image have lower spatial frequencies so we can often improve the image by low pass filtering with an averaging filter. This technique works well if the noise does not vary that much from the image. Speckle or impulsive noise is very extreme noise that looks like white or black spots on an image. This section looks at how averaging filters affect the image.

\begin{enumerate}
\item{\label{sec:nfilt1} Display the images $I_g$ and $I_{sp}$ in grayscale. These images have been corrupted by two different types of noise. You will use an averaging filter to remove the noise in this image. Filter each image with a 3x3 averaging filter and a 7x7 averaging filter. Display the results in image pairs. Discuss the results. Assess the trade-off between noise reduction and image resolution (how clear objects are in the image).

Compare the results for filtering different types of noise. Does the filtering work better for one type of noise? If so, why?
}
\item{\label{sec:nfilt2} \label{sec:imfilt} Using the images $I_g$ and $I_{sp}$, filter each image with the high pass filters (a) and (c) from Section~\ref{sec:hpf}. Display the results in image pairs. Discuss the results focusing on how clear the edges in each direction appear.

Compare the results for filtering different types of noise. Does the filtering work better for one type of noise? If so, why?
}
\end{enumerate}


\subsection{Filtering Edges in Noise}
\label{sec:nefilt}

Edge detection is a key processing step for finding objects with computer vision systems. Unfortunately, noise can disrupt finding edges, so a commonly used strategy is to first blur the image to remove noise, then to use a high pass filter to find edges. Most practical edge detection operators use multi-stage algorithms to detect a wide range of edges in images. One example is Canny Edge Detection, developed by John F. Canny in 1986. The Canny Edge Detector first smooths the image to blur the noise to prevent false detections. Then, it uses a set of high pass filters to detect edges and estimate their direction. Further processing filters and thins the edges. In this lab, we implement the first two steps.

\begin{enumerate}
\item{\label{sec:nefilt1} \textbf{Smoothing:} Apply a low pass filter to smooth the image to remove noise to prevent false detections caused by noise. The typical filter used by this function is a Gaussian filter kernel which smooths the edges of the filter. Use the filter parameters below to filter image $I_g$:
\begin{lstlisting}
>> hgauss = [2 4 5 4 2; 4 9 12 9 4; 5 12 15 12 5; ...
          4 9 12 9 4; 2 4 5 4 2]/159;
\end{lstlisting}
Display the result and comment on what you see in the image.
}
\item{\label{sec:nefilt2} \textbf{Edge Detection and Intensity Gradient:} An edge in an image may point in a variety of directions, so the Canny algorithm uses filters to detect horizontal, vertical, and diagonal edges in the smoothed image. Use the Sobel Edge Detector filters given below on the result from the previous step with each filter separately to get images filtered in the $x$ and $y$ directions, $G_x$ and $G_y$.
\begin{center}
%\begin{table}[]
\begin{tabular}{llllll}
\hline
\multicolumn{3}{|c|}{(a)} & \multicolumn{3}{|c|}{(b)}     \\ \thickhline
\multicolumn{1}{|c|}{-1} & \multicolumn{1}{c|}{-2} & \multicolumn{1}{c|}{-1} & \multicolumn{1}{|c|}{-1} & \multicolumn{1}{c|}{0} & \multicolumn{1}{c|}{1}     \\ \hline
\multicolumn{1}{|c|}{0} & \multicolumn{1}{c|}{0} & \multicolumn{1}{c|}{0} & \multicolumn{1}{|c|}{-2} & \multicolumn{1}{c|}{0} & \multicolumn{1}{c|}{2}     \\ \hline
\multicolumn{1}{|c|}{1} & \multicolumn{1}{c|}{2} & \multicolumn{1}{c|}{1} & \multicolumn{1}{|c|}{-1} & \multicolumn{1}{c|}{0} & \multicolumn{1}{c|}{1}     \\ \hline
\end{tabular}
%\end{table}
\end{center}
Finally, combine the two filtered images to compute the intensity gradient and local direction at each pixel in the image:
\begin{align*}
G &= \sqrt{G_x^2 + G_y^2} \\
\theta &= \atantwo{(G_y,G_x)}
\end{align*}
Display the intensity gradient image, $G$, using the command \texttt{imshow(G,[0,255])} and comment on your results from Section~\ref{sec:nfilt} Part~\ref{sec:imfilt} above. Are the edges clearer in one case or the other?
}
\end{enumerate}

\subsection{Thought Question}
\label{sec:thought}
What will happen if the order of processing for an LSI low pass filter followed by an LSI high-pass filter is reversed for edge detection? Explain if order matters or not.


\section{Report Checklist}

Be sure the following are included in your report.
\begin{enumerate}
\item Section~\ref{sec:prelab}: answer the three pre-lab questions
\item Section~\ref{sec:lpf}: write out 5x5 filter
\item Section~\ref{sec:lpf}: report on filter effects for $N=5$, $9$, and $13$; include images
\item Section~\ref{sec:hpf}: describe results of applying the four filters, answer related questions; include images
\item Section~\ref{sec:nfilt}, Part~\ref{sec:nfilt1}: answer questions about noise reduction vs. resolution; include images
\item Section~\ref{sec:nfilt}, Part~\ref{sec:nfilt2}: answer questions about edge enhancement; include images
\item Section~\ref{sec:nefilt}, Part~\ref{sec:nefilt1}: display filtered image and comment
\item Section~\ref{sec:nefilt}, Part~\ref{sec:nefilt2}: display gradient and direction images; compare to Section~\ref{sec:nfilt}, Part~\ref{sec:nfilt2}
\item Section~\ref{sec:thought}: answer the thought question
\end{enumerate}

\begin{thebibliography}{9}

\bibitem{gonzalez}
R.C. Gonzalez and R.E. Woods,
\textit{Digital Image Processing},
Addison Wesley, 1993

 
\bibitem{matlabIPT}
``MATLAB Image Processing Toolbox Documentation'',
\textit{The Mathworks},
\\\url{https://www.mathworks.com/help/images/index.html?s_tid=CRUX_lftnav}

\bibitem{cannyWiki}
``Canny Edge Etector'',
\textit{Wikipedia},
\\\url{https://en.wikipedia.org/wiki/Canny_edge_detector}

\bibitem{canny1986}
J. Canny,
``A Computational Approach to Edge Detection,''
\textit{IEEE Transactions on Pattern Analysis and Machine Intelligence,}
vol. PAMI-8, no. 6, pp. 679-698, Nov. 1986

\bibitem{green2002}
Bill Green,
``Canny Edge Detection Tutorial'',
\\\url{https://web.archive.org/web/20160324173252/http:/dasl.mem.drexel.edu/alumni/bGreen/www.pages.drexel.edu/_weg22/can_tut.html}

\bibitem{kalra}
Prem K. Kalra,
``Canny Edge Detection'',
\\\url{https://web.archive.org/web/20150421090938/http:/www.cse.iitd.ernet.in/~pkalra/csl783/canny.pdf}

\bibitem{nasa}
\textit{NASA Image and Video Library}
\url{https://images.nasa.gov/}

\end{thebibliography}



%\end{thebibliography}





\end{document} 