Added initial astrometry source selector, base class and tests.

python/lsst/meas/algorithms/__init__.py

@@ -35,6 +35,8 @@
 from .subtractBackground import *
 from .measureApCorr import *
 from .flaggedStarSelector import *
+from .sourceSelector import *
+from .astrometrySourceSelector import *
 
 from .version import *
 

python/lsst/meas/algorithms/astrometrySourceSelector.py

@@ -0,0 +1,121 @@
+# See COPYRIGHT file at the top of the source tree.
+
+from __future__ import absolute_import, division, print_function
+
+import numpy as np
+
+from lsst.afw import table
+import lsst.pex.config as pexConfig
+from .sourceSelector import BaseSourceSelectorConfig, BaseSourceSelectorTask, sourceSelectorRegistry
+from lsst.pipe.base import Struct
+
+
+class AstrometrySourceSelectorConfig(BaseSourceSelectorConfig):
+    sourceFluxType = pexConfig.Field(
+        doc = "Type of source flux; typically one of Ap or Psf",
+        dtype = str,
+        default = "Ap",
+    )
+    minSnr = pexConfig.Field(
+        dtype = float,
+        doc = "Minimum allowed signal-to-noise ratio for sources used for matching "
+              "(in the flux specified by sourceFluxType); <= 0 for no limit",
+        default = 10,
+    )
+
+
+class AstrometrySourceSelectorTask(BaseSourceSelectorTask):
+    """
+    !Select sources that are useful for astrometry.
+
+    Good astrometry sources have high signal/noise, are non-blended, and
+    did not have certain "bad" flags set during source extraction. They need not
+    be PSF sources, just have reliable centroids.
+    """
+    ConfigClass = AstrometrySourceSelectorConfig
+
+    def __init__(self, *args, **kwargs):
+        BaseSourceSelectorTask.__init__(self, *args, **kwargs)
+
+    def selectSources(self, sourceCat, matches=None):
+        """
+        !Return a catalog of sources: a subset of sourceCat.
+
+        @param[in] sourceCat  catalog of sources that may be sources
+                                (an lsst.afw.table.SourceCatalog)
+
+        @return a pipeBase.Struct containing:
+        - sourceCat  a catalog of sources
+        """
+        self._getSchemaKeys(sourceCat.schema)
+
+        result = table.SourceCatalog(sourceCat.table)
+        for source in sourceCat:
+            if self._isGood(source) and not self._isBad(source):
+                result.append(source)
+        return Struct(sourceCat=result)
+
+    def _getSchemaKeys(self, schema):
+        """Extract and save the necessary keys from schema with asKey."""
+        self.parentKey = schema["parent"].asKey()
+        self.centroidKey = table.Point2DKey(schema["slot_Centroid"])
+        self.centroidFlagKey = schema["slot_Centroid_flag"].asKey()
+        self.nChildKey = schema["deblend_nChild"].asKey()
+
+        self.edgeKey = schema["base_PixelFlags_flag_edge"].asKey()
+        self.interpolatedCenterKey = schema["base_PixelFlags_flag_interpolatedCenter"].asKey()
+        self.saturatedKey = schema["base_PixelFlags_flag_saturated"].asKey()
+
+        fluxPrefix = "slot_%sFlux_" % (self.config.sourceFluxType,)
+        self.fluxField = fluxPrefix + "flux"
+        self.fluxKey = schema[fluxPrefix + "flux"].asKey()
+        self.fluxFlagKey = schema[fluxPrefix + "flag"].asKey()
+        self.fluxSigmaKey = schema[fluxPrefix + "fluxSigma"].asKey()
+
+    def _isMultiple(self, source):
+        """Return True if source is likely multiple sources."""
+        if (source.get(self.parentKey) != 0) or (source.get(self.nChildKey) != 0):
+            return True
+        footprint = source.getFootprint()
+        return footprint is not None and len(footprint.getPeaks()) > 1
+
+    def _hasCentroid(self, source):
+        """Return True if the source has a valid centroid"""
+        centroid = source.get(self.centroidKey)
+        return np.all(np.isfinite(centroid)) and not source.getCentroidFlag()
+
+    def _goodSN(self, source):
+        """Return True if source has Signal/Noise > config.minSnr."""
+        return (self.config.minSnr <= 0 or
+                (source.get(self.fluxKey)/source.get(self.fluxSigmaKey) > self.config.minSnr))
+
+    def _isUsable(self, source):
+        """
+        Return True if the source is usable for matching, even if it may have a poor centroid.
+
+        For a source to be usable it must:
+        - have a valid centroid
+        - not be deblended
+        - have a valid flux (of the type specified in this object's constructor)
+        - have adequate signal-to-noise
+        """
+        return self._hasCentroid(source) \
+            and not self._isMultiple(source) \
+            and not source.get(self.fluxFlagKey) \
+            and self._goodSN(source)
+
+    def _isGood(self, source):
+        """
+        Return True if source is usable for matching and likely has a good centroid.
+
+        The additional tests for a good centroid, beyond isUsable, are:
+        - not interpolated in the center (this includes saturated sources,
+            so we don't test separately for that)
+        - not near the edge
+        """
+        return self._isUsable(source) \
+            and not source.get(self.saturatedKey) \
+            and not source.get(self.interpolatedCenterKey) \
+            and not source.get(self.edgeKey)
+
+sourceSelectorRegistry.register("astrometry", AstrometrySourceSelectorTask)

