Image Processing Workbook I.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image Processing\n",
    "\n",
    "This ipython (sorry, Jupyter) notebook contains the examples that I'll be covering in the 1-dimensional part of my image processing session.  There are no external data files to worry about.\n",
    "\n",
    "I should warn you that this is tested on python 2; I think it'll run on python 3 but I haven't actually checked."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 0\n",
    "\n",
    "Import numpy and matplotlib, and make a simple plot of a sine wave\n",
    "\n",
    "_N.b._  I usually start with a piece of boilerplate roughly like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "# A nice alternative to inline (which allows interaction with the plots, but can be confusing) is\n",
    "# %matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Answer 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A Toy 1-D Model\n",
    "\n",
    "We'll make a simple one dimensional model of a star field.  Images of real stars are complicated, but we'll assume that the profile is a Gaussian.  I write a Gaussian with mean $\\mu$ and variance $\\sigma^2$ as $N(\\mu, \\sigma^2)$.\n",
    "\n",
    "Almost all stars are essentially point sources as viewed from the Earth, so stars look like the Point Spread Function (PSF) produced by a combination of the atmosphere, the telescope, and the detector.  There's no standard notation for the PSF, but I always call it $\\phi$.  The Full Width Half Maximum (FWHM) is a measure of the PSF's width.\n",
    "\n",
    "In addition to flux from the stars that we care about there's a smooth background of extra light that comes from a number of sources (atmospheric emission, scattering of star and moonlight, dark current in the detector...); for now we'll just treat this as an annoying constant that I'll call `The Sky'.\n",
    "\n",
    "In reality CCD data is only measured where we have pixels, but under certain conditions (band-limited; satisfying the Nyquist condition) it turns out that the pixellation is unimportant, so we'll ignore it for this session."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 1\n",
    "\n",
    "Write a simulator that simulates noise-free 1-D data.  Provide a function with signature\n",
    "\n",
    "    def phi(x, xc=0.0, fwhm=2):\n",
    "        \"\"\"Return a numpy array with a star centred at the point xc and FWHM evaluated at the points x\"\"\"\n",
    "        \n",
    "Use a Gaussian PSF, and set the sky level to S; for a Gaussian $N(0, \\beta^2)$, the FWHM is $2\\sqrt{2\\ln(2)}\\,\\beta \\sim 2.3548\\beta$\n",
    "\n",
    "Plot a few of your beautiful simulations.  Once you've got it working, wrap it up in a function\n",
    "\n",
    "    simulate(x, S=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### Answer 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 2\n",
    "\n",
    "There are lots of sources of noise in astronomical data, but for now let's assume that the only thing that matters is the finite number of photons that we detect.  Detecting $n$ photons in a pixel results in a Poisson distribution; if\n",
    "$n$'s mean value is $\\mu$, then its variance is also $\\mu$.  If $\\mu \\gg 1$, $P(\\mu) \\sim N(\\mu, \\mu)$.   You can do better than this if you need to;  if $\\mu > 4$ (!), Anscombe showed that  $x \\equiv 2\\sqrt{x + 3/8}$ is approximately Gaussian, $N(2\\sqrt{\\mu + 3/8} - 1/(4\\sqrt\\mu), 1)$ -- but you probably only care if you're looking at the tails of the distribution.\n",
    "\n",
    "Add noise to your simulation.  You can get Poisson variates from numpy by saying _e.g._:\n",
    "\n",
    "    mu = 100\n",
    "    print(np.random.poisson(lam=mu))\n",
    "    print(np.random.normal(loc=mu, scale=np.sqrt(mu)))    # Here's how to get a Gaussian approximation\n",
    "    \n",
    "If you want to set the random number seed so that the noise added is always the same, say something like\n",
    "\n",
    "    np.random.seed(666)\n",
    "    \n",
    "If you can't see your stars you might want to make them brighter --- remember that we added a background noise of 10 when we set the sky level to 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mu = 100\n",
    "print(np.random.poisson(lam=mu))\n",
    "print(np.random.normal(loc=mu, scale=np.sqrt(mu)))    # Here's how to get a Gaussian approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Answer 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 3\n",
    "\n",
    "Let's investigate measuring the flux in a single isolated star.  Start my modifying your previous answer to simulate a noisy single star with $\\beta = 2$ (a FWHM of c. 5 pixels), centred at x=0, with total flux F0.  In reality estimating the sky background is not trivial, but for now let's simply subtract its known value.  Plot one of your simulations, with F0=1000 and S=100.\n",
    "\n",
    "The simplest way to measure the flux in a star is to sum the counts within an aperture.  Modify your code to estimate the flux within 5 pixels of the centre, then run a Monte-Carlo simulation to estimate the mean and variance of your estimator.\n",
    "\n",
    "Your mean should be close to F0 --- if it isn't, take another look at your code.  Your value won't be exactly F0; is this a statistical anomaly, or is it something more interesting?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Answer 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 4\n",
    "\n",
    "Package your estimator into a function and estimate the mean and variance as a function of $R$; plot your results.  What value of $R$ gives the smallest variance in $F_0$?\n",
    "\n",
    "For small $R$ are we measuring $F_0$ correctly?  If not, make appropriate corrections and remake your plots.  Does your conclusion change?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Answer 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 5\n",
    "\n",
    "The lecture notes imply that we'd do better to use the PSF to measure the flux.  Write a Monte-Carlo simulation, and compare your results with the previous task"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Answer 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "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": 0
}