%---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}

\end{comment}

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

\begin{titlepage}

\author{Steven W.\ Roberts \and Tomasz J.\ Cholewo}
\title{EE\,618 Project \#5\\
  On Image Restoration
}
\date{\today}
\maketitle
\vfil
\tableofcontents
\vfil
\thispagestyle{empty}
\end{titlepage}
%====================================================================
\chapter{Power Spectral Estimation}

%--------------------------------------------------------------------
\section{One-dimensional Power Spectral Estimation}

\subsection{Periodogram method}

  We have to adapt \Mma{} \texttt{Fourier} function which is defined as
  \[ L_s = \frac{1}{\sqrt{n}} \sum_{r=1}^n l_r e^{2 j \pi (r-1)(s-1)/n}
  \]
  to the more common conventions (real sequences assumed):
\begin{mathematica}
FFT[l_List] := 
  Sqrt[N[Apply[Times, Dimensions[l]]]] Conjugate[ Chop[Fourier[N[l]]] ]; \
IFFT[L_List] := 
  Chop[1./Sqrt[N[Apply[Times, Dimensions[L]]]] InverseFourier[ Conjugate[N[L]]] ];
\end{mathematica}
Normalized periodogram of a sequence is then defined as:
\begin{mathematica}
Periodogram[l_List] := Abs[ FFT[l] ]^2 / Apply[Times, Dimensions[l]];
\end{mathematica}

\subsection{Periodogram smoothed with Hamming window}
Hamming window for $n$ samples can be obtained from:
\begin{mathematica}
HammingWindow[n_Integer] := Table[
    0.54 + 0.46 Cos[N[Pi] (k / ((n - 1) / 2) - 1)],
    {k, 0, n - 1}
];
\end{mathematica}
Its plot for $n = 256$ is shown in Figure~\ref{hamming}.
\begin{comment}
\begin{mathematica}
PSTeX[ListPlot[HammingWindow[256], ojoin], "hamming"];
\end{mathematica}
\end{comment}
\figi{hamming}{Hamming window for $n = 256$}
  The following function calculates the smoothed version of the
  periodogram.
\begin{mathematica}
SmoothedPeriodogram[l_List] := Periodogram[ HammingWindow[Length[l]] l ];
\end{mathematica}


\subsection{Tests}

  Function \texttt{test} returns a list of $n$ samples of the given
  function obtained by sampling it with frequency $f$ and with
  added noise multiplied by coefficient $k$.