python/lsst/meas/algorithms/sourceSelector.py

@@ -0,0 +1,79 @@
+# See COPYRIGHT file at the top of the source tree.
+
+from __future__ import absolute_import, division, print_function
+
+import abc
+import numpy as np
+
+from lsst.afw import table
+import lsst.pex.config as pexConfig
+import lsst.pipe.base as pipeBase
+
+__all__ = ["BaseSourceSelectorConfig", "BaseSourceSelectorTask", "sourceSelectorRegistry"]
+
+
+class BaseSourceSelectorConfig(pexConfig.Config):
+    badFlags = pexConfig.ListField(
+        doc = "List of flags which cause a source to be rejected as bad",
+        dtype = str,
+        default = [
+            "base_PixelFlags_flag_edge",
+            "base_PixelFlags_flag_interpolatedCenter",
+            "base_PixelFlags_flag_saturatedCenter",
+            "base_PixelFlags_flag_crCenter",
+            "base_PixelFlags_flag_bad",
+            "base_PixelFlags_flag_interpolated",
+        ],
+    )
+
+
+class BaseSourceSelectorTask(pipeBase.Task):
+    """!Base class for source selectors
+
+    Register all source selectors with the sourceSelectorRegistry using:
+        sourceSelectorRegistry.register(name, class)
+    """
+    __metaclass__ = abc.ABCMeta
+
+    ConfigClass = BaseSourceSelectorConfig
+    _DefaultName = "sourceSelector"
+
+    def __init__(self, **kwargs):
+        """!Initialize a source selector."""
+        pipeBase.Task.__init__(self, **kwargs)
+
+    def run(self, sourceCat, maskedImage=None, **kwargs):
+        """!Select sources and return them.
+
+        @param[in] sourceCat  catalog of sources that may be sources (an lsst.afw.table.SourceCatalog)
+        @param[in] maskedImage  the maskedImage containing the sources, for plotting.
+
+        @return an lsst.pipe.base.Struct containing:
+        - sourceCat  catalog of sources that were selected
+        """
+        return self.selectSources(maskedImage=maskedImage, sourceCat=sourceCat, **kwargs)
+
+    @abc.abstractmethod
+    def selectSources(self, sourceCat, matches=None):
+        """!Return a catalog of sources: a subset of sourceCat.
+
+        @param[in] sourceCat  catalog of sources that may be sources (an lsst.afw.table.SourceCatalog)
+
+        @return a pipeBase.Struct containing:
+        - sourceCat  a catalog of sources
+        """
+
+        # NOTE: example implementation, returning all sources that have no bad flags set.
+        result = table.SourceCatalog(sourceCat.table)
+        for source in sourceCat:
+            if not self._isBad(source):
+                result.append(source)
+        return pipeBase.Struct(sourceCat=result)
+
+    def _isBad(self, source):
+        """Return True if any of config.badFlags are set for this source."""
+        return reduce(lambda x, y: np.logical_or(x, source.get(y)), self.config.badFlags, False)
+
+sourceSelectorRegistry = pexConfig.makeRegistry(
+    doc="A registry of source selectors (subclasses of BaseSourceSelectorTask)",
+)

tests/testAstrometrySourceSelector.py

