%---TeX-start-of-header---
\documentclass{report}
\usepackage{epsfig,mathematica,figi}

\mmanobreak

\begin{comment}
\begin{mathematica}
zero = 0.00001; \
ohide = DisplayFunction->Identity; \
oshow = DisplayFunction -> $DisplayFunction; \
ojoin = PlotJoined -> True; \
oall = PlotRange->All; \
osurf = HiddenSurface->False; 
\end{mathematica}
\begin{mathematica}
<<tek.m; \
<<ip.m; \
<<pgm.m;
\end{mathematica}
\begin{mathematica}
Periodogram[l_List] := Abs[ FFT[l] ]^2 / Apply[Times, Dimensions[l]];
\end{mathematica}

\end{comment}


\newcommand{\bigpagexv}[2]{
\begin{figure}[p!]
\centering
\lfigiii{#11_1}{#11_2}{#11_3}{}
\medskip
\lfigiii{#12_1}{#12_2}{#12_3}{}
\medskip
\lfigiii{#13_1}{#13_2}{#13_3}{}
\medskip
\lfigiii{#14_1}{#14_2}{#14_3}{}
\medskip
\lfigiii{#15_1}{#15_2}{#15_3}{}
\medskip
  \caption{#2}
  \label{#1}
\end{figure}
}

\newcommand{\bigpagexii}[2]{
\begin{figure}[p!]
\centering
\lfigiii{#11_1}{#11_2}{#11_3}{}
\medskip
\lfigiii{#12_1}{#12_2}{#12_3}{}
\medskip
\lfigiii{#13_1}{#13_2}{#13_3}{}
\medskip
\lfigiii{#14_1}{#14_2}{#14_3}{}
\medskip
  \caption{#2}
  \label{#1}
\end{figure}
}

\begin{document}
\psfiginit
%---TeX-end-of-header---

\begin{titlepage}

\author{Tomasz J.\ Cholewo}
\title{EE\,618 Project \#5\\
  On Image Restoration -- supplementary supplement
}
\date{\today}
\maketitle
\vfil
\tableofcontents
\vfil
\thispagestyle{empty}
\end{titlepage}
%====================================================================
%--------------------------------------------------------------------
\chapter{Image Restoration}

\begin{comment}
\section{Function definitions}
\subsection{Wiener filtering}

Function \texttt{WienerFilter} returns modified Wiener filter sequence
for given signal power spectral estimate $P_f$, white noise power $P_v$
and coefficients $\alpha$ and $\beta$.  It is assumed that signals are
zero-mean.
\begin{mathematica}
WienerFilter[Pf_List, 0, _, _] := 1; \
WienerFilter[Pf_List, _, 0, _] := 1; \
WienerFilter[Pf_List, Pv_, alpha_, beta_] := 
    (Pf / (alpha Pv + Pf))^beta;
\end{mathematica}
%--------------------------------------------------------------------
\subsection{Adaptive Wiener filtering}

Function \texttt{MeanList} creates a list of local averages over the
window having width $2M+1$.
\begin{mathematica}
MeanList2D[l_List, 0] := l; \
MeanList2D[l_List, M_Integer] :=
    Table[ Sum[index2d[l, i + j, k + m], {j, -M, M}, {m, -M, M}] / (2. M+1)^2,
 {i, 1, Length[l]}, {k, 1, Length[l[[1]]]}];
\end{mathematica}
Function \texttt{VarianceList} calculates the variance in a similar way.
It can be used for signal power estimation.
\begin{mathematica}
VarianceList2D[l_List, M_Integer] := MeanList2D[ (l - MeanList2D[l, M])^2, M ];
\end{mathematica}
The following function implements adaptive Wiener filter for disturbing
noise power $P_v$ and window size $2M+1$.
\begin{mathematica}
AdaptiveWienerFilter2D[g_List, _, 0] := g; \
AdaptiveWienerFilter2D[g_List, M_Integer, Pv_] := Module[{ave, var},
    ave = MeanList2D[g, M];
    var = MeanList2D[ (g - ave)^2, M ];  (* VarianceList inlined here *)
    ave + var / (var + Pv) (g - ave)
];
\end{mathematica}
%--------------------------------------------------------------------
\subsection{Inverse filtering}

The following function calculates inverse filter for the given blurring
sequence $b$ using threshold value $\gamma$.
\begin{mathematica}
InverseFilter2D[{1}, _] := 1; \
InverseFilter2D[b_List, gamma_] := 
  Map[(If[ Abs[#] gamma >= 1, 1/#, If[ # == 0, gamma, gamma Abs[#] / #]])&,
      FFT[b], {2}];
\end{mathematica}
%--------------------------------------------------------------------
\subsection{Restoration}

Function \texttt{Mean} calculates the average of list.
\begin{mathematica}
Mean[l_List] := N[Apply[Plus, l] / Length[l]];
\end{mathematica}
Function \texttt{Variance} calculates the variance of a list.
\begin{mathematica}
Variance[l_List] := Mean[ (l - Mean[l])^2 ];
\end{mathematica}
All previous efforts are summarized in the following two functions:
\texttt{ModifiedRestoration}:
\begin{mathematica}
ModifiedRestoration2D[g_List, Pv_, alpha_, beta_, b_List, gamma_] := 
    Module[{mf, G, Pf},
    mf = Mean[Flatten[g]];
    G = FFT[g - mf];
    Pf = Periodogram[g - mf];
    IFFT[ InverseFilter2D[b, gamma] WienerFilter[Pf, Pv, alpha, beta] G ] + mf 
];
\end{mathematica}
and \texttt{AdaptiveRestoration}:
\begin{mathematica}
AdaptiveRestoration2D[g_List, M_Integer, Pv_, b_List, gamma_] := 
    Module[{mf, G, Pf},
    mf = Mean[Flatten[g]];
    G = FFT[ AdaptiveWienerFilter2D[g - mf, M, Pv] ];
    IFFT[ InverseFilter2D[b, gamma] G ] + mf 
];
\end{mathematica}
%--------------------------------------------------------------------
\section{Tests}
\subsection{Image ``square''}

To be able to finish the project before the deadline I have to operate
on downscaled images what significantly reduces the computational
complexity but not at all the need for deep understanding of what is
going on.
\begin{mathematica}
sq32 = ReadPGM["square32"];
\end{mathematica}
The variance of the big image is exactly the same.
\begin{mathematica}
sqvar = Variance[Flatten[sq32]]
.
Out[14]= 585.937
\end{mathematica}
The desired values of standard deviation for noise generation purposes
can be obtained for $\mbox{SNR} = 0.1$, $1$, $5$, $10$, and $\infty$\,dB from
the formula:
\begin{mathematica}
sqdevs = Reverse[Sqrt[sqvar / 10.^({0.1, 1, 5, 10, Infinity} / 10.)]]
.
Out[15]= {0, 7.65466, 13.6121, 21.5738, 23.9291}
\end{mathematica}
Function \texttt{SNR} calculates signal-to-noise ratio for two sequences:
\begin{mathematica}
SNR[f_List, v_List] := 10. Log[10., Variance[f] / Variance[v]];
\end{mathematica}
Let's check the resulting value of SNR for noise with assumed SNR equal
$5$\,dB.
\begin{mathematica}
SNR[Flatten[sq32], Flatten[Table[GaussianRandom[0, sqdevs[[3]], 12], {32},{32}]]]
.
Out[17]= 5.06572
\end{mathematica}
Function \texttt{SNRImprovement} will be useful for comparison of restoration
results:
\begin{mathematica}
SNRImprovement[f_List, g_List, p_List] := Module[
    {fgvar = Variance[f - g], fpvar = Variance[f - p]},
    If [fpvar != 0, 10. Log[10., fgvar / fpvar], If[fgvar == 0, 0, Infinity]]
];
.
General::spell1: 
   Possible spelling error: new symbol name "fpvar"
     is similar to existing symbol "fgvar".
\end{mathematica}
Function \texttt{Pad} is used for creating shifted Gaussian kernel
sequence of arbitrary length $size$.
\begin{mathematica}
Pad[l_List, size_Integer] := RotateLeft[ 
    Join[ l, Table[0.0, {size - Length[l]}]],
    Floor[Length[l]/2]
];
\end{mathematica}
The following functions add blur and noise to the 1-D sequence.
\begin{mathematica}
Blur2D[l_List, 0] := l; \
Blur2D[l_List, b_] := ConvolveSeparable2D[l, GaussianKernel[b]];\
BlurAndNoise2D[l_List, b_, 0] := Blur2D[l, b]; \
BlurAndNoise2D[l_List, b_, n_] :=
    Map[ GaussianRandom[0, n, 12] + #&, Blur2D[l, b], {2} ];
\end{mathematica}
I will operate on three blur sequences:
\begin{mathematica}
blurs = {0, 3, 5};
\end{mathematica}

Using the above function one can create a set of disturbed images for
all noise values and three blur kernel widths which is shown in Figure~\ref{2disturbed}:
\begin{mathematica}
disturbed2d = Map[ BlurAndNoise2D[sq32, #[[2]], #[[1]]]&, 
    Outer[List, sqdevs, blurs], {2}];
\end{mathematica}
\begin{comment}
\begin{mathematica}
bigpage2D[d_List, dn_String] := MapIndexed[
WritePGM[Transpose @ Reverse @ Round[#1], dn <> ToString[#2[[1]]] <> "_" <> ToString[#2[[2]]] ]&, d, {2}];
\end{mathematica}
\begin{mathematica}
bigpage2D[disturbed2d, "2disturbed"];
\end{mathematica}
\end{comment}
\bigpagexv{2disturbed}{Set of disturbed images: in columns blur kernel
    widths are 0, 3 and 5, and in rows noise levels are $\mbox{SNR} =
     \infty$, $10$, $5$, $1$, and $0.1$\,dB, accordingly (so in left
     upper is an undisturbed image).}

Function \texttt{bigimprovement} allows for numerical analysis of the
improvement over the whole set of images presented on one page.
\begin{mathematica}
bigimprovement2D[result_List] :=
    Chop[MapThread[ SNRImprovement[Flatten[sq32], Flatten[#1], Flatten[#2]]&, {disturbed2d, result}, 2]];
\end{mathematica}
%--------------------------------------------------------------------
\subsubsection{Wiener filtering}
Classical Wiener filter uses parameters $\alpha = \beta = 1$.
\begin{mathematica}
wiener2D[gamma_] := Chop[MapIndexed[ ModifiedRestoration2D[#,
  sqdevs[[#2[[1]]]], 1., 1., 
  If[# != 0, ((Outer[Times,#,#])& @ Pad[GaussianKernel[#], 32]), {1}]& [ blurs[[#2[[2]]]] ],
  gamma] &, disturbed2d, {2}]];
\end{mathematica}
\begin{comment}
\begin{mathematica}
bigpage2D[wiener2 = wiener2D[0.5], "2wiener"];
\end{mathematica}
\begin{mathematica}
bigpage2D[wiener2a = wiener2D[1.], "2wienera"];
\end{mathematica}
\end{comment}
\bigpagexv{2wiener}{Results of Wiener filtering for $\gamma = 0.5$.}
\bigpagexv{2wienera}{Results of Wiener filtering for $\gamma = 1$.}

Below there are SNR improvements for the Figure~\ref{2wiener}:
\begin{mathematica}
bigimprovement2D[wiener2]//MatrixForm
.
Out[171]//MatrixForm= -Infinity   -46.4409    -10.5648

                      0.566196    -4.62314    -4.05122

                      0.421714    -0.288329   -0.433644

                      0.27191     2.63644     2.1426

                      0.277302    2.86994     2.76769
\end{mathematica}
and for Figure~\ref{2wienera}:
\begin{mathematica}
bigimprovement2D[wiener2a]//MatrixForm
.
Out[173]//MatrixForm= -Infinity   0           0

                      0.566196    0.623219    0.388761

                      0.421714    0.407788    0.357506

                      0.27191     0.300504    0.281664

                      0.277302    0.265323    0.264721
\end{mathematica}
Again, as in the 1-D case, we can see that the choice of the parameter
$\gamma$ should be very careful.
%--------------------------------------------------------------------
\subsubsection{Modified Wiener filtering}

Function \texttt{ModifiedImprovement} returns the value of improvement
for modified Wiener restoration.
\begin{mathematica}
ModifiedImprovement2D[f_List, g_List, alpha_, beta_, gamma_, sdev_, blur_] := 
    SNRImprovement[Flatten[f], Flatten[g], Flatten[ModifiedRestoration2D[g, sdev, alpha, beta, 
        If[# != 0, ((Outer[Times,#,#])& @ Pad[GaussianKernel[#], 32]), {1}]& [ blur ], gamma]]];
\end{mathematica}
We can check that it in fact yields results consistent with the above
table entry for kernel width 3, $\mbox{SNR} = 10$ and $\gamma = 1$:
\begin{mathematica}
ModifiedImprovement2D[sq32, disturbed2d[[2,2]], 1., 1., 1., 
  sqdevs[[2]], blurs[[2]]]
.
Out[206]= 0.623219
\end{mathematica}

Function \texttt{MaxImprovement2D} finds parameters $\alpha$ and $\beta$
guaranteeing the best restored sequence for disturbed sequence $g$.
\begin{mathematica}
MaxImprovement2D[a_Integer, b_Integer] := 
{-#1, {alpha, beta} /. #2}& @@ 
    FindMinimum[ -ModifiedImprovement2D[sq32, disturbed2d[[a,b]], 
    Abs[alpha], Abs[beta],
    1., sqdevs[[a]], blurs[[b]]], {alpha, {2.1, 2.15}}, {beta, {0.7,
        0.75}},  AccuracyGoal->2, PrecisionGoal->2]
.
General::spell1: 
   Possible spelling error: new symbol name "blurs"
     is similar to existing symbol "blur".
\end{mathematica}
Here we calculate optimal parameters for the whole page of sequences.
Please notice that the sequences without noise are omitted because their
reconstructions do not depend on $\alpha$ and $\beta$.  Answers are
returned in the format $\{\mbox{SNR improvement}, \{\alpha, \beta\}\}$:
\begin{mathematica}
(modified2d = Map[MaxImprovement2D @@ # &, 
    Outer[List, {2,3,4,5}, {1,2,3}], {2}] ) // MatrixForm
. 
    3.0361                2.99721               1.78048
    {24.7024, 0.900068}   {23.7896, 0.900669}   {3.13007, 4.07944}

    4.06602               4.08555               3.50855
    {57.4379, 0.876143}   {58.6562, 0.897299}   {20.7739, 2.02522}

    5.62798               5.81043               4.85896
    {49.8611, 1.90571}    {50.9752, 1.91215}    {126.43, 0.930217}

    6.26383               6.31742               6.03307
    {166.263, 0.915281}   {177.637, 0.909385}   {60.8015, 2.03773}
\end{mathematica}
The above values show that quite big improvements are possible, though
for rather unexpected values of $\alpha$ and $\beta$.

The following function generates the optimal restored sequences for
$\gamma = 1$.
\begin{mathematica}
modifiedwiener2 = Chop[MapIndexed[ 
  ModifiedRestoration2D[#, sqdevs[[#2[[1]]+1]], 
  modified2d[[#2[[1]], #2[[2]]]] [[2,1]], modified2d[[#2[[1]], #2[[2]]]] [[2,2]],
  If[# != 0, Outer[Times, #, #]& @ Pad[GaussianKernel[#], 32], {1}]& [ blurs[[#2[[2]]]] ],
  10.] &, disturbed2d[[{2,3,4,5}]], {2}]];
\end{mathematica}
\begin{comment}
\begin{mathematica}
bigpage2D[modifiedwiener2, "2modifiedwiener"];
\end{mathematica}
\end{comment}
\bigpagexii{2modifiedwiener}{Results of optimal modified Wiener filtering
  for $\gamma = 1$ (images with blur only were omitted).}

%--------------------------------------------------------------------
\subsubsection{Adaptive restoration}

   The following function calculates adaptive-filtered sequences for the
   adaptation area width $2n+1$ and inverse filter parameter $\gamma$.
   Results are shown in Figures~\ref{2adaptive} and~\ref{2adaptivea}.
   Values $\gamma = 10$ and 1000 were selected to show the increase of the
   high frequency components in the second case.  For $\gamma < 10$ the
   inverse filtering is not effective enough.
\begin{mathematica}
adaptive2D[n_, gamma_] := Chop[MapIndexed[ AdaptiveRestoration2D[#, n,
  sqdevs[[#2[[1]]]], 
  If[# != 0, Outer[Times, #, #]& @ Pad[GaussianKernel[#], 32], {1}]& [ blurs[[#2[[2]]]] ],
  gamma] &, disturbed2d, {2}]];
\end{mathematica}
\begin{comment}
\begin{mathematica}
bigpage2D[adaptive2 = adaptive2D[2, 0.5], "2adaptive"];
\end{mathematica}
\begin{mathematica}
bigpage2D[adaptive2a = adaptive2D[2, 1.], "2adaptivea"];
\end{mathematica}
\end{comment}
\bigpagexv{2adaptive}{Results of adaptive filtering for adaptation width 5
  and $\gamma = 0.5$.}
\bigpagexv{2adaptivea}{Results of adaptive filtering for adaptation width 11
  and $\gamma = 1$.}

Below there are SNR improvements for the Figure~\ref{2adaptive}:
\begin{mathematica}
bigimprovement2D[adaptive2]//MatrixForm
.
Out[194]//MatrixForm= -Infinity   -46.4409    -10.5648

                      0.80978     -4.58553    -4.02632

                      0.54151     -0.252478   -0.410964

                      0.328233    2.67199     2.16552

                      0.337172    2.90519     2.78797
\end{mathematica}
and for the Figure~\ref{2adaptivea}:
\begin{mathematica}
bigimprovement2D[adaptive2a]//MatrixForm
.
Out[196]//MatrixForm= -Infinity   0           0

                      0.80978     0.927417    0.603247

                      0.54151     0.525911    0.452338

                      0.328233    0.36587     0.339291

                      0.337172    0.31797     0.306699
\end{mathematica}
As one can see, adaptive filter performed worse in this case.

Summarizing: the best performace for the most disturbed image and
$\gamma = 1$ had modified Wiener filter, the next was adaptive filter
and the worst results were obtained from normal Wiener filter.

\end{comment}

\subsection{"Jet" image restoration.}

I will operate on an with lower resolution.  The original is of very low
quality so it will not hurt much.
\begin{mathematica}
jet = ReadPGM["jet128"];
\end{mathematica}
Here I create a disturbed version by adding noise with standard
deviation 10 and blurring it with the Gaussian kernel having width 5.
\begin{mathematica}
jetdisturbed = Abs @ BlurAndNoise2D[jet, 5, 10.];
\end{mathematica}
\begin{comment}
\begin{mathematica}
WritePGM[Round[jetdisturbed], "jetdisturbed"];
\end{mathematica}
\end{comment}
I assume $\gamma=1$ for all restoration methods:
\begin{mathematica}
jetwiener = ModifiedRestoration2D[jetdisturbed, 10., 1., 1.,
  (Outer[Times,#,#])& @ Pad[GaussianKernel[5], 128], .5 ];
\end{mathematica}
\begin{comment}
\begin{mathematica}
WritePGM[Round[jetwiener], "jetwiener"];
\end{mathematica}
\end{comment}
Let's find the SNR improvement:
\begin{mathematica}
SNRImprovement[Flatten[jet], Flatten[jetdisturbed], Flatten[jetwiener]]
.
Out[92]= 0.26776
\end{mathematica}
%================
First we find optimal coordinates:
\begin{mathematica}
jetmaximprovement = {-#1, {alpha, beta} /. #2}& @@ 
    FindMinimum[ -SNRImprovement[Flatten[jet], Flatten[jetdisturbed], 
    Flatten[ModifiedRestoration2D[jetdisturbed, 10., Abs[alpha], Abs[beta], 
        (Outer[Times,#,#])& @ Pad[GaussianKernel[5], 128], 1.]]],
    {alpha, {16, 17}}, {beta, {4, 5}},  AccuracyGoal->1,
    PrecisionGoal->1]
.
{0.749258, {2.54091, 4.17072}}
\end{mathematica}
And then use them for restoration:
\begin{mathematica}
jetmodified = ModifiedRestoration2D[jetdisturbed, 10., 
    jetmaximprovement[[2,1]], jetmaximprovement[[2,2]],
    (Outer[Times,#,#])& @ Pad[GaussianKernel[5], 128], 1. ];
\end{mathematica}
\begin{comment}
\begin{mathematica}
WritePGM[Abs @ Round[jetmodified], "jetmodified"];
\end{mathematica}
\end{comment}
Let's double-check the SNR improvement:
\begin{mathematica}
SNRImprovement[Flatten[jet], Flatten[jetdisturbed], Flatten[jetmodified]]
.
Out[107]= 0.749258
\end{mathematica}
%================
For adaptive restoration I will use the $7 \times 7$ region:
\begin{mathematica}
jetadaptive = AdaptiveRestoration2D[jetdisturbed, 3, 10., 
  (Outer[Times,#,#])& @ Pad[GaussianKernel[5], 128], 1. ];
\end{mathematica}
\begin{comment}
\begin{mathematica}
WritePGM[Round[jetadaptive], "jetadaptive"];
\end{mathematica}
\end{comment}
Let's find the SNR improvement:
\begin{mathematica}
SNRImprovement[Flatten[jet], Flatten[jetdisturbed], Flatten[jetadaptive]]
.
Out[28]= 0.227241
\end{mathematica}

The results are shown in the figures~\ref{jet128} and~\ref{jetwiener}.
\figii{jet128}{jetdisturbed}{Image ``jet'': (a)~original; (b)~disturbed by
  adding noise with standard deviation 20 and blurring with Gaussian
  kernel having width 5.}
\figiii{jetwiener}{jetmodified}{jetadaptive}{Restored image ``jet'':
  (a)~Wiener restoration ($\Delta\mbox{SNR}=0.27\,dB$); (b)~modified Wiener
  restoration($\Delta\mbox{SNR}=0.75\,dB$); (c)~adaptive restoration
  ($\Delta\mbox{SNR}=0.23\,dB$).}

\section{Conclusions}
As in 1-D case the optimal modified Wiener restoration was the best
method.  The two other methods probably also are capable of giving good
results but they should be optimized with respect to parameters $\gamma$
and adaptation width.  The sensitivity to $\gamma$ is much better
visible in the 2-D case, bad choice causes stronger degradation of an
image.  Results obtained from the two last filters are very similar.

\end{document}
%--------------------------------------------------------------------