cleaning_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cleaning images\n",
    "\n",
    "## Here, I will show how the limb model introduced in __[curve_fitting](https://github.com/fakahil/Projects/blob/master/curve_fitting.ipynb)__ can be used to derive the (point spread function) of the telescope, and correcting the images for instrumental smearing.\n",
    "\n",
    "### The 1D total point spread function of the telescope is given by the following expression:\n",
    "$$ 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)  $$\n",
    "\n",
    "### In Fourier space, the PSF turns into MTF with the following expression:\n",
    "$$MTF(k) =  w_1 \\times MTF_{PD} + w_2\\times G_2(k) + w_3 \\times G_3(x) $$\n",
    "$$MTF(k=0) = (1-w_2-w_2-w_3)$$\n",
    "### 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](https://github.com/fakahil/Projects/blob/master/curve_fitting.ipynb)__. So from now on, we will assume that they are already computed and can be passed to any value we want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math as m\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt\n",
    "import pyfits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Enter here the values of the PSF parameters:\n",
    "\n",
    "w_1 = 0.528\n",
    "w_2 =  0.289\n",
    "w_3 = 0.101\n",
    "w_4 = 1- w_1 -w_2 -w_3\n",
    "s_2 =  0.721\n",
    "s_3 = 2.818"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 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. \n",
    "### __STEP 1__:   to load the image to be corrected and turn it to a quadratic image.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "## loading the input image:\n",
    "path = '/home/fatima/Desktop/project_3/'\n",
    "Input = pyfits.getdata(path+'sufi300_lev3_2009.fits', ignore_missing_end=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "714 1972\n"
     ]
    }
   ],
   "source": [
    "## Dimensions of input image\n",
    "w_in = Input.shape[1]\n",
    "h_in = Input.shape[0]\n",
    "print w_in, h_in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now I will embed the initial image to have dimensions of 1972x1972, the pixels outside the input image will have a value of zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "Input_emb = np.zeros((h_in,h_in))\n",
    "Input_emb[:,(h_in-w_in)/2:(h_in+w_in)/2] = Input  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ###  __STEP 2__:   Turn the MTF_PD into a quadratic image  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "## load the MTF:\n",
    "mtf_pd = pyfits.getdata(path+'mtf_300_0013865-0014441.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "714 1972\n"
     ]
    }
   ],
   "source": [
    "w_mtf = mtf_pd.shape[1]\n",
    "h_mtf = mtf_pd.shape[0]\n",
    "print w_mtf, h_mtf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Here I will use Fourier transform to embed the 2D mtf into a quadratic image:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.fft import fft, fftfreq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = np.fft.fft2(mtf_pd)\n",
    "a = np.fft.fftshift(a)\n",
    "b = np.zeros((h_in,h_in),dtype=complex)\n",
    "b[:,(h_in-w_mtf)/2:(h_in+w_mtf)/2] = a \n",
    "b = np.fft.ifftshift(b)\n",
    "b = np.fft.ifft2(b)\n",
    "mtf_pd = np.abs(b)\n",
    "mtf_pd = mtf_pd/mtf_pd.max() ## Here, I normalized the MTF so that the maximum signal is equal to 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1972, 1972)\n"
     ]
    }
   ],
   "source": [
    "print mtf_pd.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ###  __STEP 3__:   Turn the 1D Gaussian functions into 2D functions  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "## 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 \n",
    "imscale = 0.0198 # this is the sampling in image space (arcseconds per pixel)\n",
    "\n",
    "freqscale = 1./(h_in*imscale) # this is the sampling in Fourier space (arcseconds^-1 per pixel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "G2= np.zeros((h_in,h_in))\n",
    "G3= np.zeros((h_in,h_in))\n",
    "i,j = np.indices([h_in,h_in])*freqscale\n",
    "center = np.array([(i.max()-i.min())/2.0, (j.max()-j.min())/2.0])\n",
    "xc = int(center[1])\n",
    "yc = int(center[0])\n",
    "r = np.hypot(i - yc, j - xc)\n",
    "G2 = np.exp(-2 * np.pi**2 * s_2**2* r**2 )\n",
    "G3 = np.exp(-2 * np.pi**2 * s_3**2* r**2 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ### __STEP 4__:   Building the total MTF  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "otf = np.zeros((h_in,h_in))\n",
    "otf = w_1*mtf_pd+ w_2*G2+ w_3*G3\n",
    "otf[yc,xc] = otf[yc,xc]+w_4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ###  __STEP 4__:   Deconvolution, which will be a simple division in Fourier space  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "Input_fft = np.fft.fft2(Input_emb)\n",
    "Input_deconv = np.real(np.fft.ifft2(Input_fft/otf)) ## the division\n",
    "Input_deconv = Input_deconv[:,(h_in-w_in)/2:(h_in+w_in)/2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## And voila!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## End"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}