%
%
%  I wanted to add another problem, but I'm out of time.  Do it next year.
%
%  linear blur + noise on cameraman.
%  first, use the inverse filter with just the blur.
%  then do inverse filter on blur + noise,  should be disaster!
%  then try the wiener filter.
%
\documentclass[12pt]{article} 

%\usepackage{times}
\usepackage{exscale}
\usepackage{graphicx}
\usepackage{amsmath,amssymb,eucal}
\input{/home/joebob/texfiles/Class/ece3793/homework/mydefs}

\setlength{\oddsidemargin}{0.0pt}
\setlength{\evensidemargin}{0.0pt}
\setlength{\topmargin}{-0.75in}
\setlength{\textwidth}{6.5in}
\setlength{\textheight}{9.25in}

\begin{document}
\pagestyle{empty}


{\Large\bf
\centerline{ECE 5273}
\centerline{Homework 5}
}


\bigskip
\par\noindent
Spring 2021 \hfill Dr.~Havlicek
\baselineskip=16pt
\parskip=16pt

\begin{enumerate}
%
% P1
%
\item
Obtain the image {\tt salesman.bin} from the course web site.  This is a
$256 \times 256$ grayscale image with 8-bit pixels.  In this problem you will
use three different methods to apply a $7 \times 7$ linear low-pass average filter
to the image.  The resulting filtered images should be exactly the same in all three
cases.  

The filter impulse response is a $7 \times 7$ square where each pixel is
equal to $1/49$.  In terms of the notation on page 2.54 of the course notes, the
filter window is SQUARE(49).  The reference point is located in the center, so that the
output image is not shifted relative to the input image.

{\bf Note:} this is a linear translation invariant (LTI) filter.  The output is
linear convolution, {\em not\/} circular convolution.

\begin{enumerate}
\item
Implement the $7 \times 7$ linear average filter using image domain convolution.  Set
each pixel in the output image equal to the average of the pixels in a $7 \times 7$
neighborhood about the corresponding pixel in the input image.  Handle edge effects
by zero padding (i.e., use a value of zero for pixels where the window ``hangs over"
the edges of the input image).  Show the input and output images with full scale contrast.

{\bf Note:} to obtain agreement with the DFT filtering methods that you will implement
in parts (b) and (c), you need to apply the filter at {\em every\/} pixel of the 
input image.  At all four edges of the image, this creates an issue about how to
implement the zero padding.  One way to handle the issue is to treat the first
three rows, the last three rows, the first three columns, and the last three columns
as special cases.  An alternative way that is probably easier is the following:
add three rows of all zero pixels to the top and bottom of the image and add three
columns of all zero pixels to the left and right sides of the image.  Then you will
have a new image of size 262 $\times$ 262.  You can then apply the filter to this
new enlarged image starting at pixel (4,4) (for Matlab indexing; for C indexing it is
pixel (3,3)) and ending at pixel (259,259) (for Matlab indexing; for C indexing it is
pixel (258,258)).

\item
Implement the same filter by pointwise multiplication of DFT's using the method of
Example 3 on pages 5.61-5.62 of the course notes.  In this case, your impulse response
image $\mathbf H$ will be of size $128 \times 128$ and will have all pixels equal to
zero except for a $7 \times 7$ square of pixels in the center that will all be equal to
$1/49$.  As in Example 3, we have $(p_0, q_0) = (64, 64)$, which is Matlab
$(\mbox{row}, \mbox{col}) = (r, s) = (65, 65)$.  This means that the {\em center\/} 
of the
square (the reference point shown in green on page 2.54 of the course notes) will be
for $(m, n) = (0, 0)$ and should be located at $(p, q) = (64, 64)$ in the $\mathbf H$
image (which is again $(r, s) = (65, 65)$ in Matlab).

{\bf Note:} according to Example 3 on pages 5.61-5.62 of the course notes, the
zero-padded size should be 383 $\times$ 383 
(since 256 + 128 - 1 = 383).  If you are using C and the
FFT routine {\tt fft2d()} from the course web site, this will cause a problem
because {\tt fft2d()} assumes that the size is {\em even}.  To solve this problem,
increase the zero-padded size by one to 384 $\times$ 384.

Show the following images with full scale contrast:
\begin{itemize}
\item The original input image.
\item The zero padded original image.
\item The zero padded impulse response image.
\item The centered DFT log-magnitude spectrum of the zero padded input image.
\item The centered DFT log-magnitude spectrum of the zero padded impulse response image.
\item The centered DFT log-magnitude spectrum of the zero padded output image.
\item The zero padded output image.
\item The final $256 \times 256$ output image.
\end{itemize}

