diff --git a/CEP/PyBDSM/doc/source/conf.py b/CEP/PyBDSM/doc/source/conf.py
index 5094c7165ccadcf70b68f6a7f18a6124f23ac7cd..a7508205f40def162a2ac03200bdf5f19bc5d46c 100644
--- a/CEP/PyBDSM/doc/source/conf.py
+++ b/CEP/PyBDSM/doc/source/conf.py
@@ -41,7 +41,7 @@ master_doc = 'index'
 
 # General information about the project.
 project = u'PyBDSM'
-copyright = u'2014, David Rafferty and Niruj Mohan'
+copyright = u'2016, David Rafferty and Niruj Mohan'
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
@@ -50,7 +50,7 @@ copyright = u'2014, David Rafferty and Niruj Mohan'
 # The short X.Y version.
 version = '1.8'
 # The full version, including alpha/beta/rc tags.
-release = '1.8.4'
+release = '1.8.6'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
diff --git a/CEP/PyBDSM/doc/source/process_image.rst b/CEP/PyBDSM/doc/source/process_image.rst
index 3c184fc1c6809680345d2ce05e3ea51be6d83788..0f6f2e318e857dbec1c959898db5a68681732fe6 100644
--- a/CEP/PyBDSM/doc/source/process_image.rst
+++ b/CEP/PyBDSM/doc/source/process_image.rst
@@ -906,6 +906,7 @@ The options for this module are as follows:
                                    original image
       :term:`atrous_sum` .......... True : Fit to the sum of images of the remaining wavelet
                                    scales
+      :term:`use_scipy_fft` ....... True : Use fast SciPy FFT for convolution
 
 .. glossary::
 
@@ -939,6 +940,9 @@ The options for this module are as follows:
         This parameter is a Boolean (default is ``True``). If ``True``, fitting is done on an image that is the sum of the remaining wavelet scales. Using the sum will generally result in improved signal.
         If ``False``, fitting is done on only the wavelet scale under consideration.
 
+    use_scipy_fft
+        This parameter is a Boolean (default is ``True``). If ``True``, the SciPy FFT function will be used instead of the custom version. The SciPy version is much faster but also uses much more memory.
+
 .. _psf_vary_do:
 
 PSF variation module
diff --git a/CEP/PyBDSM/doc/source/whats_new.rst b/CEP/PyBDSM/doc/source/whats_new.rst
index f15f39f4ee9a0106f7c4f33689cd0d74c3bfd29c..7d3ac6267c03e72ada221da7302ae4e20405a610 100644
--- a/CEP/PyBDSM/doc/source/whats_new.rst
+++ b/CEP/PyBDSM/doc/source/whats_new.rst
@@ -4,6 +4,25 @@
 What's New
 **********
 
+Version 1.8.6 (2016/01/20)
+
+    * Fix to bug that caused incorrect island mask when two islands are very close together.
+
+    * Fix to bug that caused crash when image is masked and the ``src_ra_dec`` option is used.
+
+
+Version 1.8.5 (2015/11/30):
+
+    * Fix to bug in ``export_image`` that resulted in incorrect output image when both ``trim_box`` and ``pad_image`` were used.
+
+    * Fix to bug in wavelet module related to merging of islands.
+
+    * Fix to bug in polarization module related to numbering of new islands.
+
+    * Added option to use much faster (but also much more memory intensive) SciPy ``fftconvolve`` function instead of custom PyBDSM one. The option (``use_scipy_fft``) defaults to ``True``.
+
+    * Increased number of digits for values in output text catalogs
+
 Version 1.8.4 (2015/08/06):
 
     * Improved speed of wavelet module.
@@ -14,7 +33,7 @@ Version 1.8.4 (2015/08/06):
 
     * Fix to bug that caused a failure to write shapelet models in FITS format.
 
-    * Fix to bug that caused a crash when both atrous_do = True and output_all = True.
+    * 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.
 
@@ -36,7 +55,7 @@ Version 1.8.1 (2014/01/14):
 
     * Added option (``bbs_patches = 'mask'``) to allow patches in an output BBS sky model to be defined using a mask image (set with the ``bbs_patches_mask`` option).
 
-    * Fix to bug that caused the ``incl_empty`` option to be ignored when ``format='fits'`` in the ``write_catalog`` task.
+    * Fix to bug that caused the ``incl_empty`` option to be ignored when ``format = 'fits'`` in the ``write_catalog`` task.
 
     * Enabled output of images in CASA format in the ``export_image`` task (``img_format = 'casa'``).
 