@@ -0,0 +1,121 @@
+# See COPYRIGHT file at the top of the source tree.
+
+from __future__ import division, absolute_import, print_function
+
+import unittest
+import numpy as np
+
+import lsst.utils
+from lsst.meas.base.tests import TestDataset
+from lsst.afw import table
+from lsst.meas.algorithms import sourceSelector
+
+
+badFlags = ["base_PixelFlags_flag_edge",
+            "base_PixelFlags_flag_interpolatedCenter",
+            "base_PixelFlags_flag_saturated",
+            "base_PixelFlags_flag_saturatedCenter",
+            "base_PixelFlags_flag_crCenter",
+            "base_PixelFlags_flag_bad",
+            "base_PixelFlags_flag_interpolated",
+            "slot_Centroid_flag",
+            "slot_ApFlux_flag",
+            ]
+
+
+def add_good_source(src, num=0):
+    """
+    Insert a likely-good source into the catalog.
+
+    num is added to various values to distinguish them in catalogs with multiple objects.
+    """
+    src.addNew()
+    src['coord_ra'][-1] = 1.+num
+    src['coord_dec'][-1] = 2.+num
+    src['slot_Centroid_x'][-1] = 10.+num
+    src['slot_Centroid_y'][-1] = 20.+num
+    src['slot_ApFlux_flux'][-1] = 100.+num
+    src['slot_ApFlux_fluxSigma'][-1] = 1.
+
+
+class TestAstrometrySourceSelector(lsst.utils.tests.TestCase):
+    def setUp(self):
+        schema = TestDataset.makeMinimalSchema()
+        schema.addField("truth_flag", type="Flag")  # for slot_Centroid_flag
+        schema.addField("slot_ApFlux_flux", type=float)
+        schema.addField("slot_ApFlux_fluxSigma", type=float)
+        for flag in badFlags:
+            schema.addField(flag, type="Flag")
+
+        self.src = table.SourceCatalog(schema)
+        self.sourceSelector = sourceSelector.sourceSelectorRegistry['astrometry']()
+
+    def tearDown(self):
+        del self.src
+        del self.sourceSelector
+
+    def testSelectSources_good(self):
+        for i in range(5):
+            add_good_source(self.src, i)
+        result = self.sourceSelector.selectSources(self.src)
+        # TODO: assertEqual doesn't work right on source catalogs.
+        # self.assertEqual(result.sourceCat, self.src)
+        for x in self.src['id']:
+            self.assertIn(x, result.sourceCat['id'])
+
+    def testSelectSources_bad(self):
+        for i, flag in enumerate(badFlags):
+            add_good_source(self.src, i)
+            self.src[i].set(flag, True)
+        result = self.sourceSelector.selectSources(self.src)
+        for i, x in enumerate(self.src['id']):
+            self.assertNotIn(x, result.sourceCat['id'], "should not have found %s"%badFlags[i])
+
+    def testSelectSources_bad_centroid(self):
+        add_good_source(self.src, 1)
+        self.src[0].set('slot_Centroid_x', np.nan)
+        result = self.sourceSelector.selectSources(self.src)
+        self.assertNotIn(self.src['id'][0], result.sourceCat['id'])
+
+    def testSelectSources_is_parent(self):
+        add_good_source(self.src, 1)
+        self.src[0].set('parent', 1)
+        result = self.sourceSelector.selectSources(self.src)
+        self.assertNotIn(self.src['id'][0], result.sourceCat['id'])
+
+    def testSelectSources_has_children(self):
+        add_good_source(self.src, 1)
+        self.src[0].set('deblend_nChild', 1)
+        result = self.sourceSelector.selectSources(self.src)
+        self.assertNotIn(self.src['id'][0], result.sourceCat['id'])
+
+    def testSelectSources_highSN_cut(self):
+        add_good_source(self.src, 1)
+        add_good_source(self.src, 2)
+        self.src['slot_ApFlux_flux'][0] = 20.
+        self.src['slot_ApFlux_flux'][1] = 1000.
+
+        self.sourceSelector.config.minSnr = 100
+        result = self.sourceSelector.selectSources(self.src)
+        self.assertNotIn(self.src[0]['id'], result.sourceCat['id'])
+        self.assertIn(self.src[1]['id'], result.sourceCat['id'])
+
+    def testSelectSources_no_SN_cut(self):
+        self.sourceSelector.config.minSnr = 0
+        add_good_source(self.src, 1)
+        self.src['slot_ApFlux_flux'][0] = 0
+        result = self.sourceSelector.selectSources(self.src)
+        self.assertIn(self.src[0]['id'], result.sourceCat['id'])
+
+
+# for MemoryTestCase
+def setup_module(module):
+    lsst.utils.tests.init()
+
+
+class MyMemoryTestCase(lsst.utils.tests.MemoryTestCase):
+    pass
+
+if __name__ == "__main__":
+    lsst.utils.tests.init()
+    unittest.main()