Add mask multiplication into algorithm.

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

@@ -117,6 +117,7 @@ def __init__(self, algoDir):
 
         # True if at least one of images has a blended object
         self.blend_exists = False
+        self.mask_growth_iter = None
 
     def reset(self):
         """Reset the calculation for the new input images with the same
@@ -812,8 +813,37 @@ def _singleItr(self, I1, I2, model, tol=1e-3):
                 boundaryT = self.getBoundaryThickness()
                 # If the compensable image has no blended centroids
                 # this function will just create a single masked donut
-                I1.makeBlendedMask(self._inst, model, boundaryT, 1, blendPadding=4, compensated=True)
-                I2.makeBlendedMask(self._inst, model, boundaryT, 1, blendPadding=4, compensated=True)
+                I1.makeBlendedMask(self._inst, model, boundaryT, 1, compensated=True)#, blendPadding=1)
+                I2.makeBlendedMask(self._inst, model, boundaryT, 1, compensated=True)#, blendPadding=1)
+
+                # Create shifted mask from non-blended mask
+                if self.blend_exists:
+                    for compIm in [I1, I2]:
+                        compIm.makeMask(
+                            self._inst, model, boundaryT, 1
+                        )
+                        if self.mask_growth_iter is None:
+                            (
+                                dilatedMask,
+                                numPaddingIter,
+                            ) = compIm.autoDilateBlendMask(
+                                compIm.mask_pupil
+                            )
+                        else:
+                            numPaddingIter = self.mask_growth_iter
+
+                        finalMask, shiftedMask = compIm.createBlendedCoadd(
+                            compIm.mask_pupil,
+                            blendPadding=numPaddingIter,
+                            returnShiftedMask=True,
+                        )
+                        # Mask only blended areas in final stamp
+                        compIm.updateImage(
+                            compIm.getImg() * np.invert(
+                                np.array(shiftedMask, dtype=bool)
+                            )
+                        )
+
                 self._makeMasterMask(I1, I2, self.getPoissonSolverName())
 
                 # Load the offAxis correction coefficients

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

@@ -1586,11 +1586,11 @@ def createBlendedCoadd(self, maskArray, blendPadding=None, returnShiftedMask=Fal
         Raises
         ------
         AssertionError
-            blendPadding must be None or an integer >= 0.
+            blendPadding must be None or an integer > 0.
         """
 
         if blendPadding is not None:
-            assert blendPadding >= 0, "blendPadding must be None or an integer >= 0."
+            assert blendPadding > 0, "blendPadding must be None or an integer > 0."
 
         # Switch x,y order because x,y locations in the catalogs and the
         # LSST Science Pipeline objects are consistent with each other but
@@ -1660,14 +1660,18 @@ def calcNumPeaks(imageArray, blendMask):
             medianPixel = np.median(
                 maskedPixelVals[np.where(maskedPixelVals > maskedPixelHist[1][1])]
             )
+            binEdges = np.linspace(0.5 * medianPixel, 1.5*medianPixel, num=11)
+            binEdges[-1] = 2.5 * medianPixel
+            binEdges = np.append(binEdges, 3.0 * medianPixel)
             maskedPixelHist = np.histogram(
-                maskedPixelVals, bins=12, range=(0, 2.5 * medianPixel)
+                maskedPixelVals, bins=binEdges, # range=(0.5 * medianPixel, 2.5 * medianPixel)
             )
+            print(maskedPixelHist)
             # Find the highest bins
             binRanking = np.argsort(maskedPixelHist[0])
             # The highest bin will be the one around 0,
             # the second highest should be the main donut
-            maxNonZeroBin = binRanking[-2]
+            maxNonZeroBin = binRanking[-1]
             # Find the peaks
             histPeaks = argrelmax(maskedPixelHist[0])  # [0]
             # Count the number of peaks beyond the main donut peak
@@ -1679,6 +1683,11 @@ def calcNumPeaks(imageArray, blendMask):
 
         blendMask = self.createBlendedCoadd(copy(shiftedMask), blendPadding=None)
         numMaskedPeaks = calcNumPeaks(self.getImg(), blendMask)
+
+        # Return if no need to pad mask
+        if numMaskedPeaks == 0:
+            return blendMask, None
+
         paddingIter = 0
         while numMaskedPeaks > 0:
             paddingIter += 1