9.4

'''
Nad slikom snimak.png simulirati efekat degradacije atmosferske turbulencije definisane modelom:

gde su  (u,v)  prostorne koordinate DFT slike sa centrom u sredini slike a  k  parametar degradacije. Koristiti vrednosti  k=0.0025 .

Nad degradiranom slikom dodati AWGN varijanse 1 i pokušati ukloniti degradaciju inverznim filtriranjem a potom koristeći Wiener-ov filtar. Prikazati kako izgledaju karakteristike degradacije i Wiener-ovog filtra.
-------------------------------------------------------
Napomena: U okviru računanja fft2 koristiti simetrično proširenje slike.
'''
import plotly.express as px
from skimage import io
import numpy as np

img = io.imread('snimak.png')
fig = px.imshow(img, zmin=0, zmax=255, color_continuous_scale='gray')
fig.show()

pad = 128
img_p = np.pad(img, pad_width = pad, mode='reflect')

IMG = np.fft.fftshift(np.fft.fft2(img_p))
fig = px.imshow(np.log(np.abs(IMG)+1), color_continuous_scale='gray')
fig.show()

k = 0.0025

u = np.arange(IMG.shape[0]).reshape(-1, 1) - np.floor(IMG.shape[0]/2)
v = np.arange(IMG.shape[1]) - np.floor(IMG.shape[1]/2)

H = np.exp(-k*(u**2+v**2)**(5/6))

fig = px.imshow(np.log(np.abs(H)+1), color_continuous_scale='gray')
fig.show()

IMG1 = IMG * H
fig = px.imshow(np.log(np.abs(IMG1)+1), color_continuous_scale='gray')
fig.show()

img1 = np.real(np.fft.ifft2(np.fft.ifftshift(IMG1)))
img1 = img1[pad:-pad,pad:-pad]

fig = px.imshow(img1,zmin=0, zmax=255, color_continuous_scale='gray')
fig.show()

rng = np.random.default_rng()
noise = rng.standard_normal(size=img1.shape)*np.sqrt(1) + 0

img_noisy = img1 + noise

fig = px.imshow(img_noisy, zmin=0, zmax=255, color_continuous_scale='gray')
fig.show()

img_noisy_p = np.pad(img_noisy, pad_width = pad, mode='reflect')

IMG_NOISY = np.fft.fftshift(np.fft.fft2(img_noisy_p))

fig = px.imshow(np.log(np.abs(IMG_NOISY)+1), color_continuous_scale='gray')
fig.show()

#INVERZNO FILTRIRANJE

IMG_INVF = IMG_NOISY/H

img_invf = np.real(np.fft.ifft2(np.fft.ifftshift(IMG_INVF)))
img_invf = img_invf[pad:-pad,pad:-pad]

fig = px.imshow(img_invf, color_continuous_scale='gray')
fig.show()

K = 1/img.var()

W = np.conj(H)/(np.abs(H)**2+K)

fig = px.imshow(np.log(np.abs(W)+1), color_continuous_scale='gray')
fig.show()

IMG_INVW = IMG_NOISY * W

img_invw = np.real(np.fft.ifft2(np.fft.ifftshift(IMG_INVW)))
img_invw = img_invw[pad:-pad,pad:-pad]

fig = px.imshow(img_invw, zmin=0, zmax=255, color_continuous_scale='gray')
fig.show()