From ee3b81a855dc4e41e29528b4de0f3ead67f06b32 Mon Sep 17 00:00:00 2001
From: Arno Schoenmakers <schoenmakers@astron.nl>
Date: Wed, 12 Nov 2014 15:29:11 +0000
Subject: [PATCH] Task #6735: Merging fix in PyBDSM (#6871, r. 30408) into
 release brach

---
 CEP/PyBDSM/src/apps/pybdsm.in           |  4 +-
 CEP/PyBDSM/src/python/__init__.py       |  2 +-
 CEP/PyBDSM/src/python/_version.py       |  4 ++
 CEP/PyBDSM/src/python/collapse.py       |  8 ++--
 CEP/PyBDSM/src/python/functions.py      | 55 ++++++++++++++-----------
 CEP/PyBDSM/src/python/gaul2srl.py       |  2 +-
 CEP/PyBDSM/src/python/gausfit.py        | 14 +++----
 CEP/PyBDSM/src/python/image.py          |  6 +--
 CEP/PyBDSM/src/python/interface.py      | 32 +++++++-------
 CEP/PyBDSM/src/python/islands.py        | 16 +++----
 CEP/PyBDSM/src/python/multi_proc.py     | 20 +++++----
 CEP/PyBDSM/src/python/opts.py           |  8 ++--
 CEP/PyBDSM/src/python/output.py         | 50 +++++++++++-----------
 CEP/PyBDSM/src/python/plotresults.py    | 16 +++----
 CEP/PyBDSM/src/python/polarisation.py   |  4 +-
 CEP/PyBDSM/src/python/preprocess.py     |  8 +++-
 CEP/PyBDSM/src/python/psf_vary.py       | 22 +++++-----
 CEP/PyBDSM/src/python/pybdsm.py         | 16 +++----
 CEP/PyBDSM/src/python/readimage.py      | 20 ++++-----
 CEP/PyBDSM/src/python/rmsimage.py       | 14 +++----
 CEP/PyBDSM/src/python/shapefit.py       |  2 +-
 CEP/PyBDSM/src/python/shapelets.py      |  5 ++-
 CEP/PyBDSM/src/python/spectralindex.py  |  8 ++--
 CEP/PyBDSM/src/python/wavelet_atrous.py | 11 +++--
 24 files changed, 184 insertions(+), 163 deletions(-)

diff --git a/CEP/PyBDSM/src/apps/pybdsm.in b/CEP/PyBDSM/src/apps/pybdsm.in
index e592c6ffc76..ff70747efd1 100644
--- a/CEP/PyBDSM/src/apps/pybdsm.in
+++ b/CEP/PyBDSM/src/apps/pybdsm.in
@@ -14,7 +14,7 @@ PYTHONPATH=${PYTHONPATH:+${PYTHONPATH}:}@PYTHON_INSTALL_DIR@
 # And execute pybdsm.py. On a Mac, use pythonw if it exists instead of python
 # to avoid problems with the matplotlib OS X backend.
 if [ `uname` == "Darwin" ] && [ -f @PYTHON_EXECUTABLE@w ]; then
-    exec @PYTHON_EXECUTABLE@w @PYTHON_INSTALL_DIR@/lofar/bdsm/pybdsm.py
+    exec @PYTHON_EXECUTABLE@w -W ignore @PYTHON_INSTALL_DIR@/lofar/bdsm/pybdsm.py
 else
-    exec @PYTHON_EXECUTABLE@ @PYTHON_INSTALL_DIR@/lofar/bdsm/pybdsm.py
+    exec @PYTHON_EXECUTABLE@ -W ignore @PYTHON_INSTALL_DIR@/lofar/bdsm/pybdsm.py
 fi
diff --git a/CEP/PyBDSM/src/python/__init__.py b/CEP/PyBDSM/src/python/__init__.py
index 0d2eb1bb95b..083aacaac4d 100644
--- a/CEP/PyBDSM/src/python/__init__.py
+++ b/CEP/PyBDSM/src/python/__init__.py
@@ -232,7 +232,7 @@ def process_image(input, **kwargs):
     # If load_pars fails (returns None), assume that input is an image file. If it's not a
     # valid image file (but is an existing file), an error will be raised
     # by img.process() during reading of the file.
-    if img == None:
+    if img is None:
         if os.path.exists(input):
             img = Image({'filename': input})
         else:
diff --git a/CEP/PyBDSM/src/python/_version.py b/CEP/PyBDSM/src/python/_version.py
index 9597d520177..234a7b5a3d9 100644
--- a/CEP/PyBDSM/src/python/_version.py
+++ b/CEP/PyBDSM/src/python/_version.py
@@ -27,6 +27,10 @@ def changelog():
     PyBDSM Changelog.
     -----------------------------------------------------------------------
 
+    2014/11/07 - Fix to bug that caused a crash when both atrous_do = True
+        and output_all = True. Fixed a bug that caused a crash on machines
+        with only one core.
+
     2014/09/26 - Version 1.8.3
 
     2014/09/26 - Fix to bug that caused a crash when using the wavelet
diff --git a/CEP/PyBDSM/src/python/collapse.py b/CEP/PyBDSM/src/python/collapse.py
index 62d78d6df2f..96e007e2701 100644
--- a/CEP/PyBDSM/src/python/collapse.py
+++ b/CEP/PyBDSM/src/python/collapse.py
@@ -223,7 +223,7 @@ def avspc_direct(c_list, image, rmsarr, c_wts, wtarr=None):
     shape2 = image.shape[1:]
     ch0 = N.zeros(shape2, dtype=N.float32)
     sumwts = 0.0
-    if wtarr == None:
+    if wtarr is None:
       wtarr = N.zeros(len(c_list))
       for i, ch in enumerate(c_list):
         im = image[ch]
@@ -253,7 +253,7 @@ def avspc_blanks(c_list, image, rmsarr, c_wts, wtarr=None):
     shape2 = image.shape[1:]
     ch0 = N.zeros(shape2, dtype=N.float32)
     sumwtim = N.zeros(shape2, dtype=N.float32)
-    if wtarr == None:
+    if wtarr is None:
       wtarr = N.zeros(len(c_list))
       for i, ch in enumerate(c_list):
         im = image[ch]
@@ -283,7 +283,7 @@ def avspc_blanks(c_list, image, rmsarr, c_wts, wtarr=None):
 def init_freq_collapse(img, wtarr):
     # Place appropriate, post-collapse frequency info in img
     # Calculate weighted average frequency
-    if img.opts.frequency_sp != None:
+    if img.opts.frequency_sp is not None:
         c_list = img.opts.collapse_av
         if c_list == []: c_list = N.arange(img.image_arr.shape[1])
         freqs = img.opts.frequency_sp
@@ -304,7 +304,7 @@ def init_freq_collapse(img, wtarr):
         sumwts = 0.0
         sumfrq = 0.0
         spec_indx = img.wcs_obj.wcs.spec
-        if spec_indx == -1 and img.opts.frequency_sp == None:
+        if spec_indx == -1 and img.opts.frequency_sp is None:
             raise RuntimeError("Frequency information not found in header and frequencies "\
                          "not specified by user")
         else:
