<jupyter_start><jupyter_text>Week7  Image Recovery (Part 2/2)The prior knowledge is incorporated here into Restoration in a different way (as a sample of a random field). We'll have knowledge of two random fields and work with those. We'll cover Wiener filters and other type of Estimation models. Also, we'll discuss the Bayesian Inference.  Weiner Restoration FilterIn order to define the Filter we only need the auto-correlation function of an image and the cross-correlation function between two images or the Fourier transfer, which we refer to as the Power Spectra.We'll look into stochastic approaches:Wiener (or Wiener-Kolmogorov) FilterBayesian Formulation: ML (Maximum likelihood), MAP, Hierarchical Bayesian.Let's start with definitions:Autocorrelation:$$R_{f f}(i,j,k,l)=E|f(i,j)f^*(k,l)|$$For a good estimationm I need a high number of realizations, i.e. images, (called ensemble). Another good definition is Wide Sense Stationarity $WWS$:$$R_{f f}(i,j,k,l)=R_{f f}(i-k,j-l) = R_{f f} (n_1,n_2)$$It is irrelevant where the axes are located. I just want to know the distance among the images. The $R_{f f}$ model can be choosen with help of literature. For example, the Isotropic Exponential Decay model:$$R(n_1,n_2) = c \cdot \gamma^{-(|n_1|+|n_2|)}$$As we usually have a single image, we'll use the notion of Ergodicity, that tells that Ensemble Average equals Spatial Average:$$R_{f f}(n_1,n_2)=\lim_{N\rightarrow \infty} \frac{1}{(2N+1)^2} \sum_{k_1=-N}^N \sum_{k_2=-N}^N f(k_1,k_2) f^*(k_1-n_1,k_2-n_2)$$Finally, the Power Spectrum (Fourier Transform ($F$) of AutoCorrelation function):$$P_{f f} (\omega_1,\omega_2) = F\{ R_{f f}(n_1,n_2) \}$$In the Non-Causal Wiener Restoration Filter, one Degradation Model is, with $w$ the additive noise:$$y(i,j) = h(i,j) * * f(i,j)+w(i,j)$$The objective is to find $\hat{f}$ that minimize the error:$$\hat{f}(i,j)=argmin_{f(i,j)} E\left[ |e(i,j)^2| \right] =argmin_{f(i,j)} E\left[ |f(i,j)-\hat{f}(i,j)|^2 \right] $$We need the expectation $E$ because $f$ is the realization of a random field of 2d random process... so we need a parameter for all the realizations. Consider an LSI restoration filter $r$ for which the error is the least as possible one:$$\hat{f}(i,j) = r(i,j)* * y(i,j) $$If $f(i,j)$ is $WSS$ with Autocorrelation $R_{f f}(i,j)$ :$$R_{yy}(i,j)=R_{f f}(i,j)* * h(i,j)* * h^* (-i,-j)$$$$P_{yy}(\omega_1,\omega_2)=|H(\omega_1,\omega_2)|^2 \cdot P_{f f}(\omega_1,\omega_2)$$Cross-Correlations input-output and output-input are: respectively:$$R_{fy}(i,j)=R_{f f}(i,j)* * h^* (-i,-j) \hspace{2cm} P_{fy}(\omega_1,\omega_2)=H^*(\omega_1,\omega_2) \cdot P_{f f}(\omega_1,\omega_2)$$$$R_{yf}(i,j)=R_{f f}(i,j)* * h(i,j) \hspace{2cm} P_{y f}(\omega_1,\omega_2)=H(\omega_1,\omega_2) \cdot P_{f f}(\omega_1,\omega_2)$$The Wiener Restoration Filter has $f(i,j)$ and $w(i,j)$ which are WSS. It is based on Orthogonally Principle: the Error is Orthogonal to the data:$$E\left[ e(i,j) y^*(k,l) \right] = 0$$which leads to:$$R_{f y}(i,j)=R_{yy}(i,j)* * r(i,j) \hspace{1cm} R(\omega_1, \omega_2)=\frac{P_{fy}(\omega_1,\omega_2)}{P_{yy}(\omega_1,\omega_2)}$$Common assumptions are that $f(i,j)$ and $w(i,j)$ are uncorrelated, therefore $E|f(i,j)w^*(i,j)|=E|f(i,j)| \cdot E|w^*(i,j)|$. The other assumption is that both $f(i,j)$ and $w(i,j)$, i.e. image and noise, are zero mean. These will give:$$P_{f y}(\omega_1,\omega_2)=H^*(\omega_1,\omega_2) \cdot P_{f f}(\omega_1,\omega_2)$$and$$P_{y y}(\omega_1,\omega_2)=|H(\omega_1,\omega_2)|^2 \cdot P_{f f}(\omega_1,\omega_2) + P_{w w} (\omega_1, \omega_2)$$Eventually:$$R(\omega_1, \omega_2) = \frac{H^*(\omega_1,\omega_2) \cdot P_{f f}(\omega_1,\omega_2)}{|H(\omega_1,\omega_2)|^2 \cdot P_{f f}(\omega_1,\omega_2) + P_{w w} (\omega_1, \omega_2)}$$  Wiener v.s. Constrained Least-Squares Restoration Filter They differ only in one term: the Stabilizing Term. If we model Power Spectra (function of stabilizing term) and use the white and gaussian model, the results of Wiener and CLS becomes identical. If noise is white, we have:$$R(\omega_1, \omega_2) = \frac{H^*(\omega_1,\omega_2) }{|H(\omega_1,\omega_2)|^2  + \frac{P_{w w} (\omega_1, \omega_2)}{P_{f f}(\omega_1,\omega_2)}} = \frac{H^*(\omega_1,\omega_2) }{|H(\omega_1,\omega_2)|^2  + \frac{\sigma^2_{w w} }{P_{f f}(\omega_1,\omega_2)}}$$The CLS was instead:$$R_{CLS}(\omega_1, \omega_2) = \frac{H^*(\omega_1,\omega_2) }{|H(\omega_1,\omega_2)|^2  + \alpha | C(\omega_1,\omega_2)|^2}$$That means, by substitution, that we get identity if:$$\alpha = \sigma^2_{ww} \hspace{2cm} |C(\omega_1,\omega_2)|^2=\frac{1}{P_{f f} (\omega_1,\omega_2)}$$$|C|^2$ and $1/P$ are magnitudes of frequency response of High Pass Filter.<jupyter_code>%matplotlib inline
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
from scipy import signal, misc, ndimage
import numpy as np
import matplotlib.pyplot as plt
import os
import PIL
from PIL import ImageFilter
import skimage
import cv2

