Merge pull request #74 from lsst-ts/issue/67/DonutDetector

Add DonutDetector code

doc/content.rst

@@ -33,6 +33,7 @@ This module is a high-level module to use other modules.
 * **ParamReader**: Parameter reader class to read the yaml configuration files used in the calculation.
 * **DonutImageCheck**: Donut image check class to judge the donut image is effective or not.
 * **CreatePhosimDonutTemplates**: Create donut templates on camera detectors using Phosim. See :doc:`here <phosimDonutTemplates>` for more information on generating Phosim donut templates.
+* **DonutDetector**: Detect donuts directly from an out of focus image by convolution with a template image.
 
 .. _lsst.ts.wep-modules_wep_bsc:
 

doc/uml/wepClass.uml

@@ -1,5 +1,6 @@
 @startuml
 class DonutImageCheck
+class DonutDetector
 WepController *-- ButlerWrapper
 WepController *-- CamDataCollector
 WepController *-- CamIsrWrapper

doc/versionHistory.rst

@@ -6,6 +6,14 @@
 Version History
 ##################
 
+.. _lsst.ts.wep-1.5.5:
+
+-------------
+1.5.5
+-------------
+
+* Add `DonutDetector` class.
+
 .. _lsst.ts.wep-1.5.4:
 
 -------------

python/lsst/ts/wep/DonutDetector.py

@@ -0,0 +1,149 @@
+# This file is part of ts_wep.
+#
+# Developed for the LSST Telescope and Site Systems.
+# This product includes software developed by the LSST Project
+# (https://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+import numpy as np
+import pandas as pd
+from copy import copy
+
+
+from lsst.ts.wep.Utility import CentroidFindType
+from lsst.ts.wep.cwfs.CentroidFindFactory import CentroidFindFactory
+from scipy.spatial.distance import cdist
+
+
+class DonutDetector(object):
+    """
+    Class to detect donuts directly from an out of focus image
+    by convolution with a template image.
+    """
+
+    def detectDonuts(
+        self, expArray, template, blendRadius, peakThreshold=0.95, dbscanEps=5
+    ):
+        """
+        Detect and categorize donut sources as blended/unblended
+
+        Parameters
+        -------
+        expArray: numpy.ndarray
+            The input image data.
+        template: numpy.ndarray
+            Donut template appropriate for the image.
+        blendRadius: float
+            Minimum distance in pixels two donut centers need to
+            be apart in order to be tagged as unblended.
+        peakThreshold: float, optional
+            This value is a specifies a number between 0 and 1 that is
+            the fraction of the highest pixel value in the convolved image.
+            The code then sets all pixels with a value below this to 0 before
+            running the K-means algorithm to find peaks that represent possible
+            donut locations. (The default is 0.95)
+        dbscanEps: float, optional
+            Maximum distance scikit-learn DBSCAN algorithm allows "between two
+            samples for one to considered in the neighborhood of the other".
+            (The default is 5.0)
+
+        Returns
+        -------
+        pandas.DataFrame
+            Dataframe identifying donut positions and if they
+            are blended with other donuts. If blended also identfies
+            which donuts are blended with which.
+        """
+
+        centroidFinder = CentroidFindFactory.createCentroidFind(
+            CentroidFindType.ConvolveTemplate
+        )
+        binaryExp = centroidFinder.getImgBinary(copy(expArray))
+        centroidX, centroidY, donutRad = centroidFinder.getCenterAndRfromTemplateConv(
+            binaryExp,
+            templateImgBinary=template,
+            nDonuts=-1,
+            peakThreshold=peakThreshold,
+            dbscanEps=dbscanEps,
+        )
+
+        donutDf = pd.DataFrame(
+            np.array([centroidX, centroidY]).T, columns=["x_center", "y_center"]
+        )
+        donutDf = self.identifyBlendedDonuts(donutDf, blendRadius)
+
+        return donutDf
+
+    def identifyBlendedDonuts(self, donutDf, blendRadius):
+        """
+        Label donuts as blended/unblended if the centroids are within
+        the blendRadius number of pixels.
+
+        Parameters
+        ----------
+        donutDf: pandas.DataFrame
+            Dataframe identifying donut positions with labels
+            'x_center' and 'y_center'.
+        blendRadius: float
+            Minimum distance in pixels two donut centers need to
+            be apart in order to be tagged as unblended.
+
+        Returns
+        -------
+        pandas.DataFrame
+            Dataframe identifying donut positions and if they
+            are blended with other donuts. If blended also identfies
+            which donuts are blended with which.
+        """
+
+        # Find distances between each pair of objects
+        donutCenters = [donutDf["x_center"].values, donutDf["y_center"].values]
+        donutCenters = np.array(donutCenters).T
+        distMatrix = cdist(donutCenters, donutCenters)
+        # Don't need repeats of each pair
+        distMatrixUpper = np.triu(distMatrix)
+
+        # Identify blended pairs of objects by distance
+        blendedPairs = np.array(
+            np.where((distMatrixUpper > 0.0) & (distMatrixUpper < blendRadius))
+        ).T
+        blendedCenters = np.unique(blendedPairs.flatten())
+
+        # Add blended information into dataframe
+        donutDf["blended"] = False
+        donutDf.loc[blendedCenters, "blended"] = True
+        donutDf["blended_with"] = None
+        for donutOne, donutTwo in blendedPairs:
+            if donutDf.loc[donutOne, "blended_with"] is None:
+                donutDf.at[donutOne, "blended_with"] = []
+            if donutDf.loc[donutTwo, "blended_with"] is None:
+                donutDf.at[donutTwo, "blended_with"] = []
+            donutDf.loc[donutOne, "blended_with"].append(donutTwo)
+            donutDf.loc[donutTwo, "blended_with"].append(donutOne)
+
+        # Count the number of other donuts overlapping
+        # each donut
+        donutDf["num_blended_neighbors"] = 0
+        for donutIdx in range(len(donutDf)):
+            if donutDf["blended_with"].iloc[donutIdx] is None:
+                continue
+
+            donutDf.at[donutIdx, "num_blended_neighbors"] = len(
+                donutDf["blended_with"].loc[donutIdx]
+            )
+
+        return donutDf

