# Image Restoration

## Cleaning images¶

### Here, I will show how the limb model introduced in curve_fitting can be used to derive the (point spread function) of the telescope, and correcting the images for instrumental smearing.¶

#### The 1D total point spread function of the telescope is given by the following expression:¶

$$PSF(x) = w_1 \times PSF_{PD} + w_2\times g_2(x) + w_3 \times g_3(x) + (1-w_2-w_2-w_3)$$

#### In Fourier space, the PSF turns into MTF with the following expression:¶

$$MTF(k) = w_1 \times MTF_{PD} + w_2\times G_2(k) + w_3 \times G_3(x)$$$$MTF(k=0) = (1-w_2-w_2-w_3)$$

#### Where $PSF_{PD}$ is the radial average of a 2D PSF retrieved from phase diversity reconstruction. The parameters $w_1, w_2, w_3$ are the weights and $g_2(x), g_3(x)$ are 2 Gaussian functions with widths of $s_2, s_3$, respectively. The latter quantities (weights and sigmas) were retrieved from the limb profile fitting explained in curve_fitting. So from now on, we will assume that they are already computed and can be passed to any value we want.¶

In [3]:
import numpy as np
import math as m
import scipy
import matplotlib.pyplot as plt
import pyfits

In [4]:
## Enter here the values of the PSF parameters:

w_1 = 0.528
w_2 =  0.289
w_3 = 0.101
w_4 = 1- w_1 -w_2 -w_3
s_2 =  0.721
s_3 = 2.818


#### STEP 1: to load the image to be corrected and turn it to a quadratic image.¶

In [6]:
## loading the input image:
path = '/home/fatima/Desktop/project_3/'
Input = pyfits.getdata(path+'sufi300_lev3_2009.fits', ignore_missing_end=True)

In [7]:
## Dimensions of input image
w_in = Input.shape[1]
h_in = Input.shape[0]
print w_in, h_in

714 1972


#### Now I will embed the initial image to have dimensions of 1972x1972, the pixels outside the input image will have a value of zero¶

In [8]:
Input_emb = np.zeros((h_in,h_in))
Input_emb[:,(h_in-w_in)/2:(h_in+w_in)/2] = Input


#### STEP 2: Turn the MTF_PD into a quadratic image¶

In [9]:
## load the MTF:
mtf_pd = pyfits.getdata(path+'mtf_300_0013865-0014441.fits')

In [10]:
w_mtf = mtf_pd.shape[1]
h_mtf = mtf_pd.shape[0]
print w_mtf, h_mtf

714 1972


#### Here I will use Fourier transform to embed the 2D mtf into a quadratic image:¶

In [11]:
from numpy.fft import fft, fftfreq

In [12]:
a = np.fft.fft2(mtf_pd)
a = np.fft.fftshift(a)
b = np.zeros((h_in,h_in),dtype=complex)
b[:,(h_in-w_mtf)/2:(h_in+w_mtf)/2] = a
b = np.fft.ifftshift(b)
b = np.fft.ifft2(b)
mtf_pd = np.abs(b)
mtf_pd = mtf_pd/mtf_pd.max() ## Here, I normalized the MTF so that the maximum signal is equal to 1

In [13]:
print mtf_pd.shape

(1972, 1972)


#### STEP 3: Turn the 1D Gaussian functions into 2D functions¶

In [14]:
## We will assume spherical symmetry. Since we will be working in Fourier space, we have to change the sampling from X to 1/N.X
imscale = 0.0198 # this is the sampling in image space (arcseconds per pixel)

freqscale = 1./(h_in*imscale) # this is the sampling in Fourier space (arcseconds^-1 per pixel)

In [15]:
G2= np.zeros((h_in,h_in))
G3= np.zeros((h_in,h_in))
i,j = np.indices([h_in,h_in])*freqscale
center = np.array([(i.max()-i.min())/2.0, (j.max()-j.min())/2.0])
xc = int(center[1])
yc = int(center[0])
r = np.hypot(i - yc, j - xc)
G2 = np.exp(-2 * np.pi**2 * s_2**2* r**2 )
G3 = np.exp(-2 * np.pi**2 * s_3**2* r**2 )


#### STEP 4: Building the total MTF¶

In [16]:
otf = np.zeros((h_in,h_in))
otf = w_1*mtf_pd+ w_2*G2+ w_3*G3
otf[yc,xc] = otf[yc,xc]+w_4


#### STEP 4: Deconvolution, which will be a simple division in Fourier space¶

In [17]:
Input_fft = np.fft.fft2(Input_emb)
Input_deconv = np.real(np.fft.ifft2(Input_fft/otf)) ## the division
Input_deconv = Input_deconv[:,(h_in-w_in)/2:(h_in+w_in)/2]


### And voila!¶

In [ ]:
## End