Merge pull request #21 from lsst-sitcom/tickets/DM-36403

DM-36403: Add utility function for summing footprint fluxes

python/lsst/summit/utils/utils.py

@@ -20,11 +20,13 @@
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 import os
+from typing import Iterable
 import numpy as np
 import logging
 from scipy.ndimage import gaussian_filter
 import lsst.afw.image as afwImage
 import lsst.afw.detection as afwDetect
+from lsst.afw.detection import Footprint, FootprintSet
 import lsst.afw.math as afwMath
 import lsst.daf.base as dafBase
 import lsst.geom as geom
@@ -308,6 +310,80 @@ def detectObjectsInExp(exp, nSigma=10, nPixMin=10, grow=0):
     return footPrintSet
 
 
+def fluxesFromFootprints(footprints, parentImage, subtractImageMedian=False):
+    """Calculate the flux from a set of footprints, given the parent image,
+    optionally subtracting the whole-image median from each pixel as a very
+    rough background subtraction.
+
+    Parameters
+    ----------
+    footprints : `lsst.afw.detection.FootprintSet` or
+                 `lsst.afw.detection.Footprint` or
+                 `iterable` of `lsst.afw.detection.Footprint`
+        The footprints to measure.
+    parentImage : `lsst.afw.image.Image`
+        The parent image.
+    subtractImageMedian : `bool`, optional
+        Subtract a whole-image median from each pixel in the footprint when
+        summing as a very crude background subtraction. Does not change the
+        original image.
+
+    Returns
+    -------
+    fluxes : `list` of `float`
+        The fluxes for each footprint.
+
+    Raises
+    ------
+    TypeError : raise for unsupported types.
+    """
+    median = 0
+    if subtractImageMedian:
+        median = np.nanmedian(parentImage.array)
+
+    # poor person's single dispatch
+    badTypeMsg = ("This function works with FootprintSets, single Footprints, and iterables of Footprints. "
+                  f"Got {type(footprints)}: {footprints}")
+    if isinstance(footprints, FootprintSet):
+        footprints = footprints.getFootprints()
+    elif isinstance(footprints, Iterable):
+        if not isinstance(footprints[0], Footprint):
+            raise TypeError(badTypeMsg)
+    elif isinstance(footprints, Footprint):
+        footprints = [footprints]
+    else:
+        raise TypeError(badTypeMsg)
+
+    return np.array([fluxFromFootprint(fp, parentImage, backgroundValue=median) for fp in footprints])
+
+
+def fluxFromFootprint(footprint, parentImage, backgroundValue=0):
+    """Calculate the flux from a footprint, given the parent image, optionally
+    subtracting a single value from each pixel as a very rough background
+    subtraction, e.g. the image median.
+
+    Parameters
+    ----------
+    footprint : `lsst.afw.detection.Footprint`
+        The footprint to measure.
+    parentImage : `lsst.afw.image.Image`
+        Image containing the footprint.
+    backgroundValue : `bool`, optional
+        The value to subtract from each pixel in the footprint when summing
+        as a very crude background subtraction. Does not change the original
+        image.
+
+    Returns
+    -------
+    flux : `float`
+        The flux in the footprint
+    """
+    if backgroundValue:  # only do the subtraction if non-zero for speed
+        xy0 = parentImage.getBBox().getMin()
+        return footprint.computeFluxFromArray(parentImage.array - backgroundValue, xy0)
+    return footprint.computeFluxFromImage(parentImage)
+
+
 def humanNameForCelestialObject(objName):
     """Returns a list of all human names for obj, or [] if none are found.
 

tests/test_utils.py

@@ -23,25 +23,31 @@
 
 import copy
 import itertools
+from typing import Iterable
 import unittest
 
 import astropy.time
 import astropy.units as u
+import lsst.afw.detection as afwDetect
 import lsst.afw.image as afwImage
+import lsst.afw.geom as afwGeom
 import lsst.geom as geom
 import lsst.utils.tests
 import numpy as np
 from astro_metadata_translator import makeObservationInfo
 from lsst.obs.base import createInitialSkyWcsFromBoresight
 from lsst.obs.base.makeRawVisitInfoViaObsInfo import MakeRawVisitInfoViaObsInfo
 from lsst.obs.lsst import Latiss
+
 from lsst.summit.utils.utils import (getExpPositionOffset,
                                      getFieldNameAndTileNumber,
                                      getAirmassSeeingCorrection,
                                      getFilterSeeingCorrection,
                                      quickSmooth,
                                      getQuantiles,
+                                     fluxesFromFootprints,
                                      )
+
 from lsst.obs.lsst.translators.latiss import AUXTEL_LOCATION
 
 
@@ -207,6 +213,62 @@ def test_quantiles(self):
             np.testing.assert_almost_equal(edges1, edges2, decimal=decimal)
 
 
+class ImageBasedTestCase(lsst.utils.tests.TestCase):
+    def test_fluxFromFootprint(self):
+        image = afwImage.Image(
+            np.arange(8100, dtype=np.int32).reshape(90, 90),
+            xy0=lsst.geom.Point2I(10, 12),
+            dtype="I"
+        )
+
+        radius = 3
+        spans = afwGeom.SpanSet.fromShape(radius, afwGeom.Stencil.CIRCLE, offset=(27, 30))
+        footprint1 = afwDetect.Footprint(spans)
+
+        # The extracted footprint should be the same as the product of the
+        # spans and the overlapped bow with the image
+        truth1 = spans.asArray() * image.array[15:22, 14:21]
+
+        radius = 3
+        spans = afwGeom.SpanSet.fromShape(radius, afwGeom.Stencil.CIRCLE, offset=(44, 49))
+        footprint2 = afwDetect.Footprint(spans)
+        truth2 = spans.asArray() * image.array[34:41, 31:38]
+
+        allFootprints = [footprint1, footprint2]
+        footprintSet = afwDetect.FootprintSet(image.getBBox())
+        footprintSet.setFootprints(allFootprints)
+
+        # check it can accept a footprintSet, and single and iterables of
+        # footprints
+        with self.assertRaises(TypeError):
+            fluxesFromFootprints(10, image)
+
+        with self.assertRaises(TypeError):
+            fluxesFromFootprints([8, 6, 7, 5, 3, 0, 9], image)
+
+        # check the footPrintSet
+        fluxes = fluxesFromFootprints(footprintSet, image)
+        expectedLength = len(footprintSet.getFootprints())
+        self.assertEqual(len(fluxes), expectedLength)  # always one flux per footprint
+        self.assertIsInstance(fluxes, Iterable)
+        self.assertAlmostEqual(fluxes[0], np.sum(truth1))
+        self.assertAlmostEqual(fluxes[1], np.sum(truth2))
+
+        # check the list of footprints
+        fluxes = fluxesFromFootprints(allFootprints, image)
+        expectedLength = 2
+        self.assertEqual(len(fluxes), expectedLength)  # always one flux per footprint
+        self.assertIsInstance(fluxes, Iterable)
+        self.assertAlmostEqual(fluxes[0], np.sum(truth1))
+        self.assertAlmostEqual(fluxes[1], np.sum(truth2))
+
+        # ensure that subtracting the image median from fluxes leave image
+        # pixels untouched
+        oldImageArray = copy.deepcopy(image.array)
+        fluxes = fluxesFromFootprints(footprintSet, image, subtractImageMedian=True)
+        np.testing.assert_array_equal(image.array, oldImageArray)
+
+
 class TestMemory(lsst.utils.tests.MemoryTestCase):
     pass