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