diff --git a/CEP/PyBDSM/src/python/functions.py b/CEP/PyBDSM/src/python/functions.py
index 67a036403e5..281f6028921 100755
--- a/CEP/PyBDSM/src/python/functions.py
+++ b/CEP/PyBDSM/src/python/functions.py
@@ -399,7 +399,7 @@ def moment(x,mask=None):
     """
     import numpy as N
 
-    if mask == None:
+    if mask is None:
         mask=N.zeros(x.shape, dtype=bool)
     m1=N.zeros(1)
     m2=N.zeros(x.ndim)
@@ -432,7 +432,7 @@ def fit_mask_1d(x, y, sig, mask, funct, do_err, order=0, p0 = None):
       if isinstance(sig, list): sig = N.array(sig)
       xfit=x[ind]; yfit=y[ind]; sigfit=sig[ind]
 
-      if p0 == None:
+      if p0 is None:
         if funct == poly:
            p0=N.array([0]*(order+1))
            p0[1]=(yfit[0]-yfit[-1])/(xfit[0]-xfit[-1])
@@ -461,7 +461,7 @@ def fit_mask_1d(x, y, sig, mask, funct, do_err, order=0, p0 = None):
         sys.stdout = original_stdout  # turn STDOUT back on
 
       if do_err:
-        if cov != None:
+        if cov is not None:
           if N.sum(sig != 1.) > 0:
             err = N.array([sqrt(abs(cov[i,i])) for i in range(len(p))])
           else:
@@ -594,17 +594,17 @@ def fit_gaus2d(data, p_ini, x, y, mask = None, err = None):
     import numpy as N
     import sys
 
-    if mask != None and mask.shape != data.shape:
+    if mask is not None and mask.shape != data.shape:
         print 'Data and mask array dont have the same shape, ignoring mask'
         mask = None
-    if err != None and err.shape != data.shape:
+    if err is not None and err.shape != data.shape:
         print 'Data and error array dont have the same shape, ignoring error'
         err = None
 
-    if mask == None: mask = N.zeros(data.shape, bool)
+    if mask is None: mask = N.zeros(data.shape, bool)
     g_ind = N.where(~N.ravel(mask))[0]
 
-    if err == None:
+    if err is None:
         errorfunction = lambda p: N.ravel(gaus_2d(p, x, y) - data)[g_ind]
     else:
         errorfunction = lambda p: N.ravel((gaus_2d(p, x, y) - data)/err)[g_ind]
@@ -780,7 +780,7 @@ def get_errors(img, p, stdav, bm_pix=None):
         errors = errors + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
       else:
         sq2 = sqrt(2.0)
-        if bm_pix == None:
+        if bm_pix is None:
             bm_pix = N.array([img.pixel_beam()[0]*fwsig, img.pixel_beam()[1]*fwsig, img.pixel_beam()[2]])
         dumr = sqrt(abs(size[0] * size[1] / (4.0 * bm_pix[0] * bm_pix[1])))
         dumrr1 = 1.0 + bm_pix[0] * bm_pix[1] / (size[0] * size[0])
@@ -889,13 +889,13 @@ def fit_mulgaus2d(image, gaus, x, y, mask = None, fitfix = None, err = None, adj
     import numpy as N
     import sys
 
-    if mask != None and mask.shape != image.shape:
+    if mask is not None and mask.shape != image.shape:
         print 'Data and mask array dont have the same shape, ignoring mask'
         mask = None
-    if err != None and err.shape != image.shape:
+    if err is not None and err.shape != image.shape:
         print 'Data and error array dont have the same shape, ignoring error'
         err = None
-    if mask == None: mask = N.zeros(image.shape, bool)
+    if mask is None: mask = N.zeros(image.shape, bool)
 
     g_ind = N.where(~N.ravel(mask))[0]
 
@@ -906,7 +906,7 @@ def fit_mulgaus2d(image, gaus, x, y, mask = None, fitfix = None, err = None, adj
         p_ini = p_ini + g2param(g, adj)
       p_ini = N.array(p_ini)
 
-      if fitfix == None: fitfix = [0]*ngaus
+      if fitfix is None: fitfix = [0]*ngaus
       ind = N.ones(6*ngaus)                                     # 1 => fit ; 0 => fix
       for i in range(ngaus):
         if fitfix[i] == 1: ind[i*6+1:i*6+6] = 0
@@ -915,7 +915,7 @@ def fit_mulgaus2d(image, gaus, x, y, mask = None, fitfix = None, err = None, adj
       ind = N.array(ind)
       p_tofit = p_ini[N.where(ind==1)[0]]
       p_tofix = p_ini[N.where(ind==0)[0]]
-      if err == None: err = N.ones(image.shape)
+      if err is None: err = N.ones(image.shape)
 
       errorfunction = lambda p, x, y, p_tofix, ind, image, err, g_ind: \
                      N.ravel((gaus_2d_itscomplicated(p, x, y, p_tofix, ind)-image)/err)[g_ind]
@@ -1006,7 +1006,7 @@ def get_maxima(im, mask, thr, shape, beam, im_pos=None):
     from copy import deepcopy as cp
     import numpy as N
 
-    if im_pos == None:
+    if im_pos is None:
         im_pos = im
     im1 = cp(im)
     ind = N.array(N.where(~mask)).transpose()
@@ -1073,7 +1073,7 @@ def read_image_from_file(filename, img, indir, quiet=False):
     from distutils.version import StrictVersion
 
     mylog = mylogger.logging.getLogger("PyBDSM."+img.log+"Readfile")
-    if indir == None or indir == './':
+    if indir is None or indir == './':
         prefix = ''
     else:
         prefix = indir + '/'
@@ -1226,7 +1226,7 @@ def read_image_from_file(filename, img, indir, quiet=False):
     # entries.
     for i in range(4):
         key_val_raw = hdr.get('CUNIT' + str(i+1))
-        if key_val_raw != None:
+        if key_val_raw is not None:
             if 'M/S' in key_val_raw or 'm/S' in key_val_raw or 'M/s' in key_val_raw:
                 hdr['CUNIT' + str(i+1)] = 'm/s'
             if 'HZ' in key_val_raw or 'hZ' in key_val_raw or 'hz' in key_val_raw:
@@ -1276,7 +1276,7 @@ def read_image_from_file(filename, img, indir, quiet=False):
     img._original_naxis = data_shape
     img._original_shape = (shape_out[2], shape_out[3])
     img._xy_hdr_shift = (0, 0)
-    if img.opts.trim_box != None:
+    if img.opts.trim_box is not None:
         img.trim_box = img.opts.trim_box
         xmin, xmax, ymin, ymax = img.trim_box
         if xmin < 0: xmin = 0
@@ -1395,6 +1395,11 @@ def write_image_to_file(use, filename, image, img, outdir=None,
         xmin = 0
         ymin = 0
 
+    if not hasattr(img, '_telescope'):
+        telescope = None
+    else:
+        telescope = img._telescope
+
     if filename == 'SAMP':
         import tempfile
         if not hasattr(img,'samp_client'):
@@ -1404,14 +1409,14 @@ def write_image_to_file(use, filename, image, img, outdir=None,
 
         # Broadcast image to SAMP Hub
         temp_im = make_fits_image(N.transpose(image), wcs_obj, img.beam,
-            img.frequency, img.equinox, img._telescope, xmin=xmin, ymin=ymin,
+            img.frequency, img.equinox, telescope, xmin=xmin, ymin=ymin,
             is_mask=is_mask)
         tfile = tempfile.NamedTemporaryFile(delete=False)
         temp_im.writeto(tfile.name, clobber=clobber)
         send_fits_image(img.samp_client, img.samp_key, 'PyBDSM image', tfile.name)
     else:
         # Write image to FITS file
-        if outdir == None:
+        if outdir is None:
             outdir = img.indir
         if not os.path.exists(outdir) and outdir != '':
             os.makedirs(outdir)
@@ -1426,7 +1431,7 @@ def write_image_to_file(use, filename, image, img, outdir=None,
             else:
                 return
         temp_im = make_fits_image(N.transpose(image), wcs_obj, img.beam,
-            img.frequency, img.equinox, img._telescope, xmin=xmin, ymin=ymin,
+            img.frequency, img.equinox, telescope, xmin=xmin, ymin=ymin,
             is_mask=is_mask, shape=(img.shape[1], img.shape[0], img.shape[2],
             img.shape[3]))
         if use == 'rap':
@@ -1729,7 +1734,7 @@ def open_isl(mask, index):
     labels, n_subisl = nd.label(open, connectivity)  # get label/rank image for open. label = 0 for masked pixels
     labels, mask = assign_leftovers(mask, open, n_subisl, labels)  # add the leftover pixels to some island
 
-    if labels != None:
+    if labels is not None:
         isl_pixs = [len(N.where(labels==i)[0]) for i in range(1,n_subisl+1)]
         isl_pixs = N.array(isl_pixs)/float(N.sum(isl_pixs))
     else:
@@ -1872,10 +1877,10 @@ def isl_tosplit(isl, opts):
 
                                 # take open 3 or 5
     open3, open5 = False, False
-    if n_subisl3 > 0 and isl_pixs3 != None:                                 # open 3 breaks up island
+    if n_subisl3 > 0 and isl_pixs3 is not None:                                 # open 3 breaks up island
       max_sub3 = N.max(isl_pixs3)
       if max_sub3 < frac_bigisl3 : open3 = True       # if biggest sub island isnt too big
-    if n_subisl5 > 0 and isl_pixs5 != None:                                 # open 5 breaks up island
+    if n_subisl5 > 0 and isl_pixs5 is not None:                                 # open 5 breaks up island
       max_sub5 = N.max(isl_pixs5)                     # if biggest subisl isnt too big OR smallest extra islands add upto 10 %
       if (max_sub5 < 0.75*max_sub3) or (N.sum(N.sort(isl_pixs5)[:len(isl_pixs5)-n_subisl3]) > size_extra5):
         open5 = True
@@ -1913,7 +1918,7 @@ def ch0_aperture_flux(img, posn_pix, aperture_pix):
     """
     import numpy as N
 
-    if aperture_pix == None:
+    if aperture_pix is None:
         return [0.0, 0.0]
 
     # Make ch0 and rms subimages
@@ -1968,7 +1973,7 @@ def make_src_mask(mask_size, posn_pix, aperture_pix):
     import numpy as N
 
     xsize, ysize = mask_size
-    if aperture_pix == None:
+    if aperture_pix is None:
         return N.zeros((xsize, ysize), dtype=N.int)
 
     # Make subimages
