Merge pull request #11 from lofar-astron/master

Merge master

bdsf/gaul2srl.py

@@ -44,6 +44,7 @@ def __call__(self, img):
         sources = []
         dsources = []
         no_gaus_islands = []
+        no_gaus_islands_flag_values = []
         for iisl, isl in enumerate(img.islands):
             isl_sources = []
             isl_dsources = []
@@ -65,6 +66,10 @@ def __call__(self, img):
                 if not img.waveletimage:
                     dg = isl.dgaul[0]
                     no_gaus_islands.append((isl.island_id, dg.centre_pix[0], dg.centre_pix[1]))
+                    flag_values = []
+                    for fg in isl.fgaul:
+                        flag_values.append(fg.flag)
+                    no_gaus_islands_flag_values.append(flag_values)
                     # Put in the dummy Source as the source and use negative IDs
                     g_list = isl.dgaul
                     dsrc_index, dsource = self.process_single_gaussian(img, g_list, dsrc_index, code = 'S')
@@ -84,18 +89,29 @@ def __call__(self, img):
                 message += ':\n'
             else:
                 message += 's:\n'
-            for isl_id in no_gaus_islands:
-                message += '    Island #%i (x=%i, y=%i)\n' % isl_id
+            for isl_id, flag_list in zip(no_gaus_islands, no_gaus_islands_flag_values):
+                message += '    Island #%i (x=%i, y=%i): ' % isl_id
+                if len(flag_list) > 0:
+                    flags_str = '{}'.format(', '.join([str(f) for f in flag_list]))
+                    if len(flag_list) == 1:
+                        pl_str = ''
+                    else:
+                        pl_str = 's'
+                    message += 'fit with {0} Gaussian{1} with flag{1} = {2}\n'.format(len(flag_list), pl_str, flags_str)
+                else:
+                    message += '\n'
             if len(no_gaus_islands) == 1:
                 message += 'Please check this island. If it is a valid island and\n'
             else:
                 message += 'Please check these islands. If they are valid islands and\n'
             if img.opts.atrous_do:
                 message += 'should be fit, try adjusting the flagging options (use\n'\
-                           'show_fit with "ch0_flagged=True" to see the flagged Gaussians).'
+                           'show_fit with "ch0_flagged=True" to see the flagged Gaussians\n'\
+                           'and "help \'flagging_opts\'" to see the meaning of the flags).'
             else:
                 message += 'should be fit, try adjusting the flagging options (use\n'\
-                           'show_fit with "ch0_flagged=True" to see the flagged Gaussians)\n'\
+                           'show_fit with "ch0_flagged=True" to see the flagged Gaussians\n'\
+                           'and "help \'flagging_opts\'" to see the meaning of the flags)\n'\
                            'or enabling the wavelet module (with "atrous_do=True").'
             message += '\nTo include empty islands in output source catalogs, set\n'\
                         'incl_empty=True in the write_catalog task.'
@@ -203,7 +219,7 @@ def in_same_island(self, pair, img, g_list, isl, subim, subn, subm, delc):
         from . import functions as func
 
         def same_island_min(pair, g_list, subim, delc, tol=0.5):