Verify that the output image 
(obtained after performing a full-scale contrast stretch and rounding each pixel
to the nearest integer)
is the same as the one in part (a).  
In Matlab, if the output image from part (a) is in the array {\tt Y1a} and the output 
image from this part is in the array {\tt Y1b}, then you can do this using the Matlab
statement

{\footnotesize
\begin{verbatim}
disp(['(b): max difference from part (a): ' num2str(max(max(abs(Y1b-Y1a))))])
\end{verbatim}
}

\item
Implement the same filter again using the zero-phase impulse response and
DFT method of Example 5 on page
5.76 of the course notes.  This method is explained on pages 5.66-5.75.  As in
Example 5, use a size of $256 \times 256$ for the impulse response image $\mathbf H$.
Your solution should be almost the same as the Matlab code in Example 5, except that 
you will replace the Gaussian impulse response in lines 3 and 4 with a
$7 \times 7$ square of pixels that are equal to $1/49$.

To do this in Matlab, initialize the original non-zero-phase impulse response 
({\tt h} in Example 5) to all zeros and then set elements $(\mbox{126:132},\mbox{126:132})$ 
equal to $1/49$.  This will give you a $7 \times 7$ square of pixels with value $1/49$
centered at $(p_0, q_0) = (128, 128)$,
which is Matlab $(\mbox{row}, \mbox{col}) = (r, s) = (129, 129)$.  Then use the
Matlab {\tt fftshift} function to obtain the zero-phase impulse response image 
exactly as in Example 5.  If you are using C or some other language
then you will need
to program the equivalent operations yourself, but you can still use the example
Matlab code on page 5.76 as a guideline.

\newpage

Show the following images with full scale contrast:
\begin{itemize}
\item The original input image.
\item The $256 \times 256$ zero-phase impulse response image (this is {\tt h2} in
Example 5).
\item The $512 \times 512$ zero padded zero-phase impulse response image (this is
{\tt h2ZP} in Example 5).
\item The final $256 \times 256$ output image.
\end{itemize}

Verify that the output image 
(obtained after performing a full-scale contrast stretch and rounding each pixel
to the nearest integer)
is the same as the one in part (a).
\end{enumerate}

%
% P2
%
\item
Obtain the images {\tt girl2.bin}, {\tt girl2Noise32Hi.bin}, and
{\tt girl2Noise32.bin} from the course web site.  All three are
$256 \times 256$ grayscale images with 8-bit pixels.
\begin{itemize}
\item {\tt girl2.bin} is the original image, also known as {\em Tiffany}.
\item {\tt girl2Noise32Hi.bin} has hi-pass Gaussian white noise added.  The DFT of the
noise is nonzero for $\sqrt{u^2 + v^2} > 64$ cpi and zero for $\sqrt{u^2 + v^2} \leq 64$.
\item {\tt girl2Noise32.bin} has broadband Gaussian white noise added.
\end{itemize}

\begin{enumerate}
\item
Read and display all three images.  Compute the MSE of {\tt girl2Noise32Hi} and
{\tt girl2Noise32} relative to {\tt girl2} using the formula given on page 5.114
of the course notes.
\item
Apply an isotropic ideal low-pass filter as on page 5.111 of the course notes 
with $U_{\rm cutoff} = 64$ cpi to all three images.  Since the cutoff frequency of the
hi-pass noise is known exactly in this case, it is best to implement the filter using
circular convolution in order to avoid the edge effects associated with zero
padding and linear convolution.
Therefore, for all three images, use the formula on page 5.111
of the notes to make a $256 \times 256$ $\widetilde{\mathbf H}$ DFT array.  In Matlab,
you can do it like this:

\begin{verbatim}
U_cutoff = 64;
[U,V] = meshgrid(-128:127,-128:127);
HLtildeCenter = double(sqrt(U.^2 + V.^2) <= U_cutoff);
HLtilde = fftshift(HLtildeCenter);
\end{verbatim}
Then, simply pointwise multiply this $\widetilde{\mathbf H}$ array times the
$256 \times 256$ DFT array of each image to obtain the three filtered images.


Display each filtered image and compute the MSE (relative to the original {\tt girl2} image)
using the formula on
page 5.114 of the course notes.  For the two noisy images, also compute
the ISNR using the formula on page 5.115 of the course notes.
You should display the images with full scale contrast.  However, for computing the
MSE and ISNR, it is {\bf important} to use the floating point images obtained directly
from the inverse
DFT without rounding and without applying any full scale contrast stretch.

You should find that this ideal LPF is very effective for removing the hi-pass noise
in {\tt girl2Noise32Hi}.  However, it will introduce ringing artifacts in the
original image {\tt girl2} and it will not work as well on the broadband noise
in {\tt girl2Noise32}.

