05_Intro_to_Source_Detection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img align=\"left\" src = https://project.lsst.org/sites/default/files/Rubin-O-Logo_0.png width=250 style=\"padding: 10px\"> \n",
    "<b>Introduction to Source Detection</b> <br>\n",
    "Last verified to run on 2021-10-12 with LSST Science Pipelines release w_2021_40 <br>\n",
    "Contact author: Alex Drlica-Wagner, Melissa Graham <br>\n",
    "Target audience: All DP0 delegates. <br>\n",
    "Container Size: medium <br>\n",
    "Questions welcome at <a href=\"https://community.lsst.org/c/support/dp0\">community.lsst.org/c/support/dp0</a> <br>\n",
    "Find DP0 documentation and resources at <a href=\"https://dp0-1.lsst.io\">dp0-1.lsst.io</a> <br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Credit:** This tutorial was originally developed by Alex Drlica-Wagner and Imran Hasan in the context of the LSST Stack Club. Please consider acknowledging Alex Drlica-Wagner and Imran Hasan in any publications or software releases that make use of this notebook's contents."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning Objectives\n",
    "This notebook provides a brief introduction to running the LSST Science Pipelines source detection, measurement, and deblending algorithms. It does not go into depth about optimizing detection, measurement, or deblending parameters for different types of sources. \n",
    "\n",
    "Some source detection and measurement details come from Robert Lupton's [Tune Detection.ipynb](https://github.com/RobertLuptonTheGood/notebooks/blob/master/Demos/Tune%20Detection.ipynb) and [Kron.ipynb](https://github.com/RobertLuptonTheGood/notebooks/blob/master/Demos/Kron.ipynb).\n",
    "Interaction with `lsst.afw.display` was also improved by studying Michael Wood-Vasey's [DC2_Postage Stamps.ipynb](https://github.com/LSSTDESC/DC2-analysis/blob/master/tutorials/dm_butler_postage_stamps.ipynb).\n",
    "More information on footprints can be found on the Stack Club notebook by Imran Hasan [here](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/SourceDetection/Footprints.ipynb).\n",
    "\n",
    "After working through this notebook you should be able to\n",
    "   1. Run the `lsst.meas.algorithm` source detection, deblending, and measurement tasks.\n",
    "   2. Plot the resulting source catalogs on an image.\n",
    "   3. Examine the `footprint` of the detected sources.\n",
    "\n",
    "Other techniques that are demonstrated, but not emphasized, in this notebook are\n",
    "   1. Use the `butler` to access a specific `calexp`.\n",
    "   2. Create an image cutout and use `lsst.afw.display` to plot it.\n",
    "\n",
    "\n",
    "### Set Up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import python packages\n",
    "import os\n",
    "import warnings\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import IFrame, display, Markdown\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from matplotlib.patches import Rectangle\n",
    "from astropy.visualization import ZScaleInterval\n",
    "\n",
    "%matplotlib inline\n",
    "# %matplotlib ipympl # currently slow, but may be a good option in the future"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import LSST Science Pipelines packages (see pipelines.lsst.io)\n",
    "import lsst.daf.base as dafBase\n",
    "from lsst.daf.butler import Butler\n",
    "import lsst.afw.image as afwImage\n",
    "import lsst.afw.display as afwDisplay\n",
    "import lsst.afw.table as afwTable\n",
    "import lsst.geom as geom\n",
    "\n",
    "# Import Pipeline tasks \n",
    "from lsst.pipe.tasks.characterizeImage import CharacterizeImageTask\n",
    "from lsst.pipe.tasks.calibrate import CalibrateTask\n",
    "from lsst.meas.algorithms.detection import SourceDetectionTask\n",
    "from lsst.meas.deblender import SourceDeblendTask\n",
    "from lsst.meas.base import SingleFrameMeasurementTask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prevent some helpful but ancillary warning messages from printing\n",
    "#   during some LSST DM Release calls\n",
    "warnings.simplefilter(\"ignore\", category=UserWarning)\n",
    "warnings.simplefilter(\"ignore\", category=FutureWarning)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use lsst.afw.display with the matplotlib backend\n",
    "afwDisplay.setDefaultBackend('matplotlib')\n",
    "plt.rcParams['figure.figsize'] = (8.0, 8.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What version of the Stack are we using?\n",
    "! echo $IMAGE_DESCRIPTION\n",
    "! eups list -s | grep lsst_distrib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.0 Data access\n",
    "\n",
    "Here we use the `butler` to access a `calexp` from the DP0.1 dataset. More information on the `butler` and `calexp`, and how to determine the `dataId` (e.g., the visit and detector numbers) are available in other tutorials. Here we start assuming that information is known, and assuming a basic understanding of `butler` use.\n",
    "\n",
    "If a warning in pink is output below, it is related to the recent gen2 to gen3 butler upgrade, and can be ignored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the dataId using just visit and detector\n",
    "dataId = {'visit': 512055, 'detector': 75}\n",
    "\n",
    "# For DC2 gen3, these are the only optoins\n",
    "repo = 's3://butler-us-central1-dp01'\n",
    "collection = '2.2i/runs/DP0.1'\n",
    "\n",
    "# Use the butler to get the calexp\n",
    "butler = Butler(repo, collections=collection)\n",
    "calexp = butler.get('calexp', **dataId)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "\n",
    "As described in other tutorials, the `calexp` object possess more than just the raw pixel data of the image. It also contains a `mask`, which stores information about various pixels in a bit mask.\n",
    "\n",
    "Here are some optional commands to explore the calexp. Uncomment one of the code lines to learn more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If you want to investigate the contents of the masked image:\n",
    "# calexp.maskedImage\n",
    "\n",
    "# If you just want one of the three components:\n",
    "# calexp.maskedImage.image\n",
    "# calexp.maskedImage.mask\n",
    "# calexp.maskedImage.variance\n",
    "\n",
    "# These also work:\n",
    "# calexp.image\n",
    "# calexp.mask\n",
    "# calexp.variance\n",
    "\n",
    "# The calexp also contains the PSF, the WCS, and the photometric calibration\n",
    "# calexp.getPsf()\n",
    "# calexp.getWcs()\n",
    "# calexp.getPhotoCalib()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we are interested in performing our own source detection and measurement, we choose to clear the existing `DETECTED` mask plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Unset the `DETECTED` bits of the mask plane\n",
    "calexp.mask.removeAndClearMaskPlane('DETECTED')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the calexp we just retrieved\n",
    "plt.figure()\n",
    "afw_display = afwDisplay.Display()\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(calexp.image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Add the subtracted sky background back into the image\n",
    "\n",
    "Retrieve the subtracted background for the same dataId. This section is optional."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkgd = butler.get('calexpBackground', **dataId)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Display the background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "afw_display = afwDisplay.Display()\n",
    "afw_display.scale('linear', 'zscale')\n",
    "afw_display.mtv(bkgd.getImage())\n",
    "plt.title(\"Local Polynomial Background\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add the background into the calexp, and re-display the calexp. Note the scale in the sidebar now goes up to thousands of counts instead of hundreds of counts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note: executing this cell multiple times will add the background multiple times\n",
    "calexp.maskedImage += bkgd.getImage()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "afw_display = afwDisplay.Display()\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(calexp.image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Generate and display a smaller cutout"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our next step is to generate a cutout image. This is done by creating a bounding box and passing it to the `Factory` method of our calexp (a `lsst.afw.image.Exposure` object). Below we explain the specific arguments that we are passing to `Factory`:\n",
    "```\n",
    "calexp : the ExposureF we are starting from\n",
    "bbox   : the bounding box of the cutout\n",
    "origin : the image pixel origin is local to the cutout array\n",
    "deep   : copy the data rather than passing by reference\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the pixel coordinates of a known low surface brightness \"galaxy\"\n",
    "x_target, y_target = 1700, 2100\n",
    "width, height = 400, 400\n",
    "xmin, ymin = x_target-width//2, y_target-height//2\n",
    "point = geom.Point2D(x_target,y_target)\n",
    "\n",
    "# Define a small region for a cutout\n",
    "bbox = geom.Box2I()\n",
    "bbox.include(geom.Point2I(xmin, ymin))\n",
    "bbox.include(geom.Point2I(xmin + width, ymin + height))\n",
    "\n",
    "# An alternative way to defined the same cutout region\n",
    "# bbox = geom.Box2I(geom.Point2I(xmin, ymin), geom.Extent2I(width, height))\n",
    "\n",
    "# Generate the cutout image\n",
    "cutout = calexp.Factory(calexp, bbox, origin=afwImage.LOCAL, deep=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Follow the same procedure as before to plot the cutout\n",
    "plt.figure()\n",
    "afw_display = afwDisplay.Display()\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(cutout.image)\n",
    "plt.gca().axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.0 Source Detection, Deblending, and Measurement\n",
    "\n",
    "We now want to run the LSST Science Pipelines' source detection, deblending, and measurement tasks. While we run all three tasks, this notebook is mostly focused on the detection of sources.\n",
    "\n",
    "Recall that these tasks were imported up at the top of this notebook, from `lsst.pipe` and `lsst.meas`. More information can be found at [pipelines.lsst.io](https://pipelines.lsst.io/) (the search bar at the top left of that page is a very handy way to find documentation for a specific task).\n",
    "\n",
    "We start by creating a minimal schema for the source table. The schema describes the output properties that will be measured for each source. This schema will be passed to all of the tasks, as we call each in turn, and each task will add columns to this schema as it measures sources in the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a basic schema to use with these tasks\n",
    "schema = afwTable.SourceTable.makeMinimalSchema()\n",
    "print(schema)\n",
    "\n",
    "# Create a container which will be used to record metadata\n",
    "#  about algorithm execution\n",
    "algMetadata = dafBase.PropertyList()\n",
    "print('algMetadata: ')\n",
    "algMetadata"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Configuring Tasks\n",
    "\n",
    "Each task possesses an associated configuration class. The properties of these configuration classes can be determined from the classes themselves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment the following line to view help\n",
    "#  for the CharacterizeImageTask configuration\n",
    "# Replace 'CharacterizeImageTask' with a different\n",
    "#  task name to view additional help information\n",
    "\n",
    "# help(CharacterizeImageTask.ConfigClass())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a starting point, like the `schema` and `algMetadata` above, here we set some basic config parameters and instantiate the tasks to get you started. In this case, we configure several different tasks:\n",
    "\n",
    "* CharacterizeImageTask: Characterizes the image properties (e.g., PSF, etc.)\n",
    "* SourceDetectionTask: Detects sources\n",
    "* SourceDeblendTask: Deblend sources into constituent children\n",
    "* SingleFrameMeasurementTask: Measures source properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Characterize the image properties\n",
    "config = CharacterizeImageTask.ConfigClass()\n",
    "config.psfIterations = 1\n",
    "charImageTask = CharacterizeImageTask(config=config)\n",
    "\n",
    "# Detect sources\n",
    "config = SourceDetectionTask.ConfigClass()\n",
    "# detection threshold in units of thresholdType\n",
    "config.thresholdValue = 10\n",
    "# units for thresholdValue\n",
    "config.thresholdType = \"stdev\"\n",
    "sourceDetectionTask = SourceDetectionTask(schema=schema, config=config)\n",
    "\n",
    "# Deblend sources\n",
    "sourceDeblendTask = SourceDeblendTask(schema=schema)\n",
    "\n",
    "# Measure source properties\n",
    "config = SingleFrameMeasurementTask.ConfigClass()\n",
    "sourceMeasurementTask = SingleFrameMeasurementTask(schema=schema,\n",
    "                                                   config=config,\n",
    "                                                   algMetadata=algMetadata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that if you want to change the value of a config parameter (e.g., `psfIterations`), do not change it in the already constructed task. Instead, change your config object and then construct a new characterize image task. Like so:\n",
    "> `config.psfIterations = 3` <br>\n",
    "> `charImageTask = CharacterizeImageTask(config=config)`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like the configs, we can use `help` to explore each task and the methods used to run it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# help(charImageTask)\n",
    "\n",
    "# Uncomment the following line, position your cursor after the period,\n",
    "#  and press tab to see a list of all methods. Then recomment the line\n",
    "#  because \"Task.\" is not executable and will cause an error.\n",
    "# charImageTask.\n",
    "\n",
    "# Use the help function on any of the methods to learn more:\n",
    "# help(charImageTask.writeSchemas)\n",
    "\n",
    "# E.g., find out what options there are for config.thresholdType\n",
    "# help(SourceDetectionTask.ConfigClass)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the each of the tasks configured, we can now move on to running the source detection, deblending, and measurement. First we create `SourceTable` for holding the output of our source analysis. The columns and characteristics of this table are defined by the `schema` that we created in our configuration step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tab = afwTable.SourceTable.make(schema)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Image Characterization\n",
    "\n",
    "Next we characterize our image. This calculates various global properties of the image, such as the full-width half-max of its point spread function (PSF FWHM)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Image characterization (this cell may take a few seconds)\n",
    "result = charImageTask.run(calexp)\n",
    "\n",
    "# Get the PSF at our point of interest\n",
    "psf = calexp.getPsf()\n",
    "sigma = psf.computeShape(point).getDeterminantRadius()\n",
    "pixelScale = calexp.getWcs().getPixelScale().asArcseconds()\n",
    "\n",
    "# The factor of 2.355 converts from std to fwhm\n",
    "#print('psf fwhm = {:.2f} arcsec'.format(sigma*pixelScale*2.355))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the image characterized, we are now interested in running the source detection, deblending, and measurement tasks. Each of these tasks is called with the `run` method. The parameters of this method can be investigated using `help`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We are specifically interested in the `SourceMeasurementTask`\n",
    "# help(sourceMeasurementTask.run)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Source detection (this cell may take a few seconds)\n",
    "result = sourceDetectionTask.run(tab, calexp)\n",
    "type(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The source detection task returns an [`lsst.pipe.base.struct.Struct`](http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/classlsst_1_1pipe_1_1base_1_1struct_1_1_struct.html). A `Struct` is just a generalized container for storing multiple output components and accessessing them as attributes. The content of this `Struct` can be investigated with the `getDict` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k, v in result.getDict().items():\n",
    "    print(k, type(v))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The members of the `Struct` can be accessed either through dictionary keys or as attributes of the `Struct`. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result.numPosPeaks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we use `config.thresholdValue = 10` we get 943 positive peaks. If we rerun with `config.thresholdValue = 3` we would find 5367 positive peaks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sources = result.sources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the columns for sources will be NaN at first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment to see that most values are NaN\n",
    "# sources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that if we desire we can save some of these processed objects to disk as FITS tables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sources.writeFits(\"outputTable.fits\")\n",
    "# calexp.writeFits(\"example1-out.fits\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Deblend and Measure Sources\n",
    "\n",
    "Next we run the `SourceDeblendTask` and `SingleFrameMeasurementTask`. A deeper investigation of these tasks is beyond the scope of this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Source deblending\n",
    "sourceDeblendTask.run(calexp, sources)\n",
    "\n",
    "# Source measurement (catch future warning about machine precision)\n",
    "sourceMeasurementTask.run(measCat=sources, exposure=calexp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get a better look at the output sources, we need to make sure that the `SourceCatalog` is contiguous in memory. Converting to an `astropy` table provides a human-readable output format. A deeper dive into `SourceCatalog` is beyond the scope of this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The copy makes sure that the sources are sequential in memory\n",
    "sources = sources.copy(True)\n",
    "\n",
    "# Investigate the output source catalog\n",
    "sources.asAstropy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Overplot Sources on Image\n",
    "\n",
    "We can now overplot our detected sources on the calexp or cutout image using `afwDisplay`.\n",
    "\n",
    "<a id='display-error'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display the cutout and sources with afw display\n",
    "image = cutout.image\n",
    "\n",
    "plt.figure()\n",
    "afw_display = afwDisplay.Display()\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(image)\n",
    "plt.gca().axis('off')\n",
    "\n",
    "# We use display buffering to avoid re-drawing the image\n",
    "#  after each source is plotted\n",
    "with afw_display.Buffering():\n",
    "    for s in sources:\n",
    "        afw_display.dot('+', s.getX(), s.getY(), ctype=afwDisplay.RED)\n",
    "        afw_display.dot('o', s.getX(), s.getY(), size=20, ctype='orange')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.0 Footprints\n",
    "\n",
    "Object footprints are an integral component of the high-level CCD processing tasks (e.g., detection, measurement, and deblending). To quote [Bosch et al. (2017)](https://arxiv.org/pdf/1705.06766.pdf), \n",
    "\n",
    "> Footprints record the exact above-threshold detection region on a CCD. These are similar to  SExtractor’s “segmentation map\", in that they identify which pixels belong to which detected objects\n",
    "\n",
    "This quote draws an analogy between footprints and segmentation maps, since they both identify pixels with values above some threshold. This is a useful similarity, since it gives us a place to start understanding the properties of footprints."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result of the `SourceDetectionTask` stores the footprints associated to detected objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grab the above-threshold footprints that were detected,\n",
    "#  and assign them to the variable `fps`\n",
    "fpset = result.positive\n",
    "fps = fpset.getFootprints()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can get a rough view of the first source's footprint from its span\n",
    "fps[0].getSpans()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can almost see the footprint by looking at the 1's and 0's here. Keep in mind that the first row of this array will be the *bottom* row of the image. Later, when we display the footprint, its general pattern will appear \"upside down\" compared to this pattern of 1s. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Heavy Footprints\n",
    "\n",
    "At the moment, our footprints indicate which pixels they consist of, but not the values of those pixels from the image. To extract the actual values of the pixels that correspond to the ones in the footprint span, we need to convert our footprint into a `HeavyFootprint`. HeavyFootprints have all of the qualities of Footprints, but additionally 'know' about pixel level data from the image, variance, and mask planes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First we demonstrate that the footprint is NOT heavy\n",
    "fps[0].isHeavy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Next, we make all the footprints heavy at the same time\n",
    "#  by operating on the footprint set\n",
    "fpset.makeHeavy(calexp.getMaskedImage())\n",
    "\n",
    "# This means we have to redefine fps:\n",
    "hfps = fpset.getFootprints()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# All of the arrays here will be flattend 1D arrays\n",
    "#  of pixels from the footprint.\n",
    "# Uncomment this line to print the array and see that\n",
    "#  now, it contains pixel values.\n",
    "\n",
    "# hfps[0].getImageArray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the span set to reassemble the image array into the footprint. Above we saw that the image array is a 1D numpy array, but the footprint itself is 2 dimensional. Fortunately, the span set has an `unflatten` method that rearranges the image array into the proper 2 dimensional shape. If you want to change the colormap, see [matplotlib colormap options](https://matplotlib.org/stable/tutorials/colors/colormaps.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.imshow(fps[0].getSpans().unflatten(hfps[0].getImageArray()),\n",
    "           cmap='bone', origin='lower')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The heavy footprint also comes with a 1d mask array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hfps[0].getMaskArray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To understand these values, lets look at the mask plane's dictionary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calexp.getMask().getMaskPlaneDict()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The values are the exponent of the bitmask. So pixels only marked detected will be 2^5 = 32. Pixels that are both on the edge of the original image and detected will be 2^5 + 2^4 = 48. We will visualize the mask plane values in a similar manner as before, except that we will be displaying the values of the mask array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "im = plt.imshow(fps[0].getSpans().unflatten(hfps[0].getMaskArray()),\n",
    "                origin='lower')\n",
    "\n",
    "# Create a new axis, \"cax\" on the right side of the image display.\n",
    "# The width of cax will be 5% of the axis \"ax\".\n",
    "# The padding between cax and ax will be 1% of the axis.\n",
    "ax = plt.gca()\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=\"1%\")\n",
    "plt.colorbar(im, cax=cax, ticks=[0, 32, 32+16])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More information on footprints can be found on the Stack Club notebook by Imran Hasan [here](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/SourceDetection/Footprints.ipynb)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "LSST",
   "language": "python",
   "name": "lsst"
  },
  "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}