diff --git a/CEP/PyBDSM/src/python/gaul2srl.py b/CEP/PyBDSM/src/python/gaul2srl.py
index 7458be7d518..8f2dcf9d729 100644
--- a/CEP/PyBDSM/src/python/gaul2srl.py
+++ b/CEP/PyBDSM/src/python/gaul2srl.py
@@ -39,7 +39,7 @@ class Op_gaul2srl(Op):
         mylog = mylogger.logging.getLogger("PyBDSM."+img.log+"Gaul2Srl")
         mylogger.userinfo(mylog, 'Grouping Gaussians into sources')
         img.aperture = img.opts.aperture
-        if img.aperture != None and img.aperture <= 0.0:
+        if img.aperture is not None and img.aperture <= 0.0:
             mylog.warn('Specified aperture is <= 0. Skipping aperture fluxes.')
             img.aperture = None
 
diff --git a/CEP/PyBDSM/src/python/gausfit.py b/CEP/PyBDSM/src/python/gausfit.py
index d71d1463eda..3fc57d14bce 100644
--- a/CEP/PyBDSM/src/python/gausfit.py
+++ b/CEP/PyBDSM/src/python/gausfit.py
@@ -111,7 +111,7 @@ class Op_gausfit(Op):
                 # These dummy Gaussians all have an ID of -1. They do not
                 # appear in any of the source or island Gaussian lists except
                 # the island dgaul list.
-                if opts.src_ra_dec != None:
+                if opts.src_ra_dec is not None:
                     # Center the dummy Gaussian on the user-specified source position
                     posn_isl = (isl.shape[0]/2.0, isl.shape[1]/2.0)
                     posn_img = (isl.shape[0]/2.0 + isl.origin[0], isl.shape[1]/2.0 + isl.origin[1])
@@ -173,7 +173,7 @@ class Op_gausfit(Op):
         # Check if there are many Gaussians with deconvolved size of 0 in one
         # axis but not in the other. Don't bother to do this for wavelet images.
         fraction_1d = self.check_for_1d_gaussians(img)
-        if fraction_1d > 0.5 and img.beam != None and img.waveletimage == False:
+        if fraction_1d > 0.5 and img.beam is not None and img.waveletimage == False:
             mylog.warn('After deconvolution, more than 50% of Gaussians are '\
                            "1-D. Unless you're fitting an extended source, "\
                            "beam may be incorrect.")
@@ -189,7 +189,7 @@ class Op_gausfit(Op):
         """
         import functions as func
 
-        if opts == None:
+        if opts is None:
             opts = img.opts
         iter_ngmax  = 10
         maxsize = opts.splitisl_maxsize
@@ -261,7 +261,7 @@ class Op_gausfit(Op):
         import functions as func
         from const import fwsig
 
-        if ffimg == None:
+        if ffimg is None:
             fit_image = isl.image-isl.islmean
         else:
             fit_image = isl.image-isl.islmean-ffimg
@@ -288,12 +288,12 @@ class Op_gausfit(Op):
         gaul = []
         iter = 0
         ng1 = 0
-        if ini_gausfit == None:
+        if ini_gausfit is None:
             ini_gausfit = opts.ini_gausfit
 
         if ini_gausfit not in ['default', 'simple', 'nobeam']:
             ini_gausfit = 'default'
-        if ini_gausfit == 'simple' and ngmax == None:
+        if ini_gausfit == 'simple' and ngmax is None:
           ngmax = 25
         if ini_gausfit == 'default' or opts.fix_to_beam:
           gaul, ng1, ngmax = self.inigaus_fbdsm(isl, thr0, beam, img)
@@ -409,7 +409,7 @@ class Op_gausfit(Op):
         import functions as func
         sgaul = []; sfgaul = []
         gaul = []; fgaul = []
-        if opts == None:
+        if opts is None:
             opts = img.opts
         thresh_isl = opts.thresh_isl
         thresh_pix = opts.thresh_pix
diff --git a/CEP/PyBDSM/src/python/image.py b/CEP/PyBDSM/src/python/image.py
index 8a4ad2f322c..cb4e1778211 100644
--- a/CEP/PyBDSM/src/python/image.py
+++ b/CEP/PyBDSM/src/python/image.py
@@ -70,7 +70,7 @@ class Image(object):
         if name.endswith("_arr"):
             if self.do_cache:
                 map_data = func.retrieve_map(self, name)
-                if map_data != None:
+                if map_data is not None:
                     return map_data
                 else:
                     return object.__getattribute__(self, name)
@@ -135,11 +135,11 @@ class Image(object):
         """Load parameter values."""
         import interface
         import os
-        if loadfile == None or loadfile == '':
+        if loadfile is None or loadfile == '':
             loadfile = self.opts.filename + '.pybdsm.sav'
         if os.path.exists(loadfile):
             timg, err = interface.load_pars(loadfile)
-            if timg != None:
+            if timg is not None:
                 orig_filename = self.opts.filename
                 self.opts = timg.opts
                 self.opts.filename = orig_filename # reset filename to original
diff --git a/CEP/PyBDSM/src/python/interface.py b/CEP/PyBDSM/src/python/interface.py
index bb2fd8bda3b..dd4d1070edd 100644
--- a/CEP/PyBDSM/src/python/interface.py
+++ b/CEP/PyBDSM/src/python/interface.py
@@ -96,7 +96,7 @@ def get_op_chain(img):
                 'outlist',
                 'cleanup']
     prev_opts = img._prev_opts
-    if prev_opts == None:
+    if prev_opts is None:
         return img, default_chain
     new_opts = img.opts.to_dict()
 
@@ -370,7 +370,7 @@ def save_pars(img, savefile=None, quiet=False):
     import tc
     import sys
 
-    if savefile == None or savefile == '':
+    if savefile is None or savefile == '':
         savefile = img.opts.filename + '.pybdsm.sav'
 
     # convert opts to dictionary
@@ -397,7 +397,7 @@ def list_pars(img, opts_list=None, banner=None, use_groups=True):
     opts = img.opts.to_list()
 
     # Filter list
-    if opts_list != None:
+    if opts_list is not None:
         opts_temp = []
         for o in opts:
             if o[0] in opts_list:
@@ -471,7 +471,7 @@ def group_opts(opts):
     gp = []
     for i in range(len(opts)):
         grp = opts[i][1].group()
-        if grp != None and grp not in groups:
+        if grp is not None and grp not in groups:
             groups.append(opts[i][1].group())
 
     groups.sort()
@@ -528,7 +528,7 @@ def print_opts(grouped_opts_list, img, banner=None):
     nc = '\033[0m'    # normal text color
     ncb = '\033[1m'    # normal text color bold
 
-    if banner != None:
+    if banner is not None:
         print banner
     spcstr = ' ' * minwidth # spaces string for second or later lines
     infix = nc + ': ' + nc # infix character used to separate values from comments
@@ -626,9 +626,9 @@ def print_opts(grouped_opts_list, img, banner=None):
                     else:
                         if isinstance(val, float):
                             val = round_float(val)
-                        if k == 'beam_spectrum' and val != None:
+                        if k == 'beam_spectrum' and val is not None:
                             val = round_list_of_tuples(val)
-                        if k == 'frequency_sp' and val != None:
+                        if k == 'frequency_sp' and val is not None:
                             val = round_list(val)
                         valstr = v1 + str(val) + v2
                         width_par_val = max(minwidth, len(k) + len(str(val)) + 6)
@@ -807,7 +807,7 @@ def export_image(img, outfile=None, img_format='fits', pad_image = False,
         print '\033[91mERROR\033[0m: img_format must be "fits" or "casa"'
         return False
     filename = outfile
-    if filename == None or filename == '':
+    if filename is None or filename == '':
         filename = img.imagename + '_' + img_type + '.' + format
     if os.path.exists(filename) and clobber == False:
         print '\033[91mERROR\033[0m: File exists and clobber = False.'
@@ -1031,7 +1031,7 @@ def write_catalog(img, outfile=None, format='bbs', srcroot=None, catalog_type='g
         filename = output.write_fits_list(img, filename=filename,
                                              incl_chan=incl_chan, incl_empty=incl_empty,
                                              clobber=clobber, objtype=catalog_type)
-        if filename == None:
+        if filename is None:
             print '\033[91mERROR\033[0m: File exists and clobber = False.'
             return False
         else:
@@ -1042,7 +1042,7 @@ def write_catalog(img, outfile=None, format='bbs', srcroot=None, catalog_type='g
                                               incl_chan=incl_chan, incl_empty=incl_empty,
                                               sort_by='index', format = format,
                                               clobber=clobber, objtype=catalog_type)
-        if filename == None:
+        if filename is None:
             print '\033[91mERROR\033[0m: File exists and clobber = False.'
             return False
         else:
@@ -1057,7 +1057,7 @@ def write_catalog(img, outfile=None, format='bbs', srcroot=None, catalog_type='g
                                             patch=patch, correct_proj=correct_proj,
                                             sort_by='flux',
                                             clobber=clobber)
-        if filename == None:
+        if filename is None:
             print '\033[91mERROR\033[0m: File exists and clobber = False.'
             return False
         else:
@@ -1072,7 +1072,7 @@ def write_catalog(img, outfile=None, format='bbs', srcroot=None, catalog_type='g
                                             patch=patch,
                                             sort_by='flux',
                                             clobber=clobber)
-        if filename == None:
+        if filename is None:
             print '\033[91mERROR\033[0m: File exists and clobber = False.'
             return False
         else:
@@ -1082,7 +1082,7 @@ def write_catalog(img, outfile=None, format='bbs', srcroot=None, catalog_type='g
         filename = output.write_ds9_list(img, filename=filename,
                                             srcroot=srcroot, incl_empty=incl_empty,
                                             clobber=clobber, objtype=catalog_type)
-        if filename == None:
+        if filename is None:
             print '\033[91mERROR\033[0m: File exists and clobber = False.'
             return False
         else:
@@ -1094,7 +1094,7 @@ def write_catalog(img, outfile=None, format='bbs', srcroot=None, catalog_type='g
             return False
         filename = output.write_star(img, filename=filename,
                                         clobber=clobber)
-        if filename == None:
+        if filename is None:
             print '\033[91mERROR\033[0m: File exists and clobber = False.'
             return False
         else:
@@ -1106,7 +1106,7 @@ def write_catalog(img, outfile=None, format='bbs', srcroot=None, catalog_type='g
             return False
         filename = output.write_kvis_ann(img, filename=filename,
                                             clobber=clobber)
-        if filename == None:
+        if filename is None:
             print '\033[91mERROR\033[0m: File exists and clobber=False.'
             return False
         else:
@@ -1115,7 +1115,7 @@ def write_catalog(img, outfile=None, format='bbs', srcroot=None, catalog_type='g
     if format == 'casabox':
         filename = output.write_casa_gaul(img, filename=filename,
                                       incl_empty=incl_empty, clobber=clobber)
-        if filename == None:
+        if filename is None:
             print '\033[91mERROR\033[0m: File exists and clobber=False.'
         else:
             print '--> Wrote CASA clean box file ' + filename
diff --git a/CEP/PyBDSM/src/python/islands.py b/CEP/PyBDSM/src/python/islands.py
index f6742e70dc2..02fe77e24ef 100644
--- a/CEP/PyBDSM/src/python/islands.py
+++ b/CEP/PyBDSM/src/python/islands.py
@@ -49,7 +49,7 @@ class Op_islands(Op):
         opts = img.opts
 
         minsize = opts.minpix_isl
-        if minsize == None:
+        if minsize is None:
             minsize = int(img.pixel_beamarea()/3.0) # 1/3 of beam area in pixels
             if minsize < 6:
                 minsize = 6 # Need at least 6 pixels to obtain good fits
@@ -98,7 +98,7 @@ class Op_islands(Op):
             img.minpix_isl = det_img.minpix_isl
             mylogger.userinfo(mylog, "\nContinuing processing using primary image")
         else:
-            if opts.src_ra_dec != None:
+            if opts.src_ra_dec is not None:
                 mylogger.userinfo(mylog, "Constructing islands at user-supplied source locations")
                 img.islands = self.coords_to_isl(img, opts)
             else:
@@ -205,7 +205,7 @@ class Op_islands(Op):
 
         coords = opts.src_ra_dec # list of RA and Dec tuples
         isl_radius_pix = opts.src_radius_pix
-        if isl_radius_pix == None:
+        if isl_radius_pix is None:
             isl_radius_pix = img.beam2pix(img.beam)[0] # twice beam major axis radius at half max (= FWHM)
 
         res = []
@@ -318,10 +318,10 @@ class Island(object):
                 noise_mask[mask[bbox]] = True
                 isl_mask[mask[bbox]] = True
         else:
-            if origin == None:
+            if origin is None:
                 origin = [b.start for b in bbox]
             isl_mask = mask
-            if noise_mask == None:
+            if noise_mask is None:
                 noise_mask = mask
             data = img
             bbox_rms_im = rms
@@ -393,11 +393,11 @@ class Island(object):
     def copy(self, pixel_beamarea, image=None, mean=None, rms=None):
         mask = self.mask_active
         noise_mask = self.mask_noisy
-        if image == None:
+        if image is None:
             image = self.image
-        if mean == None:
+        if mean is None:
             mean = N.zeros(mask.shape, dtype=N.float32) + self.mean
-        if rms == None:
+        if rms is None:
             rms =  N.zeros(mask.shape, dtype=N.float32) + self.rms
 
         bbox = self.bbox
diff --git a/CEP/PyBDSM/src/python/multi_proc.py b/CEP/PyBDSM/src/python/multi_proc.py
index f4c348d6cd7..597e8fd7c59 100644
--- a/CEP/PyBDSM/src/python/multi_proc.py
+++ b/CEP/PyBDSM/src/python/multi_proc.py
@@ -9,13 +9,13 @@ Adapted from a module by Brian Refsdal at SAO, available at AstroPython
 
 """
 import numpy