I = misc.imread(os.path.join('media','bear.png') );
# read in the original, sharp and noise-free image
original = misc.imread(os.path.join('media','bear.png') , mode='L').astype(float)/255
H = np.size(original,0)
W = np.size(original,1)
x = range(W)
y = range(H)
[X, Y] = np.meshgrid(x,y)
whitenoise = np.random.rand(H,W)

Pff = np.power(np.abs(np.fft.fftshift(np.fft.fft2(original))),2)
Pww = np.power(np.abs(np.fft.fftshift(np.fft.fft2(whitenoise))),2)

fig=plt.figure(figsize=(8, 8))
ax1=plt.subplot(2, 2, 1)
plt.imshow(original,cmap='gray')
plt.title('Original')

ax2 = fig.add_subplot(222, projection='3d')
ax2.plot_surface(Y,X,Pff, vmin=0, vmax=1e3, cmap='viridis', linewidth=0.1)
plt.title('$P_{ff}(\omega_1,\omega_2)$')
ax2.set_zlim([0, 1e4]);

ax3 = fig.add_subplot(223, projection='3d')
ax3.plot_surface(Y,X,np.log10(Pww), vmin=0, vmax=4, cmap='viridis', linewidth=0.1)
plt.title('$P_{ww}(\omega_1,\omega_2)$')
ax3.set_zlim([0, 5]);

ax4 = fig.add_subplot(224, projection='3d')
ax4.plot_surface(Y,X,np.log10(Pww/Pff), vmin=0, vmax=3, cmap='viridis', linewidth=0.1)
plt.title('$P_{ww}(\omega_1,\omega_2)/P_{ff}(\omega_1,\omega_2)$')
ax4.set_zlim([0, 5]);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1975530179445351, 'fraction_numerical': 0.032137030995106036, 'mean_word_length': 5.842633928571429, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '258893', 'license': 'MIT', 'n_tokens_mistral': 2525, 'n_tokens_neox': 2419, 'n_words': 610, 'path': 'week7.ipynb', 'repo_name': 'JulesAM/Fundamentals-of-Digital-Image-and-Video-Processing-course', 'size': '6130'}