\vspace{12pt}

{\bf Note:} {\em if\/} you are using C, and {\em if\/} you are using the routine
{\tt fft2d()} from the course web site, then you may run into an issue here
in Problem 2(b) (if you are using Matlab, then you can ignore this paragraph). 
The reason you may run into an issue is that {\tt fft2d()} is designed to input and
output {\em centered\/} DFT's.  On a forward transform, it performs an fftshift
operation after calling the 2D FFT.  On a reverse transform, it performs an
fftshift before calling the 2D FFT.  Now, for doing DFT domain filtering, it
does not really matter if you are using centered DFT's or un-centered DFT's,
{\em as long as you use\/} {\tt fft2d()} {\em to compute both transforms: the
transform of the input image and the transform of the filter}.  However, if, as
we are doing here in Problem 2(b), you are only going to use {\tt fft2d()} for
{\em one\/} of the transforms, then you have to be careful!  Here in Problem 2(b),
you are going to compute the $\widetilde{\mathbf H}$ 
image yourself {\em without\/} calling
{\tt fft2d()}.  So, if you {\em are\/} going to call {\tt fft2d()} for the DFT of the
input images, then realize that they will be {\em centered\/} DFT's.  This means
that you must also compute the $\widetilde{\mathbf H}$ image as a 
{\em centered\/} DFT.  Otherwise,
you will get the wrong answer when you multiply the two DFT's pointwise.
That is a difference between using C and {\tt fft2d()} as compared to
the Matlab code shown above.

\item
Apply the Gaussian low-pass filter on page 5.113 of the course notes to all three images
with $U_{\rm cutoff} = 64$.  Use the bottom formula with $N=256$ and $U_{\rm cutoff} = 64$
to design the space constant $\sigma$.  Then use the top formula to make a
$256 \times 256$ $\widetilde{\mathbf H}$ DFT array.  Apply {\tt fftshift} to ``un-center"
this DFT array and invert to obtain the $256 \times 256$ zero phase impulse response
(as explained in the note for part (b),
if you are using C and {\tt fft2d()} then you do {\em not\/} have to un-center
the DFT array before inverting).
Apply {\tt fftshift} again to the zero phase impulse response to center it.  Then use
the method of Example 4 on pages 5.64 and 5.65 of the course notes to obtain the three
filtered images.

Here is an example of how to use Matlab to compute the $512 \times 512$ zero padded
impulse response image:
\begin{verbatim}
U_cutoff_H = 64;
SigmaH = 0.19 * 256 / U_cutoff_H;
[U,V] = meshgrid(-128:127,-128:127);
HtildeCenter = exp((-2*pi^2*SigmaH^2)/(256.^2)*(U.^2 + V.^2));
Htilde = fftshift(HtildeCenter);
H  = ifft2(Htilde);
H2 = fftshift(H);
ZPH2 = zeros(512,512);
ZPH2(1:256,1:256) = H2;
\end{verbatim}
You can then take the $512 \times 512$ DFT of {\tt ZPH2} and multiply it pointwise with
the DFT's of the $512 \times 512$ zero padded images to obtain the DFT's of the
three zero padded filtered images.

Display each $256 \times 256$ filtered image and compute the MSE using the formula on
page 5.114 of the course notes.  For the two noisy images, also compute
the ISNR using the formula on page 5.115 of the course notes.
As in part (b), for computing the
MSE and ISNR it is {\bf important} to use the floating point images obtained directly
from the inverse
DFT without rounding and without applying any full scale contrast stretch.

Compared to the ideal LPF,
you should find that this filter produces a much better MSE on the original image
since it does not suffer from ringing artifacts.  However, this filter will not
perform as well as the ideal LPF on the noisy images.

\item
Use the method of Example 4 on pages 5.64 and 5.65 of the course notes to
apply the Gaussian low-pass filter on page 5.113 of the course notes to all three images
again, this time with $U_{\rm cutoff} = 77.5$.

Display each $256 \times 256$ filtered image and compute the MSE using the formula on
page 5.114 of the course notes.  For the two noisy images, also compute
the ISNR using the formula on page 5.115 of the course notes.
As in part (b), for computing the
MSE and ISNR it is {\bf important} to use the floating point images obtained directly
from the inverse
DFT without rounding and without applying any full scale contrast stretch.

You should find that this filter distorts the original image much less than the other
two filters.  Overall, this filter
will perform better than the other two against the broadband noise in {\tt  girl2Noise32}.
\end{enumerate}






\end{enumerate}


\bigskip
\bigskip
{\large\bf DUE: 4/20/2021, 11:59 PM}

\end{document}