-_multi=False
-_ncpus=1
+_multi = False
+_ncpus = 1
 
 try:
     # May raise ImportError
     import multiprocessing
-    _multi=True
+    _multi = True
 
     # May raise NotImplementedError
     _ncpus = min(multiprocessing.cpu_count(), 8)
@@ -55,7 +55,7 @@ def worker(f, ii, chunk, out_q, err_q, lock, bar, bar_state):
         vals.append(result)
 
         # update statusbar
-        if bar != None:
+        if bar is not None:
             if bar_state['started']:
                 bar.pos = bar_state['pos']
                 bar.spin_pos = bar_state['spin_pos']
@@ -146,16 +146,18 @@ def parallel_map(function, sequence, numcores=None, bar=None, weights=None):
 
     if not _multi or size == 1:
         results = map(function, sequence)
-        if bar != None:
+        if bar is not None:
             bar.stop()
         return results
 
 
-    # Set default number of cores to use. Leave one core free for pyplot.
+    # Set default number of cores to use. Try to leave one core free for pyplot.
     if numcores is None:
         numcores = _ncpus - 1
     if numcores > _ncpus - 1:
         numcores = _ncpus - 1
+    if numcores < 1:
+        numcores = 1
 
     # Returns a started SyncManager object which can be used for sharing
     # objects between processes. The returned manager object corresponds
@@ -170,7 +172,7 @@ def parallel_map(function, sequence, numcores=None, bar=None, weights=None):
     err_q = manager.Queue()
     lock = manager.Lock()
     bar_state = manager.dict()
-    if bar != None:
+    if bar is not None:
         bar_state['pos'] = bar.pos
         bar_state['spin_pos'] = bar.spin_pos
         bar_state['started'] = bar.started
@@ -181,7 +183,7 @@ def parallel_map(function, sequence, numcores=None, bar=None, weights=None):
         numcores = size
 
     # group sequence into numcores-worth of chunks
-    if weights == None or numcores == size:
+    if weights is None or numcores == size:
         # No grouping specified (or there are as many cores as
         # processes), so divide into equal chunks
         sequence = numpy.array_split(sequence, numcores)
@@ -209,7 +211,7 @@ def parallel_map(function, sequence, numcores=None, bar=None, weights=None):
 
     try:
         results = run_tasks(procs, err_q, out_q, len(sequence))
-        if bar != None:
+        if bar is not None:
             if bar.started:
                 bar.stop()
         return results
diff --git a/CEP/PyBDSM/src/python/opts.py b/CEP/PyBDSM/src/python/opts.py
index 823cbdd0871..8cac0a16312 100644
--- a/CEP/PyBDSM/src/python/opts.py
+++ b/CEP/PyBDSM/src/python/opts.py
@@ -1378,7 +1378,7 @@ class Opts(object):
                 # and then try to parse it
                 if hasattr(self, k):
                     if isinstance(self.__getattribute__(k), bool):
-                        if isinstance(v, bool) or v == None:
+                        if isinstance(v, bool) or v is None:
                             # just enter the bool into the parameter
                             pass
                         elif isinstance(v, basestring):