python/lsst/ts/wep/cwfs/CentroidConvolveTemplate.py

@@ -24,7 +24,7 @@
 from lsst.ts.wep.cwfs.CentroidDefault import CentroidDefault
 from lsst.ts.wep.cwfs.CentroidRandomWalk import CentroidRandomWalk
 from scipy.signal import correlate
-from sklearn.cluster import KMeans
+from sklearn.cluster import KMeans, DBSCAN
 
 
 class CentroidConvolveTemplate(CentroidDefault):
@@ -128,7 +128,12 @@ def getCenterAndRfromImgBinary(
         return x[0], y[0], radius
 
     def getCenterAndRfromTemplateConv(
-        self, imageBinary, templateImgBinary=None, nDonuts=1, peakThreshold=0.95
+        self,
+        imageBinary,
+        templateImgBinary=None,
+        nDonuts=1,
+        peakThreshold=0.95,
+        dbscanEps=5.0,
     ):
         """
         Get the centers of the donuts by convolving a binary template image
@@ -147,14 +152,22 @@ def getCenterAndRfromTemplateConv(
             Binary image of template donut. If set to None then the image
             is convolved with itself. (The default is None)
         nDonuts: int, optional
-            Number of donuts there should be in the binary image. Needs to
-            be >= 1. (The default is 1)
+            Number of donuts there should be in the binary image. If the number
+            is >= 1 then K-Means clustering will be used to return the
+            specified number of donut centers. However, this can also be set to
+            -1 if the number of donuts is unknown and it will perform DBSCAN
+            clustering to find and return a set of donut centers.
+            (The default is 1)
         peakThreshold: float, optional
             This value is a specifies a number between 0 and 1 that is
             the fraction of the highest pixel value in the convolved image.
             The code then sets all pixels with a value below this to 0 before
             running the K-means algorithm to find peaks that represent possible
             donut locations. (The default is 0.95)
+        dbscanEps: float, optional
+            Maximum distance scikit-learn DBSCAN algorithm allows "between two
+            samples for one to considered in the neighborhood of the other".
+            (The default is 5.0)
 
         Returns
         -------
@@ -169,8 +182,10 @@ def getCenterAndRfromTemplateConv(
         if templateImgBinary is None:
             templateImgBinary = copy(imageBinary)
 
-        nDonutsAssertStr = "nDonuts must be an integer >= 1"
-        assert (nDonuts >= 1) & (type(nDonuts) is int), nDonutsAssertStr
+        nDonutsAssertStr = "nDonuts must be an integer >= 1 or -1"
+        assert ((nDonuts >= 1) | (nDonuts == -1)) & (
+            type(nDonuts) is int
+        ), nDonutsAssertStr
 
         # We set the mode to be "same" because we need to return the same
         # size image to the code.
@@ -185,20 +200,36 @@ def getCenterAndRfromTemplateConv(
         rankedConvolveCutoff = rankedConvolve[:cutoff]
         nx, ny = np.unravel_index(rankedConvolveCutoff, np.shape(imageBinary))
 
-        # Then to find peaks in the image we use K-Means with the
-        # specified number of donuts
-        kmeans = KMeans(n_clusters=nDonuts)
-        labels = kmeans.fit_predict(np.array([nx, ny]).T)
-
-        # Then in each cluster we take the brightest pixel as the centroid
+        # Donut centers lists
         centX = []
         centY = []
-        for labelNum in range(nDonuts):
-            nxLabel, nyLabel = np.unravel_index(
-                rankedConvolveCutoff[labels == labelNum][0], np.shape(imageBinary)
-            )
-            centX.append(nxLabel)
-            centY.append(nyLabel)
+
+        if nDonuts >= 1:
+            # Then to find peaks in the image we use K-Means with the
+            # specified number of donuts
+            kmeans = KMeans(n_clusters=nDonuts)
+            labels = kmeans.fit_predict(np.array([nx, ny]).T)
+
+            # Then in each cluster we take the brightest pixel as the centroid
+            for labelNum in range(nDonuts):
+                nxLabel, nyLabel = np.unravel_index(
+                    rankedConvolveCutoff[labels == labelNum][0], np.shape(imageBinary)
+                )
+                centX.append(nxLabel)
+                centY.append(nyLabel)
+        elif nDonuts == -1:
+            # Use DBSCAN to find clusters of points when the
+            # number of donuts is unknown
+            labels = DBSCAN(eps=dbscanEps).fit_predict(np.array([ny, nx]).T)
+
+            # Save the centroid as the brightest pixel
+            # within each identified cluster
+            for labelNum in np.unique(labels):
+                nxLabel, nyLabel = np.unravel_index(
+                    rankedConvolveCutoff[labels == labelNum][0], np.shape(imageBinary)
+                )
+                centX.append(nxLabel)
+                centY.append(nyLabel)
 
         # Get the radius of the donut from the template image
         radius = np.sqrt(np.sum(templateImgBinary) / np.pi)

tests/cwfs/test_centroidConvolveTemplate.py

@@ -150,23 +150,23 @@ def testNDonutsAssertion(self):
 
         singleDonut, doubleDonut, eff_radius = self._createData(20, 40, 160)
 
-        nDonutsAssertMsg = "nDonuts must be an integer >= 1"
+        nDonutsAssertMsg = "nDonuts must be an integer >= 1 or -1"
         with self.assertRaises(AssertionError, msg=nDonutsAssertMsg):
             cX, cY, rad = self.centroidConv.getCenterAndRfromTemplateConv(
                 singleDonut, nDonuts=0
             )
 
         with self.assertRaises(AssertionError, msg=nDonutsAssertMsg):
             cX, cY, rad = self.centroidConv.getCenterAndRfromTemplateConv(
-                singleDonut, nDonuts=-1
+                singleDonut, nDonuts=-2
             )
 
         with self.assertRaises(AssertionError, msg=nDonutsAssertMsg):
             cX, cY, rad = self.centroidConv.getCenterAndRfromTemplateConv(
                 singleDonut, nDonuts=1.5
             )
 
-    def testGetCenterAndRFromTemplateConv(self):
+    def testGetCenterAndRFromTemplateConvKMeans(self):
 
         singleDonut, doubleDonut, eff_radius = self._createData(20, 40, 160)
 
@@ -186,6 +186,26 @@ def testGetCenterAndRFromTemplateConv(self):
         self.assertEqual(doubleCY, [80.0, 80.0])
         self.assertAlmostEqual(rad, eff_radius, delta=0.1)
 
+    def testGetCenterAndRFromTemplateConvDBSCAN(self):
+
+        singleDonut, doubleDonut, eff_radius = self._createData(20, 40, 160)
+
+        # Test recovery of single donut
+        singleCX, singleCY, rad = self.centroidConv.getCenterAndRfromTemplateConv(
+            singleDonut, nDonuts=-1
+        )
+        self.assertEqual(singleCX, [80.0])
+        self.assertEqual(singleCY, [80.0])
+        self.assertAlmostEqual(rad, eff_radius, delta=0.1)
+
+        # Test recovery of two donuts at once
+        doubleCX, doubleCY, rad = self.centroidConv.getCenterAndRfromTemplateConv(
+            doubleDonut, templateImgBinary=singleDonut, nDonuts=-1
+        )
+        self.assertCountEqual(doubleCX, [50.0, 110.0])
+        self.assertEqual(doubleCY, [80.0, 80.0])
+        self.assertAlmostEqual(rad, eff_radius, delta=0.1)
+
 
 if __name__ == "__main__":