Tickets/DM-42228: add donut quality check

doc/versionHistory.rst

@@ -6,6 +6,14 @@
 Version History
 ##################
 
+.. _lsst.ts.wep-9.6.0:
+
+-------------
+9.6.0
+-------------
+
+* Add donut image quality checking.
+
 .. _lsst.ts.wep-9.5.4:
 
 -------------

python/lsst/ts/wep/donutImageCheck.py

@@ -26,7 +26,7 @@
 
 
 class DonutImageCheck(object):
-    def __init__(self, numOfBins=256, entroThres=3.5):
+    def __init__(self, numOfBins=256, entroThres=3.5, returnEntro=False):
         """Donut image check class to judge the donut image is effective or
         not.
 
@@ -36,6 +36,9 @@ def __init__(self, numOfBins=256, entroThres=3.5):
             Number of bins in the histogram. (the default is 256.)
         entroThres : float, optional
             Threshold of entropy (the default is 3.5.)
+        returnEntro: bool, optional
+            Whether to return the calculated entropy value
+            (the default is False).
         """
 
         # Number of bins in the histogram
@@ -44,6 +47,9 @@ def __init__(self, numOfBins=256, entroThres=3.5):
         # Threshold of entropy
         self.entroThres = entroThres
 
+        # Whether to return the calculated value
+        self.returnEntro = returnEntro
+
     def isEffDonut(self, donutImg):
         """Is effective donut image or not.
 
@@ -64,9 +70,15 @@ def isEffDonut(self, donutImg):
         # Square the distribution to magnify the difference in entropy
         imgEntropy = entropy(hist**2)
         if (imgEntropy < self.entroThres) and (imgEntropy != 0):
-            return True
+            eff = True
+        else:
+            eff = False
+
+        # Return the actual entropy value if needed
+        if self.returnEntro:
+            return eff, imgEntropy
         else:
-            return False
+            return eff
 
 
 if __name__ == "__main__":

python/lsst/ts/wep/task/calcZernikesTask.py

@@ -30,6 +30,7 @@
 import lsst.pex.config as pexConfig
 import lsst.pipe.base as pipeBase
 import numpy as np
+from astropy.table import Table
 from lsst.pipe.base import connectionTypes
 from lsst.ts.wep.task.combineZernikesSigmaClipTask import CombineZernikesSigmaClipTask
 from lsst.ts.wep.task.donutStamps import DonutStamps
@@ -53,6 +54,7 @@ class CalcZernikesTaskConnections(
         storageClass="StampsBase",
         name="donutStampsIntra",
     )
+
     outputZernikesRaw = connectionTypes.Output(
         doc="Zernike Coefficients from all donuts",
         dimensions=("visit", "detector", "instrument"),
@@ -65,6 +67,18 @@ class CalcZernikesTaskConnections(
         storageClass="NumpyArray",
         name="zernikeEstimateAvg",
     )
+    donutsExtraQuality = connectionTypes.Output(
+        doc="Quality information for extra-focal donuts",
+        dimensions=("visit", "detector", "instrument"),
+        storageClass="DataFrame",
+        name="donutsExtraQuality",
+    )
+    donutsIntraQuality = connectionTypes.Output(
+        doc="Quality information for intra-focal donuts",
+        dimensions=("visit", "detector", "instrument"),
+        storageClass="DataFrame",
+        name="donutsIntraQuality",
+    )
 
 
 class CalcZernikesTaskConfig(
@@ -86,6 +100,31 @@ class CalcZernikesTaskConfig(
             + "is CombineZernikesSigmaClipTask.)"
         ),
     )
+    selectWithEntropy = pexConfig.Field(
+        dtype=bool,
+        default=False,
+        doc="Whether to use entropy in deciding to use the donut.",
+    )
+
+    selectWithSignalToNoise = pexConfig.Field(
+        dtype=bool,
+        default=False,
+        doc="Whether to use signal to noise ratio in deciding to use the donut.",
+    )
+
+    minSignalToNoise = pexConfig.Field(
+        dtype=float,
+        default=5000,
+        doc=str(
+            "The S/N threshold to use (keep donuts only above the threshold)."
+            + "The default S/N = 5000 corresponds to the default entropy threshold (3.5)"
+        ),
+    )
+    maxEntropy = pexConfig.Field(
+        dtype=float,
+        default=3.5,
+        doc=str("The entropy threshold to use (keep donuts only below the threshold)."),
+    )
 
 
 class CalcZernikesTask(pipeBase.PipelineTask, metaclass=abc.ABCMeta):
@@ -108,6 +147,12 @@ def __init__(self, **kwargs) -> None:
         self.combineZernikes = self.config.combineZernikes
         self.makeSubtask("combineZernikes")
 
+        # Inherit config parameters
+        self.selectWithEntropy = self.config.selectWithEntropy
+        self.selectWithSignalToNoise = self.config.selectWithSignalToNoise
+        self.minSignalToNoise = self.config.minSignalToNoise
+        self.maxEntropy = self.config.maxEntropy
+
     @timeMethod
     def run(
         self,
@@ -120,13 +165,75 @@ def run(
             return pipeBase.Struct(
                 outputZernikesRaw=np.full(19, np.nan),
                 outputZernikesAvg=np.full(19, np.nan),
+                donutsExtraQuality=[],
+                donutsIntraQuality=[],
             )
+        # Which donuts to use for Zernike estimation
+        # initiate these by selecting all donuts
+        entropySelectExtra = np.ones(len(donutStampsExtra), dtype="bool")
+        entropySelectIntra = np.ones(len(donutStampsIntra), dtype="bool")
+        entropyExtra = np.zeros(len(donutStampsIntra))
+        entropyIntra = np.zeros(len(donutStampsIntra))
+        if self.selectWithEntropy:
+            entropyExtra = np.asarray(donutStampsExtra.metadata.getArray("ENTROPY"))
+            entropyIntra = np.asarray(donutStampsIntra.metadata.getArray("ENTROPY"))
+            entropySelectExtra = entropyExtra < self.maxEntropy
+            entropySelectIntra = entropyIntra < self.maxEntropy
+
+        snSelectExtra = np.ones(len(donutStampsExtra), dtype="bool")
+        snSelectIntra = np.ones(len(donutStampsIntra), dtype="bool")
+        snExtra = np.zeros(len(donutStampsIntra))
+        snIntra = np.zeros(len(donutStampsIntra))
+        if self.selectWithSignalToNoise:
+            snExtra = np.asarray(donutStampsExtra.metadata.getArray("SN_VAR"))
+            snSelectExtra = self.minSignalToNoise < snExtra
+
+            snIntra = np.asarray(donutStampsIntra.metadata.getArray("SN_VAR"))
+            snSelectIntra = self.minSignalToNoise < snIntra
+
+        # AND condition : if both selectWithEntropy
+        # and selectWithSignalToNoise, then
+        # only donuts that pass with SN criterion as well
+        # as entropy criterion are selected
+        selectExtra = entropySelectExtra * snSelectExtra
+        selectIntra = entropySelectIntra * snSelectIntra
+
+        # store information about which extra-focal donuts were selected
+        donutsExtraQuality = Table(
+            data=[
+                snExtra,
+                entropyExtra,
+                entropySelectExtra,
+                snSelectExtra,
+                selectExtra,
+            ],
+            names=["SN", "ENTROPY", "ENTROPY_SELECT", "SN_SELECT", "FINAL_SELECT"],
+        ).to_pandas()
+
+        # store information about which intra-focal donuts were selected
+        donutsIntraQuality = Table(
+            data=[
+                snIntra,
+                entropyIntra,
+                entropySelectIntra,
+                snSelectIntra,
+                selectIntra,
+            ],
+            names=["SN", "ENTROPY", "ENTROPY_SELECT", "SN_SELECT", "FINAL_SELECT"],
+        ).to_pandas()
+
+        donutStampsExtraSelect = np.array(donutStampsExtra)[selectExtra]
+        donutStampsIntraSelect = np.array(donutStampsIntra)[selectIntra]
 
         # Estimate Zernikes from the collection of stamps
-        zkCoeffRaw = self.estimateZernikes.run(donutStampsExtra, donutStampsIntra)
+        zkCoeffRaw = self.estimateZernikes.run(
+            donutStampsExtraSelect, donutStampsIntraSelect
+        )
         zkCoeffCombined = self.combineZernikes.run(zkCoeffRaw.zernikes)
 
         return pipeBase.Struct(
             outputZernikesAvg=np.array(zkCoeffCombined.combinedZernikes),
             outputZernikesRaw=np.array(zkCoeffRaw.zernikes),
+            donutsExtraQuality=donutsExtraQuality,
+            donutsIntraQuality=donutsIntraQuality,
         )