@@ -1404,7 +1404,7 @@ class Opts(object):
         a string of a single opt name.
 
         If None, set all opts to default values."""
-        if opt_names == None:
+        if opt_names is None:
             TCInit(self)
         else:
             if isinstance(opt_names, str):
@@ -1437,7 +1437,7 @@ class Opts(object):
         opts_list = []
         for k, v in self.__class__.__dict__.iteritems():
             if isinstance(v, tc.TC):
-                if group != None:
+                if group is not None:
                     if v.group() == group:
                         opts_list.append((k, v))
                 else:
@@ -1464,7 +1464,7 @@ class Opts(object):
         opts_list = []
         for k, v in self.__class__.__dict__.iteritems():
             if isinstance(v, tc.TC):
-                if group != None:
+                if group is not None:
                     if v.group() == group:
                         opts_list.append(k)
                 else:
diff --git a/CEP/PyBDSM/src/python/output.py b/CEP/PyBDSM/src/python/output.py
index c636a1fb325..8aaf0d96915 100644
--- a/CEP/PyBDSM/src/python/output.py
+++ b/CEP/PyBDSM/src/python/output.py
@@ -227,7 +227,7 @@ def write_bbs_gaul(img, filename=None, srcroot=None, patch=None,
     outstr_list = make_bbs_str(img, outl, outn, patl, incl_empty=incl_empty,
                                correct_proj=correct_proj)
 
-    if filename == None:
+    if filename is None:
         filename = img.imagename + '.sky_in'
     if os.path.exists(filename) and clobber == False:
         return None
@@ -260,7 +260,7 @@ def write_lsm_gaul(img, filename=None, srcroot=None, patch=None,
                                                root=srcroot, sort_by=sort_by)
     outstr_list = make_lsm_str(img, outl, outn, incl_empty=incl_empty)
 
-    if filename == None:
+    if filename is None:
         filename = img.imagename + '.lsm'
     if os.path.exists(filename) and clobber == False:
         return None
@@ -300,7 +300,7 @@ def write_ds9_list(img, filename=None, srcroot=None, deconvolve=False,
                             str(dsrc.source_id))
         outn = [outn]
     outstr_list = make_ds9_str(img, outl, outn, deconvolve=deconvolve, objtype=objtype, incl_empty=incl_empty)
-    if filename == None:
+    if filename is None:
         filename = img.imagename + '.' + objtype + '.reg'
     if os.path.exists(filename) and clobber == False:
         return None
@@ -328,7 +328,7 @@ def write_ascii_list(img, filename=None, sort_by='indx', format = 'ascii',
             outl[0] += img.dsources
     outstr_list = make_ascii_str(img, outl, objtype=objtype,
                                  incl_empty=incl_empty, format=format)
-    if filename == None:
+    if filename is None:
         if objtype == 'gaul':
             filename = img.imagename + '.gaul'
         elif objtype == 'srl':
@@ -351,7 +351,7 @@ def write_casa_gaul(img, filename=None,  incl_empty=False, clobber=False):
     mylog = mylogger.logging.getLogger("PyBDSM."+img.log+"Output")
     outl, outn, patl = list_and_sort_gaussians(img, patch=None)
     outstr_list = make_casa_str(img, outl)
-    if filename == None:
+    if filename is None:
         filename = img.imagename + '.box'
     if os.path.exists(filename) and clobber == False:
         return None
@@ -401,7 +401,7 @@ def write_fits_list(img, filename=None, sort_by='index', objtype='gaul',
                 nmax = isl.shapelet_nmax
         nmax += 1
 
-    if img.opts.aperture != None:
+    if img.opts.aperture is not None:
         incl_aper = True
     else:
         incl_aper = False
@@ -442,7 +442,7 @@ def write_fits_list(img, filename=None, sort_by='index', objtype='gaul',
         tbhdu.header['INIMAGE'] = (img.filename, 'Filename of image')
         tbhdu.header['FREQ0'] = (float(freq), 'Reference frequency')
         tbhdu.header['EQUINOX'] = (img.equinox, 'Equinox')
-    if filename == None:
+    if filename is None:
         filename = img.imagename + '.' + objtype + '.fits'
     if os.path.exists(filename) and clobber == False:
         return None
@@ -457,7 +457,7 @@ def write_kvis_ann(img, filename=None, sort_by='indx',
     import os
 
     mylog = mylogger.logging.getLogger("PyBDSM."+img.log+"Output")
-    if filename == None:
+    if filename is None:
         filename = img.imagename + '.kvis.ann'
     if os.path.exists(filename) and clobber == False:
         return None
@@ -491,7 +491,7 @@ def write_star(img, filename=None, sort_by='indx',
     import os
 
     mylog = mylogger.logging.getLogger("PyBDSM."+img.log+"Output")
-    if filename == None:
+    if filename is None:
         filename = img.imagename + '.star'
     if os.path.exists(filename) and clobber == False:
         return None
@@ -539,7 +539,7 @@ def make_bbs_str(img, glist, gnames, patchnames, objtype='gaul',
                                "MajorAxis, MinorAxis, Orientation, "\
                                "ReferenceFrequency='"+freq+"', "\
                                "SpectralIndex='[]'\n\n")
-    elif patchnames[0] == None:
+    elif patchnames[0] is None:
         outstr_list.append("format = Name, Type, Ra, Dec, I, Q, U, V, "\
                                "MajorAxis, MinorAxis, Orientation, "\
                                "ReferenceFrequency='"+freq+"', "\
@@ -552,11 +552,11 @@ def make_bbs_str(img, glist, gnames, patchnames, objtype='gaul',
     if objtype == 'shap':
         patchname_last = ''
         for pindx, patch_name in enumerate(patchnames): # loop over patches
-          if patch_name != None and patch_name != patchname_last:
+          if patch_name is not None and patch_name != patchname_last:
               outstr_list.append(', , ' + patch_name + ', 00:00:00, +00.00.00\n')
               patchname_last = patch_name
               names_in_patch = gnames[pindx]
-              if patch_name == None:
+              if patch_name is None:
                   outstr_list.append(src_name + sep + stype + sep + sra + sep +
                                      sdec + sep + total + sep + Q_flux + sep +
                                      U_flux + sep + V_flux + sep +
@@ -571,7 +571,7 @@ def make_bbs_str(img, glist, gnames, patchnames, objtype='gaul',
     else:
         patchname_last = ''
         for pindx, patch_name in enumerate(patchnames): # loop over patches
-          if patch_name != None and patch_name != patchname_last:
+          if patch_name is not None and patch_name != patchname_last:
               outstr_list.append(', , ' + patch_name + ', 00:00:00, +00.00.00\n')
               patchname_last = patch_name
           gaussians_in_patch = glist[pindx]
@@ -603,7 +603,7 @@ def make_bbs_str(img, glist, gnames, patchnames, objtype='gaul',
                   deconvstr = deconv1 + ', ' + deconv2 + ', ' + deconv3
                   specin = '-0.8'
                   if 'spectralindex' in img.completed_Ops:
-                      if g.spec_indx != None and N.isfinite(g.spec_indx):
+                      if g.spec_indx is not None and N.isfinite(g.spec_indx):
                           specin = str("%.3e" % (g.spec_indx))
                   sep = ', '
                   if img.opts.polarisation_do:
@@ -614,7 +614,7 @@ def make_bbs_str(img, glist, gnames, patchnames, objtype='gaul',
                       Q_flux = '0.0'
                       U_flux = '0.0'
                       V_flux = '0.0'
-                  if patch_name == None:
+                  if patch_name is None:
                       outstr_list.append(src_name + sep + stype + sep + sra + sep +
                                          sdec + sep + total + sep + Q_flux + sep +
                                          U_flux + sep + V_flux + sep +
@@ -717,7 +717,7 @@ def make_lsm_str(img, glist, gnames, incl_empty=False):
             deconvstr = deconv1 + ' ' + deconv2 + ' ' + deconv3
             specin = '-0.8'
             if 'spectralindex' in img.completed_Ops:
-                if g.spec_indx != None and N.isfinite(g.spec_indx):
+                if g.spec_indx is not None and N.isfinite(g.spec_indx):
                     specin = str("%.3e" % (g.spec_indx))
             sep = ' '
             if img.opts.polarisation_do:
@@ -740,7 +740,7 @@ def make_ds9_str(img, glist, gnames, deconvolve=False, objtype='gaul', incl_empt
     """Makes a list of string entries for a ds9 region file."""
     outstr_list = []
     freq = "%.5e" % img.frequency
-    if img.equinox == None:
+    if img.equinox is None:
         equinox = 'fk5'
     else:
         if int(img.equinox) == 2000:
@@ -806,7 +806,7 @@ def make_ascii_str(img, glist, objtype='gaul', format='ascii', incl_empty=False)
     outstr_list.append('# Reference frequency of the detection ("ch0") image: %s Hz\n' % freq)
     outstr_list.append('# Equinox : %s \n\n' % img.equinox)
     val_list = []
-    if img.opts.aperture != None:
+    if img.opts.aperture is not None:
         incl_aper = True
     else:
         incl_aper = False
@@ -820,7 +820,7 @@ def make_ascii_str(img, glist, objtype='gaul', format='ascii', incl_empty=False)
                                                               incl_aper=incl_aper,
                                                               incl_empty=incl_empty,
                                                               nchan=img.nchan)
-        if cvals != None:
+        if cvals is not None:
             cformats[-1] += "\n"
             if format == 'ascii':
                 if i == 0:
@@ -837,7 +837,7 @@ def make_fits_list(img, glist, objtype='gaul', nmax=30, incl_empty=False):
     import functions as func
 
     out_list = []
-    if img.opts.aperture != None:
+    if img.opts.aperture is not None:
         incl_aper = True
     else:
         incl_aper = False
@@ -849,7 +849,7 @@ def make_fits_list(img, glist, objtype='gaul', nmax=30, incl_empty=False):
                                                       incl_aper=incl_aper,
                                                       incl_empty=incl_empty,
                                                       nmax=nmax, nchan=img.nchan)
-        if cvals != None:
+        if cvals is not None:
             out_list.append(cvals)
     out_list = func.trans_gaul(out_list)
     return out_list
@@ -948,7 +948,7 @@ def list_and_sort_gaussians(img, patch=None, root=None,
     import functions as func
 
     # Define lists
-    if root == None:
+    if root is None:
         root = img.parentname
     gauslist = []
     gausname = []
@@ -1047,7 +1047,7 @@ def list_and_sort_gaussians(img, patch=None, root=None,
                 patchindx.append(p)
 
     # Sort
-    if patch == 'single' or patch == None:
+    if patch == 'single' or patch is None:
         outlist = [list(gauslist)]
         outlist_sorted = [list(gauslist)]
         outnames = [list(gausname)]
@@ -1180,7 +1180,7 @@ def make_output_columns(obj, fits=False, objtype='gaul', incl_spin=False,
                     val = obj.__getattribute__(name)
                     colname = obj.__class__.__dict__[name]._colname
                     units = obj.__class__.__dict__[name]._units
-                    if units == None:
+                    if units is None:
                         units = ' '
                     if isinstance(val, list):
                         # This is a list, so handle it differently. We assume the next
@@ -1190,7 +1190,7 @@ def make_output_columns(obj, fits=False, objtype='gaul', incl_spin=False,
                         val_next = obj.__getattribute__(next_name)
                         colname_next = obj.__class__.__dict__[next_name]._colname
                         units_next = obj.__class__.__dict__[next_name]._units
-                        if units_next == None:
+                        if units_next is None:
                             units_next = ' '
                         for i in range(len(val)):
                             cvals.append(val[i])
diff --git a/CEP/PyBDSM/src/python/plotresults.py b/CEP/PyBDSM/src/python/plotresults.py
index 7e5172bdff3..19a8aab73be 100644
--- a/CEP/PyBDSM/src/python/plotresults.py
+++ b/CEP/PyBDSM/src/python/plotresults.py
@@ -151,7 +151,7 @@ def plotresults(img, ch0_image=True, rms_image=True, mean_image=True,
         else:
             src_list = img.sources
             sed_src = get_src(src_list, 0)
-            if sed_src == None:
+            if sed_src is None:
                 print 'No sources found. Skipping source SED plots.'
             else:
                 images.append('seds')
@@ -201,7 +201,7 @@ def plotresults(img, ch0_image=True, rms_image=True, mean_image=True,
 
     im_mean = img.clipped_mean
     im_rms = img.clipped_rms
-    if img.resid_gaus == None:
+    if img.resid_gaus is None:
         low = 1.1*abs(img.min_value)
     else:
         low = N.max([1.1*abs(img.min_value),1.1*abs(N.nanmin(img.resid_gaus))])
@@ -411,10 +411,10 @@ def on_pick(event):
                 str(round(pflux,4)) + ' Jy/beam'
 
         # Transmit src_id, gaus_id, and coordinates to SAMP Hub (if we are connected)
-        if do_broadcast and samp_key != None:
-            if samp_gaul_table_url != None:
+        if do_broadcast and samp_key is not None:
+            if samp_gaul_table_url is not None:
                 func.send_highlight_row(samp_client, samp_key, samp_gaul_table_url, gaus_id)
-            if samp_srl_table_url != None:
+            if samp_srl_table_url is not None:
                 func.send_highlight_row(samp_client, samp_key, samp_srl_table_url, src_id)
             func.send_coords(samp_client, samp_key, g.centre_sky)
 
@@ -537,7 +537,7 @@ def on_press(event):
                     return
         ax_indx = images.index('seds')
         sed_src = get_src(src_list, srcid)
-        if sed_src == None:
+        if sed_src is None:
             print 'Source not found!'
             return
         srcid_cur = srcid
@@ -572,7 +572,7 @@ def on_press(event):
         num_pix_unmasked = float(N.size(N.where(mask == False), 1))
         mean_rms = N.nansum(img_rms[xmin:xmax, ymin:ymax])/num_pix_unmasked
         mean_map_flux = N.nansum(img_mean[xmin:xmax, ymin:ymax])/pixels_per_beam
-        if img_gaus_mod == None:
+        if img_gaus_mod is None:
             gaus_mod_flux = 0.0
         else:
             gaus_mod_flux = N.nansum(img_gaus_mod[xmin:xmax, ymin:ymax])/pixels_per_beam
@@ -583,7 +583,7 @@ def on_press(event):
             % (mean_map_flux,)
         print '  Gaussian model flux density ........... : %f Jy'\
             % (gaus_mod_flux,)
-        if img_shap_mod != None:
+        if img_shap_mod is not None:
             shap_mod_flux = N.nansum(img_shap_mod[xmin:xmax, ymin:ymax])/pixels_per_beam
             print '  Shapelet model flux density ........... : %f Jy'\
                 % (shap_mod_flux,)
diff --git a/CEP/PyBDSM/src/python/polarisation.py b/CEP/PyBDSM/src/python/polarisation.py
index e1b816aaf05..5130defd1ab 100644
--- a/CEP/PyBDSM/src/python/polarisation.py
+++ b/CEP/PyBDSM/src/python/polarisation.py
@@ -527,9 +527,9 @@ class Op_polarisation(Op):
                  Op_gausfit(), Op_gaul2srl(), Op_make_residimage()]
 
         opts = img.opts.to_dict()
-        if img.opts.pi_thresh_isl != None:
+        if img.opts.pi_thresh_isl is not None:
             opts['thresh_isl'] = img.opts.pi_thresh_isl
-        if img.opts.pi_thresh_pix != None:
+        if img.opts.pi_thresh_pix is not None:
             opts['thresh_pix'] = img.opts.pi_thresh_pix
         opts['thresh'] = 'hard'
         opts['polarisation_do'] = False
diff --git a/CEP/PyBDSM/src/python/preprocess.py b/CEP/PyBDSM/src/python/preprocess.py
index b851a96b848..4b65f7b173e 100644
--- a/CEP/PyBDSM/src/python/preprocess.py
+++ b/CEP/PyBDSM/src/python/preprocess.py
@@ -108,7 +108,7 @@ class Op_preprocess(Op):
 
         ### max/min pixel value & coordinates
         shape = image.shape[0:2]
-        if mask != None:
+        if mask is not None:
             img.blankpix = N.sum(mask)
         if img.blankpix == 0:
           max_idx = image.argmax()
@@ -143,7 +143,11 @@ class Op_preprocess(Op):
         ### if image seems confused, then take background mean as zero instead
         alpha_sourcecounts = 2.5  # approx diff src count slope. 2.2?
         if opts.bmpersrc_th is None:
-          n = (image >= 5.*crms).sum()
+          if mask is not None:
+              unmasked = N.where(~img.mask_arr)
+              n = (image[unmasked] >= 5.*crms).sum()
+          else:
+              n = (image >= 5.*crms).sum()
           if n <= 0:
             n = 1
             mylog.info('No pixels in image > 5-sigma.')
diff --git a/CEP/PyBDSM/src/python/psf_vary.py b/CEP/PyBDSM/src/python/psf_vary.py
index aefcf730c56..f35c4731d5b 100644
--- a/CEP/PyBDSM/src/python/psf_vary.py
+++ b/CEP/PyBDSM/src/python/psf_vary.py
@@ -73,7 +73,7 @@ class Op_psf_vary(Op):
             else:
                 snrcut = opts.psf_snrcut
             img.psf_snrcut = snrcut
-            if opts.psf_high_snr != None:
+            if opts.psf_high_snr is not None:
                 if opts.psf_high_snr < 10.0:
                     mylogger.userinfo(mylog, "Value of psf_high_snr too low; increasing to 10")
                     high_snrcut = 10.0
@@ -258,7 +258,7 @@ class Op_psf_vary(Op):
                     bar.stop()
 
                     # Interpolate Gaussian parameters
-                    if img.aperture == None:
+                    if img.aperture is None:
                         psf_maps = [psf_maj, psf_min, psf_pa, psfratio]
                     else:
                         psf_maps = [psf_maj, psf_min, psf_pa, psfratio, psfratio_aper]
@@ -271,15 +271,15 @@ class Op_psf_vary(Op):
                         psf_maps, itertools.repeat(psfcoords),
                         itertools.repeat(image.shape)), numcores=opts.ncores,
                         bar=bar)
-                    if img.aperture == None:
+                    if img.aperture is None:
                         psf_maj_int, psf_min_int, psf_pa_int, psf_ratio_int = map_list
                     else:
                         psf_maj_int, psf_min_int, psf_pa_int, psf_ratio_int, psf_ratio_aper_int = map_list
 
                     # Smooth if desired
-                    if img.opts.psf_smooth != None:
+                    if img.opts.psf_smooth is not None:
                         sm_scale = img.opts.psf_smooth / img.pix2beam([1.0, 1.0, 0.0])[0] / 3600.0 # pixels
-                        if img.opts.aperture == None:
+                        if img.opts.aperture is None:
                             psf_maps = [psf_maj_int, psf_min_int, psf_pa_int, psf_ratio_int]
                         else:
                             psf_maps = [psf_maj_int, psf_min_int, psf_pa_int, psf_ratio_int, psf_ratio_aper_int]
@@ -291,7 +291,7 @@ class Op_psf_vary(Op):
                             itertools.izip(itertools.repeat(self.blur_image),
                             psf_maps, itertools.repeat(sm_scale)), numcores=opts.ncores,
                             bar=bar)
-                        if img.aperture == None:
+                        if img.aperture is None:
                             psf_maj_int, psf_min_int, psf_pa_int, psf_ratio_int = map_list
                         else:
                             psf_maj_int, psf_min_int, psf_pa_int, psf_ratio_int, psf_ratio_aper_int = map_list
@@ -301,7 +301,7 @@ class Op_psf_vary(Op):
                     psf_min_int = N.array(psf_min_int)
                     psf_pa_int = N.array(psf_pa_int)
                     psf_ratio_int = N.array(psf_ratio_int)
-                    if img.aperture == None:
+                    if img.aperture is None:
                         psf_ratio_aper_int = N.zeros(psf_maj_int.shape, dtype=N.float32)
                     else:
                         psf_ratio_aper_int = N.array(psf_ratio_aper_int, dtype=N.float32)
@@ -511,7 +511,7 @@ class Op_psf_vary(Op):
         f_s = f_sclip[0]*f_sclip[1]
 
         # Add bright sources
-        if bright_snr_cut != None:
+        if bright_snr_cut is not None:
             if bright_snr_cut < 20.0:
                 bright_snr_cut = 20.0
             bright_srcs = N.where(snr >= bright_snr_cut)
@@ -890,12 +890,12 @@ class Op_psf_vary(Op):
             src_wts_aper = []
             for gt in tile_gauls:
                 src = img.sources[gt[0]]
-                if img.aperture != None:
+                if img.aperture is not None:
                     src_ratio_aper.append(src.peak_flux_max / src.aperture_flux)
                     src_wts_aper.append(src.total_flux / src.aperture_fluxE)
                 src_ratio.append(src.peak_flux_max / src.total_flux)
                 src_wts.append(src.total_flux / src.total_fluxE)
-            if img.aperture != None:
+            if img.aperture is not None:
                 psfratio_aper.append(sum(N.asarray(src_ratio_aper)*src_wts_aper)/sum(src_wts_aper))
             else:
                 psfratio_aper.append(0.0)
@@ -1048,7 +1048,7 @@ class Op_psf_vary(Op):
         from scipy.ndimage import gaussian_filter
 
         sx = n
-        if ny != None:
+        if ny is not None:
             sy = ny
         else:
             sy = n
diff --git a/CEP/PyBDSM/src/python/pybdsm.py b/CEP/PyBDSM/src/python/pybdsm.py
index 87d18c4ac19..a1782adc078 100644
--- a/CEP/PyBDSM/src/python/pybdsm.py
+++ b/CEP/PyBDSM/src/python/pybdsm.py
@@ -40,7 +40,7 @@ def inp(cur_cmd=None):
     success = _set_pars_from_prompt()
     if not success:
         return
-    if cur_cmd != None:
+    if cur_cmd is not None:
         if not hasattr(cur_cmd, 'arg_list'):
             print '\033[31;1mERROR\033[0m: not a valid task'
             return
@@ -65,7 +65,7 @@ def go(cur_cmd=None):
     success = _set_pars_from_prompt()
     if not success:
         return
-    if cur_cmd == None:
+    if cur_cmd is None:
         if not hasattr(_img, '_current_cmd'):
             print '\033[31;1mERROR\033[0m: no task is set'
             return
@@ -84,7 +84,7 @@ def default(cur_cmd=None):
     given, the parameters of the current task are reset.
     """
     global _img
