SourceDetection/Footprints.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Introduction to Footprints\n",
    "<br>Developer(s): **Imran Hasan** ([@ih64](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@ih64))\n",
    "<br>Maintainer(s): **Alex Drlica-Wagner** ([@kadrlica](https://github.com/LSSTScienceCollaborations/StackClub/issues/new?body=@kadrlica))\n",
    "<br>Level: **Intermediate**\n",
    "<br>Last Verified to Run: **2021-09-10**\n",
    "<br>Verified Stack Release: **w_2021_33**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Learning Objectives:\n",
    "\n",
    "After working through this tutorial you should be able to: \n",
    "1. Understand footprints, and their relationship to detection, and measurement\n",
    "2. Use methods and attributes from the footprint class to get information like centroids, size, shape, of the footprint\n",
    "3. Perform operations on footprints, like dilation an errosion\n",
    "4. Understand the relaitonship between peaks and footprints, and their role in deblending\n",
    "4. Use footprints to access pixel level data in their parent exposures\n",
    "\n",
    "### Logistics\n",
    "This notebook is intended to be run at `lsst-lsp-stable.ncsa.illinois.edu` or `data.lsst.cloud` from a local git clone of the [StackClub](https://github.com/LSSTScienceCollaborations/StackClub) repo.\n",
    "\n",
    "### Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Site, host, and stack version\n",
    "! echo $EXTERNAL_INSTANCE_URL\n",
    "! echo $HOSTNAME\n",
    "! eups list -s | grep lsst_distrib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For this tutorial we'll need the following modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import os\n",
    "import warnings\n",
    "# Filter some warnings printed by v16.0 of the stack\n",
    "warnings.simplefilter(\"ignore\", category=FutureWarning)\n",
    "warnings.simplefilter(\"ignore\", category=UserWarning)\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import matplotlib.animation as animation\n",
    "from IPython.display import IFrame, display, Markdown, HTML\n",
    "import lsst.daf.persistence as dafPersist\n",
    "import lsst.afw.table as afwTable\n",
    "from lsst.meas.algorithms import SourceDetectionTask\n",
    "from lsst.meas.base import NoiseReplacer, NoiseReplacerConfig\n",
    "import lsst.obs.base as obsBase\n",
    "import lsst.afw.detection as afwDetection\n",
    "import lsst.afw.display as afwDisplay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "afwDisplay.setDefaultBackend('matplotlib') \n",
    "plt.rcParams['figure.figsize'] = (8.0, 8.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Footprints\n",
    "\n",
    "To paraphrase from [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",
    "As you might expect, this means footprints are integral to high-level CCD processing tasks&mdash;like detection, measurement, and deblending&mdash;which directly impact science results. Because footprints are so closely related to these very important processes, we will take a close look at them in this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The Basics: Footprint Geometry and Location\n",
    "In the quote above, an analogy was drawn between footprints and segmentation maps, as they both identify above threshold pixels. As we first introduce footprints, we will concentrate on this similarity as it gives us a place to start understanding the location and geometeric properties of footprints. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "URL = os.getenv('EXTERNAL_INSTANCE_URL')\n",
    "if URL.endswith('data.lsst.cloud'): # IDF\n",
    "    repo = \"s3://butler-us-central1-dp01\"\n",
    "elif URL.endswith('ncsa.illinois.edu'): # NCSA\n",
    "    repo = \"/repo/dc2\"\n",
    "else:\n",
    "    raise Exception(f\"Unrecognized URL: {URL}\")\n",
    "\n",
    "collection = '2.2i/runs/DP0.1'\n",
    "dataset_type = \"deepCoadd\"\n",
    "dataId = {'tract': 4226, 'band': 'i', 'patch': (0) + 7*(4) }  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lsst.daf.butler import Butler \n",
    "butler = Butler(repo,collections=collection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grab the data\n",
    "deepCoadd = butler.get(dataset_type, dataId=dataId)\n",
    "# Let's grab the source catalog too\n",
    "cat = butler.get(dataset_type+'_meas', dataId=dataId)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Take a look at the image of the patch\n",
    "plt.figure(figsize=(30,30))\n",
    "afw_display = afwDisplay.Display()\n",
    "afw_display.setMaskTransparency(100)\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(deepCoadd.getMaskedImage())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# set up minimal detection task\n",
    "config = SourceDetectionTask.ConfigClass()\n",
    "# detect sources at 10 sigma \n",
    "config.thresholdValue = 10       \n",
    "config.thresholdType = \"stdev\"   \n",
    "\n",
    "sourceDetectionTask =   SourceDetectionTask(config=config)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the `detectFootprints` method in `SourceDetectionTask` to find and store the detected footprints in the image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run detection on our frame\n",
    "detect_res = sourceDetectionTask.detectFootprints(deepCoadd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`detect_res` is a `Struct`. `Struct` objects are similar to python dictionaries, in that they have key value pairs. However, here you access the values by using the keys as an attribute. If you are and IDL user, these will remind you of `Structures` in IDL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "detect_res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at `detect_res`, we can see that there are attributes associated with sources detected positive (above threshold) and detected negative (below threshold). Let's concentrate on the positive detections for now, so we can look at astrophysical sources. The positive attribute holds a `FootprintSet`, an object that contains all the positively detected footprints. The `numPos` attribute tells us the number of footprints that were detected positive. The `numPosPeaks` attribute tells us the number of peaks that were detected positive. You'll notice that the number of peaks is larger because more than one peak can belong to a footprint. We will discuss peaks more fully in a later section, but for now think of them as what they sound like&mdash;the peaks in the surface brightness profile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets grab the above threshold footprints that were detected and assign them to a varriable\n",
    "fpset = detect_res.positive"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the variable `fpset` is a `FootprintSet` object. While it contains many footprints itself, a few more steps are required to obtain individual footprints. Namely, we will ask for a `FootprintSet` to give us a list of footprints, which we can iterate over"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fps = fpset.getFootprints()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have drilled our way down to getting our hands on footprints. Lets discuss some of their methods and attributes. To begin, we will concentrate on footprint geometry and location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the region\n",
    "# This will return a box with the pixel coordinates of the CCD that the footprint lies on\n",
    "fps[0].getRegion()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The center of the footprint, in the x,y coordinates on the CCD chip\n",
    "fps[0].getCentroid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bounding box of the footprint. This is the smallest rectangle you can enscribe around the footprint. Not all of the pixels in the bounding box are pixels in the footprint, since the footprint is not necessarily rectangular in shape. the BBox will always be inside the region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fps[0].getBBox()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the area of the footprint. **NOTE:** this is not the area of the bounding box. It is the number of pixels associated with the footprint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fps[0].getArea()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The shape of a footprint can be calculated. It is assumed to be an ellipse, and the relevent parameters are defined with the quadurpole, $Q_{i,j}$, a matrix containing the second moments, where\n",
    "\n",
    "$$Q_{i,j} = \\int x_i x_j \\mathrm{d}x^i \\mathrm{d}x^j$$\n",
    "\n",
    "Notice this definition *does not* include the surface brightness of the detected object. That is to say, the quadropole computed for footprints describes the geometry of the pixels in the footprint, *not* the intensity weighted shape of its detected object.\n",
    "\n",
    "The elements for a footprint can be computed with the `getShape()` method. Note that the matix is symmetric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fps[0].getShape()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at $Q_{i,j}$, the diagonal elements are close in magnitude, and the off diagonal elements are close to zero. We expect the footprint to be roughly circular. Now that we have all the pertinent information about where the footprint is and what pixels belong to it, we can display it and see if our expectation is true.\n",
    "\n",
    "## Visualizing Footprints & Heavy Footprints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get information about which pixels in the footprints bounding box actually belong to the footprint by using its `SpanSet`. The DM stack is object oriented, and often classes inherit or contain other classes. That is ture here too. The spanSet is a collection of `n` `spans`, where each span has `m` elements, where `n` and `m` are the height and width of the bounding box, respectively. \n",
    "\n",
    "If an element in a span is zero, this means the pixel is inside the footprint's bounding box but does not belong to the footprint. If an element is one, the pixel is both in the bounding box and associated with the footprint.\n",
    "\n",
    "Lets dump the output of the spanset from the first footprint in our list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fps[0].getSpans()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can almost see the footprint by looking at the 1's and zeros here. To extract the actual pixel values that correspond to the ones in the span, we need an additional step. At the moment, our footprints can tell you if a pixel belongs to it or not, buthave no way of accessing pixel values on the `deepCoadd`. To remedy this, we will turn our footprint into a `HeavyFootprint`. HeavyFootprints have all of the qualities of Footprints, but additionally 'know' about pixel level data from the deepCoadd image, variance, and mask planes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first we domonstrate the footprint is not heavy\n",
    "fps[0].isHeavy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we will make all the footprins heavy at the same time by operating on the footprint set\n",
    "fpset.makeHeavy(deepCoadd.getMaskedImage())\n",
    "# we have to redefine fps\n",
    "hfps = fpset.getFootprints()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prove the footprints are now indeed heavy\n",
    "hfps[0].isHeavy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all of the arrays here will be flattend 1D arrays of pixels from the footprint\n",
    "hfps[0].getImageArray()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the spanset 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 we will use, which can rearrange the image array into the proper 2 dimensional shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(fps[0].getSpans().unflatten(hfps[0].getImageArray()),\n",
    "           cmap='bone', origin='lower')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can do a simmilar exercise with the mask plane values, instead of the image plane values. We can first grab the (flattened) mask plane values of the footprint"
   ]
  },
  {
   "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": [
    "deepCoadd.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 and detected will be 2^5 + 2^4 = 48. Now we will visualize this in a similar manner to the imshow exercise we did before, only now we are *only* using data for the footprint because we are using the span."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "ax = plt.gca()\n",
    "# create an axes on the right side of ax. The width of cax will be 5%\n",
    "# of ax and the padding between cax and ax will be fixed at 0.05 inch.\n",
    "divider = make_axes_locatable(ax)\n",
    "im = plt.imshow(fps[0].getSpans().unflatten(hfps[0].getMaskArray()),\n",
    "                origin='lower')\n",
    "\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "\n",
    "plt.colorbar(im, cax=cax, ticks=[0, 32, 32+16])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Erosion and Dilation on Footprints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can grow or shrink footprints by using the `dilate()` and `erode()` methods. You may imagine doing this if you want to ensure you are getting all the flux associated with the detected source, or break apart one big footprint that contains many sources into several which each contain far fewer. Using these methods requires that we provide a scale by which to grow or shrink the footprint. This scale is often related to the PSF scale, e.g., scaling the RMS or $\\sigma$ of a gaussian PSF. \n",
    "\n",
    "In the dilation process, we can grow out footprint by including its surroundinb pixels. During erosion, we shrinking the object by excluding its outermost pixels. \n",
    "\n",
    "Let's demonstrate this below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first get the gaussian fwhm of the PSF \n",
    "psf = deepCoadd.getPsf()\n",
    "sigma = psf.computeShape().getDeterminantRadius()\n",
    "\n",
    "# print how many pixels belong to the footprint before we dilate it\n",
    "print('before dilating the footprint has an area of {}'.format(hfps[0].getArea()))\n",
    "\n",
    "# now we will grow the footprint by 2 x sigma. We need to pass an int\n",
    "hfps[0].dilate(int(2*sigma))\n",
    "\n",
    "# print how many pixels belong to the footprint ater dilation\n",
    "print('after dilating the footprint has an area of {}'.format(hfps[0].getArea()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lets chip away at the footprint by using the erode method now\n",
    "hfps[0].erode(int(2*sigma))\n",
    "\n",
    "# print how many pixels belong to the footprint ater dilation\n",
    "print('after eroding the footprint has an area of {}'.format(hfps[0].getArea()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To give us some more intuition on how dilation effects the footprints, we will show cut outs of several footprints before they have been dilated (top row) and after they have been dilated. To exagerate the effect, we will use a large scale length in the dilation. Note the change in the range of the x and y axes in these images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_fp(ax, fp, **kwargs):\n",
    "    span = fp.getSpans()\n",
    "    im_data = fp.getImageArray()\n",
    "    ax.imshow(span.unflatten(im_data), origin='lower', **kwargs)\n",
    "    return\n",
    "\n",
    "f, ax = plt.subplots(nrows=3, ncols=4, figsize=(10,6))\n",
    "\n",
    "# prepare a list with indicies of a few footprints\n",
    "fp_ids = [100,200,300,400]\n",
    "\n",
    "# make cut outs of footprints before dilation\n",
    "for i, j in zip(range(len(fp_ids)), fp_ids):\n",
    "    plot_fp(ax[0,i], hfps[j])\n",
    "\n",
    "# dilate footprints\n",
    "for j in fp_ids:\n",
    "    hfps[j].dilate(int(5*sigma))\n",
    "\n",
    "# we need to make the dilated footprints heavy again\n",
    "# updating the spans by dilation does not automatically update their associated image data\n",
    "fpset.makeHeavy(deepCoadd.getMaskedImage())\n",
    "hfps = fpset.getFootprints()\n",
    "\n",
    "# cut outs of footprints after dilation\n",
    "for i, j in zip(range(len(fp_ids)), fp_ids):\n",
    "    plot_fp(ax[1,i], hfps[j])\n",
    "    \n",
    "# now lets erode by the same scale so we can demonstrate that \n",
    "# erode and dilate are inverses of each other\n",
    "for j in fp_ids:\n",
    "    hfps[j].erode(int(5*sigma))\n",
    "\n",
    "# make them heavy again    \n",
    "fpset.makeHeavy(deepCoadd.getMaskedImage())\n",
    "hfps = fpset.getFootprints()\n",
    "\n",
    "# plot the last row\n",
    "for i, j in zip(range(len(fp_ids)), fp_ids):\n",
    "    plot_fp(ax[2,i], hfps[j])\n",
    "    \n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Footprints, Peaks and Deblending\n",
    "\n",
    "In the previous sections, we paid close attention to a fairly bright, isolated sources and their associated footprints that were identified as being above a 10$\\sigma$ threshold. In pratice, especially as we detect with lower thresholds, footprints from a first-pass threshold detection can contain several astronomical objects. This is because the wings of the intensity profiles of individual, but closely packed, sources can bleed into each other, forming one above threshold region, and thus one footprint. \n",
    "\n",
    "Formally, pixels which are local maxima in the smoothed image during detection are regarded as peaks, and each peak is assumed to be an astrophysical source. In the subsequent deblending process, the deblender looks at all the peaks in each footprints, and gives them their own deblended child footprint. In addition to identifying realestate to deblended children footprints, the deblender also allocates flux from the parent footprint to its children. As a result, the flux in a deblended heavy footprint will not necessarily equal the flux at the same pixels in its image.\n",
    "\n",
    "In this section we will examine footprints that contain more than one peak, and the interplay between footprints, peaks, and deblending. To do so, we will read footprints from data that has already gone through detection, deblending, and measurement. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following few cells, we will use a few `afwTable` tricks. Because we are interested in understanding peaks and deblending, we are going to sort the source catalog to find the source with the largest number of deblended children. You can learn more about `afwTables` in the [afwTables: A Guided Tour](https://github.com/LSSTScienceCollaborations/StackClub/blob/master/Basics/afw_table_guided_tour.ipynb) notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We grabbed the source catalogs at the start of the notebook\n",
    "print(type(cat))\n",
    "\n",
    "# First get the key for the deblend_nChild field, and id. We need these to sort on them\n",
    "nChild_key = cat.getSchema().find('deblend_nChild').key\n",
    "idKey = cat.getIdKey()\n",
    "\n",
    "# Sort the catalog on number of children\n",
    "cat.sort(nChild_key)\n",
    "# The catalog is sorted low to high, so grab the last element\n",
    "# Use the nChild key to get the number of children, and the id key to get the parent's ID\n",
    "num_children = cat[-1].get(nChild_key)\n",
    "parent_id = cat[-1].getId()\n",
    "\n",
    "# Footprints are stored in the source catalog, we can access them easily\n",
    "parent_fp = cat[-1].getFootprint()\n",
    "print('Parent source {} has {} deblended children'.format(parent_id, num_children))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The footprint of the parent object contains catalog of the peaks it contains. This catalog gives the peaks and thier locations and peak pixel value in the *smoothed* image. During the detection process, the image in question is smoothed with a spatial filter. A nice demonstration of this is avaliable in the [LowSurfaceBrightness notebook](https://nbviewer.jupyter.org/github/LSSTScienceCollaborations/StackClub/blob/rendered/SourceDetection/LowSurfaceBrightness.nbconvert.ipynb). The peaks' centers and counts are in this smoothed image, and may differ slightly from their subsequent children footprints after deblending."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parent_fp = cat[-1].getFootprint()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parent_fp.getPeaks()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will set up a visualization to see our parent source, the locations of its peaks (as red crosses), and the centers of the deblended children footprints (as orange circles)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the centroid of the parent\n",
    "parent_x = cat[-1].getX()\n",
    "parent_y = cat[-1].getY()\n",
    "\n",
    "parent_x = int(parent_x)\n",
    "parent_y = int(parent_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sort the catalog on ID so we can get a view of the \n",
    "# catalog that contains the children of the parent source\n",
    "cat.sort(idKey)\n",
    "child_cat = cat.getChildren(parent_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# credit to ADW for this afwDisplay snippet\n",
    "plt.figure()\n",
    "afw_display = afwDisplay.Display()\n",
    "afw_display.setMaskTransparency(100)\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(deepCoadd.getMaskedImage()[parent_x-200:parent_x+200, parent_y-200:parent_y+200])\n",
    "\n",
    "# We use display buffering to avoid re-drawing the image after each source is plotted\n",
    "with afw_display.Buffering():\n",
    "    # loop over child footprints, placing a red plus on their centers\n",
    "    for s in child_cat:\n",
    "        afw_display.dot('+', s.getX(), s.getY(), ctype=afwDisplay.RED)\n",
    "    # loop over peaks, placing an orange circle on their centers    \n",
    "    for p in parent_fp.getPeaks():\n",
    "        afw_display.dot('o', p.getIx(), p.getIy(), ctype=afwDisplay.ORANGE, size=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A few comments are in order about the above image.\n",
    "\n",
    "Even though you can identify by eye there is are several distinct sources here, the detection algorithm grouped them together into one big footprint. Remember, that is totally fine, as footprints can have multiple \"peaks\" that correspond to different sources. Take a look at one of the large galaxies. It has several peaks associated with it, even though by eye we might call it one large astrophysical object. The deblender shredded these sources into several children sources. This is a known failure mode of deblenders. Deblending is hard, and we will not attempt to solve it here, instead opting to restrict our attention to footprints. If you are interested in learning more about deblending beyond the scope of footprints, Fred Moolekamp has prepared [deblending notebooks](https://github.com/LSSTScienceCollaborations/StackClub/tree/master/Deblending) for the StackClub.\n",
    "\n",
    "The peaks are marked in the detection process. The peaks are formally detections of distinct astrophysical sources. Notice this is different than Source Extractor. although by eye you can tell they are likely just noise peaks. \n",
    "\n",
    "The peaks in the parent footprint (red crosses) and  the centers of the deblended children are almost all right about on top of each other. There is, however, a region where there is some miss-match. To see this more clearly, we will re-make the plot but zoom in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "afw_display = afwDisplay.Display()\n",
    "afw_display.setMaskTransparency(100)\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "#afw_display.mtv(deepCoadd.getMaskedImage()[3100:3200, 14600:14750])\n",
    "afw_display.mtv(deepCoadd.getMaskedImage()[parent_x-100:parent_x+100, parent_y-100:parent_y+100])\n",
    "\n",
    "# We use display buffering to avoid re-drawing the image after each source is plotted\n",
    "with afw_display.Buffering():\n",
    "    # loop over child footprints, placing a red plus on their centers\n",
    "    for s in child_cat:\n",
    "        afw_display.dot('+', s.getX(), s.getY(), ctype=afwDisplay.RED)\n",
    "    # loop over peaks, placing an orange circle on their centers    \n",
    "    for p in parent_fp.getPeaks():\n",
    "        afw_display.dot('o', p.getIx(), p.getIy(), ctype=afwDisplay.ORANGE, size=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Effects of Erosion on Peaks, and Contiguity\n",
    "\n",
    "Now that we have covered erosion/dilation and peaks, we will discuss the interplay between the two. Previously we saw we can change the shape of a footprint by using erosion and dilation operations. Because footprints also contain a list of peaks contained by its span, these operations will in turn effect peaks. Let's continue using the elliptical galaxy with many noise peaks from the previous section. \n",
    "\n",
    "We will fist plot the elliptical image data and overlay its peaks with red crosses. Then we will erode the eliptical and see the impact this has on peaks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a heavy footprint. This will let us grab pixel info associated with the footprint\n",
    "heavy_parent_fp = afwDetection.HeavyFootprintF(parent_fp, deepCoadd.getMaskedImage())\n",
    "\n",
    "f, ax = plt.subplots(ncols=2, nrows=1)\n",
    "plot_fp(ax[0], heavy_parent_fp,norm=LogNorm())\n",
    "# iterate over the peaks and overplot them on the image as red crosses\n",
    "for p in heavy_parent_fp.getPeaks():\n",
    "    ax[0].plot((p.getIx() - heavy_parent_fp.getBBox().getBeginX()),\n",
    "            (p.getIy() - heavy_parent_fp.getBBox().getBeginY()), 'r+')\n",
    "# the title will tell us how many peaks are in the FP before erosion\n",
    "ax[0].set_title('before erosion, {} peaks'.format(len(heavy_parent_fp.getPeaks())))\n",
    "\n",
    "# erode the footprint by a large ammount, then repeat the process.\n",
    "heavy_parent_fp.erode(10)\n",
    "# assign pixel values to the new eroded heavy footprint\n",
    "eroded_heavy_parent_fp = afwDetection.HeavyFootprintF(heavy_parent_fp, deepCoadd.getMaskedImage())\n",
    "\n",
    "plot_fp(ax[1], eroded_heavy_parent_fp,norm=LogNorm())\n",
    "\n",
    "for p in eroded_heavy_parent_fp.getPeaks():\n",
    "    ax[1].plot((p.getIx() - heavy_parent_fp.getBBox().getBeginX()),\n",
    "            (p.getIy() - heavy_parent_fp.getBBox().getBeginY()), 'r+')\n",
    "\n",
    "ax[1].set_title('after erosion, {} peaks'.format(len(eroded_heavy_parent_fp.getPeaks())))\n",
    "    \n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A few things happened to the footprint after erosion. \n",
    "\n",
    "1. Its area decreased, as expected\n",
    "\n",
    "2. The number of peaks decreased. When you use the `erode` method, any peaks that fall outside the footprint after erosion (called _orphan peaks_ in DM parlance) are discarded. If you happen to alter the span set of a footprint in antother way, you may also use the `removeOrphanPeaks` method by hand to cull any peaks that are outside the span.\n",
    "\n",
    "3. The footprint was broken up and is no longer contiguous. We can check the contiguity explicitly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eroded_heavy_parent_fp.isContiguous()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can remedy this by using the `split` method on the non-contiguous footprint. This will give us several contiguous footprints, in which all peaks from the initial peak are contained"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# split the non-contiguous fp into a vector of contiguous ones\n",
    "split_fps = eroded_heavy_parent_fp.split()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(nrows=1, ncols=2, figsize=(20, 7))\n",
    "\n",
    "# plot up the split fps\n",
    "plot_fp(ax[0],afwDetection.HeavyFootprintF(split_fps[0], deepCoadd.getMaskedImage()),norm=LogNorm())\n",
    "plot_fp(ax[1],afwDetection.HeavyFootprintF(split_fps[1], deepCoadd.getMaskedImage()),norm=LogNorm())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have established that deblending can go wrong, and give us unreliable footprints at times. But how does this impact the measurement process? How much flux actually goes to peaks? We will build our intuition for the impact of deblending on heavy footprints in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Footprints, Noise Replacer, and Deblending\n",
    "\n",
    "As we mentioned at the start of our notebook, footprints work hand in hand with detection, deblending, and measurement. In the following few cells we will loosely follow high level steps that are taken in single frame measurement, so we can understand how footprint deblending can impact measurements and science results.\n",
    "\n",
    "Refering again to [Bosch et al. (2017)](https://arxiv.org/pdf/1705.06766.pdf), the procedure for source measurement is\n",
    "\n",
    "> 1. We replace all Footprints in the image with random Gaussian noise with the same variance as the original noise in those pixels.\n",
    "> 2. We insert the deblended pixels for a particular source to be measured back into the image (replacing the noise pixels).\n",
    "> 3. We run all measurement plug-ins on the current source.\n",
    "> 4. We re-replace the Footprint of the current source with noise.\n",
    "> 5. Repeat this process for the next Footprint\n",
    "\n",
    "In our discussion below, we will concentrate on one parent footprint, and visualize its deblended children being swapped back in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will sort the source catalog and grab the object that has a bright psf flux and plenty of children\n",
    "cat.sort(cat.getPsfFluxSlot().getMeasKey())\n",
    "cat = cat.copy(deep=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the 55th brightest source should be bright enough\n",
    "# sources much brighter than this tend to be saturated stars,\n",
    "# but feel free to play around with this\n",
    "cat.sort(nChild_key)\n",
    "idx = -55\n",
    "# verify there are plenty of children\n",
    "cat[idx].get(nChild_key)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the centroid of this source so we can visualize it in a few steps\n",
    "x, y = cat[idx].getCentroid()\n",
    "ix = int(x)\n",
    "iy= int(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# grab the parent ID\n",
    "pId = cat[idx].getId()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we start playing around with swapping footprints in and out of the frame, lets take a look at the calexp as it is. We can see a nice bright source at the center, along with a few distinct fainter sources in the background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat.sort()\n",
    "children = cat.getChildren(pId)['id']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "children"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "afw_display = afwDisplay.Display()\n",
    "afw_display.setMaskTransparency(100)\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(deepCoadd.getMaskedImage()[ix - 100:ix + 100,iy - 100:iy + 100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will see how we can exchange the detection footprints with noise by using the `NoiseReplacer` object in the stack. We will emulate steps that are taken in `SingleFrameMeasurementTask`, where measurement algorithms are run on sources. The next two cells are inspired by a few lines from [sfm.py](https://github.com/lsst/meas_base/blob/master/python/lsst/meas/base/sfm.py#L305). We will first need to organize our footprints into a dictionary, where the source ID's are the keys and the associated value pairs are the footprints themselves.\n",
    "\n",
    "We will also point out that you can grab footprints directly from the source catalog."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare footprints in dictionary\n",
    "fp_dict = {measRecord.getId(): (measRecord.getParent(), measRecord.getFootprint())\n",
    "                      for measRecord in cat}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# instantiate NoiseReplacer. show it our calexp and footprints\n",
    "nr_config = NoiseReplacerConfig()\n",
    "noiseReplacer = NoiseReplacer(nr_config, deepCoadd, fp_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we create the `noiseReplacer` object, noise properties are calculated using the pixel values whose associated mask plane pixels have not been set. This noise distribution is used to replace the footprints with noise that is generated to be representative of the true noise properties of the deepCoadd. Now we will display the deepCoadd again to see the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "afw_display = afwDisplay.Display()\n",
    "afw_display.setMaskTransparency(100)\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(deepCoadd.getMaskedImage()[ix-100:ix+100,iy-100:iy+100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As advertised, we swapped out the data in our footprints with noise. As an aside, if you look closely you can see an 'edge' in the data where the parent object used to be, suggesting perhaps the faint tails of the elliptical's surface brightness profile were missed in detection and now remain in our noise image. We could imagine dilating the footprint to capture these.\n",
    "\n",
    "Now we will draw some more inspiration from [sfm.py](https://github.com/lsst/meas_base/blob/master/python/lsst/meas/base/sfm.py#L353) where we will insert our parent footprint of interest back into the image using its ID from the source catalog."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "noiseReplacer.insertSource(pId)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "afw_display = afwDisplay.Display()\n",
    "afw_display.setMaskTransparency(100)\n",
    "afw_display.scale('asinh', 'zscale')\n",
    "afw_display.mtv(deepCoadd.getMaskedImage()[ix-100:ix+100,iy-100:iy+100])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have used our `NoiseReplacer` instance to ensure this parent footprint is the only source in the image, and all other's have been replaced with noise. The sequence during the measurement process on parent and children footprints would roughly follow like this:\n",
    "\n",
    "+ The noise replacer swaps out all sources with noise\n",
    "+ A parent footprint is restored into the image\n",
    "+ Measurements are run on the parent footprint\n",
    "+ The parent footprint is replaced with noise\n",
    "+ A deblended child footprint is restored into the image\n",
    "+ Measurements are run on the first child\n",
    "\n",
    "This process would iterate over all parent and child sources. In the next cell we will make a small animation to let us visualize this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "# ims is a list of lists, each row is a list of artists to draw in the\n",
    "# current frame; here we are just animating one artist, the image, in\n",
    "# each frame\n",
    "ims = []\n",
    "\n",
    "src_ids = [pId] +  children.tolist()\n",
    "\n",
    "for i in src_ids:\n",
    "    noiseReplacer.insertSource(i)\n",
    "    data = deepCoadd.getMaskedImage()[ix-100:ix+100,iy-100:iy+100].getImage().array\n",
    "    im = plt.imshow(data, origin='lower', animated=True,\n",
    "                    vmin=0, vmax=5,\n",
    "                    visible=False)\n",
    "    ims.append([im])\n",
    "    noiseReplacer.removeSource(i)\n",
    "\n",
    "# the plt.close catches a spurious extra image\n",
    "plt.close()\n",
    "\n",
    "# makes an animation of the parent, then child footprints\n",
    "HTML(animation.ArtistAnimation(fig, ims, interval=1000, blit=True,\n",
    "                                repeat_delay=1000).to_jshtml())\n",
    "# takes a minute to get going, "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At first glance, in the animation above it seems the central bright galaxies and fainter smaller sources are successfully deblended. However, as the movie progresses we can see that the bright galaxy is shredded by the deblender. This is a common failure mode of the deblender. There are peaks in the outer wings of the galaxy; the deblender receives this list of peaks and assigns flux to them, creating these wispy looking sources that are inserted into the image. Because the deblender conserves the flux of the parent source, any flux allocated to the whispy sources is flux that cannot be allocated to the main child footprint of the main galaxy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will show a few cut outs  of the parent (top left), and first few deblended sources.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(nrows=2, ncols=3, sharey=True, sharex=True, figsize=(12,8))\n",
    "\n",
    "for i, a in zip(range(0, 6), ax.flatten()):\n",
    "    noiseReplacer.insertSource(src_ids[i])\n",
    "    data = deepCoadd.getMaskedImage()[ix-100:ix+100,iy-100:iy+100].getImage().array\n",
    "    a.imshow(data, origin='lower', animated=True)\n",
    "    noiseReplacer.removeSource(src_ids[i])\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is good pratice to use the end method for the `NoiseReplacer` after you are finished. This will restore our deepCoadd to the state it was in before we began replacing sources with noise, and so forth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "noiseReplacer.end()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "In this notebook, we learned how to make footprints from detection tasks, create heavy footprints, and access pixel level data from heavy footprints. We also learned the role footprints play in detection, deblending, and measurement."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Additional resources:\n",
    "\n",
    "1. [Documentation on footprints](https://github.com/lsst/afw/blob/b4f0cd2fac5183cac01ac5e8ac92be759afed8a3/doc/lsst.afw.detection/Footprint.rst) from the pipelines documentation in the afw package, courtesy Jonathan Sick: \n",
    "2. Footprint [Doxygen page](http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/classlsst_1_1afw_1_1detection_1_1_footprint.html)\n",
    "3. Jim Bosch's [HSC Software Pipeline paper](https://arxiv.org/pdf/1705.06766.pdf), particularly Section 4.8\n",
    "4. [Low Surface Brightness Stack Club Notebook](https://nbviewer.jupyter.org/github/LSSTScienceCollaborations/StackClub/blob/rendered/SourceDetection/LowSurfaceBrightness.nbconvert.ipynb) by Alex Drlica-Wagner"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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"
  },
  "livereveal": {
   "scroll": true,
   "start_slideshow_at": "selected"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}