{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial\n",
    "\n",
    "This tutorial shows the basic steps of using SEP to detect objects in an image and perform some basic aperture photometry.\n",
    "\n",
    "Here, we use the `fitsio` package, just to read the test image, but you can also use `astropy.io.fits` for this purpose (or any other FITS reader)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sep_pjw as sep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# additional setup for reading the test image and displaying plots\n",
    "import fitsio\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rcParams\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "rcParams['figure.figsize'] = [10., 8.]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we'll read an example image from a FITS file and display it, just to show what we're dealing with. The example image is just 256 x 256 pixels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# read image into standard 2-d numpy array\n",
    "data = fitsio.read(\"../data/image.fits\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# show the image\n",
    "m, s = np.mean(data), np.std(data)\n",
    "plt.imshow(data, interpolation='nearest', cmap='gray', vmin=m-s, vmax=m+s, origin='lower')\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background subtraction\n",
    "\n",
    "Most optical/IR data must be background subtracted before sources can be detected. In SEP, background estimation and source detection are two separate steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# measure a spatially varying background on the image\n",
    "bkg = sep.Background(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are various options for controlling the box size used in estimating the background. It is also possible to mask pixels. For example:\n",
    "```python\n",
    "bkg = sep.Background(data, mask=mask, bw=64, bh=64, fw=3, fh=3)\n",
    "```\n",
    "See the reference section for descriptions of these parameters.\n",
    "\n",
    "This returns an `Background` object that holds information on the spatially varying background and spatially varying background noise level.  We can now do various things with this `Background` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# get a \"global\" mean and noise of the image background:\n",
    "print(bkg.globalback)\n",
    "print(bkg.globalrms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# evaluate background as 2-d array, same size as original image\n",
    "bkg_image = bkg.back()\n",
    "# bkg_image = np.array(bkg) # equivalent to above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# show the background\n",
    "plt.imshow(bkg_image, interpolation='nearest', cmap='gray', origin='lower')\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# evaluate the background noise as 2-d array, same size as original image\n",
    "bkg_rms = bkg.rms()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# show the background noise\n",
    "plt.imshow(bkg_rms, interpolation='nearest', cmap='gray', origin='lower')\n",
    "plt.colorbar();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# subtract the background\n",
    "data_sub = data - bkg\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can also subtract the background from the data array in-place by doing `bkg.subfrom(data)`.\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "\n",
    "**Warning:**\n",
    "\n",
    "If the data array is not background-subtracted or the threshold is too low, you will tend to get one giant object when you run object detection using `sep.extract`. Or, more likely, an exception will be raised due to exceeding the internal memory constraints of the `sep.extract` function.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Object detection\n",
    "\n",
    "Now that we've subtracted the background, we can run object detection on the background-subtracted data. You can see the background noise level is pretty flat. So here we're setting the detection threshold to be a constant value of $1.5 \\sigma$ where $\\sigma$ is the global background RMS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "objects = sep.extract(data_sub, 1.5, err=bkg.globalrms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`sep.extract` has many options for controlling detection threshold, pixel masking, filtering, and object deblending. See the reference documentation for details.\n",
    "\n",
    "`objects` is a NumPy structured array with many fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# how many objects were detected\n",
    "len(objects)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`objects['x']` and `objects['y']` will give the centroid coordinates of the objects. Just to check where the detected objects are, we'll over-plot the object coordinates with some basic shape parameters on the image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from matplotlib.patches import Ellipse\n",
    "\n",
    "# plot background-subtracted image\n",
    "fig, ax = plt.subplots()\n",
    "m, s = np.mean(data_sub), np.std(data_sub)\n",
    "im = ax.imshow(data_sub, interpolation='nearest', cmap='gray',\n",
    "               vmin=m-s, vmax=m+s, origin='lower')\n",
    "\n",
    "# plot an ellipse for each object\n",
    "for i in range(len(objects)):\n",
    "    e = Ellipse(xy=(objects['x'][i], objects['y'][i]),\n",
    "                width=6*objects['a'][i],\n",
    "                height=6*objects['b'][i],\n",
    "                angle=objects['theta'][i] * 180. / np.pi)\n",
    "    e.set_facecolor('none')\n",
    "    e.set_edgecolor('red')\n",
    "    ax.add_artist(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`objects` has many other fields, giving information such as second moments, and peak pixel positions and values. See the reference documentation for `sep.extract` for descriptions of these fields. You can see the available fields:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# available fields\n",
    "objects.dtype.names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aperture photometry\n",
    "\n",
    "Finally, we'll perform simple circular aperture photometry with a 3 pixel radius at the locations of the objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "flux, fluxerr, flag = sep.sum_circle(data_sub, objects['x'], objects['y'],\n",
    "                                     3.0, err=bkg.globalrms, gain=1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`flux`, `fluxerr` and `flag` are all 1-d arrays with one entry per object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# show the first 10 objects results:\n",
    "for i in range(10):\n",
    "    print(\"object {:d}: flux = {:f} +/- {:f}\".format(i, flux[i], fluxerr[i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finally a brief word on byte order\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**Note:**\n",
    "\n",
    "If you are using SEP to analyze data read from FITS files with\n",
    "[astropy.io.fits](http://astropy.readthedocs.org/en/stable/io/fits/)\n",
    "you may see an error message such as:\n",
    "\n",
    "```\n",
    "ValueError: Input array with dtype '>f4' has non-native byte order.\n",
    "Only native byte order arrays are supported. To change the byte\n",
    "order of the array 'data', do 'data = data.astype(data.dtype.newbyteorder(\"=\"))'\n",
    "```\n",
    "\n",
    "It is usually easiest to do this byte-swap operation directly after\n",
    "reading the array from the FITS file. The exact procedure is slightly different,\n",
    "depending on the version of ``numpy``. For ``numpy<2.0``, the operation was:\n",
    "\n",
    "```python\n",
    "# Byte order changed in-place\n",
    ">>> data = data.byteswap(inplace=True).newbyteorder() \n",
    "# Data copied to a new array\n",
    ">>> new_data = data.byteswap().newbyteorder()\n",
    "```\n",
    "\n",
    "For ``numpy>=2.0``, the correct operation is one of the following:\n",
    "\n",
    "```python\n",
    "# Copies data to a new array, preserves the ordering of the original array\n",
    ">>> new_data = data.astype(data.dtype.newbyteorder(\"=\")) \n",
    "# The same outcome as the previous operation\n",
    ">>> new_data = data.byteswap()\n",
    ">>> new_data = new_data.view(new_data.dtype.newbyteorder(\"=\"))\n",
    "# Changes data in-place\n",
    ">>> data = data.byteswap()\n",
    ">>> data = data.view(data.dtype.newbyteorder(\"=\"))\n",
    "```\n",
    "\n",
    "If you do this in-place operation, ensure that there are no other\n",
    "references to ``data``, as they will be rendered nonsensical.\n",
    "\n",
    "For the interested reader, this byteswap operation is necessary because\n",
    "``astropy.io.fits`` always returns big-endian byte order arrays, even on\n",
    "little-endian machines. This is due to the FITS standard requiring\n",
    "big-endian arrays (see the \n",
    "[FITS Standard](https://fits.gsfc.nasa.gov/fits_standard.html)), and \n",
    "``astropy.io.fits`` aiming for backwards compatibility.\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}