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.

import numpy as np
import math as m
import scipy
import matplotlib.pyplot as plt
import pyfits

## 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

The images to be cleaned are 2D and non-quadratic, so now I will carry out several steps before convolving them with the total PSF.

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

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

## 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

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

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

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:

from numpy.fft import fft, fftfreq

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

print mtf_pd.shape

(1972, 1972)

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

## 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)

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

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

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!

## End