-    if cur_cmd == None:
+    if cur_cmd is None:
         if not hasattr(_img, '_current_cmd'):
             print '\033[31;1mERROR\033[0m: no task is set'
             return
@@ -133,7 +133,7 @@ def tget(filename=None):
     if hasattr(filename, 'arg_list'):
         filename = None
 
-    if filename == None or filename == '':
+    if filename is None or filename == '':
         if os.path.isfile('pybdsm.last'):
             filename = 'pybdsm.last'
         else:
@@ -184,7 +184,7 @@ def tput(filename=None, quiet=False):
     success = _set_pars_from_prompt()
     if not success:
         return
-    if filename == None or filename == '':
+    if filename is None or filename == '':
         filename = 'pybdsm.last'
 
     # convert opts to dictionary
@@ -248,7 +248,7 @@ def _replace_vals_in_namespace(opt_names=None):
     global _img
     f = sys._getframe(len(inspect.stack())-1)
     f_dict = f.f_locals
-    if opt_names == None:
+    if opt_names is None:
         opt_names = _img.opts.get_names()
     if isinstance(opt_names, str):
         opt_names = [opt_names]
@@ -520,7 +520,7 @@ class bdsmDocHelper(pydoc.Helper):
                 else:
                     valstr = str(default_val)
                 default_val_text = 'Default value: ' + valstr