\begin{mathematica}
test1d[f_, k_, n_] := Table[ 
    20 Cos[2 N[Pi] 0.2 i/f] + 20 Cos[2 N[Pi] i/f] + 20 Cos[2 N[Pi] 1.2 i/f] + 
        k GaussianRandom[0, 0.5, 12],
    {i, 0, n - 1}
];
\end{mathematica}
The highest frequency in the given signal is equal to $1.2$, so each
$f \geq 2.4$ will be appropriate.  Value $f_s = 128/50 = 2.56$
is chosen to put expected spikes right on grid. 
\begin{mathematica}
fs = 128/50.0; 
\end{mathematica}
Resulting plots for $k = 0$, 5 and 10 are generated below and shown in
Figure~{\ref{s0s5s10}}.
Function \texttt{labs} generates a list of general-purpose labels.
\begin{mathematica}
labs[n_, v_] := Table[ {i, i v/n }, {i, 0, n , n / 4}]
\end{mathematica}
\begin{mathematica}
PSTeX[ListPlot[test1d[fs, #, 256], ojoin, oall,
    Ticks->{labs[256, 256/fs], Automatic}, AxesLabel->{"t", "x(t)"}
], "s" <> ToString[#]]& /@ {0, 5, 10};
\end{mathematica}
\pageiii{s0}{s5}{s10}{Signal $x(t)$ sampled with frequency $f_s = 2.56$
  and 256~samples: (a)~$k=0$; (b)~$k=5$; (c)~$k=10$.}

Results of power spectral estimation are shown in Figures~\ref{p0n128}
and~\ref{p0n256}. 
\begin{mathematica}
PSTeX[ListPlot[Periodogram[test1d[fs, #[[1]], #[[2]]]], ojoin, oall,
  Ticks->{labs[#[[2]], fs], Automatic}, AxesLabel->{"f [Hz]", "S(f)"}
], "p" <> ToString[#[[1]]] <> "n" <> ToString[#[[2]]] ]& /@
    Distribute[{ {0, 5, 10}, {128, 256}}, List];
\end{mathematica}
\begin{mathematica}
PSTeX[ListPlot[SmoothedPeriodogram[test1d[fs, #[[1]], #[[2]]]], ojoin, oall,
    Ticks->{labs[#[[2]], fs], Automatic}, AxesLabel->{"f [Hz]", "S(f)"}
], "h" <> ToString[#[[1]]] <> "n" <> ToString[#[[2]]] ]& /@
    Distribute[{ {0, 5, 10}, {128, 256}}, List];
\end{mathematica}

\pagevi{p0n128}{h0n128}{p5n128}{h5n128}{p10n128}{h10n128}{Power spectral
  estimates obtained from periodogram (left column) and smoothed
  periodogram (right column) for 128~samples and for different levels of
  noise: (a),(b)~$k=0$; (c),(d)~$k=5$; (e),(f)~$k=10$.}
\pagevi{p0n256}{h0n256}{p5n256}{h5n256}{p10n256}{h10n256}{Power spectral
  estimates obtained from periodogram (left column) and smoothed
  periodogram (right column) for 256~samples and for different levels of
  noise: (a),(b)~$k=0$; (c),(d)~$k=5$; (e),(f)~$k=10$.}

%--------------------------------------------------------------------
\section{Image Power Spectral Estimation}
Functions defined so far for periodogram work also with images.  For
smoothed periodogram in the two-dimensional case it is necessary to
create a two-dimensional Hamming window:
\begin{mathematica}
HammingWindow2D[n_Integer] := 
    Outer[ Times, HammingWindow[n], HammingWindow[n]];
\end{mathematica}
Here we define 2-D smoothed periodogram.
\begin{mathematica}
SmoothedPeriodogram2D[l_List] := Periodogram[ l HammingWindow2D[Length[l]] ];
\end{mathematica}

\subsection{Tests}
The 2-D test sequence is defined as follows:
\begin{mathematica}
test2d[f_, n_] := Table[ 
    20 Cos[2 N[Pi] 0.2 j/f] + 20 Cos[2 N[Pi] i/f] + 
    20 Cos[2 N[Pi] 1.2 (i + j)/f],
    {i, 0, n - 1}, {j, 0, n - 1}
];
\end{mathematica}
Again, the best choice for sampling frequency is 2.56\,samples/sec but
it causes that plots again look a liitle bit crowdy, so the plot shown
in Figure~{test2d} is sampled with frequency 64\,samples/sec.
\begin{comment}
\begin{mathematica}
PSTeX[ListPlot3D[test2d[16, 32], AxesLabel->{"t1 [sec]", "t2 [sec]", "x(t1,t2)"}, 
    Ticks->{labs[128, 2.], labs[128, 2.], Automatic}
], "test2d"];
\end{mathematica}
\end{comment}
\figi{test2d}{Test image sampled with $f_s = 64\,\mbox{Hz}$.}
The following function creates a contour plot with all labels and ticks
used for obtaining graphical depictions of power spectral estimates.
\begin{mathematica}
cplot[l_List, fs_] := ListContourPlot[20. Log[10., l + zero], 
    ContourShading->False,
    ContourSmoothing->Automatic, FrameLabel->{"f1 [Hz]", "f2 [Hz]"}, 
    FrameTicks->{labs[Length[l], fs], labs[Length[l], fs]}
];
\end{mathematica}
\begin{comment}
\begin{mathematica}
PSTeX[cplot[Periodogram[test2d[fs, #]], fs], 
    "pp" <> ToString[#]]& /@ {64, 128}; \
PSTeX[cplot[SmoothedPeriodogram2D[test2d[fs, #]], fs], 
    "hh" <> ToString[#]]& /@ {64, 128};
\end{mathematica}
\end{comment}
The results are presented in Figure~\ref{pp64}.
\pageiv{pp64}{hh64}{pp128}{hh128}{Contour plots of power spectrum
  estimations for a test image, in the left column periodograms, in the
  right column smoothed periodograms: (a),(b)~$n = 64$; (c),(d)~$n=128$.}

Estimation is performed also for image ``square128''; in plots I assume
that it was obtained with sampling frequency 2~samples/sec.
\begin{mathematica}
sq128 = ReadPGM["square128"];
\end{mathematica}
This image is shown in Figure~\ref{sq}.
\begin{comment}
\begin{mathematica}
PSTeX[ListPlot3D[ Table[sq128[[i,j]],{i,1,128,4},{j,1,128,4}], AxesLabel->{"t1 [sec]", "t2 [sec]", "x(t1,t2)"}, 
    Ticks->{labs[128, 64.], labs[128, 64.], Automatic}
], "sq"];
\end{mathematica}
\end{comment}
\figi{sq}{Image ``square128''.}

\begin{mathematica}
PSTeX[cplot[Periodogram[sq128], 2.], "ppsq"]; \
PSTeX[cplot[SmoothedPeriodogram2D[sq128], 2.], "hhsq"];
\end{mathematica}
Results are shown in Figure~\ref{ppsq}.
\figii{ppsq}{hhsq}{Power spectral estimations of image ``square128'':
  (a)~periodogram; (b)~smoothed periodogram.}

%====================================================================
\chapter{Signal Restoration}
%--------------------------------------------------------------------
\section{One-dimensional Signal Restoration}

\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}
MeanList[l_List, 0] := l; \
MeanList[l_List, M_Integer] :=
    Table[ Sum[index[l, i + j], {j, -M, M}] / (2 M + 1), {i, 1, Length[l]}];
\end{mathematica}
Function \texttt{VarianceList} calculates the variance in a similar way.
It can be used for signal power estimation.
\begin{mathematica}
VarianceList[l_List, M_Integer] := MeanList[ (l - MeanList[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}
AdaptiveWienerFilter[g_List, _, 0] := g; \
AdaptiveWienerFilter[g_List, M_Integer, Pv_] := Module[{ave, var},
    ave = MeanList[g, M];
    var = MeanList[ (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}
InverseFilter[{1}, _] := 1; \
InverseFilter[b_List, gamma_] := 
  (If[ Abs[#] gamma >= 1, 1/#, If[ # == 0, gamma, gamma Abs[#] / #]])& /@ 
      FFT[b];
\end{mathematica}
Magnitude response of one exemplary filter is shown in
Figure~\ref{invfilt}.
\begin{mathematica}
PSTeX[ ListPlot[Abs[InverseFilter[Pad[GaussianKernel[17], 64], 10.]], 
  ojoin], "invfilt"];
\end{mathematica}
\figi{invfilt}{Magnitude of inverse filter for Gaussian kernel having
  width 17 padded to 64 samples and $\gamma = 10$.}
%--------------------------------------------------------------------
\subsection{Restoration}

Function \texttt{Mean} calculates the average of list.
\begin{mathematica}
Mean[l_List] := 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}
ModifiedRestoration[g_List, Pv_, alpha_, beta_, b_List, gamma_] := 
    Module[{mf, G, Pf},
    mf = Mean[g];
    G = FFT[g - mf];
    Pf = Periodogram[g - mf];
    IFFT[ InverseFilter[b, gamma] WienerFilter[Pf, Pv, alpha, beta] G ] + mf 
];
\end{mathematica}
and \texttt{AdaptiveRestoration}:
\begin{mathematica}
AdaptiveRestoration[g_List, M_Integer, Pv_, b_List, gamma_] := 
    Module[{mf, G, Pf},
    mf = Mean[g];
    G = FFT[ AdaptiveWienerFilter[g - mf, M, Pv] ];
    IFFT[ InverseFilter[b, gamma] G ] + mf 
];
\end{mathematica}
%--------------------------------------------------------------------
\subsection{Tests}

The given sequence can be calculated by the following function:
\begin{mathematica}
testrec = Table[ If[ n < 10, 0.0, If[n < 25, 0.9, If[n < 33, 0.0,
    If[n < 37, 0.5, If[n < 40, 0.0, If[n < 44, 0.55,
    If[n < 60, 0.0, If[n < 70, 0.1 (n - 60), If[n < 80, 0.1 (80 - n),
    If[n < 85, 0.0, If[n < 88, 0.25, If[n < 99, 0.0, 
    If[n < 102, 1.0, If[n < 122, 0.0, 0.0]]]]]]]]]]]]]],
    {n, 0, 127}
];
\end{mathematica}
It is shown in Figure~\ref{test1d}.
\begin{comment}
\begin{mathematica}
PSTeX[ListPlot[testrec, ojoin], "test1d"];
\end{mathematica}
\end{comment}
\figi{test1d}{Original test sequence for signal restoration.}
Let us calculate the variance of the test sequence:
\begin{mathematica}
testvar = Variance[testrec]
.
Out[43]= 0.129064
\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}
sdevs = Reverse[Sqrt[testvar / 10.^({0.1, 1, 5, 10, Infinity} / 10.)]]
.
Out[44]= {0, 0.113606, 0.202023, 0.320186, 0.355142}
\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 the check resulting value of SNR for noise with assumed SNR equal
$5$\,dB.
\begin{mathematica}
SNR[testrec, Table[GaussianRandom[0, sdevs[[3]], 12], {128}]]
.
Out[46]= 5.03385
\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]]
];
\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}
Blur[l_List, 0] := l; \
Blur[l_List, b_] := Convolve[l, GaussianKernel[b]];\
BlurAndNoise[l_List, b_, 0] := Blur[l, b]; \
BlurAndNoise[l_List, b_, n_] :=
    Map[ GaussianRandom[0, n, 12] + #&, Blur[l, b] ];
\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{disturbed}:
\begin{mathematica}
disturbed = Map[ BlurAndNoise[testrec, #[[2]], #[[1]]]&, 
    Outer[List, sdevs, blurs], {2}];
.
Outer::normal: Normal expression expected at position 2 in 
    Outer[List, sdevs, {0, 3, 5}].

Part::partd: Part specification 0[[2]] is longer than depth of object.

Part::partd: Part specification 0[[1]] is longer than depth of object.

Part::partd: Part specification 0[[1.]] is longer than depth of object.

General::stop: Further output of Part::partd
     will be suppressed during this calculation.

CompiledFunction::cfr: 
   Argument 0[[1]] at position 2 should be a machine-size real number.

CompiledFunction::cfr: 
   Argument 0[[1]] at position 2 should be a machine-size real number.

CompiledFunction::cfr: 
   Argument 3[[1]] at position 2 should be a machine-size real number.

General::stop: Further output of CompiledFunction::cfr
     will be suppressed during this calculation.

Outer::normal: Normal expression expected at position 2 in 
    Outer[List, sdevs, {Convolve[<<2>>], <<2>>}].
\end{mathematica}
\begin{comment}
\begin{mathematica}
bigpage[d_List, dn_String] := PSTeX[ Show[GraphicsArray[MapIndexed[
ListPlot[#1, ojoin, ohide, DefaultFont->{"Courier", 1.5}]&, d, {2}]],ohide], dn];
\end{mathematica}
\begin{mathematica}
bigpage[disturbed, "disturbed"];
\end{mathematica}
\end{comment}
\bigpage{disturbed}{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 in on one page.
\begin{mathematica}
bigimprovement[result_List] :=
    Chop[MapThread[ SNRImprovement[testrec, #1, #2]&, {disturbed, result}, 2]];
\end{mathematica}

%--------------------------------------------------------------------
\subsubsection{Wiener filtering}
Classic Wiener filter uses parameters $\alpha = \beta = 1$.
\begin{mathematica}
wiener[gamma_] := Chop[MapIndexed[ ModifiedRestoration[#,
  sdevs[[#2[[1]]]], 1., 1., 
  If[# != 0, Pad[GaussianKernel[#], 128], {1}]& [ blurs[[#2[[2]]]] ],
  gamma] &, disturbed, {2}]];
\end{mathematica}
\begin{comment}
\begin{mathematica}
bigpage[wiener1 = wiener[1.], "wiener1"];
.
Part::partd: Part specification sdevs[[3]] is longer than depth of object.

Part::partd: Part specification sdevs[[3]] is longer than depth of object.

Part::partd: Part specification sdevs[[3]] is longer than depth of object.

General::stop: Further output of Part::partd
     will be suppressed during this calculation.

Outer::normal: Normal expression expected at position 2 in Outer[List, <<2>>].

Outer::normal: Normal expression expected at position 2 in 
    Outer[List, sdevs, {<<3>>}].
\end{mathematica}
\begin{mathematica}
bigpage[wiener10 = wiener[10.], "wiener10"];
\end{mathematica}
\begin{mathematica}
bigpage[wiener20 = wiener[20], "wiener20"];
\end{mathematica}
\end{comment}
\bigpage{wiener1}{Results of Wiener filtering for $\gamma = 1$.}
\bigpage{wiener20}{Results of Wiener filtering for $\gamma = 20$.}

Below there are SNR improvements for the Figure~\ref{wiener1}:
\begin{mathematica}
bigimprovement[wiener1]//MatrixForm
.
Out[192]//MatrixForm= -Infinity   0           0

                      -0.822407   -1.54337    -0.984666

                      1.59139     1.41573     1.72932

                      3.39242     4.43484     3.48696

                      3.31186     3.6948      4.19306
\end{mathematica}
\begin{mathematica}
bigimprovement[wiener10]//MatrixForm
.
Out[299]//MatrixForm= -Infinity   252.349     288.924

                      -0.822407   -1.5385     -0.774209

                      1.59139     1.41414     0.971131

                      3.39242     4.42268     -0.083801

                      3.31186     3.68093     2.07221
\end{mathematica}
and for the Figure~\ref{wiener20}:
\begin{mathematica}
bigimprovement[wiener20]//MatrixForm
.
Out[198]//MatrixForm= -Infinity   252.349     288.924

                      -0.822407   -1.5385     -0.774209

                      1.59139     1.41414     0.971131

                      3.39242     4.42268     -0.083801

                      3.31186     3.68093     2.07221
\end{mathematica}
It is easy to observe that increasing $\gamma$ enhances the performance
of the inverse filter, because then, without noise, it is possible to
recover exactly the original sequence.  On the other hand, it magnifies
high-frequency noise.  So there is a need for a compromise.
%--------------------------------------------------------------------
\subsubsection{Modified Wiener filtering}

Function \texttt{ModifiedImprovement} returns the value of improvement
for modified Wiener restoration.
\begin{mathematica}
ModifiedImprovement[f_List, g_List, alpha_, beta_, gamma_, sdev_, blur_] := 
    SNRImprovement[f, g, ModifiedRestoration[g, sdev, alpha, beta, 
        If[# != 0, Pad[GaussianKernel[#], 128], {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 = 20$:
\begin{mathematica}
ModifiedImprovement[testrec, disturbed[[2,2]], 1., 1., 20., 
  sdevs[[2]], blurs[[2]]]
.
Out[208]= -1.5385
\end{mathematica}
The next statement produces the 3D plot of the improvement surface -- I
will concentrate on the sequence with no blur and largest noise, because
adding inverse filtering makes drawing any conlusions much more
difficult, and it was already considered.
\begin{mathematica}
PSTeX[Plot3D[ModifiedImprovement[testrec, disturbed[[5,1]], alpha, beta,
    10., sdevs[[5]], blurs[[1]]], {alpha, 0.01, 0.4}, {beta, 0.01, 50}, 
    PlotRange->All, PlotPoints->20, Evaluate[ohide], 
    AxesLabel->{"alpha", "beta", "SNR"}], "abplot"];
\end{mathematica}
\figi{abplot}{Improvement surface for the most noisy sequence.} 

As one can see, it is not easy to find the maximum point visually -- all
we can say is that it is in this case somewhere in the $[0,0.2] \times
[0,30]$ region.

Function \texttt{MaxImprovement} finds parameters $\alpha$ and $\beta$
guaranteeing the best restored sequence for disturbed sequence $g$.
\begin{mathematica}
MaxImprovement[a_Integer, b_Integer] := 
{-#1, {alpha, beta} /. #2}& @@ 
    FindMinimum[ -ModifiedImprovement[testrec, disturbed[[a,b]], 
    Abs[alpha], Abs[beta],
    10., sdevs[[a]], blurs[[b]]], {alpha, {2.1, 2.15}}, {beta, {0.7,
        0.75}},  AccuracyGoal->3, PrecisionGoal->4]
\end{mathematica}
Answers are returned in the format $\{\mbox{SNR improvement}, \{\alpha,
\beta\}\}$:
\begin{mathematica}
MaxImprovement[5,1]
.
Out[276]= {3.41718, {2.27545, 0.743329}}
\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$.
\begin{mathematica}
(modified = Map[MaxImprovement @@ # &, 
    Outer[List, {2,3,4,5}, {1,2,3}], {2}] ) // MatrixForm
.
    0.96551                0.702727               0.276047
    {2.22793, 0.176493}    {0.209424, 0.624361}   {0.0308363, 9.55338}

    1.97507                2.03174                1.0727
    {0.929686, 0.609115}   {0.681255, 0.688507}   {0.582628, 1.81705}

    3.5269                 4.51541                2.18909
    {3.57942, 0.563067}    {0.390333, 2.46967}    {2.25915, 1.6301}

    3.41718                4.39149                4.0596
    {2.27545, 0.743329}    {8.51625, 0.476593}    {0.295161, 6.49107}
\end{mathematica}
These results look quite dispersed but it is what we should expect if we
take into account the shape of the improvement surface shown in
Figure~\ref{abplot} and a random way of generating disturbed sequences.

The following function generates the optimal restored sequences for
$\gamma = 10$.
\begin{mathematica}
modifiedwiener = Chop[MapIndexed[ 
  ModifiedRestoration[#, sdevs[[#2[[1]]+1]], 
  modified[[#2[[1]], #2[[2]]]] [[2,1]], modified[[#2[[1]], #2[[2]]]] [[2,2]],
  If[# != 0, Pad[GaussianKernel[#], 128], {1}]& [ blurs[[#2[[2]]]] ],
  10.] &, disturbed[[{2,3,4,5}]], {2}]];
\end{mathematica}
\begin{comment}
\begin{mathematica}
bigpage[modifiedwiener, "modifiedwiener"];
\end{mathematica}
\end{comment}
\bigpage{modifiedwiener}{Results of optimal modified Wiener filtering
  for $\gamma = 10$ (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{adaptive2} and~\ref{adaptive5}.
   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}
adaptive[n_, gamma_] := Chop[MapIndexed[ AdaptiveRestoration[#, n,
  sdevs[[#2[[1]]]], 
  If[# != 0, Pad[GaussianKernel[#], 128], {1}]& [ blurs[[#2[[2]]]] ],
  gamma] &, disturbed, {2}]];
\end{mathematica}
\begin{comment}
\begin{mathematica}
bigpage[adaptive2 = adaptive[2, 10], "adaptive2"];
\end{mathematica}
\begin{mathematica}
bigpage[adaptive5 = adaptive[5, 1000], "adaptive5"];
\end{mathematica}
\end{comment}
\bigpage{adaptive2}{Results of adaptive filtering for adaptation width 5
  and $\gamma = 10$.}
\bigpage{adaptive5}{Results of adaptive filtering for adaptation width 11
  and $\gamma = 1000$.}

Below there are SNR improvements for the Figure~\ref{adaptive2}:
\begin{mathematica}
bigimprovement[adaptive2]//MatrixForm
.
Out[200]//MatrixForm= -Infinity   252.349     288.924

                      1.40534     0.943673    0.823419

                      3.15482     4.30491     2.14401

                      4.76648     5.04804     2.0267

                      4.18823     4.49735     3.29203
\end{mathematica}
and for the Figure~\ref{adaptive5}:
\begin{mathematica}
bigimprovement[adaptive5]//MatrixForm
.
Out[203]//MatrixForm= -Infinity   252.349     288.924

                      -0.458444   -1.19286    -1.17356

                      2.24477     2.37971     1.52011

                      4.9216      5.21165     2.4595

                      4.49917     5.63636     3.21166
\end{mathematica}
It is easy to observe that widening adaptation region can in some cases
help to recover from noise disturbance.  Further increase causes losses
because for so quickly changing sequence there is no point in adapting
to wide regions.

The given test sequence has a lot of high-frequency components resulting
from its rapid changes.  Application of the Wiener filter, which is
inherently lowpass in nature, can only worsen the SNR for small levels
of noise by smearing these edges.  For higher levels of noise it brings
more significant improvement.

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

%--------------------------------------------------------------------
\section{Image Restoration}

\end{document}
%--------------------------------------------------------------------
%--------------------------------------------------------------------
%--------------------------------------------------------------------
\begin{mathematica}
PSTeX[ListPlot[disturbed[[5,3]], ojoin], "xdis10"];
\end{mathematica}