Merge pull request #138 from lsst-camera-dh/LSSTTD-1580_fe55_gain_fitting_improvements

Lssttd 1580 fe55 gain fitting improvements

Author: jchiang87

python/lsst/eotest/sensor/MaskedCCD.py

@@ -6,6 +6,7 @@
 @author J. Chiang <jchiang@slac.stanford.edu>
 """
 import warnings
+import numpy as np
 import astropy.io.fits as fits
 from astropy.utils.exceptions import AstropyWarning, AstropyUserWarning
 from lsst.eotest.fitsTools import fitsWriteto
@@ -227,7 +228,9 @@ def bias_image(self, amp, overscan=None, **kwargs):
             return my_image.Factory(my_image, deep=True)
         return self.bias_image_using_overscan(amp, overscan=overscan, **kwargs)
 
-    def bias_subtracted_image(self, amp, overscan=None, **kwargs):
+    def bias_subtracted_image(self, amp, overscan=None,
+                              max_pca_reduced_chisq=None,
+                              **kwargs):
         """
         Subtract a bias image to correct for the offset. A bias correction is also
         applied if a base_frame is passed. The bias image with the offset values is
@@ -252,18 +255,24 @@ def bias_subtracted_image(self, amp, overscan=None, **kwargs):
             # Make a deep copy of the bias frame.
             bias = self.bias_frame[amp].Factory(self.bias_frame[amp], deep=True)
             # Subtract x-independent component using overscan.
-            bias -= \
-                self.bias_frame.bias_image_using_overscan(amp,
-                                                          overscan=overscan, **kwargs)
+            bias -= self.bias_frame.bias_image_using_overscan(
+                amp, overscan=overscan, **kwargs)
             # Subtract x-independent component of image for this amp
             # using overscan.
             my_image -= \
                 self.bias_image_using_overscan(amp, overscan=overscan, **kwargs)
             # Subtract structured, x-dependent part.
             my_image -= bias
         elif self.ccd_pcas is not None:
-            my_image.getImage().array[:] \
-                -= self.ccd_pcas.pca_bias_correction(amp, my_image.getImage().array)
+            bias_model = self.ccd_pcas.pca_bias_correction(amp, my_image.getImage().array)
+            bias_image = afwImage.ImageF(np.array(bias_model, dtype=np.float32))
+            bbox = self.amp_geom.serial_overscan
+            chisq = self.ccd_pcas.bbox_chisq(my_image.getImage(), bias_image,
+                                             bbox)
+            if (max_pca_reduced_chisq is not None and
+                chisq/bbox.area > max_pca_reduced_chisq):
+                raise MaskedCCDBiasImageException
+            my_image.getImage().array[:] -= bias_model
         else:
             my_image -= self.bias_image(amp, overscan, **kwargs)
 

python/lsst/eotest/sensor/ccd_bias_pca.py

@@ -490,6 +490,18 @@ def pca_bias_correction(self, amp, image_array):
 
         return bias_model
 
+    @staticmethod
+    def bbox_chisq(amp_image, bias_model, bbox):
+        """
+        Compute the chi-square for the pixels in the provided bounding box,
+        using the bias_model value as the variance.
+        """
+        data = amp_image.Factory(amp_image, bbox, deep=True)
+        model = bias_model.Factory(bias_model, bbox)
+        data -= model
+        chisq = np.sum(data.array*data.array/model.array)
+        return chisq
+
     def make_bias_frame(self, raw_file, outfile, residuals_file=None,
                         amps=None):
         """

python/lsst/eotest/sensor/fe55_psf.py

@@ -53,8 +53,14 @@ def cluster_moments(dn, pos):
     sum_y2 = np.sum(vy*vy*dn)
     mean_x = sum_x / sum_0
     mean_y = sum_y / sum_0
-    std_x = np.sqrt((sum_x2 / sum_0) - (mean_x * mean_x))
-    std_y = np.sqrt((sum_y2 / sum_0) - (mean_y * mean_y))
+    var_x = (sum_x2 / sum_0) - (mean_x * mean_x)
+    if var_x < 0:
+        raise RuntimeError
+    std_x = np.sqrt(var_x)
+    var_y = (sum_y2 / sum_0) - (mean_y * mean_y)
+    if var_y < 0:
+        raise RuntimeError
+    std_y = np.sqrt(var_y)
     return (mean_x, mean_y, std_x, std_y, sum_0)
 
 
@@ -183,14 +189,14 @@ def _bg_image(self, image, ccd, nx, ny):
         return bg.getImageF()
 
     def process_image(self, ccd, amp, sigma0=0.36, dn0=1590./5.,
-                      bg_reg=(10, 10), logger=None, seqnum=0):
+                      bg_reg=(10, 40), logger=None, seqnum=0):
         """
         Process a segment and accumulate the fit results for each
         charge cluster.  The dn0 and sigma0 parameters are the
         starting values used for each fit.
         """
         try:
-            image = ccd.bias_subtracted_image(amp)
+            image = ccd.bias_subtracted_image(amp, max_pca_reduced_chisq=1)
         except MaskedCCDBiasImageException:
             print("DM stack error encountered when generating bias image ")
             print("from inferred overscan region.")