-                if opt.group() != None and opt.group() != 'hidden':
+                if opt.group() is not None and opt.group() != 'hidden':
                     group_text = '\nBelongs to group: ' + opt.group()
                 else:
                     group_text = ''
@@ -666,7 +666,7 @@ from lofar.bdsm._version import __version__, __revision__, changelog
 import os
 aps_local_val = os.environ.get('APS_LOCAL')
 check_for_newer = True
-if aps_local_val == None and check_for_newer:
+if aps_local_val is None and check_for_newer:
     try:
         import ftplib
         from distutils.version import StrictVersion
diff --git a/CEP/PyBDSM/src/python/readimage.py b/CEP/PyBDSM/src/python/readimage.py
index 2c5d56d416f..0e6378460f6 100644
--- a/CEP/PyBDSM/src/python/readimage.py
+++ b/CEP/PyBDSM/src/python/readimage.py
@@ -52,7 +52,7 @@ class Op_readimage(Op):
             img.opts.filename = img.opts.filename[:-1]
 
         # Determine indir if not explicitly given by user (in img.opts.indir)
-        if img.opts.indir == None:
+        if img.opts.indir is None:
             indir = os.path.dirname(img.opts.filename)
             if indir == '':
                 indir = './'
@@ -81,7 +81,7 @@ class Op_readimage(Op):
         # Read in data and header
         image_file = os.path.basename(img.opts.filename)
         result = read_image_from_file(image_file, img, img.indir)
-        if result == None:
+        if result is None:
             raise RuntimeError("Cannot open file " + repr(image_file) + ". " + img._reason)
         else:
             data, hdr = result
@@ -131,7 +131,7 @@ class Op_readimage(Op):
         self.init_beam(img)
         self.init_freq(img)
         year, code = self.get_equinox(img)
-        if year == None:
+        if year is None:
             mylog.info('Equinox not found in image header. Assuming J2000.')
             img.equinox = 2000.0
         else:
@@ -153,7 +153,7 @@ class Op_readimage(Op):
                 os.makedirs(img.basedir)
 
             # Now add solname (if any) and time to basedir
-            if img.opts.solnname != None:
+            if img.opts.solnname is not None:
                 img.basedir += img.opts.solnname + '_'
             img.basedir += time.strftime("%d%b%Y_%H.%M.%S")
 
@@ -432,13 +432,13 @@ class Op_readimage(Op):
                 from pywcs import WCS
 
         mylog = mylogger.logging.getLogger("PyBDSM.InitFreq")
-        if img.opts.frequency_sp != None and img.image_arr.shape[1] > 1:
+        if img.opts.frequency_sp is not None and img.image_arr.shape[1] > 1:
             # If user specifies multiple frequencies, then let
             # collapse.py do the initialization
             img.frequency = img.opts.frequency_sp[0]
             img.freq_pars = (0.0, 0.0, 0.0)
             mylog.info('Using user-specified frequencies.')
-        elif img.opts.frequency != None and img.image_arr.shape[1] == 1:
+        elif img.opts.frequency is not None and img.image_arr.shape[1] == 1:
             img.frequency = img.opts.frequency
             img.freq_pars = (img.frequency, 0.0, 0.0)
             mylog.info('Using user-specified frequency.')
@@ -482,7 +482,7 @@ class Op_readimage(Op):
                     instancemethod = type(img.wcs_obj.wcs_sky2pix)
                 img.wcs_obj.f2p = instancemethod(f2p, img.wcs_obj, WCS)
 
-                if img.opts.frequency != None:
+                if img.opts.frequency is not None:
                     img.frequency = img.opts.frequency
                 else:
                     img.frequency = img.wcs_obj.p2f(0)
