island mask checks during island characterisation

AegeanTools/source_finder.py

@@ -123,25 +123,30 @@ def find_islands(im, bkg, rms, seed_clip=5.0, flood_clip=4.0, log=log):
         ymin, ymax = f[i][1].start, f[i][1].stop
         # obey seed clip constraint
         if np.any(snr[xmin:xmax, ymin:ymax] > seed_clip):
+
             data_box = copy.copy(
                 im[xmin:xmax, ymin:ymax]
             )  # copy so that we don't blank the master data
-            data_box[
-                np.where(snr[xmin:xmax, ymin:ymax] < flood_clip)
-            ] = np.nan  # blank pixels that are outside the outerclip
-            data_box[
-                np.where(l[xmin:xmax, ymin:ymax] != i + 1)
-            ] = np.nan  # blank out other summits
+
+
+            # make mask and blank out pixels with below the noise level or
+            # are pixels that are of another island in the FoV
+            island_mask = (snr[xmin:xmax, ymin:ymax] < flood_clip) | \
+                          (l[xmin:xmax, ymin:ymax] != i + 1)
+            data_box[island_mask] = np.nan 
+
             # check if there are any pixels left unmasked
             if not np.any(np.isfinite(data_box)):
                 # self.log.info("{1} Island {0} has no non-masked pixels"
                 #               .format(i,data.shape))
                 continue
+
             island = PixelIsland()
             island.calc_bounding_box(
                 np.array(np.nan_to_num(data_box), dtype=bool),
                 offsets=[xmin, ymin]
             )
+            island.set_mask(island_mask)
             islands.append(island)
 
     return islands
@@ -195,6 +200,10 @@ def estimate_parinfo_image(islands, im, rms, wcshelper,
         i_data = im[rmin:rmax, cmin:cmax]
         i_rms = rms[rmin:rmax, cmin:cmax]
 
+        island_mask = island.mask 
+        i_data[island_mask] = np.nan
+        i_rms[island_mask] = np.nan
+
         # the curvature needs a buffer of 1 pixel to correctly identify
         # the local min/max on the edge of the region.
         # We need a 1 pix buffer (if available)
@@ -237,15 +246,12 @@ def estimate_parinfo_image(islands, im, rms, wcshelper,
         ]
         del peaks, pmask, troughs, tmask, buffx, buffy
 
-        # apply the island mask
-        i_data[np.where(np.bitwise_not(island.mask))] = np.nan
-
-        isnegative = max(
-            i_data[np.where(np.isfinite(i_data) & island.mask)]) < 0
+        # have already applied the island mask at start of loop
+        isnegative = np.nanmax(i_data[np.isfinite(i_data)]) < 0
 
         # For small islands we can't do a 6 param fit
         # Don't count the NaN values as part of the island
-        non_nan_pix = len(i_data[np.where(np.isfinite(i_data))].ravel())
+        non_nan_pix = np.sum(np.isfinite(i_data))
         if 4 <= non_nan_pix <= 6:
             log.debug("FIXED2PSF")
             is_flag |= flags.FIXED2PSF