diff --git a/CEP/PyBDSM/src/python/gaul2srl.py b/CEP/PyBDSM/src/python/gaul2srl.py
index 4a812d590619e5ed2b482f77e0735d2a5dbe07e9..3c42de63332dd58f03c2ebc6722c4387d58e67ff 100644
--- a/CEP/PyBDSM/src/python/gaul2srl.py
+++ b/CEP/PyBDSM/src/python/gaul2srl.py
@@ -365,17 +365,23 @@ class Op_gaul2srl(Op):
         if N.isnan(mompara[1]):
             mompara[1] = posn[0] - delc[0]
         x1 = N.int(N.floor(mompara[1]))
+        if x1 < 0:
+            x1 = 0
         if N.isnan(mompara[2]):
             mompara[2] = posn[1] - delc[1]
         y1 = N.int(N.floor(mompara[2]))
+        if y1 < 0:
+            y1 = 0
         xind = slice(x1, x1+2, 1); yind = slice(y1, y1+2, 1)
         if img.opts.flag_smallsrc and (N.sum(mask[xind, yind]==N.ones((2,2))*isrc) != 4):
             mylog.debug('Island = '+str(isl.island_id))
             mylog.debug('Mask = '+repr(mask[xind, yind])+'xind, yind, x1, y1 = '+repr(xind)+' '+repr(yind)+' '+repr(x1)+' '+repr(y1))
         t=(mompara[1]-x1)/(x1+1-x1)  # in case u change it later
         u=(mompara[2]-y1)/(y1+1-y1)
-        s_peak=(1.0-t)*(1.0-u)*subim_src[x1,y1]+t*(1.0-u)*subim_src[x1+1,y1]+ \
-               t*u*subim_src[x1+1,y1+1]+(1.0-t)*u*subim_src[x1,y1+1]
+        xmax = N.min([x1+1, subim_src.shape[0] - 1])
+        ymax = N.min([y1+1, subim_src.shape[1] - 1])
+        s_peak=(1.0-t)*(1.0-u)*subim_src[x1,y1]+t*(1.0-u)*subim_src[xmax,y1]+ \
+               t*u*subim_src[xmax,ymax]+(1.0-t)*u*subim_src[x1,ymax]
         if (not img.opts.flag_smallsrc) and (N.sum(mask[xind, yind]==N.ones((2,2))*isrc) != 4):
             mylog.debug('Speak '+repr(s_peak)+'Mompara = '+repr(mompara))
             mylog.debug('x1, y1 : '+repr(x1)+', '+repr(y1))
diff --git a/CEP/PyBDSM/src/python/gausfit.py b/CEP/PyBDSM/src/python/gausfit.py
index 715bf578783eb23dd34db24558485314eb9d8ad7..6060f307c88e2d07d127df1781afddc3ff264315 100644
--- a/CEP/PyBDSM/src/python/gausfit.py
+++ b/CEP/PyBDSM/src/python/gausfit.py
@@ -266,8 +266,8 @@ class Op_gausfit(Op):
         if abs(beam[0]/beam[1]) < 1.1:
             beam = (1.1*beam[0], beam[1], beam[2])
 
-        thr1 = isl.mean + opts.thresh_isl*isl.rms
-        thr2 = isl.mean + img.thresh_pix*isl.rms
+        thr1 = opts.thresh_isl*isl.rms #isl.mean + opts.thresh_isl*isl.rms
+        thr2 = img.thresh_pix*isl.rms #isl.mean + img.thresh_pix*isl.rms
         thr0 = thr1
         verbose = opts.verbose_fitting
         peak = fcn.find_peak()[0]
@@ -411,7 +411,7 @@ class Op_gausfit(Op):
         thresh_isl = opts.thresh_isl
         thresh_pix = opts.thresh_pix
         thresh = opts.fittedimage_clip
-        thr = isl.mean + thresh_isl * isl.rms
+        thr = thresh_isl * isl.rms #isl.mean + thresh_isl * isl.rms
         rms = isl.rms
 
         if opts.verbose_fitting:
diff --git a/CEP/PyBDSM/src/python/islands.py b/CEP/PyBDSM/src/python/islands.py
index b3e4c94b08d20d438552438e575fe80b32d35cca..17af33cccc64b54758eec32c344761db37706bf8 100644
--- a/CEP/PyBDSM/src/python/islands.py
+++ b/CEP/PyBDSM/src/python/islands.py
@@ -289,6 +289,7 @@ class Island(object):
         self.rms  = bbox_rms_im[in_bbox_and_unmasked].mean()
         in_bbox_and_unmasked = N.where(~N.isnan(bbox_mean_im))
         self.mean  = bbox_mean_im[in_bbox_and_unmasked].mean()
+        self.islmean  = bbox_mean_im[in_bbox_and_unmasked].mean()
         self.total_flux = N.nansum(self.image[in_bbox_and_unmasked])/beamarea
         pixels_in_isl = N.sum(~N.isnan(self.image[self.mask_active])) # number of unmasked pixels assigned to current island
         self.total_fluxE = func.nanmean(bbox_rms_im[in_bbox_and_unmasked]) * N.sqrt(pixels_in_isl/beamarea) # Jy