@@ -530,7 +530,7 @@ class Op_readimage(Op):
 
         location specifies the location in pixels (x, y) for which angle is desired
         """
-        if location == None:
+        if location is None:
             x1 = img.image_arr.shape[2] / 2.0
             y1 = img.image_arr.shape[3] / 2.0
         else:
@@ -557,7 +557,7 @@ class Op_readimage(Op):
         """
         import functions as func
 
-        if location == None:
+        if location is None:
             x1 = int(img.image_arr.shape[2] / 2.0)
             y1 = int(img.image_arr.shape[3] / 2.0)
         else:
@@ -589,7 +589,7 @@ class Op_readimage(Op):
         """
         import functions as func
 
-        if location == None:
+        if location is None:
             x1 = int(img.image_arr.shape[2] / 2.0)
             y1 = int(img.image_arr.shape[3] / 2.0)
         else:
diff --git a/CEP/PyBDSM/src/python/rmsimage.py b/CEP/PyBDSM/src/python/rmsimage.py
index dcc610b89e8..a382a4cc86e 100644
--- a/CEP/PyBDSM/src/python/rmsimage.py
+++ b/CEP/PyBDSM/src/python/rmsimage.py
@@ -58,7 +58,7 @@ class Op_rmsimage(Op):
         # scales.
         fwsig = const.fwsig
         min_adapt_threshold = 10.0
-        if opts.adaptive_thresh == None:
+        if opts.adaptive_thresh is None:
             adapt_thresh = 50.0
             start_thresh = 500.0
         else:
@@ -259,7 +259,7 @@ class Op_rmsimage(Op):
                         img.use_rms_map = True
                     else:
                         self.check_rmsmap(img, rms)
-                elif opts.rms_map != None:
+                elif opts.rms_map is not None:
                     img.use_rms_map = opts.rms_map
                 if img.use_rms_map is False:
                     mylogger.userinfo(mylog, 'Using constant background rms')
@@ -277,7 +277,7 @@ class Op_rmsimage(Op):
 
           ## if rms map is insignificant, or rms_map==False use const value
           if img.use_rms_map is False:
-            if opts.rms_value == None:
+            if opts.rms_value is None:
               rms[:]  = crmss[ipol]
             else:
               rms[:]  = opts.rms_value
@@ -639,7 +639,7 @@ class Op_rmsimage(Op):
             new_shape = (arr.shape[0] + 2*pad_border_size, arr.shape[1]
                          + 2*pad_border_size)
             arr_pad = self.pad_array(arr, new_shape)
-            if mask == None:
+            if mask is None:
                 mask_pad = None
             else:
                 mask_pad = self.pad_array(mask, new_shape)
@@ -680,7 +680,7 @@ class Op_rmsimage(Op):
             rms_map[co] = cr
 
         # Check if all regions have too few unmasked pixels
-        if mask != None and N.size(N.where(mean_map != N.inf)) == 0:
+        if mask is not None and N.size(N.where(mean_map != N.inf)) == 0:
             raise RuntimeError("No unmasked regions from which to determine "\
                          "mean and rms maps")
 
@@ -870,7 +870,7 @@ class Op_rmsimage(Op):
 
         bstat = func.bstat#_cbdsm.bstat
         a, b, c, d = ind; i, j = co
-        if mask == None:
+        if mask is None:
           m, r, cm, cr, cnt = bstat(arr[a:b, c:d], mask, kappa)
           if cnt > 198: cm = m; cr = r
           mean_map[i, j], rms_map[i, j] = cm, cr
@@ -894,7 +894,7 @@ class Op_rmsimage(Op):
 
         bstat = func.bstat #_cbdsm.bstat
         a, b, c, d = ind
-        if mask == None:
+        if mask is None:
           m, r, cm, cr, cnt = bstat(arr[a:b, c:d], mask, kappa)
           if cnt > 198: cm = m; cr = r
         else:
diff --git a/CEP/PyBDSM/src/python/shapefit.py b/CEP/PyBDSM/src/python/shapefit.py
index 7b05327f38f..92b9657fce9 100755
--- a/CEP/PyBDSM/src/python/shapefit.py
+++ b/CEP/PyBDSM/src/python/shapefit.py
@@ -68,7 +68,7 @@ class Op_shapelets(Op):
 
         Returns shapelet parameters.
         """
-        if opts == None:
+        if opts is None:
             opts = img.opts
         if opts.shapelet_gresid:
             shape = img.shape
diff --git a/CEP/PyBDSM/src/python/shapelets.py b/CEP/PyBDSM/src/python/shapelets.py
index b695949d686..bba73770d61 100755
--- a/CEP/PyBDSM/src/python/shapelets.py
+++ b/CEP/PyBDSM/src/python/shapelets.py
@@ -94,7 +94,10 @@ def shapelet_image(basis, beta, centre, hc, nx, ny, size):
     try:
         from scipy import factorial
     except ImportError:
-        from scipy.misc.common import factorial
+        try:
+            from scipy.misc.common import factorial
+        except ImportError:
+            from scipy.misc import factorial
 
     hcx = hc[nx,:]
     hcy = hc[ny,:]
diff --git a/CEP/PyBDSM/src/python/spectralindex.py b/CEP/PyBDSM/src/python/spectralindex.py
index 355e076c489..b6414261ebf 100644
--- a/CEP/PyBDSM/src/python/spectralindex.py
+++ b/CEP/PyBDSM/src/python/spectralindex.py
@@ -302,19 +302,19 @@ class Op_spectralindex(Op):
 
         shp = img.image_arr.shape
         sbeam = img.opts.beam_spectrum
-        if sbeam != None and len(sbeam) != shp[1]: sbeam = None  # sanity check
-        if sbeam == None:
+        if sbeam is not None and len(sbeam) != shp[1]: sbeam = None  # sanity check
+        if sbeam is None:
             sbeam = [img.beam]*shp[1]
 
         img.beam_spectrum = sbeam
         img.freq = N.zeros(shp[1])
         crval, cdelt, crpix = img.freq_pars
         if img.wcs_obj.wcs.spec == -1 and \
-                img.opts.frequency_sp == None:
+                img.opts.frequency_sp is None:
             raise RuntimeError("Frequency info not found in header "\
                                    "and frequencies not specified by user")
         else:
-            if img.opts.frequency_sp == None:
+            if img.opts.frequency_sp is None:
                 for ichan in range(shp[1]):
                     img.freq[ichan] = img.wcs_obj.p2f(ichan)
             else:
diff --git a/CEP/PyBDSM/src/python/wavelet_atrous.py b/CEP/PyBDSM/src/python/wavelet_atrous.py
index edaeff4e348..908c9b1a931 100644
--- a/CEP/PyBDSM/src/python/wavelet_atrous.py
+++ b/CEP/PyBDSM/src/python/wavelet_atrous.py
@@ -117,7 +117,7 @@ class Op_wavelet_atrous(Op):
             suffix = 'w' + `j`
             filename = img.imagename + '.atrous.' + suffix + '.fits'
             if img.opts.output_all:
-                func.write_image_to_file(img.use_io, filename, w, img, bdir)
+                func.write_image_to_file('fits', filename, w, img, bdir)
                 mylog.info('%s %s' % ('Wrote ', img.imagename + '.atrous.' + suffix + '.fits'))
 
             # now do bdsm on each wavelet image.
@@ -302,13 +302,13 @@ class Op_wavelet_atrous(Op):
           img.total_flux_gaus += total_flux
           mylogger.userinfo(mylog, "Total flux density in model on all scales" , '%.3f Jy' % img.total_flux_gaus)
           if img.opts.output_all:
-              func.write_image_to_file(img.use_io, img.imagename + '.atrous.cJ.fits',
+              func.write_image_to_file('fits', img.imagename + '.atrous.cJ.fits',
                                        im_new, img, bdir)
               mylog.info('%s %s' % ('Wrote ', img.imagename + '.atrous.cJ.fits'))
-              func.write_image_to_file(img.use_io, img.imagename + '.resid_wavelets.fits',
+              func.write_image_to_file('fits', img.imagename + '.resid_wavelets.fits',
                                        (img.ch0_arr - img.resid_gaus_arr + img.resid_wavelets_arr), img, bdir + '/residual/')
               mylog.info('%s %s' % ('Wrote ', img.imagename + '.resid_wavelets.fits'))
-              func.write_image_to_file(img.use_io, img.imagename + '.model_wavelets.fits',
+              func.write_image_to_file('fits', img.imagename + '.model_wavelets.fits',
                                        (img.resid_gaus_arr - img.resid_wavelets_arr), img, bdir + '/model/')
               mylog.info('%s %s' % ('Wrote ', img.imagename + '.model_wavelets.fits'))
           img.completed_Ops.append('wavelet_atrous')
@@ -401,6 +401,9 @@ class Op_wavelet_atrous(Op):
         wimg.use_io = img.use_io
         wimg.do_cache = img.do_cache
         wimg.tempdir = img.tempdir
+        wimg.shape = img.shape
+        wimg.use_io = 'fits'
+
 
 ######################################################################################################
     def subtract_wvgaus(self, opts, residim, gaussians, islands):
-- 
GitLab