diff --git a/CEP/PyBDSM/src/python/_version.py b/CEP/PyBDSM/src/python/_version.py
index d6ba134a24512463223a09f2323adcde85124b95..6b965579414c14023b5d24f5408d2c8e6c15fb72 100644
--- a/CEP/PyBDSM/src/python/_version.py
+++ b/CEP/PyBDSM/src/python/_version.py
@@ -9,7 +9,7 @@ adding to the changelog will naturally do this.
 """
 
 # Version number
-__version__ = '1.8.5'
+__version__ = '1.8.6'
 
 # Store svn Revision number. For this to work, one also needs to do:
 #
@@ -26,9 +26,16 @@ def changelog():
     """
     PyBDSM Changelog.
     -----------------------------------------------------------------------
-    2015/08/06 - Version 1.8.4
 
-    2015/10/06 - Version 1.8.5
+    2016/01/20 - Version 1.8.6
+
+    2016/01/15 - Fix to bug that caused incorrect island mask when two
+        islands are very close together.
+
+    2015/12/07 - Fix to bug that caused crash when image is masked and
+        the src_ra_dec option is used.
+
+    2015/11/30 - Version 1.8.5
 
     2015/11/25 - Fix to bug in export_image that resulted in incorrect
         output image when both trim_box and pad_image were used.
@@ -38,6 +45,9 @@ def changelog():
     2015/11/20 - Fix to bug in polarization module related to numbering of
         new islands.
 
+    2015/11/20 - Fix to bug in spectral index module related to rms map
+        calculation.
+
     2015/11/20 - Added option to use much faster (but also much more memory
         intensive) SciPy fftconvolve function instead of custom PyBDSM one.
         The option (use_scipy_fft) defaults to True.
diff --git a/CEP/PyBDSM/src/python/islands.py b/CEP/PyBDSM/src/python/islands.py
index 02fe77e24eff3a9bd1c81ff56b1e743fd8594498..5f845ee6924717e0283c0bda9d44ff7fb9c7d9ce 100644
--- a/CEP/PyBDSM/src/python/islands.py
+++ b/CEP/PyBDSM/src/python/islands.py
@@ -110,13 +110,11 @@ class Op_islands(Op):
 
             ch0_map = img.ch0_arr
             ch0_shape = ch0_map.shape
-            pyrank = N.zeros(ch0_shape, dtype=N.int32) - 1
+            pyrank = N.zeros(ch0_shape, dtype=N.int32)
             for i, isl in enumerate(img.islands):
                 isl.island_id = i
-                if i == 0:
-                    pyrank[isl.bbox] = N.invert(isl.mask_active) - 1
-                else:
-                    pyrank[isl.bbox] = N.invert(isl.mask_active) * i - isl.mask_active
+                pyrank[isl.bbox] += N.invert(isl.mask_active) * (i + 1)
+            pyrank -= 1 # align pyrank values with island ids and set regions outside of islands to -1
 
             if opts.output_all: write_islands(img)
             if opts.savefits_rankim:
@@ -219,7 +217,7 @@ class Op_islands(Op):
             labels = func.make_src_mask(image.shape,
                         isl_posn_pix, isl_radius_pix)
             if img.masked:
-                aper_mask = N.where(labels.astype(bool) & ~img.mask)
+                aper_mask = N.where(labels.astype(bool) & ~mask)
             else:
                 aper_mask = N.where(labels.astype(bool))
             if N.size(aper_mask) > img.minpix_isl:
diff --git a/CEP/PyBDSM/src/python/output.py b/CEP/PyBDSM/src/python/output.py
index c7c154a0b70c83e3e7d890b1d3545a4d933fa3b1..fdef24829e4b3ab678b41fc14fa7a61a671c0649 100644
--- a/CEP/PyBDSM/src/python/output.py
+++ b/CEP/PyBDSM/src/python/output.py
@@ -326,7 +326,7 @@ def write_ascii_list(img, filename=None, sort_by='indx', format = 'ascii',
         if incl_empty:
             # Append the dummy sources for islands without any unflagged Gaussians
             outl[0] += img.dsources
-    outstr_list = make_ascii_str(img, outl, objtype=objtype,
+    outstr_list = make_ascii_str(img, outl, objtype=objtype, incl_chan=incl_chan,
                                  incl_empty=incl_empty, format=format)
     if filename is None:
         if objtype == 'gaul':
@@ -409,7 +409,7 @@ def write_fits_list(img, filename=None, sort_by='index', objtype='gaul',
         cvals, cnames, cformats, cunits = make_output_columns(outl[0][0], fits=True,
                                                               objtype=objtype,
                                                               incl_spin=img.opts.spectralindex_do,
-                                                              incl_chan=img.opts.incl_chan,
+                                                              incl_chan=incl_chan,
                                                               incl_pol=img.opts.polarisation_do,
                                                               incl_aper=incl_aper,
                                                               incl_empty=incl_empty,
@@ -791,7 +791,8 @@ def make_ds9_str(img, glist, gnames, deconvolve=False, objtype='gaul', incl_empt
     return outstr_list
 
 
-def make_ascii_str(img, glist, objtype='gaul', format='ascii', incl_empty=False):
+def make_ascii_str(img, glist, objtype='gaul', format='ascii', incl_empty=False,
+    incl_chan=False):
     """Makes a list of string entries for an ascii region file."""
     from _version import __version__, __revision__
     outstr_list = []
@@ -815,7 +816,7 @@ def make_ascii_str(img, glist, objtype='gaul', format='ascii', incl_empty=False)
         cvals, cnames, cformats, cunits = make_output_columns(g, fits=False,
                                                               objtype=objtype,
                                                               incl_spin=img.opts.spectralindex_do,
-                                                              incl_chan=img.opts.incl_chan,
+                                                              incl_chan=incl_chan,
                                                               incl_pol=img.opts.polarisation_do,
                                                               incl_aper=incl_aper,
                                                               incl_empty=incl_empty,
@@ -833,7 +834,8 @@ def make_ascii_str(img, glist, objtype='gaul', format='ascii', incl_empty=False)
     return outstr_list
 
 
-def make_fits_list(img, glist, objtype='gaul', nmax=30, incl_empty=False):
+def make_fits_list(img, glist, objtype='gaul', nmax=30, incl_empty=False,
+    incl_chan=False):
     import functions as func
 
     out_list = []
@@ -844,7 +846,7 @@ def make_fits_list(img, glist, objtype='gaul', nmax=30, incl_empty=False):
     for g in glist[0]:
         cvals, ext1, ext2, ext3 = make_output_columns(g, fits=True, objtype=objtype,
                                                       incl_spin=img.opts.spectralindex_do,
-                                                      incl_chan=img.opts.incl_chan,
+                                                      incl_chan=incl_chan,
                                                       incl_pol=img.opts.polarisation_do,
                                                       incl_aper=incl_aper,
                                                       incl_empty=incl_empty,
diff --git a/CEP/PyBDSM/src/python/wavelet_atrous.py b/CEP/PyBDSM/src/python/wavelet_atrous.py
index 2d62be65ce400abdd2d44b541c10a5785d693edd..89253a22c94e472b49766215ece800250d7cd979 100644
--- a/CEP/PyBDSM/src/python/wavelet_atrous.py
+++ b/CEP/PyBDSM/src/python/wavelet_atrous.py
@@ -300,16 +300,16 @@ class Op_wavelet_atrous(Op):
               if stop_wav == True:
                   break
 
+          pyrank = N.zeros(img.pyrank.shape, dtype=N.int32)
           for i, isl in enumerate(img.islands):
               isl.island_id = i
               for g in isl.gaul:
                   g.island_id = i
               for dg in isl.dgaul:
                   dg.island_id = i
-              if i == 0:
-                  img.pyrank[isl.bbox] = N.invert(isl.mask_active) - 1
-              else:
-                  img.pyrank[isl.bbox] = N.invert(isl.mask_active) * isl.island_id - isl.mask_active
+              pyrank[isl.bbox] += N.invert(isl.mask_active) * (i + 1)
+          pyrank -= 1 # align pyrank values with island ids and set regions outside of islands to -1
+          img.pyrank = pyrank
 
           pdir = img.basedir + '/misc/'
           img.ngaus += ntot_wvgaus
@@ -632,16 +632,16 @@ def renumber_islands(img):
 
     Also renumbers the pyrank image.
     """
+    pyrank = N.zeros(img.pyrank.shape, dtype=N.int32)
     for i, isl in enumerate(img.islands):
         isl.island_id = i
         for g in isl.gaul:
             g.island_id = i
         for dg in isl.dgaul:
             dg.island_id = i
-        if i == 0:
-            img.pyrank[isl.bbox] = N.invert(isl.mask_active) - 1
-        else:
-            img.pyrank[isl.bbox] = N.invert(isl.mask_active) * isl.island_id - isl.mask_active
+            pyrank[isl.bbox] += N.invert(isl.mask_active) * (i + 1)
+    pyrank -= 1 # align pyrank values with island ids and set regions outside of islands to -1
+    img.pyrank = pyrank
     gaussian_list = [g for isl in img.islands for g in isl.gaul]
     img.gaussians = gaussian_list