diff --git a/CEP/PyBDSM/src/python/_version.py b/CEP/PyBDSM/src/python/_version.py
index e172c49f04c381d2ea9f413720b58e58998c194e..7f5c8a565a45f010a9355dfda32e4f4e72f108b7 100644
--- a/CEP/PyBDSM/src/python/_version.py
+++ b/CEP/PyBDSM/src/python/_version.py
@@ -27,6 +27,11 @@ def changelog():
PyBDSM Changelog.
-----------------------------------------------------------------------
+ 2012/07/05 - Fixed a bug that caused images written when output_all=
+ True to be transposed. Added frequency information to all output
+ images. Improved fitting robustness to prevent rare cases in
+ which no valid Gaussians could be fit to an island.
+
2012/07/03 - Version 1.3
2012/07/03 - Fixed a bug in calculation of the positional errors of
diff --git a/CEP/PyBDSM/src/python/collapse.py b/CEP/PyBDSM/src/python/collapse.py
index 8455e3656de5ab6b5a4f887c5835aad96dcbdce1..c10d8cba75f715f5101d272810a6348155db5fa8 100644
--- a/CEP/PyBDSM/src/python/collapse.py
+++ b/CEP/PyBDSM/src/python/collapse.py
@@ -99,7 +99,7 @@ class Op_collapse(Op):
mylog.debug('%s %s %s' % ('Channel weights : ', str1, '; unity=zero if c_wts="rms"'))
if img.opts.output_all:
- func.write_image_to_file(img.use_io, img.imagename+'.ch0_'+pol+'.fits', N.transpose(ch0), img)
+ func.write_image_to_file(img.use_io, img.imagename+'.ch0_'+pol+'.fits', ch0, img)
mylog.debug('%s %s ' % ('Writing file ', img.imagename+'.ch0_'+pol+'.fits'))
else:
@@ -243,6 +243,7 @@ def init_freq_collapse(img, wtarr):
sumwts += wtarr[i]
sumfrq += freqs[ch]*wtarr[i]
img.cfreq = sumfrq / sumwts
+ img.freq_pars = (img.cfreq, 0.0, 0.0)
else:
# Calculate from header info
c_list = img.opts.collapse_av
diff --git a/CEP/PyBDSM/src/python/functions.py b/CEP/PyBDSM/src/python/functions.py
index f1b668eb1773474bbe43d731fde1bbe071081adc..66cd4ad509c73ffddb7f57492344a5a0444e6fef 100755
--- a/CEP/PyBDSM/src/python/functions.py
+++ b/CEP/PyBDSM/src/python/functions.py
@@ -1198,18 +1198,19 @@ def write_image_to_file(use, filename, image, img, outdir=None,
# im.saveas(indir+filename)
def make_fits_image(imagedata, wcsobj, beam, freq):
- """Makes a simple FITS hdulist appropriate for images"""
+ """Makes a simple FITS hdulist appropriate for single-channel images"""
import pyfits
- hdu = pyfits.PrimaryHDU(imagedata)
+ shape_out = [1, imagedata.shape[0], imagedata.shape[1]]
+ hdu = pyfits.PrimaryHDU(imagedata.reshape(shape_out))
hdulist = pyfits.HDUList([hdu])
header = hdulist[0].header
header.update('CTYPE1', wcsobj.ctype[0])
- header.update('CTYPE2', wcsobj.ctype[1])
header.update('CRVAL1', wcsobj.crval[0])
- header.update('CRVAL2', wcsobj.crval[1])
header.update('CDELT1', wcsobj.cdelt[0])
- header.update('CDELT2', wcsobj.cdelt[1])
header.update('CRPIX1', wcsobj.crpix[0])
+ header.update('CTYPE2', wcsobj.ctype[1])
+ header.update('CRVAL2', wcsobj.crval[1])
+ header.update('CDELT2', wcsobj.cdelt[1])
header.update('CRPIX2', wcsobj.crpix[1])
if hasattr(wcsobj, 'crota'):
header.update('CROTA1', wcsobj.crota[0])
@@ -1217,6 +1218,7 @@ def make_fits_image(imagedata, wcsobj, beam, freq):
header.update('BMAJ', beam[0])
header.update('BMIN', beam[1])
header.update('BPA', beam[2])
+ header.update('CTYPE3', 'FREQ')
header.update('CRVAL3', freq[0])
header.update('CDELT3', freq[1])
header.update('CRPIX3', freq[2])
diff --git a/CEP/PyBDSM/src/python/gausfit.py b/CEP/PyBDSM/src/python/gausfit.py
index bf5642e9796375fa744df17dbee282c085e2a117..70ddcc5909ea52f4a5a2944778984701c57f234a 100644
--- a/CEP/PyBDSM/src/python/gausfit.py
+++ b/CEP/PyBDSM/src/python/gausfit.py
@@ -260,11 +260,53 @@ class Op_gausfit(Op):
if fitok and len(fgaul) == 0:
break
if not fitok:
- # If all else fails, try to fit just a single Gaussian to the island
- fitok = self.fit_iter([], 0, fcn, dof, beam, thr0, 1, 'simple', 1, verbose)
- gaul, fgaul = self.flag_gaussians(fcn.parameters, opts,
+ # If normal fitting fails, try to fit 5 or fewer Gaussians to the island
+ ngmax = 5
+ while not fitok and ngmax > 1:
+ fitok = self.fit_iter([], 0, fcn, dof, beam, thr0, 1, 'simple', ngmax, verbose)
+ ngmax -= 1
+ gaul, fgaul = self.flag_gaussians(fcn.parameters, opts,
beam, thr0, peak, shape, isl.mask_active,
isl.image, size)
+ if not fitok:
+ # If all else fails, shrink the island a little and try one last time
+ if ffimg == None:
+ fcn = MGFunction(isl.image-isl.islmean, nd.binary_dilation(isl.mask_active), 1)
+ else:
+ fcn = MGFunction(isl.image-isl.islmean-ffimg, nd.binary_dilation(isl.mask_active), 1)
+ gaul = []
+ iter = 0
+ ng1 = 0
+ ngmax = 25
+ while iter < 5:
+ iter += 1
+ fitok = self.fit_iter(gaul, ng1, fcn, dof, beam, thr0, iter, 'simple', ngmax, verbose)
+ gaul, fgaul = self.flag_gaussians(fcn.parameters, opts,
+ beam, thr0, peak, shape, isl.mask_active,
+ isl.image, size)
+ ng1 = len(gaul)
+ if fitok and len(fgaul) == 0:
+ break
+ if not fitok:
+ # If all else fails, expand the island a little and try one last time
+ if ffimg == None:
+ fcn = MGFunction(isl.image-isl.islmean, nd.binary_erosion(isl.mask_active), 1)
+ else:
+ fcn = MGFunction(isl.image-isl.islmean-ffimg, nd.binary_erosion(isl.mask_active), 1)
+ gaul = []
+ iter = 0
+ ng1 = 0
+ ngmax = 25
+ while iter < 5:
+ iter += 1
+ fitok = self.fit_iter(gaul, ng1, fcn, dof, beam, thr0, iter, 'simple', ngmax, verbose)
+ gaul, fgaul = self.flag_gaussians(fcn.parameters, opts,
+ beam, thr0, peak, shape, isl.mask_active,
+ isl.image, size)
+ ng1 = len(gaul)
+ if fitok and len(fgaul) == 0:
+ break
+
### return whatever we got
isl.mg_fcn = fcn
diff --git a/CEP/PyBDSM/src/python/make_residimage.py b/CEP/PyBDSM/src/python/make_residimage.py
index 5363815c6da4510f031ee74dac7bd1a653af09e1..7d324dc09eeee196e13e067c0bc9954ee24ce429 100644
--- a/CEP/PyBDSM/src/python/make_residimage.py
+++ b/CEP/PyBDSM/src/python/make_residimage.py
@@ -77,9 +77,9 @@ class Op_make_residimage(Op):
moddir = img.basedir + '/model/'
if not os.path.exists(resdir): os.mkdir(resdir)
if not os.path.exists(moddir): os.mkdir(moddir)
- func.write_image_to_file(img.use_io, img.imagename + '.resid_gaus.fits', N.transpose(img.resid_gaus), img, resdir)
+ func.write_image_to_file(img.use_io, img.imagename + '.resid_gaus.fits', img.resid_gaus, img, resdir)
mylog.info('%s %s' % ('Writing', resdir+img.imagename+'.resid_gaus.fits'))
- func.write_image_to_file(img.use_io, img.imagename + '.model.fits', N.transpose(img.ch0 - img.resid_gaus), img, moddir)
+ func.write_image_to_file(img.use_io, img.imagename + '.model.fits', (img.ch0 - img.resid_gaus), img, moddir)
mylog.info('%s %s' % ('Writing', moddir+img.imagename+'.model_gaus.fits'))
### residual rms and mean per island
@@ -137,7 +137,7 @@ class Op_make_residimage(Op):
img.resid_shap[pix_masked] = N.nan
if img.opts.output_all:
- func.write_image_to_file(img.use_io, img.imagename + '.resid_shap.fits', N.transpose(img.resid_shap), img, resdir)
+ func.write_image_to_file(img.use_io, img.imagename + '.resid_shap.fits', img.resid_shap, img, resdir)
mylog.info('%s %s' % ('Writing ', resdir+img.imagename+'.resid_shap.fits'))
### shapelet residual rms and mean per island
diff --git a/CEP/PyBDSM/src/python/psf_vary.py b/CEP/PyBDSM/src/python/psf_vary.py
index c772ab5a41b37b147f4a721b669ef701d125c901..f4379cfeb2956c909682f642d2189fa16b9fa1cf 100644
--- a/CEP/PyBDSM/src/python/psf_vary.py
+++ b/CEP/PyBDSM/src/python/psf_vary.py
@@ -132,7 +132,7 @@ class Op_psf_vary(Op):
volrank, vorowts = self.tesselate(vorogenP, vorogenS, tile_prop, tess_method, tess_sc, tess_fuzzy, \
generators, gencode, image.shape)
if opts.output_all:
- func.write_image_to_file(img.use_io, img.imagename + '.volrank.fits', N.transpose(volrank), img, dir)
+ func.write_image_to_file(img.use_io, img.imagename + '.volrank.fits', volrank, img, dir)
tile_list, tile_coord, tile_snr = tile_prop
ntile = len(tile_list)
@@ -252,11 +252,11 @@ class Op_psf_vary(Op):
img.psf_vary_ratio_aper = psf_ratio_aper_int
if opts.output_all:
- func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_maj.fits', N.transpose(psf_maj_int)*fwsig, img, dir)
- func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_min.fits', N.transpose(psf_min_int)*fwsig, img, dir)
- func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_pa.fits', N.transpose(psf_pa_int), img, dir)
- func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_ratio.fits', N.transpose(psf_ratio_int), img, dir)
- func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_ratio_aper.fits', N.transpose(psf_ratio_aper_int), img, dir)
+ func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_maj.fits', psf_maj_int*fwsig, img, dir)
+ func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_min.fits', psf_min_int*fwsig, img, dir)
+ func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_pa.fits', psf_pa_int, img, dir)
+ func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_ratio.fits', psf_ratio_int, img, dir)
+ func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_ratio_aper.fits', psf_ratio_aper_int, img, dir)
# Loop through source and Gaussian lists and deconvolve the sizes using appropriate beam
bar2 = statusbar.StatusBar('Correcting deconvolved source sizes ..... : ', 0, img.nsrc)
diff --git a/CEP/PyBDSM/src/python/readimage.py b/CEP/PyBDSM/src/python/readimage.py
index 7585268146ea448336d357c76634411d4ce306e0..cd4c6f2392323f565d4dbce3dce4da8fe9e15a6d 100644
--- a/CEP/PyBDSM/src/python/readimage.py
+++ b/CEP/PyBDSM/src/python/readimage.py
@@ -379,7 +379,7 @@ class Op_readimage(Op):
mylog.info('Using user-specified frequencies.')
elif img.opts.frequency != None and img.image.shape[1] == 1:
img.cfreq = img.opts.frequency
- img.freq_pars = (0.0, 0.0, 0.0)
+ img.freq_pars = (img.cfreq, 0.0, 0.0)
mylog.info('Using user-specified frequency.')
else:
found = False
diff --git a/CEP/PyBDSM/src/python/rmsimage.py b/CEP/PyBDSM/src/python/rmsimage.py
index 60004381bf96a9c2b401a55323aeb0dbfb9e3da0..184ab4dedd18109f815913b6f588230b5dba41f4 100644
--- a/CEP/PyBDSM/src/python/rmsimage.py
+++ b/CEP/PyBDSM/src/python/rmsimage.py
@@ -331,7 +331,7 @@ class Op_rmsimage(Op):
else:
resdir = img.basedir + '/background/'
if not os.path.exists(resdir): os.mkdir(resdir)
- func.write_image_to_file(img.use_io, img.imagename + '.rmsd_I.fits', N.transpose(rms), img, resdir)
+ func.write_image_to_file(img.use_io, img.imagename + '.rmsd_I.fits', rms, img, resdir)
mylog.info('%s %s' % ('Writing ', resdir+img.imagename+'.rmsd_I.fits'))
if opts.savefits_meanim or opts.output_all:
if img.waveletimage:
@@ -339,7 +339,7 @@ class Op_rmsimage(Op):
else:
resdir = img.basedir + '/background/'
if not os.path.exists(resdir): os.mkdir(resdir)
- func.write_image_to_file(img.use_io, img.imagename + '.mean_I.fits', N.transpose(mean), img, resdir)
+ func.write_image_to_file(img.use_io, img.imagename + '.mean_I.fits', mean, img, resdir)
mylog.info('%s %s' % ('Writing ', resdir+img.imagename+'.mean_I.fits'))
if opts.savefits_normim or opts.output_all:
if img.waveletimage:
@@ -347,7 +347,7 @@ class Op_rmsimage(Op):
else:
resdir = img.basedir + '/background/'
if not os.path.exists(resdir): os.mkdir(resdir)
- func.write_image_to_file(img.use_io, img.imagename + '.norm_I.fits', N.transpose((img.ch0-mean)/rms), img, resdir)
+ func.write_image_to_file(img.use_io, img.imagename + '.norm_I.fits', (img.ch0-mean)/rms, img, resdir)
mylog.info('%s %s' % ('Writing ', resdir+img.imagename+'.norm_I.fits'))
else:
img.mean_QUV.append(mean); img.rms_QUV.append(rms)
diff --git a/CEP/PyBDSM/src/python/wavelet_atrous.py b/CEP/PyBDSM/src/python/wavelet_atrous.py
index 01a5acd90e0a243a6fd0a0cb4d817fc009769af5..8af43910cb9b7a3c589a16b4453c1892ce0ee811 100644
--- a/CEP/PyBDSM/src/python/wavelet_atrous.py
+++ b/CEP/PyBDSM/src/python/wavelet_atrous.py
@@ -107,7 +107,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.transpose(), img, bdir)
+ func.write_image_to_file(img.use_io, filename, w, img, bdir)
mylog.info('%s %s' % ('Wrote ', img.imagename + '.atrous.' + suffix + '.fits'))
# now do bdsm on each wavelet image
if dobdsm:
@@ -255,13 +255,13 @@ class Op_wavelet_atrous(Op):
mylogger.userinfo(mylog, "Total flux density in model over 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',
- im_new.transpose(), img, bdir)
+ 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',
- N.transpose(img.ch0 - img.resid_gaus + img.resid_wavelets), img, bdir + '/residual/')
+ (img.ch0 - img.resid_gaus + img.resid_wavelets), 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',
- N.transpose(img.resid_gaus - img.resid_wavelets), img, bdir + '/model/')
+ (img.resid_gaus - img.resid_wavelets), img, bdir + '/model/')
mylog.info('%s %s' % ('Wrote ', img.imagename + '.model_wavelets.fits'))
img.completed_Ops.append('wavelet_atrous')