-            """ If the minimum of the reconstructed fluxes along the line joining the peak positions
+            """ If the difference between the lower peak and the minimum of the reconstructed fluxes along the line joining the peak positions
                 is greater than thresh_isl times the rms_clip, they belong to different islands. """
 
             g1 = g_list[pair[0]]

bdsf/islands.py

@@ -57,6 +57,10 @@ def __call__(self, img):
             mylogger.userinfo(mylog, "Minimum number of pixels per island", '%i' %
                           minsize)
         img.minpix_isl = minsize
+        maxsize = opts.maxpix_isl
+        if maxsize is None:
+            maxsize = N.inf
+        img.maxpix_isl = maxsize
 
         if opts.detection_image != '':
             # Use a different image for island detection. The detection
@@ -85,12 +89,11 @@ def __call__(self, img):
                 raise RuntimeError("Detection image shape does not match that of input image.")
 
             # Run through islands and correct the image and rms, mean and max values
-            img.island_labels = det_img.island_labels
             corr_islands = []
             mean_map = img.mean_arr
             rms_map = img.rms_arr
             for i, isl in enumerate(det_img.islands):
-                islcp = isl.copy(img.pixel_beamarea(), image=ch0_map[isl.bbox], mean=mean_map[isl.bbox], rms=rms_map[isl.bbox])
+                islcp = isl.copy(img.pixel_beamarea(), image=ch0_map[tuple(isl.bbox)], mean=mean_map[tuple(isl.bbox)], rms=rms_map[tuple(isl.bbox)])
                 islcp.island_id = i
                 corr_islands.append(islcp)
             img.islands = corr_islands
@@ -189,7 +192,7 @@ def ndimage_alg(self, img, opts):
             isl_peak = nd.maximum(image[s], labels[s], idx)
             isl_maxposn = tuple(N.array(N.unravel_index(N.nanargmax(image[s]), image[s].shape))+\
                           N.array((s[0].start, s[1].start)))
-            if (isl_size >= img.minpix_isl) and (isl_peak - mean[isl_maxposn])/thresh_pix > rms[isl_maxposn]:
+            if (isl_size >= img.minpix_isl) and (isl_size <= img.maxpix_isl) and (isl_peak - mean[isl_maxposn])/thresh_pix > rms[isl_maxposn]:
                 isl = Island(image, mask, mean, rms, labels, s, idx, img.pixel_beamarea())
                 res.append(isl)
                 pyrank[tuple(isl.bbox)] += N.invert(isl.mask_active)*idx // idx
@@ -223,7 +226,7 @@ def coords_to_isl(self, img, opts):
                 aper_mask = N.where(labels.astype(bool) & ~mask)
             else:
                 aper_mask = N.where(labels.astype(bool))
-            if N.size(aper_mask) > img.minpix_isl:
+            if N.size(aper_mask) >= img.minpix_isl and N.size(aper_mask) <= img.maxpix_isl:
                 labels[aper_mask] = idx
                 s = [slice(max(0, isl_posn_pix[0] - isl_radius_pix - 1),
                      min(image.shape[0], isl_posn_pix[0] + isl_radius_pix + 1)),
@@ -341,6 +344,7 @@ def __init__(self, img, mask, mean, rms, labels, bbox, idx,
         self.total_fluxE = func.nanmean(bbox_rms_im[in_bbox_and_unmasked]) * N.sqrt(pixels_in_isl/beamarea) # Jy
         self.border = self.get_border()
         self.gaul = []
+        self.fgaul = []
         self.sources = []
         self.gresid_mean = 0.0
         self.gresid_rms = 0.0

bdsf/opts.py

@@ -223,10 +223,6 @@ class Opts(object):
                                  "significant extended emission in the image, "\
                                  "it is often necessary to force the use of a "\
                                  "constant rms map by setting rms_map = False.")
-
-    rmsmean_map_filename = List(None,doc = "Name of the fits file for the mean and rms maps. Should be a list [<mean_map.fits>,<rms_map.fits>]",
-                             group = 'advanced_opts')
-
     shapelet_do = Bool(False,
                              doc = "Decompose islands into shapelets\n"\
                                  "If True, then each island is decomposed using shapelets, "\
@@ -416,6 +412,12 @@ class Opts(object):
                                  "image header. It is set to 6 pixels for all "\
                                  "wavelet images.",
                              group = "advanced_opts")
+    maxpix_isl = Option(None, Int(),
+                             doc = "Maximum number of pixels with emission per island. "\
+                                 "None -> no limit\n"\
+                                 "This is an integer and is the maximum number of pixels "\
+                                 "in an island for the island to be included.",
+                             group = "advanced_opts")
     rms_value = Option(None, Float(),
                              doc = "Value of constant rms in "\
                                  "Jy/beam to use if rms_map = False. "\
@@ -548,8 +550,10 @@ class Opts(object):
                              doc = "Tolerance for grouping of Gaussians into sources: "\
                                  "larger values will result in larger sources\n"\
                                  "Sources are created by "\
-                                 "grouping nearby Gaussians as follows: (1) If the minimum "\
-                                 "value between two Gaussians in an island is more than "\
+                                 "grouping nearby Gaussians as follows: (1) If the "\
+                                 "difference between the minimum value between two "\
+                                 "Gaussians and the lower of the peak flux densities of "\
+                                 "the Gaussians in an island is less than "\
                                  "group_tol * thresh_isl * rms_clip, "\
                                  "and (2) if the centres are seperated by a distance less "\
                                  "than 0.5*group_tol of the sum of their fwhms along the "\
@@ -572,6 +576,12 @@ class Opts(object):
                                  "3-, or 4-D cube. The detection image and the main"\
                                  "image must have the same size and be registered.",
                              group = "advanced_opts")
+    rmsmean_map_filename = List(None,
+                             doc = "Filenames of FITS files to use as the mean and rms maps, "\
+                                 "given as a list [<mean_map.fits>, <rms_map.fits>]. If "\
+                                 "supplied, the internally generated mean and rms maps "\
+                                 "are not used.",
+                             group = 'advanced_opts')
     do_mc_errors = Bool(False,
                              doc = "Estimate uncertainties for 'M'-type sources using Monte "\
                                 "Carlo method\n"\

bdsf/wavelet_atrous.py

@@ -40,8 +40,9 @@
     N.fft.ifftn = pyfftw.interfaces.numpy_fft.ifftn
     scipy.signal.signaltools.fftn = pyfftw.interfaces.scipy_fftpack.fftn
     scipy.signal.signaltools.ifftn = pyfftw.interfaces.scipy_fftpack.ifftn
+    has_pyfftw = True
 except ImportError:
-    pass
+    has_pyfftw = False
 
 
 class Op_wavelet_atrous(Op):

doc/source/algorithms.rst

@@ -8,7 +8,7 @@ image, :math:`N_s`, as:
 .. math::
 
     N_s = (\text{No. pixels} > 5\sigma)/(<\text{pix/src}>),
- 
+
 where the average number of pixels per source, :math:`<pix/src>`, is given by:
 
 .. math::
@@ -52,8 +52,8 @@ Grouping of Gaussians into sources
 ----------------------------------
 Inside each island, groups of Gaussians are deemed to be a part of the same source if:
 
-    1. no pixel on the line joining the centers of any pair of Gaussians has a (Gaussian-reconstructed) value less than the island threshold, and 
-    2. the centers are separated by a distance less than half the sum of their FWHMs along the line joining them. 
+    1. the difference between the minimum value along the line joining the centers of any pair of Gaussians and the peak value of the lower Gaussian is less than the product of the island threshold and the island rms, and
+    2. the centers are separated by a distance less than half the sum of their FWHMs along the line joining them.
 
 Once the Gaussians that belong to a source are identified, fluxes for the grouped Gaussians are summed to obtain the total flux of the source. The uncertainty on the total flux is calculated by summing the uncertainties on the total fluxes of the individual Gaussians in quadrature. The source RA and Dec position is set to the source centroid determined from moment analysis (the position of the maximum of the source is also calculated). The total source size is also measured using moment analysis (see http://en.wikipedia.org/wiki/Image_moment for an overview of moment analysis).
 
@@ -89,4 +89,4 @@ No color correction is performed when averaging channels. However, as is shown b
 
 
 
-.. [#f1] Katgert, P., Oort, M. J. A., & Windhorst, R. A. 1988, A&A, 195, 21 
+.. [#f1] Katgert, P., Oort, M. J. A., & Windhorst, R. A. 1988, A&A, 195, 21

doc/source/process_image.rst

@@ -346,7 +346,9 @@ The advanced options are:
       :term:`ini_method` .... 'intensity': Method by which inital guess for fitting of Gaussians is chosen:
                                    'intensity' or 'curvature'
       :term:`kappa_clip` ........... 3.0 : Kappa for clipped mean and rms
-      :term:`minpix_isl` .......... None : Minimal number of pixels with emission per island.
+      :term:`maxpix_isl` .......... None : Maximum number of pixels with emission per island.
+                                   None -> no limit
+      :term:`minpix_isl` .......... None : Minimum number of pixels with emission per island.
                                    None -> calculate inside program
       :term:`ncores` .............. None : Number of cores to use during fitting, None => use
                                    all
@@ -355,6 +357,10 @@ The advanced options are:
       :term:`peak_maxsize` ........ 30.0 : If island size in beam area is more than this,
                                    attempt to fit peaks separately (if
                                    peak_fit=True). Min value is 30
+      :term:`rmsmean_map_filename`  None : Filenames of FITS files to use as the mean and rms maps,
+                                   given as a list [<mean_map.fits>, <rms_map.fits>]. If
+                                   supplied, the internally generated mean and rms maps
+                                   are not used
       :term:`rms_value` ........... None : Value of constant rms in Jy/beam to use if rms_map
                                    = False. None => calculate inside program
       :term:`spline_rank` ............ 3 : Rank of the interpolating function for rms/mean
@@ -495,11 +501,12 @@ The advanced options are:
         This parameter is a float (default is 1.0) that sets the tolerance for
         grouping of Gaussians into sources: larger values will result in larger
         sources. Sources are created by grouping nearby Gaussians as follows:
-        (1) If the minimum value between two Gaussians in an island is more than
-        ``group_tol * thresh_isl * rms_clip``\, and (2) if the centers are
-        separated by a distance less than ``0.5 * group_tol`` of the sum of their
-        FWHMs along the PA of the line joining them, they belong to the same
-        island.
+        (1) If the difference between the minimum value between two Gaussians
+        and the lower of the peak flux densities of the Gaussians in an island
+        is less than ``group_tol * thresh_isl * rms_clip``\, and (2) if the
+        centers are separated by a distance less than ``0.5 * group_tol`` of the
+        sum of their FWHMs along the PA of the line joining them, they belong to
+        the same island.
 
     ini_gausfit
         This parameter is a string (default is ``'default'``). These are three
@@ -537,6 +544,11 @@ The advanced options are:
         source pixels, less number of pixels in total, or significant
         non-Gaussianity of the underlying noise will all lead to non-convergence.
 
+    maxpix_isl
+        This parameter is an integer (default is ``None``) that sets the maximum
+        number of pixels in an island for the island to be included. If
+        ``None``, there is no limit.
+
     minpix_isl
         This parameter is an integer (default is ``None``) that sets the minimum
         number of pixels in an island for the island to be included. If
@@ -562,6 +574,12 @@ The advanced options are:
         is more than this value, attempt to fit peaks iteratively (if ``peak_fit
         = True``). The minimum value is 30.
 
+    rmsmean_map_filename
+        This parameter is a list (default is ``None``) that sets the filenames of
+        FITS files to use as the mean and rms maps, given as a list
+        [<mean_map.fits>, <rms_map.fits>]. If supplied, the internally generated
+        mean and rms maps are not used.
+
     rms_value
         This parameter is a float (default is ``None``) that sets the value of
         constant rms in Jy/beam to use if ``rms_map = False``. If ``None``, the

doc/source/write_catalog.rst

@@ -232,13 +232,13 @@ The information included in the Gaussian and source catalogs varies by format an
 
 * **E_Isl_Total_flux:** the 1-:math:`\sigma` error on the total flux density of the island in which the source is located, in Jy
 
-* **Isl_rms:** the average background rms value of the island, in Jy/beam
+* **Isl_rms:** the average background rms value of the island (derived from the rms map), in Jy/beam
 
-* **Isl_mean:** the averge background mean value of the island, in Jy/beam
+* **Isl_mean:** the averge background mean value of the island  (derived from the mean map), in Jy/beam
 
-* **Resid_Isl_rms:** the average residual background rms value of the island, in Jy/beam
+* **Resid_Isl_rms:** the average residual background rms value of the island (derived from the residual map, after subtraction of fitted Gaussians), in Jy/beam
 
-* **Resid_Isl_mean:** the averge residual background mean value of the island, in Jy/beam
+* **Resid_Isl_mean:** the averge residual background mean value of the island (derived from the residual map, after subtraction of fitted Gaussians), in Jy/beam
 
 * **S_Code:** a code that defines the source structure.
     * 'S' = a single-Gaussian source that is the only source in the island

pyproject.toml

@@ -0,0 +1,4 @@
+[build-system]
+requires = ["setuptools", 
+            "wheel", 
+            "numpy"]