diff --git a/bdsf/_version.py b/bdsf/_version.py index 459f4139..7c1408c2 100644 --- a/bdsf/_version.py +++ b/bdsf/_version.py @@ -205,14 +205,14 @@ def changelog(): ignored when format='fits' in the write_catalog task. 2013/12/05 - Enabled output of images in CASA format in the export_image - task (img_format = 'casa'). Added an option to export_image task to - export an island-mask image, with ones where there is emission and - zeros elsewhere (image_type = 'island_mask'). Features in the island - mask may be optionally dilated by specifying the number of dilation - iterations with the mask_dilation parameter. Added an option to - write a CASA region file to the write_catalog task (format = - 'casabox'). Added an option to write a CSV catalog to the - write_catalog task (format = 'csv'). + task (img_format = 'casa'). Added an option to export_image task to + export an island-mask image, with ones where there is emission and + zeros elsewhere (image_type = 'island_mask'). Features in the island + mask may be optionally dilated by specifying the number of dilation + iterations with the mask_dilation parameter. Added an option to + write a CASA region file to the write_catalog task (format = + 'casabox'). Added an option to write a CSV catalog to the + write_catalog task (format = 'csv'). 2013/11/04 - Added error message when the rms is zero in some part of the rms map. @@ -220,11 +220,11 @@ def changelog(): 2013/10/16 - Version 1.8.0 2013/10/16 - Improved wavelet fitting. Added option so that wavelet - fitting can be done to the sum of images on the remaining wavelet - scales, improving the signal for fitting (controlled with the - atrous_sum option). Added option so that user can choose whether to - include new islands found only in the wavelet images in the final - fit or not (controlled with the atrous_orig_isl option). + fitting can be done to the sum of images on the remaining wavelet + scales, improving the signal for fitting (controlled with the + atrous_sum option). Added option so that user can choose whether to + include new islands found only in the wavelet images in the final + fit or not (controlled with the atrous_orig_isl option). 2013/10/10 - Fixed a bug that could lead to incomplete fitting of some islands. Improved overall convergence of fits. @@ -291,8 +291,8 @@ def changelog(): (replaced scipy.signal.fftconvolve with a custom function). 2013/08/11 - Added option to use disk caching for internally derived - images (do_cache). Caching can reduce memory usage and is - therefore useful when processing large images. + images (do_cache). Caching can reduce memory usage and is + therefore useful when processing large images. 2013/07/11 - Implemented a variable operation chain for process_image (and img.process()) that allows unneeded steps to be skipped if @@ -302,9 +302,9 @@ def changelog(): during iterative fitting of large islands. 2013/06/24 - Added option (fix_to_beam) to fix the size and position - angle of Gaussians to the restoring beam during fitting. Fix to - bug that caused the position angle used to initialize fitting to - be incorrect. + angle of Gaussians to the restoring beam during fitting. Fix to + bug that caused the position angle used to initialize fitting to + be incorrect. 2013/03/22 - Version 1.6.1 @@ -427,9 +427,9 @@ def changelog(): no Gaussians were fit to the ch0 image. 2012/09/19 - Added option (broadcast) to show_fit task to send - coordinates and row highlight request to a SAMP hub when a Gaussian - is clicked. Fixed bug in aperture flux masking that sometimes caused - the mask to be the wrong shape. + coordinates and row highlight request to a SAMP hub when a Gaussian + is clicked. Fixed bug in aperture flux masking that sometimes caused + the mask to be the wrong shape. 2012/09/18 - Added option to send images and catalogs to a SAMP hub (activated by setting outfile = 'SAMP' in the export_image and @@ -514,196 +514,196 @@ def changelog(): Changed code that determines terminal width to be more robust. 2012/05/07 - Removed dependencies on matplotlib -- if matplotlib is - not available, plotting is disabled. Corrected inconsistencies, - spelling mistakes, etc. in help text and documentation. Cleaned - up unneeded modules and files. + not available, plotting is disabled. Corrected inconsistencies, + spelling mistakes, etc. in help text and documentation. Cleaned + up unneeded modules and files. 2012/05/02 - Added option to output flux densities for every channel - found by the spectral index module. Added option to spectral index - module to allow use of flux densities that do not meet the desired - SNR. Changed flag_maxsnr criterion to also flag if the peak flux - density per beam of the Gaussian exceeds the value at its center. - Removed incl_wavelet option. + found by the spectral index module. Added option to spectral index + module to allow use of flux densities that do not meet the desired + SNR. Changed flag_maxsnr criterion to also flag if the peak flux + density per beam of the Gaussian exceeds the value at its center. + Removed incl_wavelet option. 2012/04/20 - Promoted the adaptive_rms_box parameter to the main options - listing and added the rms_box_bright option so that the user can - specify either (or both) of the rms_boxes. Fixed bug in wavelet - module so that invalid Gaussians (i.e., those that lie outside of - islands in the ch0 image) are not used when making the residual - images at each scale. Improved speed of Gaussian fitting to wavelet - images. Fixed bug that resulted in pixels found to be outside the - universe (check is enabled with the check_outsideuniv option) not - being masked properly. + listing and added the rms_box_bright option so that the user can + specify either (or both) of the rms_boxes. Fixed bug in wavelet + module so that invalid Gaussians (i.e., those that lie outside of + islands in the ch0 image) are not used when making the residual + images at each scale. Improved speed of Gaussian fitting to wavelet + images. Fixed bug that resulted in pixels found to be outside the + universe (check is enabled with the check_outsideuniv option) not + being masked properly. 2012/04/17 - Fixed bug in psf_vary module that resulted in PSF major and - minor axis maps in terms of sigma instead of FWHM. Added option - (psf_stype_only) to allow PSF fitting to non- S-type sources - (useful if sources are very distorted). + minor axis maps in terms of sigma instead of FWHM. Added option + (psf_stype_only) to allow PSF fitting to non- S-type sources + (useful if sources are very distorted). 2012/04/12 - Fixed bug in adaptive scaling code that could cause - incorrect small-scale rms_box size. Added a parameter - (adaptive_thresh) that controls the minimum threshold for sources - used to set the small-scale rms_box size. + incorrect small-scale rms_box size. Added a parameter + (adaptive_thresh) that controls the minimum threshold for sources + used to set the small-scale rms_box size. 2012/04/02 - Implemented an adaptive scaling scheme for the rms_box - parameter that shrinks the box size near bright sources and expands - it far from them (enabled with the adaptive_rms_box option when - rms_box=None). This scheme generally results in improved rms and - mean maps when both strong artifacts and extended sources are - present. Fixed bug that prevented plotting of results during wavelet - decomposition when interactive = True. + parameter that shrinks the box size near bright sources and expands + it far from them (enabled with the adaptive_rms_box option when + rms_box=None). This scheme generally results in improved rms and + mean maps when both strong artifacts and extended sources are + present. Fixed bug that prevented plotting of results during wavelet + decomposition when interactive = True. 2012/03/29 - Fixed bug in wavelet module that could cause incorrect - associations of Gaussians. Fixed bug in show_fit that displayed - incorrect model and residual images when wavelets were used. + associations of Gaussians. Fixed bug in show_fit that displayed + incorrect model and residual images when wavelets were used. 2012/03/28 - Version 1.1 2012/03/28 - Fixed bug that caused mask to be ignored when determining - whether variations in rms and mean maps is significant. Fixed bug - that caused internally derived rms_box value to be ignored. + whether variations in rms and mean maps is significant. Fixed bug + that caused internally derived rms_box value to be ignored. 2012/03/27 - Modified calculation of rms_box parameter (when rms_box - option is None) to work better with fields composed mainly of point - sources when strong artifacts are present. Tweaked flagging on FWHM - to prevent over-flagging of Gaussians in small islands. Changed - wavelet module to flag Gaussians whose centers fall outside of - islands found in the original image and removed atrous_orig_isl - option (as redundant). + option is None) to work better with fields composed mainly of point + sources when strong artifacts are present. Tweaked flagging on FWHM + to prevent over-flagging of Gaussians in small islands. Changed + wavelet module to flag Gaussians whose centers fall outside of + islands found in the original image and removed atrous_orig_isl + option (as redundant). 2012/03/26 - Modified fitting of large islands to adopt an iterative - fitting scheme that limits the number of Gaussians fit - simultaneously per iteration to 10. This change speeds up fitting of - large islands considerably. The options peak_fit and peak_maxsize - control whether iterative fitting is done. Added new Gaussian - flagging condition (flag_maxsize_fwhm) that flags Gaussians whose - sigma contour times factor extends beyond the island boundary. This - flag prevents fitting of Gaussians that extend far beyond the island - boundary. + fitting scheme that limits the number of Gaussians fit + simultaneously per iteration to 10. This change speeds up fitting of + large islands considerably. The options peak_fit and peak_maxsize + control whether iterative fitting is done. Added new Gaussian + flagging condition (flag_maxsize_fwhm) that flags Gaussians whose + sigma contour times factor extends beyond the island boundary. This + flag prevents fitting of Gaussians that extend far beyond the island + boundary. 2012/03/23 - Tweaked settings that affect fitting of Gaussians to - improve fitting in general. + improve fitting in general. 2012/03/19 - Added output of shapelet parameters to FITS tables. Fixed - issue with resizing of sources in spectral index module. + issue with resizing of sources in spectral index module. 2012/03/16 - Fixed bugs in polarisation module that caused incorrect - polarization fractions. + polarization fractions. 2012/03/13 - Improved plotting speed (by factor of ~ 4) in show_fit when - there is a large number of islands. Simplified the spectral index - module to make it more user friendly and stable. Added the option to - use a "detection" image for island detection (the detection_image - option); source properties are still measured from the main input - image. + there is a large number of islands. Simplified the spectral index + module to make it more user friendly and stable. Added the option to + use a "detection" image for island detection (the detection_image + option); source properties are still measured from the main input + image. 2012/03/01 - Fixed a bug in the polarisation module that could result in - incorrect flux densities. Changed logging module to suppress output - of ANSI color codes to the log file. + incorrect flux densities. Changed logging module to suppress output + of ANSI color codes to the log file. 2012/02/27 - Implemented fitting of Gaussians in polarisation module, - instead of simple summation of pixel values, to determine polarized - flux densities. + instead of simple summation of pixel values, to determine polarized + flux densities. 2012/02/17 - In scripts, process_image() will now accept a dictionary of - parameters as input. + parameters as input. 2012/02/10 - Sources that appear only in Stokes Q or U (and hence not in - Stokes I) are now identified and included in the polarisation - module. This identification is done using the polarized intensity - (PI) image. show_fit() and export_image() were updated to allow - display and export of the PI image. + Stokes I) are now identified and included in the polarisation + module. This identification is done using the polarized intensity + (PI) image. show_fit() and export_image() were updated to allow + display and export of the PI image. 2012/02/06 - Fixed bug in island splitting code that could result in - duplicate Gaussians. + duplicate Gaussians. 2012/02/02 - Improved polarisation module. Polarization quantities are - now calculated for Gaussians as well as sources. + now calculated for Gaussians as well as sources. 2012/01/27 - Fixed bug in psf_vary module that affected tesselation. - Fixed many small typos in parameter descriptions. + Fixed many small typos in parameter descriptions. 2012/01/18 - Fixed a bug that resulted in incorrect coordinates when the - trim_box option was used with a CASA image. Added option - (blank_zeros) to blank pixels in the input image that are exactly - zero. + trim_box option was used with a CASA image. Added option + (blank_zeros) to blank pixels in the input image that are exactly + zero. 2012/01/13 - Fixed minor bugs in the interactive shell and updated - pybdsm.py to support IPython 0.12. + pybdsm.py to support IPython 0.12. 2011/12/21 - Fixed bug in gaul2srl module due to rare cases in which an - island has a negative rms value. Fixed a memory issue in which - memory was not released after using show_fit. + island has a negative rms value. Fixed a memory issue in which + memory was not released after using show_fit. 2011/11/28 - Added option to have minpix_isl estimated automatically as - 1/3 of the beam area. This estimate should help exclude false - islands that are much smaller than the beam. This estimate is not - let to fall below 6 pixels. + 1/3 of the beam area. This estimate should help exclude false + islands that are much smaller than the beam. This estimate is not + let to fall below 6 pixels. 2011/11/11 - Fixed bugs in source generation that would lead to masking - of all pixels for certain sources during moment analysis. Adjusted - calculation of jmax in wavelet module to use island sizes (instead - of image size) if opts.atrous_orig_isl is True. + of all pixels for certain sources during moment analysis. Adjusted + calculation of jmax in wavelet module to use island sizes (instead + of image size) if opts.atrous_orig_isl is True. 2011/11/04 - Implemented new island fitting routine (enabled with the - peak_fit option) that can speed up fitting of large islands. Changed - plotting of Gaussians in show_fit to use Ellipse artists to improve - plotting speed. + peak_fit option) that can speed up fitting of large islands. Changed + plotting of Gaussians in show_fit to use Ellipse artists to improve + plotting speed. 2011/11/03 - Altered reading of images to correctly handle 4D cubes. - Fixed bug in readimage that affected filenames. + Fixed bug in readimage that affected filenames. 2011/10/26 - Extended psf_vary module to include fitting of stacked PSFs - with Gaussians, interpolation of the resulting parameters across the - image, and correction of the de- convolved source sizes using the - interpolated PSFs. Changed plotting of Gaussians in show_fit() to - use the FWHM instead of sigma. Modified error calculation of M - sources to be more robust when sources are small. Fixed spelling of - "gaussian" in bbs_patches option list. + with Gaussians, interpolation of the resulting parameters across the + image, and correction of the de- convolved source sizes using the + interpolated PSFs. Changed plotting of Gaussians in show_fit() to + use the FWHM instead of sigma. Modified error calculation of M + sources to be more robust when sources are small. Fixed spelling of + "gaussian" in bbs_patches option list. 2011/10/24 - Many small bug fixes to the psf_vary module. Fixed use of - input directory so that input files not in the current directory are - handled correctly. + input directory so that input files not in the current directory are + handled correctly. 2011/10/14 - Added residual rms and mean values to sources and source - list catalogs. These values can be compared to background rms and - mean values as a quick check of fit quality. + list catalogs. These values can be compared to background rms and + mean values as a quick check of fit quality. 2011/10/13 - Modified deconvolution to allow 1-D Gaussians and sources. - Added FREQ0, EQUINOX, INIMAGE keywords to output fits catalogs. - Fixed bug in source position angles. Adjusted column names of output - catalogs slightly to be more descriptive. + Added FREQ0, EQUINOX, INIMAGE keywords to output fits catalogs. + Fixed bug in source position angles. Adjusted column names of output + catalogs slightly to be more descriptive. 2011/10/12 - Added errors to source properties (using a Monte Carlo - method for M sources). Fixed bug in output column names. + method for M sources). Fixed bug in output column names. 2011/10/11 - Tweaked autocomplete to support IPython shell commands - (e.g., "!more file.txt"). Fixed bug in gaul2srl that resulted in - some very nearby Gaussians being placed into different sources. - Added group_tol option so that user can adjust the tolerance of how - Gaussians are grouped into sources. + (e.g., "!more file.txt"). Fixed bug in gaul2srl that resulted in + some very nearby Gaussians being placed into different sources. + Added group_tol option so that user can adjust the tolerance of how + Gaussians are grouped into sources. 2011/10/05 - Added output of source lists. Changed name of write_gaul - method to write_catalog (more general). + method to write_catalog (more general). 2011/10/04 - Added option to force source grouping by island - (group_by_isl). Added saving of parameters to a PyBDSM save file to - Op_output. + (group_by_isl). Added saving of parameters to a PyBDSM save file to + Op_output. 2011/09/21 - Fixed issue with shapelet centering failing: it now falls - back to simple moment when this happens. Fixed issue with - plotresults when shapelets are fit. + back to simple moment when this happens. Fixed issue with + plotresults when shapelets are fit. 2011/09/14 - Placed output column names and units in TC properties of - Gaussians. This allows easy standardization of the column names and - units. + Gaussians. This allows easy standardization of the column names and + units. 2011/09/13 - Fixes to trim_box and resetting of Image objects in - interface.process(). Changed thr1 --> thr2 in fit_iter in - guasfit.py, as bright sources are often "overfit" when using thr1, - leading to large negative residuals. Restricted fitting of Gaussians - to wavelet images to be only in islands found in the original image - if opts.atrous_orig_isl is True. + interface.process(). Changed thr1 --> thr2 in fit_iter in + guasfit.py, as bright sources are often "overfit" when using thr1, + leading to large negative residuals. Restricted fitting of Gaussians + to wavelet images to be only in islands found in the original image + if opts.atrous_orig_isl is True. 2011/09/08 - Version 1.0 diff --git a/bdsf/cleanup.py b/bdsf/cleanup.py index 63f0b590..e7cbb511 100644 --- a/bdsf/cleanup.py +++ b/bdsf/cleanup.py @@ -25,12 +25,12 @@ def __call__(self, img): pl.title('All gaussians including wavelet images') allgaus = img.gaussians if hasattr(img, 'atrous_gaussians'): - for gg in img.atrous_gaussians: - allgaus += gg + for gg in img.atrous_gaussians: + allgaus += gg for g in allgaus: - ellx, elly = func.drawellipse(g) - pl.plot(ellx, elly, 'r') + ellx, elly = func.drawellipse(g) + pl.plot(ellx, elly, 'r') from math import log10 bdir = img.basedir + '/misc/' @@ -50,6 +50,3 @@ def __call__(self, img): pl.close() img.completed_Ops.append('cleanup') - - - diff --git a/bdsf/collapse.py b/bdsf/collapse.py index 42647fa5..927886e5 100644 --- a/bdsf/collapse.py +++ b/bdsf/collapse.py @@ -18,176 +18,176 @@ class Op_collapse(Op): """Collapse 3D image""" def __call__(self, img): - mylog = mylogger.logging.getLogger("PyBDSM."+img.log+"Collapse") - if img.opts.polarisation_do: - pols = ['I', 'Q', 'U', 'V'] # make sure I is done first - else: - pols = ['I'] # assume I is always present - img.ch0_Q_arr = None - img.ch0_U_arr = None - img.ch0_V_arr = None - - if img.shape[1] > 1: - c_mode = img.opts.collapse_mode - chan0 = img.opts.collapse_ch0 - c_list = img.opts.collapse_av - c_wts = img.opts.collapse_wt - if not c_list: - c_list = N.arange(img.shape[1]) - if len(c_list) == 1 and c_mode=='average': - c_mode = 'single' - chan0 = c_list[0] - img.collapse_ch0 = chan0 - ch0sh = img.image_arr.shape[2:] + mylog = mylogger.logging.getLogger("PyBDSM."+img.log+"Collapse") if img.opts.polarisation_do: - ch0images = ['ch0_arr', 'ch0_Q_arr', 'ch0_U_arr', 'ch0_V_arr'] - else: - ch0images = ['ch0_arr'] - - # Check whether the collapse channel index is sensible - if chan0 < 0 or chan0 >= len(c_list): - raise RuntimeError('The channel index (set with the "collapse_ch0" option) ' - 'must be greater than zero and less than the number of ' - 'channels ({}).'.format(len(c_list))) - - # assume all Stokes images have the same blank pixels as I: - blank = N.isnan(img.image_arr[0]) - hasblanks = blank.any() - if img.opts.kappa_clip is None: - kappa = -img.pixel_beamarea() + pols = ['I', 'Q', 'U', 'V'] # make sure I is done first else: - kappa = img.opts.kappa_clip - - mean, rms, cmean, crms = chan_stats(img, kappa) - img.channel_mean = mean; img.channel_rms = rms - img.channel_clippedmean = cmean; img.channel_clippedrms = crms - - for ipol, pol in enumerate(pols): - if c_mode == 'single': - if pol == 'I': - ch0 = img.image_arr[0, chan0] - img.ch0_arr = ch0 - - # Construct weights so that desired channel has weight of 1 and all - # others have weight of 0. The init_freq_collapse function will then - # select the intended frequency - wtarr = N.zeros(len(c_list)) - wtarr[chan0] = 1.0 - init_freq_collapse(img, wtarr) - mylogger.userinfo(mylog, 'Source extraction will be ' \ - 'done on channel', '%i (%.3f MHz)' % \ - (chan0, img.frequency/1e6)) + pols = ['I'] # assume I is always present + img.ch0_Q_arr = None + img.ch0_U_arr = None + img.ch0_V_arr = None + + if img.shape[1] > 1: + c_mode = img.opts.collapse_mode + chan0 = img.opts.collapse_ch0 + c_list = img.opts.collapse_av + c_wts = img.opts.collapse_wt + if not c_list: + c_list = N.arange(img.shape[1]) + if len(c_list) == 1 and c_mode=='average': + c_mode = 'single' + chan0 = c_list[0] + img.collapse_ch0 = chan0 + ch0sh = img.image_arr.shape[2:] + if img.opts.polarisation_do: + ch0images = ['ch0_arr', 'ch0_Q_arr', 'ch0_U_arr', 'ch0_V_arr'] else: - ch0[:] = img.image_arr[ipol, chan0][:] - img.__setattr__(ch0images[ipol][:], ch0) - - elif c_mode == 'average': - if not hasblanks: - if pol == 'I': - ch0, wtarr = avspc_direct(c_list, img.image_arr[0], img.channel_clippedrms, c_wts) - else: - # use wtarr from the I image, which is always collapsed first - ch0, wtarr = avspc_direct(c_list, img.image_arr[ipol], img.channel_clippedrms, c_wts, wtarr=wtarr) + ch0images = ['ch0_arr'] + + # Check whether the collapse channel index is sensible + if chan0 < 0 or chan0 >= len(c_list): + raise RuntimeError('The channel index (set with the "collapse_ch0" option) ' + 'must be greater than zero and less than the number of ' + 'channels ({}).'.format(len(c_list))) + + # assume all Stokes images have the same blank pixels as I: + blank = N.isnan(img.image_arr[0]) + hasblanks = blank.any() + if img.opts.kappa_clip is None: + kappa = -img.pixel_beamarea() else: - if pol == 'I': - ch0, wtarr = avspc_blanks(c_list, img.image_arr[0], img.channel_clippedrms, c_wts) - else: - # use wtarr from the I image, which is always collapsed first - ch0, wtarr = avspc_blanks(c_list, img.image_arr[ipol], img.channel_clippedrms, c_wts, wtarr=wtarr) - img.__setattr__(ch0images[ipol][:], ch0) - - if pol == 'I': - img.avspc_wtarr = wtarr - init_freq_collapse(img, wtarr) - if c_wts == 'unity': - mylogger.userinfo(mylog, 'Channels averaged with '\ - 'uniform weights') - else: - mylogger.userinfo(mylog, 'Channels averaged with '\ - 'weights=(1/rms)^2') - mylogger.userinfo(mylog, 'Source extraction will be '\ - 'done on averaged ("ch0") image') - mylogger.userinfo(mylog, 'Frequency of averaged '\ - 'image', '%.3f MHz' % \ - (img.frequency/1e6,)) - str1 = " ".join(str(n) for n in c_list) - mylog.debug('%s %s' % ('Channels averaged : ', str1)) - str1 = " ".join(["%9.4e" % n for n in wtarr]) - mylog.debug('%s %s %s' % ('Channel weights : ', str1, '; unity=zero if c_wts="rms"')) - elif c_mode=='file': - mylogger.userinfo(mylog, 'Reading ch0 image from file %s' % (img.opts.collapse_file)) - image,hdr=func.read_image_from_file(img.opts.collapse_file, img, None, quiet=False) - if pol == 'I': - ch0 = image[0,0] - img.ch0_arr = ch0 - - else: - raise NotImplementedError('Polarization cubes not allowed in file mode') - else: - raise NotImplementedError('Mode supplied not implemented') # should never happen! - if img.opts.output_all: - func.write_image_to_file(img.use_io, img.imagename+'.ch0_'+pol+'.fits', ch0, - img, outdir=img.basedir) - mylog.debug('%s %s ' % ('Writing file ', img.imagename+'.ch0_'+pol+'.fits')) - - else: - # Only one channel in image - image = img.image_arr - img.ch0_arr = image[0, 0] - mylogger.userinfo(mylog, 'Frequency of image', - '%.3f MHz' % (img.frequency/1e6,)) - if img.opts.polarisation_do: - for pol in pols[1:]: - if pol == 'Q': - img.ch0_Q_arr = image[1, 0][:] - if pol == 'U': - img.ch0_U_arr = image[2, 0][:] - if pol == 'V': - img.ch0_V_arr = image[3, 0][:] - - # create mask if needed (assume all pols have the same mask as I) - image = img.ch0_arr - mask = N.isnan(image) - img.blankpix = N.sum(mask) - frac_blank = round( - float(img.blankpix) / float(image.shape[0] * image.shape[1]), - 3) - mylogger.userinfo(mylog, "Number of blank pixels", str(img.blankpix) - + ' (' + str(frac_blank * 100.0) + '%)') - - if img.opts.blank_limit is not None: - import scipy - import sys - threshold = img.opts.blank_limit - mylogger.userinfo(mylog, "Blanking pixels with values " - "below %.1e Jy/beam" % (threshold,)) - bad = (abs(image) < threshold) - original_stdout = sys.stdout # keep a reference to STDOUT - sys.stdout = func.NullDevice() # redirect the real STDOUT - count = scipy.signal.convolve2d(bad, N.ones((3, 3)), mode='same') - sys.stdout = original_stdout # turn STDOUT back on - mask_low = (count >= 5) - image[N.where(mask_low)] = N.nan - mask = N.isnan(image) - img.blankpix = N.sum(mask) - frac_blank = round( - float(img.blankpix) / float(image.shape[0] * - image.shape[1]), 3) - mylogger.userinfo(mylog, "Total number of blanked pixels", - str(img.blankpix) + ' (' + str(frac_blank * 100.0) + '%)') - - masked = mask.any() - img.masked = masked - if masked: - img.mask_arr = mask - else: - img.mask_arr = None - - if img.blankpix == image.shape[0] * image.shape[1]: - # ALL pixels are blanked! - raise RuntimeError('All pixels in the image are blanked.') - - img.completed_Ops.append('collapse') + kappa = img.opts.kappa_clip + + mean, rms, cmean, crms = chan_stats(img, kappa) + img.channel_mean = mean; img.channel_rms = rms + img.channel_clippedmean = cmean; img.channel_clippedrms = crms + + for ipol, pol in enumerate(pols): + if c_mode == 'single': + if pol == 'I': + ch0 = img.image_arr[0, chan0] + img.ch0_arr = ch0 + + # Construct weights so that desired channel has weight of 1 and all + # others have weight of 0. The init_freq_collapse function will then + # select the intended frequency + wtarr = N.zeros(len(c_list)) + wtarr[chan0] = 1.0 + init_freq_collapse(img, wtarr) + mylogger.userinfo(mylog, 'Source extraction will be ' \ + 'done on channel', '%i (%.3f MHz)' % \ + (chan0, img.frequency/1e6)) + else: + ch0[:] = img.image_arr[ipol, chan0][:] + img.__setattr__(ch0images[ipol][:], ch0) + + elif c_mode == 'average': + if not hasblanks: + if pol == 'I': + ch0, wtarr = avspc_direct(c_list, img.image_arr[0], img.channel_clippedrms, c_wts) + else: + # use wtarr from the I image, which is always collapsed first + ch0, wtarr = avspc_direct(c_list, img.image_arr[ipol], img.channel_clippedrms, c_wts, wtarr=wtarr) + else: + if pol == 'I': + ch0, wtarr = avspc_blanks(c_list, img.image_arr[0], img.channel_clippedrms, c_wts) + else: + # use wtarr from the I image, which is always collapsed first + ch0, wtarr = avspc_blanks(c_list, img.image_arr[ipol], img.channel_clippedrms, c_wts, wtarr=wtarr) + img.__setattr__(ch0images[ipol][:], ch0) + + if pol == 'I': + img.avspc_wtarr = wtarr + init_freq_collapse(img, wtarr) + if c_wts == 'unity': + mylogger.userinfo(mylog, 'Channels averaged with '\ + 'uniform weights') + else: + mylogger.userinfo(mylog, 'Channels averaged with '\ + 'weights=(1/rms)^2') + mylogger.userinfo(mylog, 'Source extraction will be '\ + 'done on averaged ("ch0") image') + mylogger.userinfo(mylog, 'Frequency of averaged '\ + 'image', '%.3f MHz' % \ + (img.frequency/1e6,)) + str1 = " ".join(str(n) for n in c_list) + mylog.debug('%s %s' % ('Channels averaged : ', str1)) + str1 = " ".join(["%9.4e" % n for n in wtarr]) + mylog.debug('%s %s %s' % ('Channel weights : ', str1, '; unity=zero if c_wts="rms"')) + elif c_mode=='file': + mylogger.userinfo(mylog, 'Reading ch0 image from file %s' % (img.opts.collapse_file)) + image,hdr=func.read_image_from_file(img.opts.collapse_file, img, None, quiet=False) + if pol == 'I': + ch0 = image[0,0] + img.ch0_arr = ch0 + + else: + raise NotImplementedError('Polarization cubes not allowed in file mode') + else: + raise NotImplementedError('Mode supplied not implemented') # should never happen! + if img.opts.output_all: + func.write_image_to_file(img.use_io, img.imagename+'.ch0_'+pol+'.fits', ch0, + img, outdir=img.basedir) + mylog.debug('%s %s ' % ('Writing file ', img.imagename+'.ch0_'+pol+'.fits')) + + else: + # Only one channel in image + image = img.image_arr + img.ch0_arr = image[0, 0] + mylogger.userinfo(mylog, 'Frequency of image', + '%.3f MHz' % (img.frequency/1e6,)) + if img.opts.polarisation_do: + for pol in pols[1:]: + if pol == 'Q': + img.ch0_Q_arr = image[1, 0][:] + if pol == 'U': + img.ch0_U_arr = image[2, 0][:] + if pol == 'V': + img.ch0_V_arr = image[3, 0][:] + + # create mask if needed (assume all pols have the same mask as I) + image = img.ch0_arr + mask = N.isnan(image) + img.blankpix = N.sum(mask) + frac_blank = round( + float(img.blankpix) / float(image.shape[0] * image.shape[1]), + 3) + mylogger.userinfo(mylog, "Number of blank pixels", str(img.blankpix) + + ' (' + str(frac_blank * 100.0) + '%)') + + if img.opts.blank_limit is not None: + import scipy + import sys + threshold = img.opts.blank_limit + mylogger.userinfo(mylog, "Blanking pixels with values " + "below %.1e Jy/beam" % (threshold,)) + bad = (abs(image) < threshold) + original_stdout = sys.stdout # keep a reference to STDOUT + sys.stdout = func.NullDevice() # redirect the real STDOUT + count = scipy.signal.convolve2d(bad, N.ones((3, 3)), mode='same') + sys.stdout = original_stdout # turn STDOUT back on + mask_low = (count >= 5) + image[N.where(mask_low)] = N.nan + mask = N.isnan(image) + img.blankpix = N.sum(mask) + frac_blank = round( + float(img.blankpix) / float(image.shape[0] * + image.shape[1]), 3) + mylogger.userinfo(mylog, "Total number of blanked pixels", + str(img.blankpix) + ' (' + str(frac_blank * 100.0) + '%)') + + masked = mask.any() + img.masked = masked + if masked: + img.mask_arr = mask + else: + img.mask_arr = None + + if img.blankpix == image.shape[0] * image.shape[1]: + # ALL pixels are blanked! + raise RuntimeError('All pixels in the image are blanked.') + + img.completed_Ops.append('collapse') ######################################################################################## @@ -198,23 +198,23 @@ def chan_stats(img, kappa): nchan = img.shape[1] mean = []; rms = []; cmean = []; crms = [] for ichan in range(nchan): - if isinstance(img, Image): # check if img is an Image or just an ndarray - im = img.image_arr[0, ichan] - else: - im = img[0, ichan] - - if N.any(im): - immask = N.isnan(im) - if immask.all(): - m, r, cm, cr = 0, 0, 0, 0 + if isinstance(img, Image): # check if img is an Image or just an ndarray + im = img.image_arr[0, ichan] + else: + im = img[0, ichan] + + if N.any(im): + immask = N.isnan(im) + if immask.all(): + m, r, cm, cr = 0, 0, 0, 0 + else: + if immask.any(): + m, r, cm, cr, cnt = bstat(im, immask, kappa) + else: + m, r, cm, cr, cnt = bstat(im, None, kappa) else: - if immask.any(): - m, r, cm, cr, cnt = bstat(im, immask, kappa) - else: - m, r, cm, cr, cnt = bstat(im, None, kappa) - else: - m, r, cm, cr = 0, 0, 0, 0 - mean.append(m); rms.append(r); cmean.append(cm); crms.append(cr) + m, r, cm, cr = 0, 0, 0, 0 + mean.append(m); rms.append(r); cmean.append(cm); crms.append(cr) return N.array(mean), N.array(rms), N.array(cmean), N.array(crms) @@ -227,24 +227,24 @@ def avspc_direct(c_list, image, rmsarr, c_wts, wtarr=None): ch0 = N.zeros(shape2, dtype=N.float32) sumwts = 0.0 if wtarr is None: - wtarr = N.zeros(len(c_list)) - for i, ch in enumerate(c_list): - im = image[ch] - r = rmsarr[ch] - if c_wts == 'unity': wt = 1.0 - if c_wts == 'rms': wt = r - if r != 0: - wt = 1.0/(wt*wt) - else: - wt = 0 - sumwts += wt - ch0 += im*wt - wtarr[i] = wt + wtarr = N.zeros(len(c_list)) + for i, ch in enumerate(c_list): + im = image[ch] + r = rmsarr[ch] + if c_wts == 'unity': wt = 1.0 + if c_wts == 'rms': wt = r + if r != 0: + wt = 1.0/(wt*wt) + else: + wt = 0 + sumwts += wt + ch0 += im*wt + wtarr[i] = wt else: - for i, ch in enumerate(c_list): - im = image[ch] - sumwts += wtarr[i] - ch0 += im*wtarr[i] + for i, ch in enumerate(c_list): + im = image[ch] + sumwts += wtarr[i] + ch0 += im*wtarr[i] ch0 = ch0/sumwts return ch0, wtarr @@ -257,27 +257,27 @@ def avspc_blanks(c_list, image, rmsarr, c_wts, wtarr=None): ch0 = N.zeros(shape2, dtype=N.float32) sumwtim = N.zeros(shape2, dtype=N.float32) if wtarr is None: - wtarr = N.zeros(len(c_list)) - for i, ch in enumerate(c_list): - im = image[ch] - r = rmsarr[ch] - if c_wts == 'unity': wt = 1.0 - if c_wts == 'rms': wt = r - if r > 1e-18 and r < 1e18: - # Set reasonable limits to avoid overflow of float32 - wt = 1.0/(wt*wt) - else: - wt = 0 - wtim = N.ones(shape2, dtype=N.float32)*wt*(~N.isnan(im)) - sumwtim += wtim - ch0 += N.nan_to_num(im)*wtim - wtarr[i] = wt + wtarr = N.zeros(len(c_list)) + for i, ch in enumerate(c_list): + im = image[ch] + r = rmsarr[ch] + if c_wts == 'unity': wt = 1.0 + if c_wts == 'rms': wt = r + if r > 1e-18 and r < 1e18: + # Set reasonable limits to avoid overflow of float32 + wt = 1.0/(wt*wt) + else: + wt = 0 + wtim = N.ones(shape2, dtype=N.float32)*wt*(~N.isnan(im)) + sumwtim += wtim + ch0 += N.nan_to_num(im)*wtim + wtarr[i] = wt else: - for i, ch in enumerate(c_list): - im = image[ch] - wtim = N.ones(shape2)*wtarr[i]*(~N.isnan(im)) - sumwtim += wtim - ch0 += N.nan_to_num(im)*wtim + for i, ch in enumerate(c_list): + im = image[ch] + wtim = N.ones(shape2)*wtarr[i]*(~N.isnan(im)) + sumwtim += wtim + ch0 += N.nan_to_num(im)*wtim ch0 = ch0/sumwtim return ch0, wtarr diff --git a/bdsf/const.py b/bdsf/const.py index 8ca3b458..b2510ff1 100644 --- a/bdsf/const.py +++ b/bdsf/const.py @@ -11,4 +11,3 @@ c=2.99792458e8 bolt=1.3806505e-23 sq2=math.sqrt(2) - diff --git a/bdsf/functions.py b/bdsf/functions.py index 39fb8b4a..212f4858 100644 --- a/bdsf/functions.py +++ b/bdsf/functions.py @@ -22,12 +22,12 @@ def sp_in(c, x): order = len(c)-1 if order == 1: - y = c[0]*N.power(x, c[1]) + y = c[0]*N.power(x, c[1]) else: - if order == 2: - y = c[0]*N.power(x, c[1])*N.power(x, c[2]*N.log(x)) - else: - print('Not yet implemented') + if order == 2: + y = c[0]*N.power(x, c[1])*N.power(x, c[2]*N.log(x)) + else: + print('Not yet implemented') return y @@ -45,9 +45,9 @@ def nanmean(x): n = N.sum(~N.isnan(x)) if n > 0: - mean = sum/n + mean = sum/n else: - mean = float("NaN") + mean = float("NaN") return mean @@ -223,15 +223,15 @@ def gaus_2d_itscomplicated(c, x, y, p_tofix, ind): val = N.zeros(x.shape) indx = N.array(ind) if len(indx) % 6 != 0: - print(" Something wrong with the parameters passed - need multiples of 6 !") + print(" Something wrong with the parameters passed - need multiples of 6 !") else: - ngaus = int(len(indx)/6) - params = N.zeros(6*ngaus) - params[N.where(indx==1)[0]] = c - params[N.where(indx==0)[0]] = p_tofix - for i in range(ngaus): - gau = params[i*6:i*6+6] - val = val + gaus_2d(gau, x, y) + ngaus = int(len(indx)/6) + params = N.zeros(6*ngaus) + params[N.where(indx==1)[0]] = c + params[N.where(indx==0)[0]] = p_tofix + for i in range(ngaus): + gau = params[i*6:i*6+6] + val = val + gaus_2d(gau, x, y) return val @@ -290,12 +290,12 @@ def drawellipse(g): rad = 180.0/N.pi if isinstance(g, Gaussian): - param = g2param(g) + param = g2param(g) else: - if isinstance(g, list) and len(g)>=6: - param = g - else: - raise RuntimeError("Input to drawellipse neither Gaussian nor list") + if isinstance(g, list) and len(g)>=6: + param = g + else: + raise RuntimeError("Input to drawellipse neither Gaussian nor list") size = [param[3], param[4], param[5]] size_fwhm = corrected_size(size) @@ -431,54 +431,54 @@ def fit_mask_1d(x, y, sig, mask, funct, do_err, order=0, p0 = None): ind=N.where(~N.array(mask))[0] if len(ind) > 1: - n=sum(mask) - - if isinstance(x, list): x = N.array(x) - if isinstance(y, list): y = N.array(y) - if isinstance(sig, list): sig = N.array(sig) - xfit=x[ind]; yfit=y[ind]; sigfit=sig[ind] - - if p0 is None: - if funct == poly: - p0=N.array([0]*(order+1)) - p0[1]=(yfit[0]-yfit[-1])/(xfit[0]-xfit[-1]) - p0[0]=yfit[0]-p0[1]*xfit[0] - if funct == wenss_fit: - p0=N.array([yfit[N.argmax(xfit)]] + [1.]) - if funct == sp_in: - ind1 = N.where(yfit > 0.)[0] - if len(ind1) >= 2: - low = ind1[0]; hi = ind1[-1] - sp = N.log(yfit[low]/yfit[hi])/N.log(xfit[low]/xfit[hi]) - p0=N.array([yfit[low]/pow(xfit[low], sp), sp] + [0.]*(order-1)) - elif len(ind1) == 1: - p0=N.array([ind1[0], -0.8] + [0.]*(order-1)) - else: - return [0, 0], [0, 0] - res=lambda p, xfit, yfit, sigfit: (yfit-funct(p, xfit))/sigfit - try: - (p, cov, info, mesg, flag)=leastsq(res, p0, args=(xfit, yfit, sigfit), full_output=True, warning=False) - except TypeError: - # This error means no warning argument is available, so redirect stdout to a null device - # to suppress printing of (unnecessary) warning messages - original_stdout = sys.stdout # keep a reference to STDOUT - sys.stdout = NullDevice() # redirect the real STDOUT - (p, cov, info, mesg, flag)=leastsq(res, p0, args=(xfit, yfit, sigfit), full_output=True) - sys.stdout = original_stdout # turn STDOUT back on - - if do_err: - 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: - chisq=sum(info["fvec"]*info["fvec"]) - dof=len(info["fvec"])-len(p) - err = N.array([sqrt(abs(cov[i,i])*chisq/dof) for i in range(len(p))]) - else: - p, err = [0, 0], [0, 0] - else: err = [0] + n=sum(mask) + + if isinstance(x, list): x = N.array(x) + if isinstance(y, list): y = N.array(y) + if isinstance(sig, list): sig = N.array(sig) + xfit=x[ind]; yfit=y[ind]; sigfit=sig[ind] + + if p0 is None: + if funct == poly: + p0=N.array([0]*(order+1)) + p0[1]=(yfit[0]-yfit[-1])/(xfit[0]-xfit[-1]) + p0[0]=yfit[0]-p0[1]*xfit[0] + if funct == wenss_fit: + p0=N.array([yfit[N.argmax(xfit)]] + [1.]) + if funct == sp_in: + ind1 = N.where(yfit > 0.)[0] + if len(ind1) >= 2: + low = ind1[0]; hi = ind1[-1] + sp = N.log(yfit[low]/yfit[hi])/N.log(xfit[low]/xfit[hi]) + p0=N.array([yfit[low]/pow(xfit[low], sp), sp] + [0.]*(order-1)) + elif len(ind1) == 1: + p0=N.array([ind1[0], -0.8] + [0.]*(order-1)) + else: + return [0, 0], [0, 0] + res=lambda p, xfit, yfit, sigfit: (yfit-funct(p, xfit))/sigfit + try: + (p, cov, info, mesg, flag)=leastsq(res, p0, args=(xfit, yfit, sigfit), full_output=True, warning=False) + except TypeError: + # This error means no warning argument is available, so redirect stdout to a null device + # to suppress printing of (unnecessary) warning messages + original_stdout = sys.stdout # keep a reference to STDOUT + sys.stdout = NullDevice() # redirect the real STDOUT + (p, cov, info, mesg, flag)=leastsq(res, p0, args=(xfit, yfit, sigfit), full_output=True) + sys.stdout = original_stdout # turn STDOUT back on + + if do_err: + 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: + chisq=sum(info["fvec"]*info["fvec"]) + dof=len(info["fvec"])-len(p) + err = N.array([sqrt(abs(cov[i,i])*chisq/dof) for i in range(len(p))]) + else: + p, err = [0, 0], [0, 0] + else: err = [0] else: - p, err = [0, 0], [0, 0] + p, err = [0, 0], [0, 0] return p, err @@ -580,17 +580,17 @@ def momanalmask_gaus(subim, mask, isrc, bmar_p, allpara=True): mompara[1:3] = m1/tot if allpara: - for coord in index: - co = N.array(coord) - m2 += (co - mompara[1:3])*(co - mompara[1:3])*subim[coord] - m11 += N.product(co - mompara[1:3])*subim[coord] - - mompara[3] = sqrt((m2[0]+m2[1]+sqrt((m2[0]-m2[1])*(m2[0]-m2[1])+4.0*m11*m11))/(2.0*tot))*fwsig - mompara[4] = sqrt((m2[0]+m2[1]-sqrt((m2[0]-m2[1])*(m2[0]-m2[1])+4.0*m11*m11))/(2.0*tot))*fwsig - dumr = atan(abs(2.0*m11/(m2[0]-m2[1]))) - dumr = atanproper(dumr, m2[0]-m2[1], 2.0*m11) - mompara[5] = 0.5*dumr*180.0/pi - 90.0 - if mompara[5] < 0.0: mompara[5] += 180.0 + for coord in index: + co = N.array(coord) + m2 += (co - mompara[1:3])*(co - mompara[1:3])*subim[coord] + m11 += N.product(co - mompara[1:3])*subim[coord] + + mompara[3] = sqrt((m2[0]+m2[1]+sqrt((m2[0]-m2[1])*(m2[0]-m2[1])+4.0*m11*m11))/(2.0*tot))*fwsig + mompara[4] = sqrt((m2[0]+m2[1]-sqrt((m2[0]-m2[1])*(m2[0]-m2[1])+4.0*m11*m11))/(2.0*tot))*fwsig + dumr = atan(abs(2.0*m11/(m2[0]-m2[1]))) + dumr = atanproper(dumr, m2[0]-m2[1], 2.0*m11) + mompara[5] = 0.5*dumr*180.0/pi - 90.0 + if mompara[5] < 0.0: mompara[5] += 180.0 return mompara def fit_gaus2d(data, p_ini, x, y, mask = None, err = None): @@ -646,11 +646,11 @@ def deconv(gaus_bm, gaus_c): rhoc = (maj2_c-min2_c)*cost-(maj2_bm-min2_bm) if rhoc == 0.0: - sigic = 0.0 - rhoa = 0.0 + sigic = 0.0 + rhoa = 0.0 else: - sigic = atan((maj2_c-min2_c)*sint/rhoc) # in radians - rhoa = ((maj2_bm-min2_bm)-(maj2_c-min2_c)*cost)/(2.0*cos(sigic)) + sigic = atan((maj2_c-min2_c)*sint/rhoc) # in radians + rhoa = ((maj2_bm-min2_bm)-(maj2_c-min2_c)*cost)/(2.0*cos(sigic)) gaus_d[2] = sigic*rad/2.0+phi_bm dumr = ((maj2_c+min2_c)-(maj2_bm+min2_bm))/2.0 @@ -665,21 +665,21 @@ def deconv(gaus_bm, gaus_c): gaus_d[0] = sqrt(abs(gaus_d[0])) gaus_d[1] = sqrt(abs(gaus_d[1])) if gaus_d[0] < gaus_d[1]: - sint = gaus_d[0] - gaus_d[0] = gaus_d[1] - gaus_d[1] = sint - gaus_d[2] = gaus_d[2]+90.0 + sint = gaus_d[0] + gaus_d[0] = gaus_d[1] + gaus_d[1] = sint + gaus_d[2] = gaus_d[2]+90.0 gaus_d[2] = gaus_d[2]+900.0 % 180 if gaus_d[0] == 0.0: - gaus_d[2] = 0.0 + gaus_d[2] = 0.0 else: - if gaus_d[1] == 0.0: - if (abs(gaus_d[2]-phi_c) > 45.0) and (abs(gaus_d[2]-phi_c) < 135.0): - gaus_d[2] = gaus_d[2]+450.0 % 180 + if gaus_d[1] == 0.0: + if (abs(gaus_d[2]-phi_c) > 45.0) and (abs(gaus_d[2]-phi_c) < 135.0): + gaus_d[2] = gaus_d[2]+450.0 % 180 # errors - #if rhoc == 0.0: + #if rhoc == 0.0: #if gaus_d[0] != 0.0: # ed_1 = gaus_c[0]/gaus_d[0]*e_1 #else: @@ -784,65 +784,65 @@ def get_errors(img, p, stdav, bm_pix=None, fixed_to_beam=False): mylog = mylogger.logging.getLogger("PyBDSM.Compute") if len(p) % 7 > 0: - mylog.error("Gaussian parameters passed have to have 7n numbers") + mylog.error("Gaussian parameters passed have to have 7n numbers") ngaus = int(len(p)/7) errors = [] for i in range(ngaus): - pp = p[i*7:i*7+7] - ### Now do error analysis as in Condon (and fBDSM) - size = pp[3:6] - size = corrected_size(size) # angle is now degrees CCW from +y-axis - if size[0] == 0.0 or size[1] == 0.0: - errors = errors + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - else: - sq2 = sqrt(2.0) - 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]))) # first term of Eq. 26 of Condon+ (1998) - dumrr1 = 1.0 + bm_pix[0] * bm_pix[1] / (size[0] * size[0]) # second term of Eq. 26 of Condon+ (1998) - dumrr2 = 1.0 + bm_pix[0] * bm_pix[1] / (size[1] * size[1]) # third term of Eq. 26 of Condon+ (1998) - dumrr3 = dumr * pp[0] / stdav # product of first and fourth terms of Eq. 26 of Condon+ (1998) - d1 = sqrt(8.0 * log(2.0)) - d2 = (size[0] * size[0] - size[1] * size[1]) / (size[0] * size[1]) # last term of Eq. 30 of Condon+ (1998) - try: - # The following three errors are calculated using Eq. 21 of Condon (1997), - # using Eq. 26 of Condon+ (1998) for rho - e_peak = pp[0] * sq2 / (dumrr3 * pow(dumrr1, 0.75) * pow(dumrr2, 0.75)) - e_maj = size[0] * sq2 / (dumrr3 * pow(dumrr1, 1.25) * pow(dumrr2, 0.25)) - e_min = size[1] * sq2 / (dumrr3 * pow(dumrr1, 0.25) * pow(dumrr2, 1.25)) - - # The following two errors are calculated using Eq. 27 of Condon+ (1998) - pa_rad = size[2] * pi / 180.0 - e_x0 = sqrt( (e_maj * N.sin(pa_rad))**2 + (e_min * N.cos(pa_rad))**2 ) / d1 - e_y0 = sqrt( (e_maj * N.cos(pa_rad))**2 + (e_min * N.sin(pa_rad))**2 ) / d1 - - # The following error is calculated using Eq. 30 of Condon+ (1998) - e_pa = 2.0 / (d2 * dumrr3 * pow(dumrr1, 0.25) * pow(dumrr2, 1.25)) - e_pa = e_pa * 180.0/pi - - # The following error is calculated using Eq. 36 of Condon+ (1998) - e_tot = pp[6] * sqrt(e_peak * e_peak / (pp[0] * pp[0]) + (0.25 / dumr / dumr) * - (e_maj * e_maj / (size[0] * size[0]) + e_min * e_min / (size[1] * size[1]))) - except: - e_peak = 0.0 - e_x0 = 0.0 - e_y0 = 0.0 - e_maj = 0.0 - e_min = 0.0 - e_pa = 0.0 - e_tot = 0.0 - if abs(e_pa) > 180.0: - e_pa = 180.0 - if fixed_to_beam: - # When the size was fixed to that of the beam during the fit, set - # uncertainties on the size to zero and reduce the error in the fluxes - # by sqrt(2) (see Eq. 25 of Condon 1997) - e_maj = 0.0 - e_min = 0.0 - e_pa = 0.0 - e_peak /= sq2 - e_tot /= sq2 - errors = errors + [e_peak, e_x0, e_y0, e_maj, e_min, e_pa, e_tot] + pp = p[i*7:i*7+7] + ### Now do error analysis as in Condon (and fBDSM) + size = pp[3:6] + size = corrected_size(size) # angle is now degrees CCW from +y-axis + if size[0] == 0.0 or size[1] == 0.0: + errors = errors + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + else: + sq2 = sqrt(2.0) + 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]))) # first term of Eq. 26 of Condon+ (1998) + dumrr1 = 1.0 + bm_pix[0] * bm_pix[1] / (size[0] * size[0]) # second term of Eq. 26 of Condon+ (1998) + dumrr2 = 1.0 + bm_pix[0] * bm_pix[1] / (size[1] * size[1]) # third term of Eq. 26 of Condon+ (1998) + dumrr3 = dumr * pp[0] / stdav # product of first and fourth terms of Eq. 26 of Condon+ (1998) + d1 = sqrt(8.0 * log(2.0)) + d2 = (size[0] * size[0] - size[1] * size[1]) / (size[0] * size[1]) # last term of Eq. 30 of Condon+ (1998) + try: + # The following three errors are calculated using Eq. 21 of Condon (1997), + # using Eq. 26 of Condon+ (1998) for rho + e_peak = pp[0] * sq2 / (dumrr3 * pow(dumrr1, 0.75) * pow(dumrr2, 0.75)) + e_maj = size[0] * sq2 / (dumrr3 * pow(dumrr1, 1.25) * pow(dumrr2, 0.25)) + e_min = size[1] * sq2 / (dumrr3 * pow(dumrr1, 0.25) * pow(dumrr2, 1.25)) + + # The following two errors are calculated using Eq. 27 of Condon+ (1998) + pa_rad = size[2] * pi / 180.0 + e_x0 = sqrt( (e_maj * N.sin(pa_rad))**2 + (e_min * N.cos(pa_rad))**2 ) / d1 + e_y0 = sqrt( (e_maj * N.cos(pa_rad))**2 + (e_min * N.sin(pa_rad))**2 ) / d1 + + # The following error is calculated using Eq. 30 of Condon+ (1998) + e_pa = 2.0 / (d2 * dumrr3 * pow(dumrr1, 0.25) * pow(dumrr2, 1.25)) + e_pa = e_pa * 180.0/pi + + # The following error is calculated using Eq. 36 of Condon+ (1998) + e_tot = pp[6] * sqrt(e_peak * e_peak / (pp[0] * pp[0]) + (0.25 / dumr / dumr) * + (e_maj * e_maj / (size[0] * size[0]) + e_min * e_min / (size[1] * size[1]))) + except: + e_peak = 0.0 + e_x0 = 0.0 + e_y0 = 0.0 + e_maj = 0.0 + e_min = 0.0 + e_pa = 0.0 + e_tot = 0.0 + if abs(e_pa) > 180.0: + e_pa = 180.0 + if fixed_to_beam: + # When the size was fixed to that of the beam during the fit, set + # uncertainties on the size to zero and reduce the error in the fluxes + # by sqrt(2) (see Eq. 25 of Condon 1997) + e_maj = 0.0 + e_min = 0.0 + e_pa = 0.0 + e_peak /= sq2 + e_tot /= sq2 + errors = errors + [e_peak, e_x0, e_y0, e_maj, e_min, e_pa, e_tot] return errors @@ -851,10 +851,10 @@ def fit_chisq(x, p, ep, mask, funct, order): ind = N.where(N.array(mask)==False)[0] if order == 0: - fit = [funct(p)]*len(p) + fit = [funct(p)]*len(p) else: - fitpara, efit = fit_mask_1d(x, p, ep, mask, funct, True, order) - fit = funct(fitpara, x) + fitpara, efit = fit_mask_1d(x, p, ep, mask, funct, True, order) + fit = funct(fitpara, x) dev = (p-fit)*(p-fit)/(ep*ep) num = order+1 @@ -866,9 +866,9 @@ def calc_chisq(x, y, ey, p, mask, funct, order): import numpy as N if order == 0: - fit = [funct(y)]*len(y) + fit = [funct(y)]*len(y) else: - fit = funct(p, x) + fit = funct(p, x) dev = (y-fit)*(y-fit)/(ey*ey) ind = N.where(~N.array(mask)) @@ -883,11 +883,11 @@ def get_windowsize_av(S_i, rms_i, chanmask, K, minchan): av_window = N.arange(2, int(len(S_i)/minchan)+1) win_size = 0 for window in av_window: - fluxes, vars, mask = variance_of_wted_windowedmean(S_i, rms_i, chanmask, window) - minsnr = N.min(fluxes[~mask]/vars[~mask]) - if minsnr > K*1.1: ### K*1.1 since fitted peak can be less than wted peak - win_size = window # is the size of averaging window - break + fluxes, vars, mask = variance_of_wted_windowedmean(S_i, rms_i, chanmask, window) + minsnr = N.min(fluxes[~mask]/vars[~mask]) + if minsnr > K*1.1: ### K*1.1 since fitted peak can be less than wted peak + win_size = window # is the size of averaging window + break return win_size @@ -901,20 +901,20 @@ def variance_of_wted_windowedmean(S_i, rms_i, chanmask, window_size): wt = wt/N.median(wt) fluxes = N.zeros(nwin); vars = N.zeros(nwin); mask = N.zeros(nwin, bool) for i in range(nwin): - strt = i*window_size; stp = (i+1)*window_size - if i == nwin-1: stp = nchan - ind = N.arange(strt,stp) - m = chanmask[ind] - index = [arg for ii,arg in enumerate(ind) if not m[ii]] - if len(index) > 0: - s = S_i[index]; r = rms_i[index]; w = wt[index] - fluxes[i] = N.sum(s*w)/N.sum(w) - vars[i] = 1.0/sqrt(N.sum(1.0/r/r)) - mask[i] = N.product(m) - else: - fluxes[i] = 0 - vars[i] = 0 - mask[i] = True + strt = i*window_size; stp = (i+1)*window_size + if i == nwin-1: stp = nchan + ind = N.arange(strt,stp) + m = chanmask[ind] + index = [arg for ii,arg in enumerate(ind) if not m[ii]] + if len(index) > 0: + s = S_i[index]; r = rms_i[index]; w = wt[index] + fluxes[i] = N.sum(s*w)/N.sum(w) + vars[i] = 1.0/sqrt(N.sum(1.0/r/r)) + mask[i] = N.product(m) + else: + fluxes[i] = 0 + vars[i] = 0 + mask[i] = True return fluxes, vars, mask @@ -936,43 +936,43 @@ def fit_mulgaus2d(image, gaus, x, y, mask = None, fitfix = None, err = None, adj ngaus = len(gaus) if ngaus > 0: - p_ini = [] - for g in gaus: - p_ini = p_ini + g2param(g, adj) - p_ini = N.array(p_ini) - - 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 - if fitfix[i] == 2: ind[i*6+3:i*6+6] = 0 - if fitfix[i] == 3: ind[i*6+1:i*6+3] = 0 - ind = N.array(ind) - p_tofit = p_ini[N.where(ind==1)[0]] - p_tofix = p_ini[N.where(ind==0)[0]] - 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] - try: - p, success = leastsq(errorfunction, p_tofit, args=(x, y, p_tofix, ind, image, err, g_ind)) - except TypeError: - # This error means no warning argument is available, so redirect stdout to a null device - # to suppress printing of warning messages - original_stdout = sys.stdout # keep a reference to STDOUT - sys.stdout = NullDevice() # redirect the real STDOUT - p, success = leastsq(errorfunction, p_tofit, args=(x, y, p_tofix, ind, image, err, g_ind)) - sys.stdout = original_stdout # turn STDOUT back on + p_ini = [] + for g in gaus: + p_ini = p_ini + g2param(g, adj) + p_ini = N.array(p_ini) + + 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 + if fitfix[i] == 2: ind[i*6+3:i*6+6] = 0 + if fitfix[i] == 3: ind[i*6+1:i*6+3] = 0 + ind = N.array(ind) + p_tofit = p_ini[N.where(ind==1)[0]] + p_tofix = p_ini[N.where(ind==0)[0]] + 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] + try: + p, success = leastsq(errorfunction, p_tofit, args=(x, y, p_tofix, ind, image, err, g_ind)) + except TypeError: + # This error means no warning argument is available, so redirect stdout to a null device + # to suppress printing of warning messages + original_stdout = sys.stdout # keep a reference to STDOUT + sys.stdout = NullDevice() # redirect the real STDOUT + p, success = leastsq(errorfunction, p_tofit, args=(x, y, p_tofix, ind, image, err, g_ind)) + sys.stdout = original_stdout # turn STDOUT back on else: - p, sucess = None, 1 + p, sucess = None, 1 para = N.zeros(6*ngaus) para[N.where(ind==1)[0]] = p para[N.where(ind==0)[0]] = p_tofix for igaus in range(ngaus): - para[igaus*6+3] = abs(para[igaus*6+3]) - para[igaus*6+4] = abs(para[igaus*6+4]) + para[igaus*6+3] = abs(para[igaus*6+3]) + para[igaus*6+4] = abs(para[igaus*6+4]) return para, success @@ -1050,45 +1050,45 @@ def get_maxima(im, mask, thr, shape, beam, im_pos=None): iniposn = [] inipeak = [] for c in ind: - goodlist = [im_pos[i,j] for i in range(c[0]-1,c[0]+2) for j in range(c[1]-1,c[1]+2) \ - if i>=0 and i=0 and j goodlist) == len(goodlist) - if peak: - iniposn.append(c) - inipeak.append(im[c]) - im1 = mclean(im1, c, beam) + goodlist = [im_pos[i,j] for i in range(c[0]-1,c[0]+2) for j in range(c[1]-1,c[1]+2) \ + if i>=0 and i=0 and j goodlist) == len(goodlist) + if peak: + iniposn.append(c) + inipeak.append(im[c]) + im1 = mclean(im1, c, beam) return inipeak, iniposn, im1 def watershed(image, mask=None, markers=None, beam=None, thr=None): - import numpy as N - from copy import deepcopy as cp - import scipy.ndimage as nd - #import matplotlib.pyplot as pl - #import pylab as pl - - if thr is None: thr = -1e9 - if mask is None: mask = N.zeros(image.shape, bool) - if beam is None: beam = (2.0, 2.0, 0.0) - if markers is None: + import numpy as N + from copy import deepcopy as cp + import scipy.ndimage as nd + #import matplotlib.pyplot as pl + #import pylab as pl + + if thr is None: thr = -1e9 + if mask is None: mask = N.zeros(image.shape, bool) + if beam is None: beam = (2.0, 2.0, 0.0) + if markers is None: inipeak, iniposn, im1 = get_maxima(image, mask, thr, image.shape, beam) ng = len(iniposn); markers = N.zeros(image.shape, int) for i in range(ng): markers[iniposn[i]] = i+2 markers[N.unravel_index(N.argmin(image), image.shape)] = 1 - im1 = cp(image) - if im1.min() < 0.: im1 = im1-im1.min() - im1 = 255 - im1/im1.max()*255 - opw = nd.watershed_ift(N.array(im1, N.uint16), markers) + im1 = cp(image) + if im1.min() < 0.: im1 = im1-im1.min() + im1 = 255 - im1/im1.max()*255 + opw = nd.watershed_ift(N.array(im1, N.uint16), markers) - return opw, markers + return opw, markers def get_kwargs(kwargs, key, typ, default): obj = True if key in kwargs: - obj = kwargs[key] + obj = kwargs[key] if not isinstance(obj, typ): - obj = default + obj = default return obj @@ -1682,9 +1682,9 @@ def connect(mask): connectivity = nd.generate_binary_structure(2,2) labels, count = nd.label(mask, connectivity) if count > 1 : - connected = 'multiple' + connected = 'multiple' else: - connected = 'single' + connected = 'single' return connected, count @@ -1699,9 +1699,9 @@ def area_polygon(points): area = 0.0 for i in range(n_tri): - p1, p2, p3 = N.array([cenx, ceny]), N.array([x[i], y[i]]), N.array([x[i+1], y[i+1]]) - t_area= N.linalg.norm(N.cross((p2 - p1), (p3 - p1)))/2. - area += t_area + p1, p2, p3 = N.array([cenx, ceny]), N.array([x[i], y[i]]), N.array([x[i+1], y[i+1]]) + t_area= N.linalg.norm(N.cross((p2 - p1), (p3 - p1)))/2. + area += t_area return area @@ -1803,12 +1803,12 @@ def check_1pixcontacts(open): connectivity = nd.generate_binary_structure(2,2) ind = N.transpose(N.where(open[1:-1,1:-1] > 0)) + [1,1] # exclude boundary to make it easier for pixel in ind: - x, y = pixel - grid = cp(open[x-1:x+2, y-1:y+2]); grid[1,1] = 0 - grid = N.where(grid == open[tuple(pixel)], 1, 0) - ll, nn = nd.label(grid, connectivity) - if nn > 1: - open[tuple(pixel)] = 0 + x, y = pixel + grid = cp(open[x-1:x+2, y-1:y+2]); grid[1,1] = 0 + grid = N.where(grid == open[tuple(pixel)], 1, 0) + ll, nn = nd.label(grid, connectivity) + if nn > 1: + open[tuple(pixel)] = 0 return open @@ -1831,29 +1831,29 @@ def assign_leftovers(mask, open, nisl, labels): npix = [len(N.where(labels==b)[0]) for b in range(1,nisl+1)] for i_subisl in range(count): - c_list = [] # is list of all bordering pixels of the sub island - ii = i_subisl+1 - coords = N.transpose(N.where(mlabels==ii)) # the coordinates of island i of left-out pixels - for co in coords: - co8 = [[x,y] for x in range(co[0]-1,co[0]+2) for y in range(co[1]-1,co[1]+2) if x >=0 and y >=0 and x mask pixels - for cc in coords: - mask = (mlabels == ii) + c_list = [] # is list of all bordering pixels of the sub island + ii = i_subisl+1 + coords = N.transpose(N.where(mlabels==ii)) # the coordinates of island i of left-out pixels + for co in coords: + co8 = [[x,y] for x in range(co[0]-1,co[0]+2) for y in range(co[1]-1,co[1]+2) if x >=0 and y >=0 and x mask pixels + for cc in coords: + mask = (mlabels == ii) # mask[cc] = True - return None, mask - if len(belongs) == 1: - for cc in coords: - labels[tuple(cc)] = belongs[0] - else: # get the border pixels of the islands - nn = [npix[b-1] for b in belongs] - addto = belongs[N.argmin(nn)] - for cc in coords: - labels[tuple(cc)] = addto + return None, mask + if len(belongs) == 1: + for cc in coords: + labels[tuple(cc)] = belongs[0] + else: # get the border pixels of the islands + nn = [npix[b-1] for b in belongs] + addto = belongs[N.argmin(nn)] + for cc in coords: + labels[tuple(cc)] = addto return labels, mask @@ -1931,17 +1931,17 @@ def isl_tosplit(isl, opts): # take open 3 or 5 open3, open5 = False, False 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 + 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 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 + 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 # index=0 => dont split if open5: index = 5; n_subisl = n_subisl5; labels = labels5 else: - if open3: index = 3; n_subisl = n_subisl3; labels = labels3 - else: index = 0 + if open3: index = 3; n_subisl = n_subisl3; labels = labels3 + else: index = 0 convex_def = convexhull_deficiency(isl) #print 'CONVEX = ',convex_def diff --git a/bdsf/gaul2srl.py b/bdsf/gaul2srl.py index d4f35138..a9dd1aa0 100644 --- a/bdsf/gaul2srl.py +++ b/bdsf/gaul2srl.py @@ -54,14 +54,14 @@ def __call__(self, img): g_list.append(g) if len(g_list) > 0: - if len(g_list) == 1: - src_index, source = self.process_single_gaussian(img, g_list, src_index, code = 'S') - sources.append(source) - isl_sources.append(source) - else: - src_index, source = self.process_CM(img, g_list, isl, src_index) - sources.extend(source) - isl_sources.extend(source) + if len(g_list) == 1: + src_index, source = self.process_single_gaussian(img, g_list, src_index, code = 'S') + sources.append(source) + isl_sources.append(source) + else: + src_index, source = self.process_CM(img, g_list, isl, src_index) + sources.extend(source) + isl_sources.extend(source) else: if not img.waveletimage: dg = isl.dgaul[0] @@ -179,36 +179,36 @@ def process_CM(self, img, g_list, isl, src_index): for pair in index: same_island = self.in_same_island(pair, img, g_list, isl, subim, subn, subm, delc) if same_island: - nsrc -= 1 - mmax, mmin = max(src_id[pair[0]],src_id[pair[1]]), min(src_id[pair[0]],src_id[pair[1]]) - arr = N.where(src_id == mmax)[0]; src_id[arr] = mmin - # now reorder src_id so that it is contiguous + nsrc -= 1 + mmax, mmin = max(src_id[pair[0]],src_id[pair[1]]), min(src_id[pair[0]],src_id[pair[1]]) + arr = N.where(src_id == mmax)[0]; src_id[arr] = mmin + # now reorder src_id so that it is contiguous for i in range(ngau): ind1 = N.where(src_id==i)[0] if len(ind1) == 0: arr = N.where(src_id > i)[0] if len(arr) > 0: - decr = N.min(src_id[arr])-i - for j in arr: src_id[j] -= decr + decr = N.min(src_id[arr])-i + for j in arr: src_id[j] -= decr nsrc = N.max(src_id)+1 # now do whats in sub_calc_para_source source_list = [] for isrc in range(nsrc): - posn = N.where(src_id == isrc)[0] - g_sublist=[] - for i in posn: - g_sublist.append(g_list[i]) - ngau_insrc = len(posn) - # Do source type C - if ngau_insrc == 1: - src_index, source = self.process_single_gaussian(img, g_sublist, src_index, code = 'C') - else: - # make mask and subim. Invalid mask value is -1 since 0 is valid srcid - mask = self.make_mask(isl, subn, subm, 1, isrc, g_sublist, delc) - src_index, source = self.process_Multiple(img, g_sublist, mask, src_index, isrc, subim, \ - isl, delc, subn, subm) - source_list.append(source) + posn = N.where(src_id == isrc)[0] + g_sublist=[] + for i in posn: + g_sublist.append(g_list[i]) + ngau_insrc = len(posn) + # Do source type C + if ngau_insrc == 1: + src_index, source = self.process_single_gaussian(img, g_sublist, src_index, code = 'C') + else: + # make mask and subim. Invalid mask value is -1 since 0 is valid srcid + mask = self.make_mask(isl, subn, subm, 1, isrc, g_sublist, delc) + src_index, source = self.process_Multiple(img, g_sublist, mask, src_index, isrc, subim, \ + isl, delc, subn, subm) + source_list.append(source) return src_index, source_list @@ -243,35 +243,35 @@ def same_island_min(pair, g_list, subim, delc, tol=0.5): same_island_min = False same_island_cont = False if maxline == 1: - same_island_min = True - same_island_cont = True - else: - if abs(pixdif[0]) > abs(pixdif[1]): - xline = N.round(min(pix1[0],pix2[0])+N.arange(maxline)) - yline = N.round((pix1[1]-pix2[1])/(pix1[0]-pix2[0])* \ - (min(pix1[0],pix2[0])+N.arange(maxline)-pix1[0])+pix1[1]) - else: - yline = N.round(min(pix1[1],pix2[1])+N.arange(maxline)) - xline = N.round((pix1[0]-pix2[0])/(pix1[1]-pix2[1])* \ - (min(pix1[1],pix2[1])+N.arange(maxline)-pix1[1])+pix1[0]) - rpixval = N.zeros(maxline, dtype=N.float32) - xbig = N.where(xline >= N.size(subim,0)) - xline[xbig] = N.size(subim,0) - 1 - ybig = N.where(yline >= N.size(subim,1)) - yline[ybig] = N.size(subim,1) - 1 - for i in range(maxline): - pixval = subim[int(xline[i]), int(yline[i])] - rpixval[i] = pixval - min_pixval = N.min(rpixval) - minind_p = N.argmin(rpixval) - maxind_p = N.argmax(rpixval) - - if minind_p in (0, maxline-1) and maxind_p in (0, maxline-1): - same_island_cont = True - if min_pixval >= min(flux1, flux2): - same_island_min = True - elif abs(min_pixval-min(flux1,flux2)) <= tol*isl.rms*img.opts.thresh_isl: same_island_min = True + same_island_cont = True + else: + if abs(pixdif[0]) > abs(pixdif[1]): + xline = N.round(min(pix1[0],pix2[0])+N.arange(maxline)) + yline = N.round((pix1[1]-pix2[1])/(pix1[0]-pix2[0])* \ + (min(pix1[0],pix2[0])+N.arange(maxline)-pix1[0])+pix1[1]) + else: + yline = N.round(min(pix1[1],pix2[1])+N.arange(maxline)) + xline = N.round((pix1[0]-pix2[0])/(pix1[1]-pix2[1])* \ + (min(pix1[1],pix2[1])+N.arange(maxline)-pix1[1])+pix1[0]) + rpixval = N.zeros(maxline, dtype=N.float32) + xbig = N.where(xline >= N.size(subim,0)) + xline[xbig] = N.size(subim,0) - 1 + ybig = N.where(yline >= N.size(subim,1)) + yline[ybig] = N.size(subim,1) - 1 + for i in range(maxline): + pixval = subim[int(xline[i]), int(yline[i])] + rpixval[i] = pixval + min_pixval = N.min(rpixval) + minind_p = N.argmin(rpixval) + maxind_p = N.argmax(rpixval) + + if minind_p in (0, maxline-1) and maxind_p in (0, maxline-1): + same_island_cont = True + if min_pixval >= min(flux1, flux2): + same_island_min = True + elif abs(min_pixval-min(flux1,flux2)) <= tol*isl.rms*img.opts.thresh_isl: + same_island_min = True return same_island_min, same_island_cont @@ -293,9 +293,9 @@ def same_island_dist(pair, g_list, tol=0.5): dist = sqrt(dy*dy + dx*dx) if dist <= tol*(fwhm1+fwhm2): - same_island = True + same_island = True else: - same_island = False + same_island = False return same_island @@ -353,16 +353,16 @@ def process_Multiple(self, img, g_sublist, mask, src_index, isrc, subim, isl, de x_ax, y_ax = N.indices(data.shape) if N.sum(~rmask) >=6: - para, ierr = func.fit_gaus2d(data, p_ini, x_ax, y_ax, rmask) - if (0.0 2.0*mompara[5]: mompara5E = 2.0*mompara[5] # Don't let errors get too large else: - mompara1E = 0.0 - mompara2E = 0.0 - mompara3E = 0.0 - mompara4E = 0.0 - mompara5E = 0.0 + mompara1E = 0.0 + mompara2E = 0.0 + mompara3E = 0.0 + mompara4E = 0.0 + mompara5E = 0.0 # Now add MC errors in quadrature with Condon (1997) errors size_skyE = [sqrt(mompara3E**2 + errors[3]**2) * sqrt(cdeltsq), diff --git a/bdsf/gausfit.py b/bdsf/gausfit.py index 9a2d3236..c1554490 100644 --- a/bdsf/gausfit.py +++ b/bdsf/gausfit.py @@ -324,8 +324,8 @@ def fit_island(self, isl, opts, img, ngmax=None, ffimg=None, ini_gausfit=None): ng1=0 ngmax=25 if ini_gausfit == 'nobeam' and not opts.fix_to_beam: - gaul = self.inigaus_nobeam(isl, thr0, beam, img) - ng1 = len(gaul); ngmax = ng1+2 + gaul = self.inigaus_nobeam(isl, thr0, beam, img) + ng1 = len(gaul); ngmax = ng1+2 if verbose: print('Initializing, ini_gausfit is',ini_gausfit,'gaul =',gaul,'ngmax =',ngmax) while iter < 5: @@ -349,14 +349,14 @@ def fit_island(self, isl, opts, img, ngmax=None, ffimg=None, ini_gausfit=None): ng1 = 0 ngmax = 25 while iter < 5: - iter += 1 - fitok = self.fit_iter(gaul, ng1, fcn, dof, beam, thr0, iter, 'simple', ngmax, verbose, g3_only) - 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 + iter += 1 + fitok = self.fit_iter(gaul, ng1, fcn, dof, beam, thr0, iter, 'simple', ngmax, verbose, g3_only) + 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 sm_isl = nd.binary_dilation(isl.mask_active) if (not fitok or len(gaul) == 0) and N.sum(~sm_isl) >= img.minpix_isl: if verbose: print('Fit still not OK, shrinking') @@ -367,14 +367,14 @@ def fit_island(self, isl, opts, img, ngmax=None, ffimg=None, ini_gausfit=None): ng1 = 0 ngmax = 25 while iter < 5: - iter += 1 - fitok = self.fit_iter(gaul, ng1, fcn, dof, beam, thr0, iter, 'simple', ngmax, verbose, g3_only) - 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 + iter += 1 + fitok = self.fit_iter(gaul, ng1, fcn, dof, beam, thr0, iter, 'simple', ngmax, verbose, g3_only) + 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 lg_isl = nd.binary_erosion(isl.mask_active) if (not fitok or len(gaul) == 0) and N.sum(~lg_isl) >= img.minpix_isl: if verbose: print('Fit still not OK, expanding') @@ -385,14 +385,14 @@ def fit_island(self, isl, opts, img, ngmax=None, ffimg=None, ini_gausfit=None): ng1 = 0 ngmax = 25 while iter < 5: - iter += 1 - fitok = self.fit_iter(gaul, ng1, fcn, dof, beam, thr0, iter, 'simple', ngmax, verbose, g3_only) - 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 + iter += 1 + fitok = self.fit_iter(gaul, ng1, fcn, dof, beam, thr0, iter, 'simple', ngmax, verbose, g3_only) + 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 or len(gaul) == 0: # If all else fails, try to use moment analysis @@ -520,24 +520,24 @@ def inigaus_fbdsm(self, isl, thr, beam, img): av = img.clipped_mean inipeak, iniposn, im1 = func.get_maxima(im, mask, thr_pos, isl.shape, beam, im_pos=im_pos) if len(inipeak) == 0: - av, stdnew, maxv, maxp, minv, minp = func.arrstatmask(im, mask) - inipeak = [maxv]; iniposn = [maxp] + av, stdnew, maxv, maxp, minv, minp = func.arrstatmask(im, mask) + inipeak = [maxv]; iniposn = [maxp] nmulsrc1 = len(iniposn) domore = True while domore: - domore = False - av, stdnew, maxv, maxp, minv, minp = func.arrstatmask(im1, mask) - if stdnew > isl.rms and maxv >= thr and maxv >= isl.mean+2.0*isl.rms: - domore = True - x1, y1 = N.array(iniposn).transpose() - dumr = N.sqrt((maxp[0]-x1)*(maxp[0]-x1)+(maxp[1]-y1)*(maxp[1]-y1)) - distbm = dumr/sqrt(beam[0]*beam[1]*fwsig*fwsig) - if N.any((distbm < 0.5) + (dumr < 2.2)): - domore = False - if domore: - iniposn.append(N.array(maxp)); inipeak.append(maxv) - im1 = func.mclean(im1, maxp, beam) + domore = False + av, stdnew, maxv, maxp, minv, minp = func.arrstatmask(im1, mask) + if stdnew > isl.rms and maxv >= thr and maxv >= isl.mean+2.0*isl.rms: + domore = True + x1, y1 = N.array(iniposn).transpose() + dumr = N.sqrt((maxp[0]-x1)*(maxp[0]-x1)+(maxp[1]-y1)*(maxp[1]-y1)) + distbm = dumr/sqrt(beam[0]*beam[1]*fwsig*fwsig) + if N.any((distbm < 0.5) + (dumr < 2.2)): + domore = False + if domore: + iniposn.append(N.array(maxp)); inipeak.append(maxv) + im1 = func.mclean(im1, maxp, beam) inipeak = N.array(inipeak); iniposn = N.array(iniposn) ind = list(N.argsort(inipeak)); ind.reverse() @@ -545,8 +545,8 @@ def inigaus_fbdsm(self, isl, thr, beam, img): iniposn = iniposn[ind] gaul = [] for i in range(len(inipeak)): - g = (float(inipeak[i]), int(iniposn[i][0]), int(iniposn[i][1])) + beam - gaul.append(g) + g = (float(inipeak[i]), int(iniposn[i][0]), int(iniposn[i][1])) + beam + gaul.append(g) return gaul, nmulsrc1, len(inipeak) @@ -586,69 +586,69 @@ def inigaus_nobeam(self, isl, thr, beam, img): av, stdnew, maxv, maxp, minv, minp = func.arrstatmask(im, mask) mom = func.momanalmask_gaus(isl.image-isl.islmean, isl.mask_active, 0, 1.0, True) if npeak <= 1: - g = (float(maxv), int(round(mom[1])), int(round(mom[2])), mom[3]/fwsig, \ - mom[4]/fwsig, mom[5]) - gaul.append(g) + g = (float(maxv), int(round(mom[1])), int(round(mom[2])), mom[3]/fwsig, \ + mom[4]/fwsig, mom[5]) + gaul.append(g) if npeak > 1: # markers start from 1=background, watershed starts from 1=background - watershed, markers = func.watershed(im, mask=isl.mask_active) - nshed = N.max(markers)-1 # excluding background - xm, ym = N.transpose([N.where(markers==i) for i in range(1,nshed+2)])[0] - coords = [c for c in N.transpose([xm,ym])[1:]] - alldists = [func.dist_2pt(c1, c2) for c1 in coords for c2 in coords if N.any(c1!=c2)] # has double - meandist = N.mean(alldists) # mean dist between all pairs of markers - # find at least some 'compact' sources - cscale = 3.0 - while True: - compact = []; invmask = [] - for ished in range(nshed): - shedmask = N.where(watershed==ished+2, False, True) + isl.mask_active # good unmasked pixels = 0 - imm = nd.binary_dilation(~shedmask, N.ones((3,3), int)) - xbad, ybad = N.where((imm==1)*(im>im[xm[ished+1], ym[ished+1]])) - imm[xbad, ybad] = 0 - invmask.append(imm); x, y = N.where(imm); xcen, ycen = N.mean(x), N.mean(y) # good pixels are now = 1 - dist = func.dist_2pt([xcen, ycen], [xm[ished+1], ym[ished+1]]) - if dist < max(cscale, meandist/4.0): - compact.append(True) # if not compact, break source + diffuse + watershed, markers = func.watershed(im, mask=isl.mask_active) + nshed = N.max(markers)-1 # excluding background + xm, ym = N.transpose([N.where(markers==i) for i in range(1,nshed+2)])[0] + coords = [c for c in N.transpose([xm,ym])[1:]] + alldists = [func.dist_2pt(c1, c2) for c1 in coords for c2 in coords if N.any(c1!=c2)] # has double + meandist = N.mean(alldists) # mean dist between all pairs of markers + # find at least some 'compact' sources + cscale = 3.0 + while True: + compact = []; invmask = [] + for ished in range(nshed): + shedmask = N.where(watershed==ished+2, False, True) + isl.mask_active # good unmasked pixels = 0 + imm = nd.binary_dilation(~shedmask, N.ones((3,3), int)) + xbad, ybad = N.where((imm==1)*(im>im[xm[ished+1], ym[ished+1]])) + imm[xbad, ybad] = 0 + invmask.append(imm); x, y = N.where(imm); xcen, ycen = N.mean(x), N.mean(y) # good pixels are now = 1 + dist = func.dist_2pt([xcen, ycen], [xm[ished+1], ym[ished+1]]) + if dist < max(cscale, meandist/4.0): + compact.append(True) # if not compact, break source + diffuse + else: + compact.append(False) + if N.any(compact): + break else: - compact.append(False) - if N.any(compact): - break - else: - # rescale to search for more compact sources - cscale*=1.5 - - if not N.all(compact): - o_avsize = [] - ind = N.where(compact)[0] - for i in ind: o_avsize.append(N.sum(invmask[i])) - avsize = sqrt(N.mean(N.array(o_avsize))) - for i in range(len(compact)): - if not compact[i]: # make them all compact - newmask = N.zeros(imm.shape, bool) - newmask[max(0,int(xm[i+1]-avsize/2)):min(im.shape[0],int(xm[i+1]+avsize/2)), \ - max(0,int(ym[i+1]-avsize/2)):min(im.shape[1],int(ym[i+1]+avsize/2))] = True - invmask[i] = invmask[i]*newmask - resid = N.zeros(im.shape, dtype=N.float32) # approx fit all compact ones - for i in range(nshed): - mask1 = ~invmask[i] - size = sqrt(N.sum(invmask))/fwsig - xf, yf = coords[i][0], coords[i][1] - p_ini = [im[xf, yf], xf, yf, size, size, 0.0] - x, y = N.indices(im.shape) - p, success = func.fit_gaus2d(im*invmask[i], p_ini, x, y) - resid = resid + func.gaus_2d(p, x, y) - gaul.append(p) - resid = im - resid - if not N.all(compact): # just add one gaussian to fit whole unmasked island - maxv = N.max(resid) # assuming resid has only diffuse emission. can be false - x, y = N.where(~isl.mask_active); xcen = N.mean(x); ycen = N.mean(y) - invm = ~isl.mask_active - #bound = invm - nd.grey_erosion(invm, footprint = N.ones((3,3), int)) # better to use bound for ellipse fitting - mom = func.momanalmask_gaus(invm, N.zeros(invm.shape, dtype=N.int16), 0, 1.0, True) - g = (maxv, xcen, ycen, mom[3]/fwsig, mom[4]/fwsig, mom[5]-90.) - gaul.append(g) - coords.append([xcen, ycen]) + # rescale to search for more compact sources + cscale*=1.5 + + if not N.all(compact): + o_avsize = [] + ind = N.where(compact)[0] + for i in ind: o_avsize.append(N.sum(invmask[i])) + avsize = sqrt(N.mean(N.array(o_avsize))) + for i in range(len(compact)): + if not compact[i]: # make them all compact + newmask = N.zeros(imm.shape, bool) + newmask[max(0,int(xm[i+1]-avsize/2)):min(im.shape[0],int(xm[i+1]+avsize/2)), \ + max(0,int(ym[i+1]-avsize/2)):min(im.shape[1],int(ym[i+1]+avsize/2))] = True + invmask[i] = invmask[i]*newmask + resid = N.zeros(im.shape, dtype=N.float32) # approx fit all compact ones + for i in range(nshed): + mask1 = ~invmask[i] + size = sqrt(N.sum(invmask))/fwsig + xf, yf = coords[i][0], coords[i][1] + p_ini = [im[xf, yf], xf, yf, size, size, 0.0] + x, y = N.indices(im.shape) + p, success = func.fit_gaus2d(im*invmask[i], p_ini, x, y) + resid = resid + func.gaus_2d(p, x, y) + gaul.append(p) + resid = im - resid + if not N.all(compact): # just add one gaussian to fit whole unmasked island + maxv = N.max(resid) # assuming resid has only diffuse emission. can be false + x, y = N.where(~isl.mask_active); xcen = N.mean(x); ycen = N.mean(y) + invm = ~isl.mask_active + #bound = invm - nd.grey_erosion(invm, footprint = N.ones((3,3), int)) # better to use bound for ellipse fitting + mom = func.momanalmask_gaus(invm, N.zeros(invm.shape, dtype=N.int16), 0, 1.0, True) + g = (maxv, xcen, ycen, mom[3]/fwsig, mom[4]/fwsig, mom[5]-90.) + gaul.append(g) + coords.append([xcen, ycen]) return gaul @@ -675,33 +675,33 @@ def fit_iter(self, gaul, ng1, fcn, dof, beam, thr, iter, inifit, ngmax, verbose= ### no error-checking here, they MUST fit fcn.reset() for ig in range(ng1): - g = gaul[ig] - self.add_gaussian(fcn, g, dof, g3_only) + g = gaul[ig] + self.add_gaussian(fcn, g, dof, g3_only) ### do a round of fitting if any initials were provided if verbose: print('About to call C++ wrapper') fitok = True if len(gaul) != 0: - fitok = fit(fcn, final=0, verbose=verbose) + fitok = fit(fcn, final=0, verbose=verbose) if verbose: print('Returned from the fit') ### iteratively add gaussians while there are high peaks ### in the image and fitting converges while fitok: - peak, coords = fcn.find_peak() - if peak < thr: ### no good peaks left - break - if len(fcn.parameters) < ngmax and iter == 1 and inifit == 'default' and len(gaul) >= ng1+1: - ng1 = ng1 + 1 - g = gaul[ng1-1] - else: - if len(fcn.parameters) < ngmax: - g = [peak, coords[0], coords[1]] + beam + peak, coords = fcn.find_peak() + if peak < thr: ### no good peaks left + break + if len(fcn.parameters) < ngmax and iter == 1 and inifit == 'default' and len(gaul) >= ng1+1: + ng1 = ng1 + 1 + g = gaul[ng1-1] else: - break - fitok &= self.add_gaussian(fcn, g, dof, g3_only) + if len(fcn.parameters) < ngmax: + g = [peak, coords[0], coords[1]] + beam + else: + break + fitok &= self.add_gaussian(fcn, g, dof, g3_only) - fitok &= fit(fcn, final=0, verbose=verbose) + fitok &= fit(fcn, final=0, verbose=verbose) ### and one last fit with higher precision ### make sure we return False when fitok==False due to lack @@ -790,12 +790,12 @@ def _flag_gaussian(self, g, beam, thr, peak, shape, opts, mask, image, size_bms) th1 = divmod(th, 180)[1] th1 = th1/180.0*pi if ss1 > 1e4 and ss2 > 1e4: - xbox = 1e9; ybox = 1e9 + xbox = 1e9; ybox = 1e9 else: - xbox = 2.0*(abs(ss1*cos(th1)*cos(th1))+abs(ss2*ss2/ss1*sin(th1)*sin(th1)))/ \ - sqrt(cos(th1)*cos(th1)+ss2*ss2/ss1/ss1*sin(th1)*sin(th1)) - ybox = 2.0*(abs(ss1*sin(th1)*sin(th1))+abs(ss2*ss2/ss1*cos(th1)*cos(th1)))/ \ - sqrt(sin(th1)*sin(th1)+ss2*ss2/ss1/ss1*cos(th1)*cos(th1)) + xbox = 2.0*(abs(ss1*cos(th1)*cos(th1))+abs(ss2*ss2/ss1*sin(th1)*sin(th1)))/ \ + sqrt(cos(th1)*cos(th1)+ss2*ss2/ss1/ss1*sin(th1)*sin(th1)) + ybox = 2.0*(abs(ss1*sin(th1)*sin(th1))+abs(ss2*ss2/ss1*cos(th1)*cos(th1)))/ \ + sqrt(sin(th1)*sin(th1)+ss2*ss2/ss1/ss1*cos(th1)*cos(th1)) ### now check all conditions border = opts.flag_bordersize @@ -817,7 +817,7 @@ def _flag_gaussian(self, g, beam, thr, peak, shape, opts, mask, image, size_bms) # Check image value at Gaussian center im_val_at_cen = nd.map_coordinates(image, [N.array([x1]), N.array([x2])]) if A > opts.flag_maxsnr*im_val_at_cen: - flag += 2 + flag += 2 borx1_1 = x1 - border if borx1_1 < 0: borx1_1 = 0 borx1_2 = x1 + border + 1 @@ -834,9 +834,9 @@ def _flag_gaussian(self, g, beam, thr, peak, shape, opts, mask, image, size_bms) if ybox > opts.flag_maxsize_isl*shape[1]: flag += 32 if s1*s2 > opts.flag_maxsize_bm*beam[0]*beam[1]: flag += 64 if opts.flag_smallsrc: - if s1*s2 < opts.flag_minsize_bm*beam[0]*beam[1]: flag += 128 + if s1*s2 < opts.flag_minsize_bm*beam[0]*beam[1]: flag += 128 if not opts.flag_smallsrc: - if s1*s2 == 0.: flag += 128 + if s1*s2 == 0.: flag += 128 if ss1/ss2 > 2.0: # Only check for fairly elliptical Gaussians, as this condition diff --git a/bdsf/interface.py b/bdsf/interface.py index bb27b9dd..b032f123 100644 --- a/bdsf/interface.py +++ b/bdsf/interface.py @@ -606,7 +606,7 @@ def print_opts(grouped_opts_list, img, banner=None): else: # Since this is the last entry in the options list and is a # tuple, it cannot be an expandable option, so make it nc color - parvalstr = nc + k + nc + ' ..' + parvalstr = nc + k + nc + ' ..' if "'" in valstr: len_without_formatting = len(k) + len(str(val)) + 5 else: @@ -991,11 +991,11 @@ def write_catalog(img, outfile=None, format='bbs', srcroot=None, catalog_type='g print('\033[91mERROR\033[0m: Image has not been fit. Please run process_image first.') return False if catalog_type == 'shap' and not 'shapelets' in img.completed_Ops: - print('\033[91mERROR\033[0m: Image has not been decomposed into shapelets. Please run process_image first.') - return False + print('\033[91mERROR\033[0m: Image has not been decomposed into shapelets. Please run process_image first.') + return False if catalog_type == 'srl' and not 'gaul2srl' in img.completed_Ops: - print('\033[91mERROR\033[0m: Gaussians have not been grouped into sources. Please run process_image first.') - return False + print('\033[91mERROR\033[0m: Gaussians have not been grouped into sources. Please run process_image first.') + return False format = format.lower() patch = bbs_patches filename = outfile diff --git a/bdsf/islands.py b/bdsf/islands.py index bc5d68df..f6f7f4a5 100644 --- a/bdsf/islands.py +++ b/bdsf/islands.py @@ -268,10 +268,10 @@ def setpara_bdsm(self, img, det_file): ops = [] for op in chain: - if isinstance(op, type): - ops.append(op()) - else: - ops.append(op) + if isinstance(op, type): + ops.append(op()) + else: + ops.append(op) return ops, opts @@ -306,7 +306,7 @@ def __init__(self, img, mask, mean, rms, labels, bbox, idx, colname='Coeff_matrix', units=None) if not copy: - ### we make bbox slightly bigger + ### we make bbox slightly bigger self.oldbbox = bbox self.oldidx = idx bbox = self.__expand_bbox(bbox, img.shape) diff --git a/bdsf/make_residimage.py b/bdsf/make_residimage.py index 04e09ee1..9ff08c03 100644 --- a/bdsf/make_residimage.py +++ b/bdsf/make_residimage.py @@ -37,9 +37,9 @@ def __call__(self, img): for g in img.gaussians: C1, C2 = g.centre_pix if hasattr(g, 'wisland_id') and img.waveletimage: - isl = img.islands[g.wisland_id] + isl = img.islands[g.wisland_id] else: - isl = img.islands[g.island_id] + isl = img.islands[g.island_id] b = self.find_bbox(thresh*isl.rms, g) bbox = N.s_[max(0, int(C1-b)):min(shape[0], int(C1+b+1)), @@ -105,14 +105,14 @@ def __call__(self, img): fimg = N.zeros(shape, dtype=N.float32) for isl in img.islands: - if hasattr(isl, 'shapelet_beta'): - if isl.shapelet_beta > 0: # make sure shapelet has nonzero scale for this island - mask=isl.mask_active - cen=isl.shapelet_centre-N.array(isl.origin) - basis, beta, nmax, cf = isl.shapelet_basis, isl.shapelet_beta, \ - isl.shapelet_nmax, isl.shapelet_cf - image_recons=reconstruct_shapelets(isl.shape, mask, basis, beta, cen, nmax, cf) - fimg[tuple(isl.bbox)] += image_recons + if hasattr(isl, 'shapelet_beta'): + if isl.shapelet_beta > 0: # make sure shapelet has nonzero scale for this island + mask=isl.mask_active + cen=isl.shapelet_centre-N.array(isl.origin) + basis, beta, nmax, cf = isl.shapelet_basis, isl.shapelet_beta, \ + isl.shapelet_nmax, isl.shapelet_cf + image_recons=reconstruct_shapelets(isl.shape, mask, basis, beta, cen, nmax, cf) + fimg[tuple(isl.bbox)] += image_recons model_shap = fimg resid_shap = img.ch0_arr - fimg diff --git a/bdsf/multi_proc.py b/bdsf/multi_proc.py index 336c5eac..4fb1eabf 100644 --- a/bdsf/multi_proc.py +++ b/bdsf/multi_proc.py @@ -130,7 +130,7 @@ def run_tasks(procs, err_q, out_q, num): # Remove extra dimension added by array_split result_list = [] for result in results: - result_list += result + result_list += result return result_list diff --git a/bdsf/mylogger.py b/bdsf/mylogger.py index ba9bfbee..124c9008 100644 --- a/bdsf/mylogger.py +++ b/bdsf/mylogger.py @@ -15,60 +15,60 @@ import copy def init_logger(logfilename, quiet=False, debug=False): - logging.USERINFO = logging.INFO + 1 - logging.addLevelName(logging.USERINFO, 'USERINFO') - logger = logging.root - logger.setLevel(logging.DEBUG) - - # First remove any existing handlers (in case PyBDSM has been run - # before in this session but the quiet or debug options have changed - while len(logger.handlers) > 0: - logger.removeHandler(logger.handlers[0]) - - # File handlers - fh = ColorStripperHandler(logfilename) - if debug: - # For log file and debug on, print name and levelname - fh.setLevel(logging.DEBUG) - fmt1 = MultiLineFormatter('%(asctime)s %(name)-20s:: %(levelname)-8s: '\ - '%(message)s', - datefmt='%a %d-%m-%Y %H:%M:%S') - else: - # For log file and debug off, don't print name and levelname as - # they have no meaning to the user. - fh.setLevel(logging.INFO) - fmt1 = MultiLineFormatter('%(asctime)s:: %(levelname)-8s: %(message)s', - datefmt='%a %d-%m-%Y %H:%M:%S') - fh.setFormatter(fmt1) - logger.addHandler(fh) - - # Console handler for warning, error, and critical: format includes levelname - # ANSI colors are used - ch = logging.StreamHandler() - ch.setLevel(logging.WARNING) - fmt2 = logging.Formatter('\033[31;1m%(levelname)s\033[0m: %(message)s') - ch.setFormatter(fmt2) - logger.addHandler(ch) - - # Console handler for USERINFO only: format does not include levelname - # (the user does not need to see the levelname, as it has no meaning to them) - # ANSI colors are allowed - chi = logging.StreamHandler() - chi.addFilter(InfoFilter()) - if quiet: - # prints nothing, since filter lets only USERINFO through - chi.setLevel(logging.WARNING) - else: - # prints only USERINFO - chi.setLevel(logging.USERINFO) - fmt3 = logging.Formatter('%(message)s') - chi.setFormatter(fmt3) - logger.addHandler(chi) + logging.USERINFO = logging.INFO + 1 + logging.addLevelName(logging.USERINFO, 'USERINFO') + logger = logging.root + logger.setLevel(logging.DEBUG) + + # First remove any existing handlers (in case PyBDSM has been run + # before in this session but the quiet or debug options have changed + while len(logger.handlers) > 0: + logger.removeHandler(logger.handlers[0]) + + # File handlers + fh = ColorStripperHandler(logfilename) + if debug: + # For log file and debug on, print name and levelname + fh.setLevel(logging.DEBUG) + fmt1 = MultiLineFormatter('%(asctime)s %(name)-20s:: %(levelname)-8s: '\ + '%(message)s', + datefmt='%a %d-%m-%Y %H:%M:%S') + else: + # For log file and debug off, don't print name and levelname as + # they have no meaning to the user. + fh.setLevel(logging.INFO) + fmt1 = MultiLineFormatter('%(asctime)s:: %(levelname)-8s: %(message)s', + datefmt='%a %d-%m-%Y %H:%M:%S') + fh.setFormatter(fmt1) + logger.addHandler(fh) + + # Console handler for warning, error, and critical: format includes levelname + # ANSI colors are used + ch = logging.StreamHandler() + ch.setLevel(logging.WARNING) + fmt2 = logging.Formatter('\033[31;1m%(levelname)s\033[0m: %(message)s') + ch.setFormatter(fmt2) + logger.addHandler(ch) + + # Console handler for USERINFO only: format does not include levelname + # (the user does not need to see the levelname, as it has no meaning to them) + # ANSI colors are allowed + chi = logging.StreamHandler() + chi.addFilter(InfoFilter()) + if quiet: + # prints nothing, since filter lets only USERINFO through + chi.setLevel(logging.WARNING) + else: + # prints only USERINFO + chi.setLevel(logging.USERINFO) + fmt3 = logging.Formatter('%(message)s') + chi.setFormatter(fmt3) + logger.addHandler(chi) class InfoFilter(logging.Filter): - # Lets only USERINFO through - def filter(self, rec): - return rec.levelno == logging.USERINFO + # Lets only USERINFO through + def filter(self, rec): + return rec.levelno == logging.USERINFO class MultiLineFormatter(logging.Formatter): def format(self, record): @@ -79,44 +79,44 @@ def format(self, record): return str def userinfo(mylog, desc_str, val_str=''): - """Writes a nicely formatted string to the log file and console - - mylog = logger - desc_str = description string / message - val_str = value string - - Message is constructed as: - 'desc_str .... : val_str' - """ - bc = '\033[1;34m' # Blue - nc = '\033[0m' # Normal text color - - if val_str == '': - sep = '' - if desc_str[:1] == '\n': - bc += '\n' - desc_str = desc_str[1:] - desc_str = bc + '--> ' + desc_str + nc - else: - sep = ' : ' - if len(desc_str) < 40: - desc_str += ' ' - if len(desc_str) < 40: - while len(desc_str) < 41: - desc_str += '.' + """Writes a nicely formatted string to the log file and console + + mylog = logger + desc_str = description string / message + val_str = value string + + Message is constructed as: + 'desc_str .... : val_str' + """ + bc = '\033[1;34m' # Blue + nc = '\033[0m' # Normal text color + + if val_str == '': + sep = '' + if desc_str[:1] == '\n': + bc += '\n' + desc_str = desc_str[1:] + desc_str = bc + '--> ' + desc_str + nc else: - while len(desc_str) < 41: - desc_str += ' ' - mylog.log(logging.USERINFO, desc_str+sep+val_str) + sep = ' : ' + if len(desc_str) < 40: + desc_str += ' ' + if len(desc_str) < 40: + while len(desc_str) < 41: + desc_str += '.' + else: + while len(desc_str) < 41: + desc_str += ' ' + mylog.log(logging.USERINFO, desc_str+sep+val_str) class ColorStripperHandler(logging.FileHandler): - def emit(self, record): - """Strips ANSI color codes from file stream""" - myrecord = copy.copy(record) - nocolor_msg = strip_color(myrecord.msg) - myrecord.msg = nocolor_msg - logging.FileHandler.emit(self, myrecord) + def emit(self, record): + """Strips ANSI color codes from file stream""" + myrecord = copy.copy(record) + nocolor_msg = strip_color(myrecord.msg) + myrecord.msg = nocolor_msg + logging.FileHandler.emit(self, myrecord) def strip_color(msg): """Strips specific ANSI color codes from an input string @@ -127,9 +127,9 @@ def strip_color(msg): nocolor_msg = '' a = msg.split('\033[1;34m') for b in a: - c = b.split('\033[0m') - for d in c: - e = d.split('\033[31;1m') - for f in e: - nocolor_msg += f + c = b.split('\033[0m') + for d in c: + e = d.split('\033[31;1m') + for f in e: + nocolor_msg += f return nocolor_msg diff --git a/bdsf/nat/__init__.py b/bdsf/nat/__init__.py index 7e803a98..a69542f4 100644 --- a/bdsf/nat/__init__.py +++ b/bdsf/nat/__init__.py @@ -1838,4 +1838,3 @@ def help(choice = None): print('Natgrid, parameters, table, regrid, aspectSlope, singlePoint, natgrids, seti, geti, setr, getr, setc, getc, getaspects, getslopes, pntinits, pnts, pntend, natgridd, setrd, getrd, getaspectd, getsloped, pntinitd, pntd, pntendd') return None - diff --git a/bdsf/output.py b/bdsf/output.py index 979ffba5..96d4ed29 100644 --- a/bdsf/output.py +++ b/bdsf/output.py @@ -1028,9 +1028,9 @@ def list_and_sort_gaussians(img, patch=None, root=None, # Unrecognized property --> Don't sort indx = list(range(len(gausindx))) for i, si in enumerate(indx): - outlist_sorted[0][i] = outlist[0][si] - outnames_sorted[0][i] = outnames[0][si] - patchnames_sorted = list(patchnames) + outlist_sorted[0][i] = outlist[0][si] + outnames_sorted[0][i] = outnames[0][si] + patchnames_sorted = list(patchnames) else: outlist_sorted = list(outlist) outnames_sorted = list(outnames) @@ -1116,8 +1116,8 @@ def make_output_columns(obj, fits=False, objtype='gaul', incl_spin=False, for n, name in enumerate(names): if hasattr(obj, name): if name in ['specin_flux', 'specin_fluxE', 'specin_freq']: - # As these are variable length lists, they must - # (unfortunately) be treated differently. + # As these are variable length lists, they must + # (unfortunately) be treated differently. val = obj.__getattribute__(name) colname = obj.__dict__[name+'_def']._colname units = obj.__dict__[name+'_def']._units diff --git a/bdsf/plotresults.py b/bdsf/plotresults.py index 19c109bd..9fd7fa46 100644 --- a/bdsf/plotresults.py +++ b/bdsf/plotresults.py @@ -378,15 +378,15 @@ def plotresults(img, ch0_image=True, rms_image=True, mean_image=True, pl.title('Pyramidal Sources for\nWavelet Scale J = ' + str(j_with_gaus[j])) for pyr in img.pyrsrcs: - for iisl, isl in enumerate(pyr.islands): - jj = pyr.jlevels[iisl] - jindx = j_with_gaus.index(jj) - col = colours[pyr.pyr_id % 6] - ind = N.where(~isl.mask_active) - cmd = "ax" + str(jindx + index_first_waveplot + 1) + \ - ".plot(ind[0]+isl.origin[0], "\ - "ind[1]+isl.origin[1], '.', color=col)" - exec(cmd) + for iisl, isl in enumerate(pyr.islands): + jj = pyr.jlevels[iisl] + jindx = j_with_gaus.index(jj) + col = colours[pyr.pyr_id % 6] + ind = N.where(~isl.mask_active) + cmd = "ax" + str(jindx + index_first_waveplot + 1) + \ + ".plot(ind[0]+isl.origin[0], "\ + "ind[1]+isl.origin[1], '.', color=col)" + exec(cmd) fig.canvas.mpl_connect('key_press_event', on_press) fig.canvas.mpl_connect('pick_event', on_pick) diff --git a/bdsf/polarisation.py b/bdsf/polarisation.py index eaeb0f98..a35e200d 100644 --- a/bdsf/polarisation.py +++ b/bdsf/polarisation.py @@ -56,271 +56,271 @@ class Op_polarisation(Op): def __call__(self, img): mylog = mylogger.logging.getLogger("PyBDSM."+img.log+"Polarisatn") if img.opts.polarisation_do: - mylog.info('Extracting polarisation properties for all sources') - pols = ['I', 'Q', 'U', 'V'] - - # Run gausfit and gual2srl on PI image to look for polarized sources - # undetected in I - fit_PI = img.opts.pi_fit - n_new = 0 - ch0_pi = N.sqrt(img.ch0_Q_arr**2 + img.ch0_U_arr**2) - img.ch0_pi_arr = ch0_pi - - if fit_PI: - from . import _run_op_list - mylogger.userinfo(mylog, "\nChecking PI image for new sources") - - mask = img.mask_arr - - # Set up image object for PI image. - pi_chain, pi_opts = self.setpara_bdsm(img) - pimg = Image(pi_opts) - pimg.beam = img.beam - pimg.pixel_beam = img.pixel_beam - pimg.pixel_beamarea = img.pixel_beamarea - pimg.log = 'PI.' - pimg.basedir = img.basedir - pimg.imagename = img.imagename - pimg.frequency = img.frequency - pimg.equinox = img.equinox - pimg.shape = img.shape - pimg.pix2beam = img.pix2beam - pimg.beam2pix = img.beam2pix - pimg.pix2gaus = img.pix2gaus - pimg.gaus2pix = img.gaus2pix - pimg.pix2sky = img.pix2sky - pimg.sky2pix = img.sky2pix - pimg.pix2coord = img.pix2coord - pimg.wcs_obj = img.wcs_obj - pimg.mask_arr = mask - pimg.masked = img.masked - pimg.ch0_arr = ch0_pi - pimg._pi = True - - success = _run_op_list(pimg, pi_chain) - if not success: - return - - img.pi_islands = pimg.islands - img.pi_gaussians = pimg.gaussians - img.pi_sources = pimg.sources - - # Now check for new sources in the PI image that are not - # found in the Stokes I image. If any new sources are found, - # adjust their IDs to follow after those found in I. - new_isl = [] - new_src = [] - new_gaus = [] - n_new_src = 0 - if len(img.islands) == 0: - isl_id = 0 - src_id = 0 - gaus_id = 0 - else: - isl_id = img.islands[-1].island_id - src_id = img.sources[-1].source_id - gaus_id = img.gaussians[-1].gaus_num - for pi_isl in pimg.islands: - new_sources = [] - for pi_src in pi_isl.sources: - if img.pyrank[int(img.sky2pix(pi_src.posn_sky_max)[0]), - int(img.sky2pix(pi_src.posn_sky_max)[1])] == -1: - src_id += 1 - pi_src._pi = True - pi_src.island_id = isl_id - pi_src.source_id = src_id - pi_src.spec_indx = N.NaN - pi_src.e_spec_indx = N.NaN - pi_src.spec_norm = N.NaN - pi_src.specin_flux = [N.NaN] - pi_src.specin_fluxE = [N.NaN] - pi_src.specin_freq = [N.NaN] - pi_src.specin_freq0 = N.NaN - for gaus in pi_src.gaussians: - gaus.island_id = isl_id - gaus.source_id = src_id - gaus.spec_indx = N.NaN - gaus.e_spec_indx = N.NaN - gaus.spec_norm = N.NaN - gaus.specin_flux = [N.NaN] - gaus.specin_fluxE = [N.NaN] - gaus.specin_freq = [N.NaN] - gaus.specin_freq0 = N.NaN - new_sources.append(pi_src) - new_src.append(pi_src) - n_new_src += 1 - for g in pi_src.gaussians: - gaus_id += 1 - new_gaus.append(g) - g.gaus_num = gaus_id - if len(new_sources) > 0: - isl_id += 1 - pi_isl.sources = new_sources - pi_isl.island_id = isl_id - pi_isl._pi = True - new_isl.append(pi_isl) - - n_new = len(new_isl) - mylogger.userinfo(mylog, "New sources found in PI image", '%i (%i total)' % - (n_new, img.nsrc+n_new)) - - if n_new > 0: - img.islands += new_isl - img.sources += new_src - img.gaussians += new_gaus - img.nsrc += n_new_src - renumber_islands(img) - - bar = statusbar.StatusBar('Calculating polarisation properties .... : ', 0, img.nsrc) - if img.opts.quiet == False: - bar.start() - - for isl in img.islands: - isl_bbox = isl.bbox - ch0_I = img.ch0_arr[tuple(isl_bbox)] - ch0_Q = img.ch0_Q_arr[tuple(isl_bbox)] - ch0_U = img.ch0_U_arr[tuple(isl_bbox)] - ch0_V = img.ch0_V_arr[tuple(isl_bbox)] - ch0_images = [ch0_I, ch0_Q, ch0_U, ch0_V] - - for i, src in enumerate(isl.sources): - # For each source, assume the morphology does not change - # across the Stokes cube. This assumption allows us to fit - # the Gaussians of each source to each Stokes image by - # simply fitting only the overall normalizations of the - # individual Gaussians. - # - # First, fit all source Gaussians to each Stokes image: - x, y = N.mgrid[isl_bbox] - gg = src.gaussians - fitfix = N.ones(len(gg)) # fit only normalization - srcmask = isl.mask_active - total_flux = N.zeros((4, len(fitfix)), dtype=N.float32) # array of fluxes: N_Stokes x N_Gaussians - errors = N.zeros((4, len(fitfix)), dtype=N.float32) # array of fluxes: N_Stokes x N_Gaussians - - for sind, image in enumerate(ch0_images): - if (sind==0 and hasattr(src, '_pi')) or sind > 0: # Fit I only for PI sources - p, ep = func.fit_mulgaus2d(image, gg, x, y, srcmask, fitfix) - for ig in range(len(fitfix)): - bm_pix = N.array([img.pixel_beam()[0], img.pixel_beam()[1], img.pixel_beam()[2]]) - total_flux[sind, ig] = p[ig*6]*p[ig*6+3]*p[ig*6+4]/(bm_pix[0]*bm_pix[1]) - p = N.insert(p, N.arange(len(fitfix))*6+6, total_flux[sind]) - if sind > 0: - rms_img = img.__getattribute__('rms_'+pols[sind]+'_arr') - else: - rms_img = img.rms_arr - if len(rms_img.shape) > 1: - rms_isl = rms_img[tuple(isl.bbox)].mean() - else: - rms_isl = rms_img - errors[sind] = func.get_errors(img, p, rms_isl)[6] - - # Now, assign fluxes to each Gaussian. - src_flux_I = 0.0 - src_flux_Q = 0.0 - src_flux_U = 0.0 - src_flux_V = 0.0 - src_flux_I_err_sq = 0.0 - src_flux_Q_err_sq = 0.0 - src_flux_U_err_sq = 0.0 - src_flux_V_err_sq = 0.0 - - for ig, gaussian in enumerate(src.gaussians): - init_gaus_attr(gaussian) - flux_I = total_flux[0, ig] - flux_I_err = abs(errors[0, ig]) - flux_Q = total_flux[1, ig] - flux_Q_err = abs(errors[1, ig]) - flux_U = total_flux[2, ig] - flux_U_err = abs(errors[2, ig]) - flux_V = total_flux[3, ig] - flux_V_err = abs(errors[3, ig]) - - if hasattr(src, '_pi'): - gaussian.total_flux = flux_I - gaussian.total_fluxE = flux_I_err - gaussian.total_flux_Q = flux_Q - gaussian.total_flux_U = flux_U - gaussian.total_flux_V = flux_V - gaussian.total_fluxE_Q = flux_Q_err - gaussian.total_fluxE_U = flux_U_err - gaussian.total_fluxE_V = flux_V_err - + mylog.info('Extracting polarisation properties for all sources') + pols = ['I', 'Q', 'U', 'V'] + + # Run gausfit and gual2srl on PI image to look for polarized sources + # undetected in I + fit_PI = img.opts.pi_fit + n_new = 0 + ch0_pi = N.sqrt(img.ch0_Q_arr**2 + img.ch0_U_arr**2) + img.ch0_pi_arr = ch0_pi + + if fit_PI: + from . import _run_op_list + mylogger.userinfo(mylog, "\nChecking PI image for new sources") + + mask = img.mask_arr + + # Set up image object for PI image. + pi_chain, pi_opts = self.setpara_bdsm(img) + pimg = Image(pi_opts) + pimg.beam = img.beam + pimg.pixel_beam = img.pixel_beam + pimg.pixel_beamarea = img.pixel_beamarea + pimg.log = 'PI.' + pimg.basedir = img.basedir + pimg.imagename = img.imagename + pimg.frequency = img.frequency + pimg.equinox = img.equinox + pimg.shape = img.shape + pimg.pix2beam = img.pix2beam + pimg.beam2pix = img.beam2pix + pimg.pix2gaus = img.pix2gaus + pimg.gaus2pix = img.gaus2pix + pimg.pix2sky = img.pix2sky + pimg.sky2pix = img.sky2pix + pimg.pix2coord = img.pix2coord + pimg.wcs_obj = img.wcs_obj + pimg.mask_arr = mask + pimg.masked = img.masked + pimg.ch0_arr = ch0_pi + pimg._pi = True + + success = _run_op_list(pimg, pi_chain) + if not success: + return + + img.pi_islands = pimg.islands + img.pi_gaussians = pimg.gaussians + img.pi_sources = pimg.sources + + # Now check for new sources in the PI image that are not + # found in the Stokes I image. If any new sources are found, + # adjust their IDs to follow after those found in I. + new_isl = [] + new_src = [] + new_gaus = [] + n_new_src = 0 + if len(img.islands) == 0: + isl_id = 0 + src_id = 0 + gaus_id = 0 + else: + isl_id = img.islands[-1].island_id + src_id = img.sources[-1].source_id + gaus_id = img.gaussians[-1].gaus_num + for pi_isl in pimg.islands: + new_sources = [] + for pi_src in pi_isl.sources: + if img.pyrank[int(img.sky2pix(pi_src.posn_sky_max)[0]), + int(img.sky2pix(pi_src.posn_sky_max)[1])] == -1: + src_id += 1 + pi_src._pi = True + pi_src.island_id = isl_id + pi_src.source_id = src_id + pi_src.spec_indx = N.NaN + pi_src.e_spec_indx = N.NaN + pi_src.spec_norm = N.NaN + pi_src.specin_flux = [N.NaN] + pi_src.specin_fluxE = [N.NaN] + pi_src.specin_freq = [N.NaN] + pi_src.specin_freq0 = N.NaN + for gaus in pi_src.gaussians: + gaus.island_id = isl_id + gaus.source_id = src_id + gaus.spec_indx = N.NaN + gaus.e_spec_indx = N.NaN + gaus.spec_norm = N.NaN + gaus.specin_flux = [N.NaN] + gaus.specin_fluxE = [N.NaN] + gaus.specin_freq = [N.NaN] + gaus.specin_freq0 = N.NaN + new_sources.append(pi_src) + new_src.append(pi_src) + n_new_src += 1 + for g in pi_src.gaussians: + gaus_id += 1 + new_gaus.append(g) + g.gaus_num = gaus_id + if len(new_sources) > 0: + isl_id += 1 + pi_isl.sources = new_sources + pi_isl.island_id = isl_id + pi_isl._pi = True + new_isl.append(pi_isl) + + n_new = len(new_isl) + mylogger.userinfo(mylog, "New sources found in PI image", '%i (%i total)' % + (n_new, img.nsrc+n_new)) + + if n_new > 0: + img.islands += new_isl + img.sources += new_src + img.gaussians += new_gaus + img.nsrc += n_new_src + renumber_islands(img) + + bar = statusbar.StatusBar('Calculating polarisation properties .... : ', 0, img.nsrc) + if img.opts.quiet == False: + bar.start() + + for isl in img.islands: + isl_bbox = isl.bbox + ch0_I = img.ch0_arr[tuple(isl_bbox)] + ch0_Q = img.ch0_Q_arr[tuple(isl_bbox)] + ch0_U = img.ch0_U_arr[tuple(isl_bbox)] + ch0_V = img.ch0_V_arr[tuple(isl_bbox)] + ch0_images = [ch0_I, ch0_Q, ch0_U, ch0_V] + + for i, src in enumerate(isl.sources): + # For each source, assume the morphology does not change + # across the Stokes cube. This assumption allows us to fit + # the Gaussians of each source to each Stokes image by + # simply fitting only the overall normalizations of the + # individual Gaussians. + # + # First, fit all source Gaussians to each Stokes image: + x, y = N.mgrid[isl_bbox] + gg = src.gaussians + fitfix = N.ones(len(gg)) # fit only normalization + srcmask = isl.mask_active + total_flux = N.zeros((4, len(fitfix)), dtype=N.float32) # array of fluxes: N_Stokes x N_Gaussians + errors = N.zeros((4, len(fitfix)), dtype=N.float32) # array of fluxes: N_Stokes x N_Gaussians + + for sind, image in enumerate(ch0_images): + if (sind==0 and hasattr(src, '_pi')) or sind > 0: # Fit I only for PI sources + p, ep = func.fit_mulgaus2d(image, gg, x, y, srcmask, fitfix) + for ig in range(len(fitfix)): + bm_pix = N.array([img.pixel_beam()[0], img.pixel_beam()[1], img.pixel_beam()[2]]) + total_flux[sind, ig] = p[ig*6]*p[ig*6+3]*p[ig*6+4]/(bm_pix[0]*bm_pix[1]) + p = N.insert(p, N.arange(len(fitfix))*6+6, total_flux[sind]) + if sind > 0: + rms_img = img.__getattribute__('rms_'+pols[sind]+'_arr') + else: + rms_img = img.rms_arr + if len(rms_img.shape) > 1: + rms_isl = rms_img[tuple(isl.bbox)].mean() + else: + rms_isl = rms_img + errors[sind] = func.get_errors(img, p, rms_isl)[6] + + # Now, assign fluxes to each Gaussian. + src_flux_I = 0.0 + src_flux_Q = 0.0 + src_flux_U = 0.0 + src_flux_V = 0.0 + src_flux_I_err_sq = 0.0 + src_flux_Q_err_sq = 0.0 + src_flux_U_err_sq = 0.0 + src_flux_V_err_sq = 0.0 + + for ig, gaussian in enumerate(src.gaussians): + init_gaus_attr(gaussian) + flux_I = total_flux[0, ig] + flux_I_err = abs(errors[0, ig]) + flux_Q = total_flux[1, ig] + flux_Q_err = abs(errors[1, ig]) + flux_U = total_flux[2, ig] + flux_U_err = abs(errors[2, ig]) + flux_V = total_flux[3, ig] + flux_V_err = abs(errors[3, ig]) + + if hasattr(src, '_pi'): + gaussian.total_flux = flux_I + gaussian.total_fluxE = flux_I_err + gaussian.total_flux_Q = flux_Q + gaussian.total_flux_U = flux_U + gaussian.total_flux_V = flux_V + gaussian.total_fluxE_Q = flux_Q_err + gaussian.total_fluxE_U = flux_U_err + gaussian.total_fluxE_V = flux_V_err + + if hasattr(src, '_pi'): + src_flux_I += flux_I + src_flux_I_err_sq += flux_I_err**2 + src_flux_Q += flux_Q + src_flux_U += flux_U + src_flux_V += flux_V + src_flux_Q_err_sq += flux_Q_err**2 + src_flux_U_err_sq += flux_U_err**2 + src_flux_V_err_sq += flux_V_err**2 + + # Calculate and store polarisation fractions and angle for each Gaussian in the island + # For this we need the I flux, which we can just take from g.total_flux and src.total_flux + flux_I = gaussian.total_flux + flux_I_err = gaussian.total_fluxE + stokes = [flux_I, flux_Q, flux_U, flux_V] + stokes_err = [flux_I_err, flux_Q_err, flux_U_err, flux_V_err] + + lpol_frac, lpol_frac_loerr, lpol_frac_hierr = self.calc_lpol_fraction(stokes, stokes_err) # linear pol fraction + lpol_ang, lpol_ang_err = self.calc_lpol_angle(stokes, stokes_err) # linear pol angle + cpol_frac, cpol_frac_loerr, cpol_frac_hierr = self.calc_cpol_fraction(stokes, stokes_err) # circular pol fraction + tpol_frac, tpol_frac_loerr, tpol_frac_hierr = self.calc_tpol_fraction(stokes, stokes_err) # total pol fraction + + gaussian.lpol_fraction = lpol_frac + gaussian.lpol_fraction_loerr = lpol_frac_loerr + gaussian.lpol_fraction_hierr = lpol_frac_hierr + gaussian.cpol_fraction = cpol_frac + gaussian.cpol_fraction_loerr = cpol_frac_loerr + gaussian.cpol_fraction_hierr = cpol_frac_hierr + gaussian.tpol_fraction = tpol_frac + gaussian.tpol_fraction_loerr = tpol_frac_loerr + gaussian.tpol_fraction_hierr = tpol_frac_hierr + gaussian.lpol_angle = lpol_ang + gaussian.lpol_angle_err = lpol_ang_err + + # Store fluxes for each source in the island + init_src_attr(src) if hasattr(src, '_pi'): - src_flux_I += flux_I - src_flux_I_err_sq += flux_I_err**2 - src_flux_Q += flux_Q - src_flux_U += flux_U - src_flux_V += flux_V - src_flux_Q_err_sq += flux_Q_err**2 - src_flux_U_err_sq += flux_U_err**2 - src_flux_V_err_sq += flux_V_err**2 - - # Calculate and store polarisation fractions and angle for each Gaussian in the island + src.total_flux = src_flux_I + src.total_fluxE = N.sqrt(src_flux_I_err_sq) + src.total_flux_Q = src_flux_Q + src.total_flux_U = src_flux_U + src.total_flux_V = src_flux_V + src.total_fluxE_Q = N.sqrt(src_flux_Q_err_sq) + src.total_fluxE_U = N.sqrt(src_flux_U_err_sq) + src.total_fluxE_V = N.sqrt(src_flux_V_err_sq) + + # Calculate and store polarisation fractions and angle for each source in the island # For this we need the I flux, which we can just take from g.total_flux and src.total_flux - flux_I = gaussian.total_flux - flux_I_err = gaussian.total_fluxE - stokes = [flux_I, flux_Q, flux_U, flux_V] - stokes_err = [flux_I_err, flux_Q_err, flux_U_err, flux_V_err] + src_flux_I = src.total_flux + src_flux_I_err = src.total_fluxE + stokes = [src_flux_I, src_flux_Q, src_flux_U, src_flux_V] + stokes_err = [src_flux_I_err, N.sqrt(src_flux_Q_err_sq), N.sqrt(src_flux_U_err_sq), N.sqrt(src_flux_V_err_sq)] lpol_frac, lpol_frac_loerr, lpol_frac_hierr = self.calc_lpol_fraction(stokes, stokes_err) # linear pol fraction lpol_ang, lpol_ang_err = self.calc_lpol_angle(stokes, stokes_err) # linear pol angle cpol_frac, cpol_frac_loerr, cpol_frac_hierr = self.calc_cpol_fraction(stokes, stokes_err) # circular pol fraction tpol_frac, tpol_frac_loerr, tpol_frac_hierr = self.calc_tpol_fraction(stokes, stokes_err) # total pol fraction - gaussian.lpol_fraction = lpol_frac - gaussian.lpol_fraction_loerr = lpol_frac_loerr - gaussian.lpol_fraction_hierr = lpol_frac_hierr - gaussian.cpol_fraction = cpol_frac - gaussian.cpol_fraction_loerr = cpol_frac_loerr - gaussian.cpol_fraction_hierr = cpol_frac_hierr - gaussian.tpol_fraction = tpol_frac - gaussian.tpol_fraction_loerr = tpol_frac_loerr - gaussian.tpol_fraction_hierr = tpol_frac_hierr - gaussian.lpol_angle = lpol_ang - gaussian.lpol_angle_err = lpol_ang_err - - # Store fluxes for each source in the island - init_src_attr(src) - if hasattr(src, '_pi'): - src.total_flux = src_flux_I - src.total_fluxE = N.sqrt(src_flux_I_err_sq) - src.total_flux_Q = src_flux_Q - src.total_flux_U = src_flux_U - src.total_flux_V = src_flux_V - src.total_fluxE_Q = N.sqrt(src_flux_Q_err_sq) - src.total_fluxE_U = N.sqrt(src_flux_U_err_sq) - src.total_fluxE_V = N.sqrt(src_flux_V_err_sq) - - # Calculate and store polarisation fractions and angle for each source in the island - # For this we need the I flux, which we can just take from g.total_flux and src.total_flux - src_flux_I = src.total_flux - src_flux_I_err = src.total_fluxE - stokes = [src_flux_I, src_flux_Q, src_flux_U, src_flux_V] - stokes_err = [src_flux_I_err, N.sqrt(src_flux_Q_err_sq), N.sqrt(src_flux_U_err_sq), N.sqrt(src_flux_V_err_sq)] - - lpol_frac, lpol_frac_loerr, lpol_frac_hierr = self.calc_lpol_fraction(stokes, stokes_err) # linear pol fraction - lpol_ang, lpol_ang_err = self.calc_lpol_angle(stokes, stokes_err) # linear pol angle - cpol_frac, cpol_frac_loerr, cpol_frac_hierr = self.calc_cpol_fraction(stokes, stokes_err) # circular pol fraction - tpol_frac, tpol_frac_loerr, tpol_frac_hierr = self.calc_tpol_fraction(stokes, stokes_err) # total pol fraction - - src.lpol_fraction = lpol_frac - src.lpol_fraction_loerr = lpol_frac_loerr - src.lpol_fraction_hierr = lpol_frac_hierr - src.cpol_fraction = cpol_frac - src.cpol_fraction_loerr = cpol_frac_loerr - src.cpol_fraction_hierr = cpol_frac_hierr - src.tpol_fraction = tpol_frac - src.tpol_fraction_loerr = tpol_frac_loerr - src.tpol_fraction_hierr = tpol_frac_hierr - src.lpol_angle = lpol_ang - src.lpol_angle_err = lpol_ang_err - if bar.started: - bar.increment() - bar.stop() - img.completed_Ops.append('polarisation') - - #################################################################################### + src.lpol_fraction = lpol_frac + src.lpol_fraction_loerr = lpol_frac_loerr + src.lpol_fraction_hierr = lpol_frac_hierr + src.cpol_fraction = cpol_frac + src.cpol_fraction_loerr = cpol_frac_loerr + src.cpol_fraction_hierr = cpol_frac_hierr + src.tpol_fraction = tpol_frac + src.tpol_fraction_loerr = tpol_frac_loerr + src.tpol_fraction_hierr = tpol_frac_hierr + src.lpol_angle = lpol_ang + src.lpol_angle_err = lpol_ang_err + if bar.started: + bar.increment() + bar.stop() + img.completed_Ops.append('polarisation') + + #################################################################################### def calc_lpol_fraction(self, stokes, err): """ Calculate linear polarisation fraction and error from: stokes = [I, Q, U, V] and err = [Ierr, Qerr, Uerr, Verr] @@ -353,7 +353,7 @@ def calc_lpol_fraction(self, stokes, err): return lfrac, loerr, uperr - #################################################################################### + #################################################################################### def calc_cpol_fraction(self, stokes, err): """ Calculate circular polarisation fraction and error from: stokes = [I, Q, U, V] and err = [Ierr, Qerr, Uerr, Verr] @@ -377,7 +377,7 @@ def calc_cpol_fraction(self, stokes, err): return cfrac, loerr, uperr - #################################################################################### + #################################################################################### def calc_tpol_fraction(self, stokes, err): """ Calculate total polarisation fraction and error from: stokes = [I, Q, U, V] and err = [Ierr, Qerr, Uerr, Verr] @@ -409,7 +409,7 @@ def calc_tpol_fraction(self, stokes, err): return tfrac, loerr, uperr - #################################################################################### + #################################################################################### def calc_lpol_angle(self, stokes, err, sig=3.0): """ Calculate linear polarisation angle and error (in degrees) from: stokes = [I, Q, U, V] and err = [Ierr, Qerr, Uerr, Verr] @@ -426,7 +426,7 @@ def calc_lpol_angle(self, stokes, err, sig=3.0): return ang, dang - #################################################################################### + #################################################################################### def debias(self, pflux, QUerr): """ Debiases the linearly polarised flux using the same method used for the NVSS catalog (see ftp://ftp.cv.nrao.edu/pub/nvss/catalog.ps). @@ -469,7 +469,7 @@ def check_frac(self, frac, loerr, uperr): uperr = 1.0 - frac return frac, loerr, uperr - #################################################################################### + #################################################################################### def setpara_bdsm(self, img): chain = [Op_preprocess, Op_rmsimage(), Op_threshold(), Op_islands(), Op_gausfit(), Op_gaul2srl(), Op_make_residimage()] @@ -487,10 +487,10 @@ def setpara_bdsm(self, img): ops = [] for op in chain: - if isinstance(op, type): - ops.append(op()) - else: - ops.append(op) + if isinstance(op, type): + ops.append(op()) + else: + ops.append(op) return ops, opts diff --git a/bdsf/preprocess.py b/bdsf/preprocess.py index 63648ce2..19e1b10d 100644 --- a/bdsf/preprocess.py +++ b/bdsf/preprocess.py @@ -28,13 +28,13 @@ def __call__(self, img): kappa = img.opts.kappa_clip if img.opts.polarisation_do: - pols = ['I', 'Q', 'U', 'V'] - ch0images = [img.ch0_arr, img.ch0_Q_arr, img.ch0_U_arr, img.ch0_V_arr] - img.clipped_mean_QUV = [] - img.clipped_rms_QUV = [] + pols = ['I', 'Q', 'U', 'V'] + ch0images = [img.ch0_arr, img.ch0_Q_arr, img.ch0_U_arr, img.ch0_V_arr] + img.clipped_mean_QUV = [] + img.clipped_rms_QUV = [] else: - pols = ['I'] # assume I is always present - ch0images = [img.ch0_arr] + pols = ['I'] # assume I is always present + ch0images = [img.ch0_arr] if hasattr(img, 'rms_mask'): mask = img.rms_mask @@ -43,29 +43,29 @@ def __call__(self, img): opts = img.opts for ipol, pol in enumerate(pols): - image = ch0images[ipol] - - ### basic stats - mean, rms, cmean, crms, cnt = bstat(image, mask, kappa) - if cnt > 198: cmean = mean; crms = rms - if pol == 'I': - if func.approx_equal(crms, 0.0, rel=None): - raise RuntimeError('Clipped rms appears to be zero. Check for regions '\ - 'with values of 0 and\nblank them (with NaNs) '\ - 'or use trim_box to exclude them.') - img.raw_mean = mean - img.raw_rms = rms - img.clipped_mean= cmean - img.clipped_rms = crms - mylog.info('%s %.4f %s %.4f %s ' % ("Raw mean (Stokes I) = ", mean*1000.0, \ - 'mJy and raw rms = ',rms*1000.0, 'mJy')) - mylog.info('%s %.4f %s %s %.4f %s ' % ("sigma clipped mean (Stokes I) = ", cmean*1000.0, \ - 'mJy and ','sigma clipped rms = ',crms*1000.0, 'mJy')) - else: - img.clipped_mean_QUV.append(cmean) - img.clipped_rms_QUV.append(crms) - mylog.info('%s %s %s %.4f %s %s %.4f %s ' % ("sigma clipped mean (Stokes ", pol, ") = ", cmean*1000.0, \ - 'mJy and ','sigma clipped rms = ',crms*1000.0, 'mJy')) + image = ch0images[ipol] + + ### basic stats + mean, rms, cmean, crms, cnt = bstat(image, mask, kappa) + if cnt > 198: cmean = mean; crms = rms + if pol == 'I': + if func.approx_equal(crms, 0.0, rel=None): + raise RuntimeError('Clipped rms appears to be zero. Check for regions '\ + 'with values of 0 and\nblank them (with NaNs) '\ + 'or use trim_box to exclude them.') + img.raw_mean = mean + img.raw_rms = rms + img.clipped_mean= cmean + img.clipped_rms = crms + mylog.info('%s %.4f %s %.4f %s ' % ("Raw mean (Stokes I) = ", mean*1000.0, \ + 'mJy and raw rms = ',rms*1000.0, 'mJy')) + mylog.info('%s %.4f %s %s %.4f %s ' % ("sigma clipped mean (Stokes I) = ", cmean*1000.0, \ + 'mJy and ','sigma clipped rms = ',crms*1000.0, 'mJy')) + else: + img.clipped_mean_QUV.append(cmean) + img.clipped_rms_QUV.append(crms) + mylog.info('%s %s %s %.4f %s %s %.4f %s ' % ("sigma clipped mean (Stokes ", pol, ") = ", cmean*1000.0, \ + 'mJy and ','sigma clipped rms = ',crms*1000.0, 'mJy')) image = img.ch0_arr # Check if pixels are outside the universe @@ -94,11 +94,11 @@ def __call__(self, img): if mask is not None: img.blankpix = N.sum(mask) if img.blankpix == 0: - max_idx = image.argmax() - min_idx = image.argmin() + max_idx = image.argmax() + min_idx = image.argmin() else: - max_idx = N.nanargmax(image) - min_idx = N.nanargmin(image) + max_idx = N.nanargmax(image) + min_idx = N.nanargmin(image) img.maxpix_coord = N.unravel_index(max_idx, shape) img.minpix_coord = N.unravel_index(min_idx, shape) @@ -126,25 +126,25 @@ def __call__(self, img): ### 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: - 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.') - mylog.info('Taking number of pixels above 5-sigma as 1.') - img.bmpersrc_th = N.product(shape)/((alpha_sourcecounts-1.)*n) - mylog.info('%s %6.2f' % ('Estimated bmpersrc_th = ', img.bmpersrc_th)) + 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.') + mylog.info('Taking number of pixels above 5-sigma as 1.') + img.bmpersrc_th = N.product(shape)/((alpha_sourcecounts-1.)*n) + mylog.info('%s %6.2f' % ('Estimated bmpersrc_th = ', img.bmpersrc_th)) else: - img.bmpersrc_th = opts.bmpersrc_th - mylog.info('%s %6.2f' % ('Taking default bmpersrc_th = ', img.bmpersrc_th)) + img.bmpersrc_th = opts.bmpersrc_th + mylog.info('%s %6.2f' % ('Taking default bmpersrc_th = ', img.bmpersrc_th)) confused = False if opts.mean_map == 'default': - if img.bmpersrc_th <= 25. or cmean/crms >= 0.1: - confused = True + if img.bmpersrc_th <= 25. or cmean/crms >= 0.1: + confused = True img.confused = confused mylog.info('Parameter confused is '+str(img.confused)) @@ -158,19 +158,19 @@ def outside_univ(self,img): noutside = 0 n, m = img.ch0_arr.shape for i in range(n): - for j in range(m): - out = False - err = '' - pix1 = (i,j) - try: - skyc = img.pix2sky(pix1) - pix2 = img.sky2pix(skyc) - if abs(pix1[0]-pix2[0]) > 0.5 or abs(pix1[1]-pix2[1]) > 0.5: out=True - except RuntimeError as err: - pass - if out or ("8" in str(err)): - noutside += 1 - ch0 = img.ch0_arr - ch0[pix1] = float("NaN") - img.ch0_arr = ch0 + for j in range(m): + out = False + err = '' + pix1 = (i,j) + try: + skyc = img.pix2sky(pix1) + pix2 = img.sky2pix(skyc) + if abs(pix1[0]-pix2[0]) > 0.5 or abs(pix1[1]-pix2[1]) > 0.5: out=True + except RuntimeError as err: + pass + if out or ("8" in str(err)): + noutside += 1 + ch0 = img.ch0_arr + ch0[pix1] = float("NaN") + img.ch0_arr = ch0 return noutside diff --git a/bdsf/psf_vary.py b/bdsf/psf_vary.py index 8fb38051..fad307ab 100644 --- a/bdsf/psf_vary.py +++ b/bdsf/psf_vary.py @@ -31,349 +31,349 @@ class Op_psf_vary(Op): def __call__(self, img): - if img.opts.psf_vary_do: - mylog = mylogger.logging.getLogger("PyBDSM."+img.log+"Psf_Vary") - mylogger.userinfo(mylog, '\nEstimating PSF variations') - opts = img.opts - dir = img.basedir + '/misc/' - plot = False # debug figures - image = img.ch0_arr - - try: - from astropy.io import fits as pyfits - old_pyfits = False - except ImportError as err: + if img.opts.psf_vary_do: + mylog = mylogger.logging.getLogger("PyBDSM."+img.log+"Psf_Vary") + mylogger.userinfo(mylog, '\nEstimating PSF variations') + opts = img.opts + dir = img.basedir + '/misc/' + plot = False # debug figures + image = img.ch0_arr + try: - from packaging.version import Version - except ImportError: - from distutils.version import StrictVersion as Version - import pyfits - if Version(pyfits.__version__) < Version('2.2'): - old_pyfits = True - else: + from astropy.io import fits as pyfits old_pyfits = False + except ImportError as err: + try: + from packaging.version import Version + except ImportError: + from distutils.version import StrictVersion as Version + import pyfits + if Version(pyfits.__version__) < Version('2.2'): + old_pyfits = True + else: + old_pyfits = False - if old_pyfits: - mylog.warning('PyFITS version is too old: psf_vary module skipped') - return + if old_pyfits: + mylog.warning('PyFITS version is too old: psf_vary module skipped') + return - if img.nisl == 0: - mylog.warning("No islands found. Skipping PSF variation estimation.") - img.completed_Ops.append('psf_vary') - return - - if opts.psf_fwhm is not None: - # User has specified a constant PSF to use, so skip PSF fitting/etc. - psf_maj = opts.psf_fwhm[0] # FWHM in deg - psf_min = opts.psf_fwhm[1] # FWHM in deg - psf_pa = opts.psf_fwhm[2] # PA in deg - mylogger.userinfo(mylog, 'Using constant PSF (major, minor, pos angle)', - '(%.5e, %.5e, %s) degrees' % (psf_maj, psf_maj, - round(psf_pa, 1))) - else: - # Use did not specify a constant PSF to use, so estimate it - over = 2 - generators = opts.psf_generators; nsig = opts.psf_nsig; kappa2 = opts.psf_kappa2 - snrtop = opts.psf_snrtop; snrbot = opts.psf_snrbot; snrcutstack = opts.psf_snrcutstack - gencode = opts.psf_gencode; primarygen = opts.psf_primarygen; itess_method = opts.psf_itess_method - tess_sc = opts.psf_tess_sc; tess_fuzzy= opts.psf_tess_fuzzy - bright_snr_cut = opts.psf_high_snr - s_only = opts.psf_stype_only - if opts.psf_snrcut < 5.0: - mylogger.userinfo(mylog, "Value of psf_snrcut too low; increasing to 5") - snrcut = 5.0 + if img.nisl == 0: + mylog.warning("No islands found. Skipping PSF variation estimation.") + img.completed_Ops.append('psf_vary') + return + + if opts.psf_fwhm is not None: + # User has specified a constant PSF to use, so skip PSF fitting/etc. + psf_maj = opts.psf_fwhm[0] # FWHM in deg + psf_min = opts.psf_fwhm[1] # FWHM in deg + psf_pa = opts.psf_fwhm[2] # PA in deg + mylogger.userinfo(mylog, 'Using constant PSF (major, minor, pos angle)', + '(%.5e, %.5e, %s) degrees' % (psf_maj, psf_maj, + round(psf_pa, 1))) else: - snrcut = opts.psf_snrcut - img.psf_snrcut = snrcut - 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 + # Use did not specify a constant PSF to use, so estimate it + over = 2 + generators = opts.psf_generators; nsig = opts.psf_nsig; kappa2 = opts.psf_kappa2 + snrtop = opts.psf_snrtop; snrbot = opts.psf_snrbot; snrcutstack = opts.psf_snrcutstack + gencode = opts.psf_gencode; primarygen = opts.psf_primarygen; itess_method = opts.psf_itess_method + tess_sc = opts.psf_tess_sc; tess_fuzzy= opts.psf_tess_fuzzy + bright_snr_cut = opts.psf_high_snr + s_only = opts.psf_stype_only + if opts.psf_snrcut < 5.0: + mylogger.userinfo(mylog, "Value of psf_snrcut too low; increasing to 5") + snrcut = 5.0 + else: + snrcut = opts.psf_snrcut + img.psf_snrcut = snrcut + 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 + else: + high_snrcut = opts.psf_high_snr else: high_snrcut = opts.psf_high_snr - else: - high_snrcut = opts.psf_high_snr - img.psf_high_snr = high_snrcut - - wtfns=['unity', 'roundness', 'log10', 'sqrtlog10'] - if 0 <= itess_method < 4: tess_method=wtfns[itess_method] - else: tess_method='unity' - - ### now put all relevant gaussian parameters into a list - ngaus = img.ngaus - nsrc = img.nsrc - num = N.zeros(nsrc, dtype=N.int32) - peak = N.zeros(nsrc) - xc = N.zeros(nsrc) - yc = N.zeros(nsrc) - bmaj = N.zeros(nsrc) - bmin = N.zeros(nsrc) - bpa = N.zeros(nsrc) - code = N.array(['']*nsrc); - rms = N.zeros(nsrc) - src_id_list = [] - for i, src in enumerate(img.sources): - src_max = 0.0 - for gmax in src.gaussians: - # Take only brightest Gaussian per source - if gmax.peak_flux > src_max: - src_max = gmax.peak_flux - g = gmax - num[i] = i - peak[i] = g.peak_flux - xc[i] = g.centre_pix[0] - yc[i] = g.centre_pix[1] - bmaj[i] = g.size_pix[0] - bmin[i] = g.size_pix[1] - bpa[i] = g.size_pix[2] - code[i] = img.sources[g.source_id].code - rms[i] = img.islands[g.island_id].rms - gauls = (num, peak, xc, yc, bmaj, bmin, bpa, code, rms) - tr_gauls = self.trans_gaul(gauls) - - # takes gaussians with code=S and snr > snrcut. - if s_only: - tr = [n for n in tr_gauls if n[1]/n[8]>snrcut and n[7] == 'S'] - else: - tr = [n for n in tr_gauls if n[1]/n[8]>snrcut] - g_gauls = self.trans_gaul(tr) - - # computes statistics of fitted sizes. Same as psfvary_fullstat.f in fBDSM. - bmaj_a, bmaj_r, bmaj_ca, bmaj_cr, ni = _cbdsm.bstat(bmaj, None, nsig) - bmin_a, bmin_r, bmin_ca, bmin_cr, ni = _cbdsm.bstat(bmin, None, nsig) - bpa_a, bpa_r, bpa_ca, bpa_cr, ni = _cbdsm.bstat(bpa, None, nsig) - - # get subset of sources deemed to be unresolved. Same as size_ksclip_wenss.f in fBDSM. - flag_unresolved = self.get_unresolved(g_gauls, img.beam, nsig, kappa2, over, img.psf_high_snr, plot) - if len(flag_unresolved) == 0: - mylog.warning('Insufficient number of sources to determine PSF variation.\nTry changing the PSF options or specify a (constant) PSF with the "psf_fwhm" option') - return + img.psf_high_snr = high_snrcut + + wtfns=['unity', 'roundness', 'log10', 'sqrtlog10'] + if 0 <= itess_method < 4: tess_method=wtfns[itess_method] + else: tess_method='unity' + + ### now put all relevant gaussian parameters into a list + ngaus = img.ngaus + nsrc = img.nsrc + num = N.zeros(nsrc, dtype=N.int32) + peak = N.zeros(nsrc) + xc = N.zeros(nsrc) + yc = N.zeros(nsrc) + bmaj = N.zeros(nsrc) + bmin = N.zeros(nsrc) + bpa = N.zeros(nsrc) + code = N.array(['']*nsrc); + rms = N.zeros(nsrc) + src_id_list = [] + for i, src in enumerate(img.sources): + src_max = 0.0 + for gmax in src.gaussians: + # Take only brightest Gaussian per source + if gmax.peak_flux > src_max: + src_max = gmax.peak_flux + g = gmax + num[i] = i + peak[i] = g.peak_flux + xc[i] = g.centre_pix[0] + yc[i] = g.centre_pix[1] + bmaj[i] = g.size_pix[0] + bmin[i] = g.size_pix[1] + bpa[i] = g.size_pix[2] + code[i] = img.sources[g.source_id].code + rms[i] = img.islands[g.island_id].rms + gauls = (num, peak, xc, yc, bmaj, bmin, bpa, code, rms) + tr_gauls = self.trans_gaul(gauls) + + # takes gaussians with code=S and snr > snrcut. + if s_only: + tr = [n for n in tr_gauls if n[1]/n[8]>snrcut and n[7] == 'S'] + else: + tr = [n for n in tr_gauls if n[1]/n[8]>snrcut] + g_gauls = self.trans_gaul(tr) + + # computes statistics of fitted sizes. Same as psfvary_fullstat.f in fBDSM. + bmaj_a, bmaj_r, bmaj_ca, bmaj_cr, ni = _cbdsm.bstat(bmaj, None, nsig) + bmin_a, bmin_r, bmin_ca, bmin_cr, ni = _cbdsm.bstat(bmin, None, nsig) + bpa_a, bpa_r, bpa_ca, bpa_cr, ni = _cbdsm.bstat(bpa, None, nsig) + + # get subset of sources deemed to be unresolved. Same as size_ksclip_wenss.f in fBDSM. + flag_unresolved = self.get_unresolved(g_gauls, img.beam, nsig, kappa2, over, img.psf_high_snr, plot) + if len(flag_unresolved) == 0: + mylog.warning('Insufficient number of sources to determine PSF variation.\nTry changing the PSF options or specify a (constant) PSF with the "psf_fwhm" option') + return - # see how much the SNR-weighted sizes of unresolved sources differ from the synthesized beam. - wtsize_beam_snr = self.av_psf(g_gauls, img.beam, flag_unresolved) + # see how much the SNR-weighted sizes of unresolved sources differ from the synthesized beam. + wtsize_beam_snr = self.av_psf(g_gauls, img.beam, flag_unresolved) - # filter out resolved sources - tr_gaul = self.trans_gaul(g_gauls) - tr = [n for i, n in enumerate(tr_gaul) if flag_unresolved[i]] - g_gauls = self.trans_gaul(tr) - mylogger.userinfo(mylog, 'Number of unresolved sources', str(len(g_gauls[0]))) + # filter out resolved sources + tr_gaul = self.trans_gaul(g_gauls) + tr = [n for i, n in enumerate(tr_gaul) if flag_unresolved[i]] + g_gauls = self.trans_gaul(tr) + mylogger.userinfo(mylog, 'Number of unresolved sources', str(len(g_gauls[0]))) - # get a list of voronoi generators. vorogenS has values (and not None) if generators='field'. - vorogenP, vorogenS = self.get_voronoi_generators(g_gauls, generators, gencode, snrcut, snrtop, snrbot, snrcutstack) - mylogger.userinfo(mylog, 'Number of generators for PSF variation', str(len(vorogenP[0]))) - if len(vorogenP[0]) < 3: - mylog.warning('Insufficient number of generators') - return - - mylogger.userinfo(mylog, 'Tesselating image') - # group generators into tiles - tile_prop = self.edit_vorogenlist(vorogenP, frac=0.9) - - # tesselate the image - 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', volrank, img, dir) - - tile_list, tile_coord, tile_snr = tile_prop - ntile = len(tile_list) - bar = statusbar.StatusBar('Determining PSF variation ............... : ', 0, ntile) - mylogger.userinfo(mylog, 'Number of tiles for PSF variation', str(ntile)) - - # For each tile, calculate the weighted averaged psf image. Also for all the sources in the image. - cdelt = list(img.wcs_obj.acdelt[0:2]) - factor=3. - psfimages, psfcoords, totpsfimage, psfratio, psfratio_aper = self.psf_in_tile(image, img.beam, g_gauls, \ - cdelt, factor, snrcutstack, volrank, tile_prop, plot, img) - npsf = len(psfimages) - - if opts.psf_use_shap: - if opts.psf_fwhm is None: - # use totpsfimage to get beta, centre and nmax for shapelet decomposition. Use nmax=5 or 6 - mask=N.zeros(totpsfimage.shape, dtype=bool) - (m1, m2, m3)=func.moment(totpsfimage, mask) - betainit=sqrt(m3[0]*m3[1])*2.0 * 1.4 - tshape = totpsfimage.shape - cen = N.array(N.unravel_index(N.argmax(totpsfimage), tshape))+[1,1] - cen = tuple(cen) - nmax = 12 - basis = 'cartesian' - betarange = [0.5,sqrt(betainit*max(tshape))] - beta, error = sh.shape_varybeta(totpsfimage, mask, basis, betainit, cen, nmax, betarange, plot) - if error == 1: print(' Unable to find minimum in beta') - - # decompose all the psf images using the beta from above - nmax=12; psf_cf=[] - for i in range(npsf): - psfim = psfimages[i] - cf = sh.decompose_shapelets(psfim, mask, basis, beta, cen, nmax, mode='') - psf_cf.append(cf) - if img.opts.quiet == False: - bar.increment() - bar.stop() - - # transpose the psf image list - xt, yt = N.transpose(tile_coord) - tr_psf_cf = N.transpose(N.array(psf_cf)) - - # interpolate the coefficients across the image. Ok, interpolate in scipy for - # irregular grids is crap. doesnt even pass through some of the points. - # for now, fit polynomial. - compress = 100.0 - x, y = N.transpose(psfcoords) - if len(x) < 3: - mylog.warning('Insufficient number of tiles to do interpolation of PSF variation') + # get a list of voronoi generators. vorogenS has values (and not None) if generators='field'. + vorogenP, vorogenS = self.get_voronoi_generators(g_gauls, generators, gencode, snrcut, snrtop, snrbot, snrcutstack) + mylogger.userinfo(mylog, 'Number of generators for PSF variation', str(len(vorogenP[0]))) + if len(vorogenP[0]) < 3: + mylog.warning('Insufficient number of generators') return - psf_coeff_interp, xgrid, ygrid = self.interp_shapcoefs(nmax, tr_psf_cf, psfcoords, image.shape, \ - compress, plot) - - psfshape = psfimages[0].shape - skip = 5 - aa = self.create_psf_grid(psf_coeff_interp, image.shape, xgrid, ygrid, skip, nmax, psfshape, \ - basis, beta, cen, totpsfimage, plot) - img.psf_images = aa - else: - if opts.psf_fwhm is None: - if ntile < 4: - mylog.warning('Insufficient number of tiles to do interpolation of PSF variation') - return - else: - # Fit stacked PSFs with Gaussians and measure aperture fluxes - bm_pix = N.array([img.pixel_beam()[0]*fwsig, img.pixel_beam()[1]*fwsig, img.pixel_beam()[2]]) - psf_maj = N.zeros(npsf) - psf_min = N.zeros(npsf) - psf_pa = N.zeros(npsf) - if img.opts.quiet == False: - bar.start() - for i in range(ntile): + mylogger.userinfo(mylog, 'Tesselating image') + # group generators into tiles + tile_prop = self.edit_vorogenlist(vorogenP, frac=0.9) + + # tesselate the image + 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', volrank, img, dir) + + tile_list, tile_coord, tile_snr = tile_prop + ntile = len(tile_list) + bar = statusbar.StatusBar('Determining PSF variation ............... : ', 0, ntile) + mylogger.userinfo(mylog, 'Number of tiles for PSF variation', str(ntile)) + + # For each tile, calculate the weighted averaged psf image. Also for all the sources in the image. + cdelt = list(img.wcs_obj.acdelt[0:2]) + factor=3. + psfimages, psfcoords, totpsfimage, psfratio, psfratio_aper = self.psf_in_tile(image, img.beam, g_gauls, \ + cdelt, factor, snrcutstack, volrank, tile_prop, plot, img) + npsf = len(psfimages) + + if opts.psf_use_shap: + if opts.psf_fwhm is None: + # use totpsfimage to get beta, centre and nmax for shapelet decomposition. Use nmax=5 or 6 + mask=N.zeros(totpsfimage.shape, dtype=bool) + (m1, m2, m3)=func.moment(totpsfimage, mask) + betainit=sqrt(m3[0]*m3[1])*2.0 * 1.4 + tshape = totpsfimage.shape + cen = N.array(N.unravel_index(N.argmax(totpsfimage), tshape))+[1,1] + cen = tuple(cen) + nmax = 12 + basis = 'cartesian' + betarange = [0.5,sqrt(betainit*max(tshape))] + beta, error = sh.shape_varybeta(totpsfimage, mask, basis, betainit, cen, nmax, betarange, plot) + if error == 1: print(' Unable to find minimum in beta') + + # decompose all the psf images using the beta from above + nmax=12; psf_cf=[] + for i in range(npsf): psfim = psfimages[i] - mask = N.zeros(psfim.shape, dtype=bool) - x_ax, y_ax = N.indices(psfim.shape) - maxv = N.max(psfim) - p_ini = [maxv, (psfim.shape[0]-1)/2.0*1.1, (psfim.shape[1]-1)/2.0*1.1, bm_pix[0]/fwsig*1.3, - bm_pix[1]/fwsig*1.1, bm_pix[2]*2] - para, ierr = func.fit_gaus2d(psfim, p_ini, x_ax, y_ax, mask) - ### first extent is major - if para[3] < para[4]: - para[3:5] = para[4:2:-1] - para[5] += 90 - ### clip position angle - para[5] = divmod(para[5], 180)[1] - - psf_maj[i] = para[3] - psf_min[i] = para[4] - posang = para[5] - while posang >= 180.0: - posang -= 180.0 - psf_pa[i] = posang - + cf = sh.decompose_shapelets(psfim, mask, basis, beta, cen, nmax, mode='') + psf_cf.append(cf) if img.opts.quiet == False: bar.increment() bar.stop() - # Interpolate Gaussian parameters - 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] - nimgs = len(psf_maps) - bar = statusbar.StatusBar('Interpolating PSF images ................ : ', 0, nimgs) - if img.opts.quiet == False: - bar.start() - map_list = mp.parallel_map(func.eval_func_tuple, - zip(itertools.repeat(self.interp_prop), - psf_maps, itertools.repeat(psfcoords), - itertools.repeat(image.shape)), numcores=opts.ncores, - bar=bar) - if img.aperture is None: - psf_maj_int, psf_min_int, psf_pa_int, psf_ratio_int = map_list + # transpose the psf image list + xt, yt = N.transpose(tile_coord) + tr_psf_cf = N.transpose(N.array(psf_cf)) + + # interpolate the coefficients across the image. Ok, interpolate in scipy for + # irregular grids is crap. doesnt even pass through some of the points. + # for now, fit polynomial. + compress = 100.0 + x, y = N.transpose(psfcoords) + if len(x) < 3: + mylog.warning('Insufficient number of tiles to do interpolation of PSF variation') + return + + psf_coeff_interp, xgrid, ygrid = self.interp_shapcoefs(nmax, tr_psf_cf, psfcoords, image.shape, \ + compress, plot) + + psfshape = psfimages[0].shape + skip = 5 + aa = self.create_psf_grid(psf_coeff_interp, image.shape, xgrid, ygrid, skip, nmax, psfshape, \ + basis, beta, cen, totpsfimage, plot) + img.psf_images = aa + else: + if opts.psf_fwhm is None: + if ntile < 4: + mylog.warning('Insufficient number of tiles to do interpolation of PSF variation') + return 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 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 is None: - psf_maps = [psf_maj_int, psf_min_int, psf_pa_int, psf_ratio_int] + # Fit stacked PSFs with Gaussians and measure aperture fluxes + bm_pix = N.array([img.pixel_beam()[0]*fwsig, img.pixel_beam()[1]*fwsig, img.pixel_beam()[2]]) + psf_maj = N.zeros(npsf) + psf_min = N.zeros(npsf) + psf_pa = N.zeros(npsf) + if img.opts.quiet == False: + bar.start() + for i in range(ntile): + psfim = psfimages[i] + mask = N.zeros(psfim.shape, dtype=bool) + x_ax, y_ax = N.indices(psfim.shape) + maxv = N.max(psfim) + p_ini = [maxv, (psfim.shape[0]-1)/2.0*1.1, (psfim.shape[1]-1)/2.0*1.1, bm_pix[0]/fwsig*1.3, + bm_pix[1]/fwsig*1.1, bm_pix[2]*2] + para, ierr = func.fit_gaus2d(psfim, p_ini, x_ax, y_ax, mask) + ### first extent is major + if para[3] < para[4]: + para[3:5] = para[4:2:-1] + para[5] += 90 + ### clip position angle + para[5] = divmod(para[5], 180)[1] + + psf_maj[i] = para[3] + psf_min[i] = para[4] + posang = para[5] + while posang >= 180.0: + posang -= 180.0 + psf_pa[i] = posang + + if img.opts.quiet == False: + bar.increment() + bar.stop() + + # Interpolate Gaussian parameters + if img.aperture is None: + psf_maps = [psf_maj, psf_min, psf_pa, psfratio] else: - psf_maps = [psf_maj_int, psf_min_int, psf_pa_int, psf_ratio_int, psf_ratio_aper_int] + psf_maps = [psf_maj, psf_min, psf_pa, psfratio, psfratio_aper] nimgs = len(psf_maps) - bar = statusbar.StatusBar('Smoothing PSF images .................... : ', 0, nimgs) + bar = statusbar.StatusBar('Interpolating PSF images ................ : ', 0, nimgs) if img.opts.quiet == False: bar.start() map_list = mp.parallel_map(func.eval_func_tuple, - zip(itertools.repeat(self.blur_image), - psf_maps, itertools.repeat(sm_scale)), numcores=opts.ncores, + zip(itertools.repeat(self.interp_prop), + psf_maps, itertools.repeat(psfcoords), + itertools.repeat(image.shape)), numcores=opts.ncores, bar=bar) 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 - # Make sure all smoothed, interpolated images are ndarrays - psf_maj_int = N.array(psf_maj_int) - 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 is None: - psf_ratio_aper_int = N.zeros(psf_maj_int.shape, dtype=N.float32) + # Smooth if desired + 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 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] + nimgs = len(psf_maps) + bar = statusbar.StatusBar('Smoothing PSF images .................... : ', 0, nimgs) + if img.opts.quiet == False: + bar.start() + map_list = mp.parallel_map(func.eval_func_tuple, + zip(itertools.repeat(self.blur_image), + psf_maps, itertools.repeat(sm_scale)), numcores=opts.ncores, + bar=bar) + 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 + + # Make sure all smoothed, interpolated images are ndarrays + psf_maj_int = N.array(psf_maj_int) + 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 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) + + # Blank with NaNs if needed + mask = img.mask_arr + if isinstance(mask, N.ndarray): + pix_masked = N.where(mask == True) + psf_maj_int[pix_masked] = N.nan + psf_min_int[pix_masked] = N.nan + psf_pa_int[pix_masked] = N.nan + psf_ratio_int[pix_masked] = N.nan + psf_ratio_aper_int[pix_masked] = N.nan + + # Store interpolated images. The major and minor axis images are + # the sigma in units of arcsec, the PA image in units of degrees east of + # north, the ratio images in units of 1/beam. + img.psf_vary_maj_arr = psf_maj_int * img.pix2beam([1.0, 1.0, 0.0])[0] * 3600.0 # sigma in arcsec + img.psf_vary_min_arr = psf_min_int * img.pix2beam([1.0, 1.0, 0.0])[0] * 3600.0 # sigma in arcsec + img.psf_vary_pa_arr = psf_pa_int + img.psf_vary_ratio_arr = psf_ratio_int # in 1/beam + img.psf_vary_ratio_aper_arr = psf_ratio_aper_int # in 1/beam + + if opts.output_all: + func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_maj.fits', img.psf_vary_maj_arr*fwsig, img, dir) + func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_min.fits', img.psf_vary_min_arr*fwsig, img, dir) + func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_pa.fits', img.psf_vary_pa_arr, img, dir) + func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_ratio.fits', img.psf_vary_ratio_arr, img, dir) + func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_ratio_aper.fits', img.psf_vary_ratio_aper_arr, 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) + if img.opts.quiet == False: + bar2.start() + for src in img.sources: + src_pos = img.sky2pix(src.posn_sky_centroid) + src_pos_int = (int(src_pos[0]), int(src_pos[1])) + gaus_c = img.gaus2pix(src.size_sky, src.posn_sky_centroid) + if opts.psf_fwhm is None: + gaus_bm = [psf_maj_int[src_pos_int]*fwsig, psf_min_int[src_pos_int]*fwsig, psf_pa_int[src_pos_int]] else: - psf_ratio_aper_int = N.array(psf_ratio_aper_int, dtype=N.float32) - - # Blank with NaNs if needed - mask = img.mask_arr - if isinstance(mask, N.ndarray): - pix_masked = N.where(mask == True) - psf_maj_int[pix_masked] = N.nan - psf_min_int[pix_masked] = N.nan - psf_pa_int[pix_masked] = N.nan - psf_ratio_int[pix_masked] = N.nan - psf_ratio_aper_int[pix_masked] = N.nan - - # Store interpolated images. The major and minor axis images are - # the sigma in units of arcsec, the PA image in units of degrees east of - # north, the ratio images in units of 1/beam. - img.psf_vary_maj_arr = psf_maj_int * img.pix2beam([1.0, 1.0, 0.0])[0] * 3600.0 # sigma in arcsec - img.psf_vary_min_arr = psf_min_int * img.pix2beam([1.0, 1.0, 0.0])[0] * 3600.0 # sigma in arcsec - img.psf_vary_pa_arr = psf_pa_int - img.psf_vary_ratio_arr = psf_ratio_int # in 1/beam - img.psf_vary_ratio_aper_arr = psf_ratio_aper_int # in 1/beam - - if opts.output_all: - func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_maj.fits', img.psf_vary_maj_arr*fwsig, img, dir) - func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_min.fits', img.psf_vary_min_arr*fwsig, img, dir) - func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_pa.fits', img.psf_vary_pa_arr, img, dir) - func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_ratio.fits', img.psf_vary_ratio_arr, img, dir) - func.write_image_to_file(img.use_io, img.imagename + '.psf_vary_ratio_aper.fits', img.psf_vary_ratio_aper_arr, 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) - if img.opts.quiet == False: - bar2.start() - for src in img.sources: - src_pos = img.sky2pix(src.posn_sky_centroid) - src_pos_int = (int(src_pos[0]), int(src_pos[1])) - gaus_c = img.gaus2pix(src.size_sky, src.posn_sky_centroid) - if opts.psf_fwhm is None: - gaus_bm = [psf_maj_int[src_pos_int]*fwsig, psf_min_int[src_pos_int]*fwsig, psf_pa_int[src_pos_int]] - else: - # Use user-specified constant PSF instead - gaus_bm = img.beam2pix(opts.psf_fwhm) - gaus_dc, err = func.deconv2(gaus_bm, gaus_c) - src.deconv_size_sky = img.pix2gaus(gaus_dc, src_pos) - src.deconv_size_skyE = [0.0, 0.0, 0.0] - for g in src.gaussians: - gaus_c = img.gaus2pix(g.size_sky, src.posn_sky_centroid) + # Use user-specified constant PSF instead + gaus_bm = img.beam2pix(opts.psf_fwhm) gaus_dc, err = func.deconv2(gaus_bm, gaus_c) - g.deconv_size_sky = img.pix2gaus(gaus_dc, g.centre_pix) - g.deconv_size_skyE = [0.0, 0.0, 0.0] + src.deconv_size_sky = img.pix2gaus(gaus_dc, src_pos) + src.deconv_size_skyE = [0.0, 0.0, 0.0] + for g in src.gaussians: + gaus_c = img.gaus2pix(g.size_sky, src.posn_sky_centroid) + gaus_dc, err = func.deconv2(gaus_bm, gaus_c) + g.deconv_size_sky = img.pix2gaus(gaus_dc, g.centre_pix) + g.deconv_size_skyE = [0.0, 0.0, 0.0] + if img.opts.quiet == False: + bar2.spin() if img.opts.quiet == False: - bar2.spin() - if img.opts.quiet == False: - bar2.increment() - bar2.stop() - img.completed_Ops.append('psf_vary') + bar2.increment() + bar2.stop() + img.completed_Ops.append('psf_vary') ################################################################################################## @@ -421,32 +421,32 @@ def bin_and_stats_ny(self, x,y,over,ptpbin,nbin,ptplastbin,nsig): if len(x1) > 0 and len(y1) > 0: nout=100; niter=0 while nout>0 and niter<6: - med1=N.median(y1[:]) - med2=10.**(N.median(N.log10(x1[:]))) - medstd=0 # calcmedianstd.f - for j in y1: medstd += (j-med1)*(j-med1) - medstd=math.sqrt(medstd/len(y1)) # - av1=N.mean(y1); std1=func.std(y1) - av2=N.mean(x1); std2=func.std(x1) - # get_medianclip_vec2 - z=N.transpose([x1, y1]) - z1=N.transpose([n for n in z if abs(n[1]-med1)<=nsig*medstd]) - nout=len(x1)-len(z1[0]) - x1=z1[0]; y1=z1[1]; - niter+=1 + med1=N.median(y1[:]) + med2=10.**(N.median(N.log10(x1[:]))) + medstd=0 # calcmedianstd.f + for j in y1: medstd += (j-med1)*(j-med1) + medstd=math.sqrt(medstd/len(y1)) # + av1=N.mean(y1); std1=func.std(y1) + av2=N.mean(x1); std2=func.std(x1) + # get_medianclip_vec2 + z=N.transpose([x1, y1]) + z1=N.transpose([n for n in z if abs(n[1]-med1)<=nsig*medstd]) + nout=len(x1)-len(z1[0]) + x1=z1[0]; y1=z1[1]; + niter+=1 xval[i]=med2; meany[i]=av1; stdy[i]=std1; mediany[i]=med1 if stdy[nbin-1]/mediany[nbin-1] > stdy[nbin-2]/mediany[nbin-2]: - stdy[nbin-1]=stdy[nbin-2]/mediany[nbin-2]*mediany[nbin-1] + stdy[nbin-1]=stdy[nbin-2]/mediany[nbin-2]*mediany[nbin-1] return xval, meany, stdy, mediany ################################################################################################## def LM_fit(self, x, y, err, funct, order=0): if funct == func.poly: - p0=N.array([y[N.argmax(x)]] + [0]*order) + p0=N.array([y[N.argmax(x)]] + [0]*order) if funct == func.wenss_fit: - p0=N.array([y[N.argmax(x)]] + [1.]) + p0=N.array([y[N.argmax(x)]] + [1.]) res=lambda p, x, y, err: (y-funct(p, x))/err (p, flag)=leastsq(res, p0, args=(x, y, err)) return p @@ -465,7 +465,7 @@ def fit_bins_func(self, x,y,over,ptpbin,nbin,ptplastbin,nsig): # sub_size_kscli i=0 while i0.75*num: - (ptpbin, nbin, ptplastbin)=self.bindata(over,nout) # get_bins in fBDSM - (s_c,s_dm) = self.fit_bins_func(xarr,yarr,over,ptpbin,nbin,ptplastbin,nsig) # size_ksclip_wenss in fBDSM - noutold = len(xarr) - z = N.transpose([xarr, yarr, s_dm[0]+s_dm[1]*N.log10(xarr)+s_dm[2]*(N.log10(xarr)**2.), \ - N.sqrt(s_c[0]*s_c[0]+s_c[1]*s_c[1]/(xarr*xarr)) ]) - z1 = N.transpose([n for n in z if abs(n[1]-n[2])/(n[2]*n[3]) unresolved - logsnr=N.log10(snr) - dumr = N.sqrt(s_c[0]*s_c[0]+s_c[1]*s_c[1]/(snr*snr)) - med = s_dm[0]+s_dm[1]*logsnr+s_dm[2]*(logsnr*logsnr) - f_sclip[idx] = N.abs((nbeam-med)/(med*dumr)) < N.array([kappa2]*num) + xarr=N.copy(snr) + yarr=N.copy(nbeam) + niter=0; nout=num; noutold=nout*2 + while niter<10 and nout >0.75*num: + (ptpbin, nbin, ptplastbin)=self.bindata(over,nout) # get_bins in fBDSM + (s_c,s_dm) = self.fit_bins_func(xarr,yarr,over,ptpbin,nbin,ptplastbin,nsig) # size_ksclip_wenss in fBDSM + noutold = len(xarr) + z = N.transpose([xarr, yarr, s_dm[0]+s_dm[1]*N.log10(xarr)+s_dm[2]*(N.log10(xarr)**2.), \ + N.sqrt(s_c[0]*s_c[0]+s_c[1]*s_c[1]/(xarr*xarr)) ]) + z1 = N.transpose([n for n in z if abs(n[1]-n[2])/(n[2]*n[3]) unresolved + logsnr=N.log10(snr) + dumr = N.sqrt(s_c[0]*s_c[0]+s_c[1]*s_c[1]/(snr*snr)) + med = s_dm[0]+s_dm[1]*logsnr+s_dm[2]*(logsnr*logsnr) + f_sclip[idx] = N.abs((nbeam-med)/(med*dumr)) < N.array([kappa2]*num) f_s = f_sclip[0]*f_sclip[1] # Add bright sources @@ -542,7 +542,7 @@ def get_unresolved(self, g_gauls, beam, nsig, kappa2, over, bright_snr_cut=20.0, if len(bright_srcs[0]) > 0: f_s[bright_srcs] = True - # now make plots + # now make plots # if plot: # bb=[b1, b2] # pl.subplot(211+idx) @@ -648,29 +648,29 @@ def edit_vorogenlist(self, vorogenP, frac): indi = N.argsort(dist) sortdist = dist[indi] if sortdist[1] < frac * sortdist[2]: # first is the element itself - if flag[indi[1]] + flag[i] == 0: # not already deleted from other pair - tile_list.append([i, indi[1]]) - tile_coord.append((coord[i]*snrgen[i]+coord[indi[1]]*snrgen[indi[1]])/(snrgen[i]+snrgen[indi[1]])) - tile_snr.append(snrgen[i]+snrgen[indi[1]]) - flag[i] = 1 - flag[indi[1]] = 1 - else: - if len(dist) > 3: - if sortdist[1]+sortdist[2] < 2.0*frac*sortdist[3]: # for 3 close-by sources - in1=indi[1] - in2=indi[2] - if flag[in1]+flag[in2]+flag[i] == 0: # not already deleted from others - tile_list.append([i, in1, in2]) - tile_coord.append((coord[i]*snrgen[i]+coord[in1]*snrgen[in1]+coord[in2]*snrgen[in2]) \ - /(snrgen[i]+snrgen[in1]+snrgen[in2])) - tile_snr.append(snrgen[i]+snrgen[in1]+snrgen[in2]) + if flag[indi[1]] + flag[i] == 0: # not already deleted from other pair + tile_list.append([i, indi[1]]) + tile_coord.append((coord[i]*snrgen[i]+coord[indi[1]]*snrgen[indi[1]])/(snrgen[i]+snrgen[indi[1]])) + tile_snr.append(snrgen[i]+snrgen[indi[1]]) flag[i] = 1 - flag[in1] = 1 - flag[in2] = 1 - else: - tile_list.append([i]) - tile_coord.append(coord[i]) - tile_snr.append(snrgen[i]) + flag[indi[1]] = 1 + else: + if len(dist) > 3: + if sortdist[1]+sortdist[2] < 2.0*frac*sortdist[3]: # for 3 close-by sources + in1=indi[1] + in2=indi[2] + if flag[in1]+flag[in2]+flag[i] == 0: # not already deleted from others + tile_list.append([i, in1, in2]) + tile_coord.append((coord[i]*snrgen[i]+coord[in1]*snrgen[in1]+coord[in2]*snrgen[in2]) \ + /(snrgen[i]+snrgen[in1]+snrgen[in2])) + tile_snr.append(snrgen[i]+snrgen[in1]+snrgen[in2]) + flag[i] = 1 + flag[in1] = 1 + flag[in2] = 1 + else: + tile_list.append([i]) + tile_coord.append(coord[i]) + tile_snr.append(snrgen[i]) # Assign any leftover generators for i in range(len(xgen)): @@ -706,22 +706,22 @@ def pixintile(self, tilecoord, pixel, tess_method, wts, tess_sc, tess_fuzzy): """ This has routines to find out which tile a given pixel belongs to. """ if tess_method == 'roundness': - #tilenum = pytess_roundness(tilecoord, pixel, wts, tess_sc, tess_fuzzy) - print(" Not yet implemented !!!! ") - return 0 - else: - xgen, ygen = tilecoord - xgen = N.asarray(xgen) - ygen = N.asarray(ygen) - ngen = len(xgen) - i,j = pixel - dist = N.sqrt((i-xgen)*(i-xgen)+(j-ygen)*(j-ygen))/wts - minind = dist.argmin() - - if tess_sc == 's': - tilenum=minind - else: + #tilenum = pytess_roundness(tilecoord, pixel, wts, tess_sc, tess_fuzzy) print(" Not yet implemented !!!! ") + return 0 + else: + xgen, ygen = tilecoord + xgen = N.asarray(xgen) + ygen = N.asarray(ygen) + ngen = len(xgen) + i,j = pixel + dist = N.sqrt((i-xgen)*(i-xgen)+(j-ygen)*(j-ygen))/wts + minind = dist.argmin() + + if tess_sc == 's': + tilenum=minind + else: + print(" Not yet implemented !!!! ") return tilenum @@ -779,30 +779,30 @@ def edit_tile(self, ltnum, g_gauls, flag_unresolved, snrcutstack, volrank, tile_ r2t = N.zeros(ntile, dtype=int) entry = -1 for itile in range(ntile): - if ngenpertile[itile] >= ltnum: - r2t[itile] = itile - else: - pixel = tile_coord[itile] - tilenum = self.pixintile(self.trans_gaul(tile_coord_n), pixel, tess_method, wts_n, tess_sc, tess_fuzzy) - r2t[itile] = tilenum + if ngenpertile[itile] >= ltnum: + r2t[itile] = itile + else: + pixel = tile_coord[itile] + tilenum = self.pixintile(self.trans_gaul(tile_coord_n), pixel, tess_method, wts_n, tess_sc, tess_fuzzy) + r2t[itile] = tilenum for itile in range(new_n): num = N.sum(r2t == itile) if num == 0: - minarr = -999 - while minarr != itile: - arr = N.where(r2t > itile)[0] - minarr = r2t[arr].min()-1 - for i in arr: r2t[i]=r2t[i]-1 + minarr = -999 + while minarr != itile: + arr = N.where(r2t > itile)[0] + minarr = r2t[arr].min()-1 + for i in arr: r2t[i]=r2t[i]-1 n_tile_list = []; n_tile_coord = []; n_tile_snr = [] for itile in range(new_n): - ind = N.where(r2t == itile)[0]; ind1 = [] - for i in ind: ind1 = ind1 + tile_list[i] - n_tile_list.append(ind1) - snrs = N.array([tile_snr[i] for i in ind]) - coords = N.array([tile_coord[i] for i in ind]) - n_tile_snr.append(N.sum(snrs)) - n_tile_coord.append(N.sum([snrs[i]*coords[i] for i in range(len(snrs))], 0)/N.sum(snrs)) + ind = N.where(r2t == itile)[0]; ind1 = [] + for i in ind: ind1 = ind1 + tile_list[i] + n_tile_list.append(ind1) + snrs = N.array([tile_snr[i] for i in ind]) + coords = N.array([tile_coord[i] for i in ind]) + n_tile_snr.append(N.sum(snrs)) + n_tile_coord.append(N.sum([snrs[i]*coords[i] for i in range(len(snrs))], 0)/N.sum(snrs)) ngenpertile=N.zeros(new_n) for itile in range(new_n): @@ -846,13 +846,13 @@ def stackpsf(self, image, beam, g_gauls, wts, cdelt, factor): if len(cutimage.shape) == 3: cutimage=cutimage[:,:,0] if 0 not in cutimage.shape: if sum(sum(N.isnan(cutimage))) == 0: - im_shift = func.imageshift(cutimage, shift) - im_shift = im_shift/peak[isrc]*wts[isrc] - subim_shift = im_shift[cc[0]-cpsf[0]:cc[0]-cpsf[0]+psfimsize,cc[1]-cpsf[1]:cc[1]-cpsf[1]+psfimsize] - if subim_shift.shape == psfimage.shape: - # Check shapes, as they can differ if source is near edge of image. - # If they do differ, don't use that source (may be distorted). - psfimage += subim_shift + im_shift = func.imageshift(cutimage, shift) + im_shift = im_shift/peak[isrc]*wts[isrc] + subim_shift = im_shift[cc[0]-cpsf[0]:cc[0]-cpsf[0]+psfimsize,cc[1]-cpsf[1]:cc[1]-cpsf[1]+psfimsize] + if subim_shift.shape == psfimage.shape: + # Check shapes, as they can differ if source is near edge of image. + # If they do differ, don't use that source (may be distorted). + psfimage += subim_shift psfimage = psfimage/wt return psfimage @@ -879,13 +879,13 @@ def psf_in_tile(self, image, beam, g_gauls, cdelt, factor, snrcutstack, volrank, xt, yt = N.transpose(tile_coord) if plot: - pl.figure(None) - colours=['b','g','r','c','m','y','k']*(len(xt)/7+1) - pl.axis([0.0, image.shape[0], 0.0, image.shape[1]]) - pl.title('Tesselated image with tile centres and unresolved sources') - for i in range(ntile): - pl.plot([xt[i]], [yt[i]], 'D'+colours[i]) - pl.text(xt[i], yt[i], str(i)) + pl.figure(None) + colours=['b','g','r','c','m','y','k']*(len(xt)/7+1) + pl.axis([0.0, image.shape[0], 0.0, image.shape[1]]) + pl.title('Tesselated image with tile centres and unresolved sources') + for i in range(ntile): + pl.plot([xt[i]], [yt[i]], 'D'+colours[i]) + pl.text(xt[i], yt[i], str(i)) for itile in range(ntile): tile_gauls = [n for n in tr if volrank[int(round(n[2])),int(round(n[3]))]-1 \ @@ -893,11 +893,11 @@ def psf_in_tile(self, image, beam, g_gauls, cdelt, factor, snrcutstack, volrank, t_gauls = self.trans_gaul(tile_gauls) srcpertile[itile] = len(tile_gauls) if plot: - pl.plot(t_gauls[2], t_gauls[3], 'x'+'k', mew=1.3)#colours[itile]) - for i, ig in enumerate(t_gauls[2]): - xx=[xt[itile], ig] - yy=[yt[itile], t_gauls[3][i]] - pl.plot(xx,yy,'-'+colours[itile]) + pl.plot(t_gauls[2], t_gauls[3], 'x'+'k', mew=1.3)#colours[itile]) + for i, ig in enumerate(t_gauls[2]): + xx=[xt[itile], ig] + yy=[yt[itile], t_gauls[3][i]] + pl.plot(xx,yy,'-'+colours[itile]) wts = N.asarray(t_gauls[1])/N.asarray(t_gauls[8]) # wt is SNR snrpertile[itile] = sum(wts) mylog.info('PSF tile #%i (center = %i, %i): %i unresolved sources, SNR = %.1f' % @@ -932,29 +932,29 @@ def psf_in_tile(self, image, beam, g_gauls, cdelt, factor, snrcutstack, volrank, totpsfimage = totpsfimage/sum(snrpertile) if plot: - pl.imshow(N.transpose(volrank), origin='lower', interpolation='nearest'); pl.colorbar() + pl.imshow(N.transpose(volrank), origin='lower', interpolation='nearest'); pl.colorbar() if plot: - pl.figure(None) - pl.clf() - ax = pl.subplot(1,1,1) - pax = ax.get_position() - start = N.array((pax.xmin, pax.ymin)) - stop = N.array((pax.xmax, pax.ymax)) - plaxis = pl.axis([0, image.shape[0], 0, image.shape[1]]) - pl.title('Stacked psf for each tile') - for itile in range(ntile): - im=psfimages[itile] - sz=0.07 - spt = int(round(snrpertile[itile]*10))/10. - titl='n='+str(int(round(srcpertile[itile])))+'; SNR='+str(spt) - posn=[psfcoords[itile][0], psfcoords[itile][1]] - normposn=N.array(stop-start, dtype=float)/N.array(image.shape[0:2])*posn+start - a=pl.axes([normposn[0]-sz/2., normposn[1]-sz/2., sz, sz]) - pl.contour(im,15) - pl.title(titl, fontsize='small') - pl.setp(a, xticks=[], yticks=[]) - pl.show() + pl.figure(None) + pl.clf() + ax = pl.subplot(1,1,1) + pax = ax.get_position() + start = N.array((pax.xmin, pax.ymin)) + stop = N.array((pax.xmax, pax.ymax)) + plaxis = pl.axis([0, image.shape[0], 0, image.shape[1]]) + pl.title('Stacked psf for each tile') + for itile in range(ntile): + im=psfimages[itile] + sz=0.07 + spt = int(round(snrpertile[itile]*10))/10. + titl='n='+str(int(round(srcpertile[itile])))+'; SNR='+str(spt) + posn=[psfcoords[itile][0], psfcoords[itile][1]] + normposn=N.array(stop-start, dtype=float)/N.array(image.shape[0:2])*posn+start + a=pl.axes([normposn[0]-sz/2., normposn[1]-sz/2., sz, sz]) + pl.contour(im,15) + pl.title(titl, fontsize='small') + pl.setp(a, xticks=[], yticks=[]) + pl.show() return psfimages, psfcoords, totpsfimage, psfratio, psfratio_aper diff --git a/bdsf/pybdsf.py b/bdsf/pybdsf.py index fe68befc..9fda3641 100755 --- a/bdsf/pybdsf.py +++ b/bdsf/pybdsf.py @@ -716,16 +716,16 @@ class CustomPrompt(Prompts): def in_prompt_tokens(self, cli=None): - return [ - (Token.Prompt, 'BDSF ['), - (Token.PromptNum, str(self.shell.execution_count)), - (Token.Prompt, ']: '), - ] + return [ + (Token.Prompt, 'BDSF ['), + (Token.PromptNum, str(self.shell.execution_count)), + (Token.Prompt, ']: '), + ] def out_prompt_tokens(self): - return [ - (Token.OutPrompt, ''), - ] + return [ + (Token.OutPrompt, ''), + ] cfg.TerminalInteractiveShell.prompts_class = CustomPrompt except ImportError: diff --git a/bdsf/readimage.py b/bdsf/readimage.py index c70652e4..47ae9e5f 100644 --- a/bdsf/readimage.py +++ b/bdsf/readimage.py @@ -366,23 +366,23 @@ def pixel_beamarea(): found = True except: ### try see if AIPS as put the beam in HISTORY as usual - for h in hdr['HISTORY']: - # Check if h is a string or a FITS Card object (long headers are - # split into Cards as of PyFITS 3.0.4) - if not isinstance(h, str): - hstr = h.value - else: - hstr = h - if N.all(['BMAJ' in hstr, 'BMIN' in hstr, 'BPA' in hstr, 'CLEAN' in hstr]): - try: - dum, dum, dum, bmaj, dum, bmin, dum, bpa = hstr.split() - except ValueError: + for h in hdr['HISTORY']: + # Check if h is a string or a FITS Card object (long headers are + # split into Cards as of PyFITS 3.0.4) + if not isinstance(h, str): + hstr = h.value + else: + hstr = h + if N.all(['BMAJ' in hstr, 'BMIN' in hstr, 'BPA' in hstr, 'CLEAN' in hstr]): try: - dum, dum, bmaj, dum, bmin, dum, bpa, dum, dum = hstr.split() + dum, dum, dum, bmaj, dum, bmin, dum, bpa = hstr.split() except ValueError: - break - beam = (float(bmaj), float(bmin), float(bpa)) - found = True + try: + dum, dum, bmaj, dum, bmin, dum, bpa, dum, dum = hstr.split() + except ValueError: + break + beam = (float(bmaj), float(bmin), float(bpa)) + found = True if not found: raise RuntimeError("No beam information found in image header.") ### convert beam into pixels (at image center) diff --git a/bdsf/rmsimage.py b/bdsf/rmsimage.py index 20e03ed1..f08bc6cd 100644 --- a/bdsf/rmsimage.py +++ b/bdsf/rmsimage.py @@ -219,162 +219,162 @@ def __call__(self, img): map_opts = (kappa, img.rms_box, opts.spline_rank) for ipol, pol in enumerate(pols): - data = ch0_images[ipol] - mean = N.zeros(data.shape, dtype=N.float32) - rms = N.zeros(data.shape, dtype=N.float32) - if len(pols) > 1: - pol_txt = ' (' + pol + ')' - else: - pol_txt = '' - - ## calculate rms/mean maps if needed - if ((opts.rms_map is not False) or (opts.mean_map not in ['zero', 'const'])) and img.rms_box[0] > min(image.shape)/4.0: + data = ch0_images[ipol] + mean = N.zeros(data.shape, dtype=N.float32) + rms = N.zeros(data.shape, dtype=N.float32) + if len(pols) > 1: + pol_txt = ' (' + pol + ')' + else: + pol_txt = '' + + ## calculate rms/mean maps if needed + if ((opts.rms_map is not False) or (opts.mean_map not in ['zero', 'const'])) and img.rms_box[0] > min(image.shape)/4.0: # rms box is too large - just use constant rms and mean - self.output_rmsbox_size(img) - mylogger.userinfo(mylog, 'Size of rms_box larger than 1/4 of image size') - mylogger.userinfo(mylog, 'Using constant background rms and mean') - img.use_rms_map = False - img.mean_map_type = 'const' - else: - if opts.rmsmean_map_filename is not None and len(opts.rmsmean_map_filename)!=0: + self.output_rmsbox_size(img) + mylogger.userinfo(mylog, 'Size of rms_box larger than 1/4 of image size') + mylogger.userinfo(mylog, 'Using constant background rms and mean') + img.use_rms_map = False + img.mean_map_type = 'const' + else: + if opts.rmsmean_map_filename is not None and len(opts.rmsmean_map_filename)!=0: # from astropy.io import fits as pyfits - def CheckShape(A): - if len(A.shape)!=4: - raise RuntimeError("Array shape should be len 4 (nch,npol,nx,ny)") - if A.shape[0]!=1: - raise RuntimeError("Array should be single channel") - if A.shape[1]!=1: - raise RuntimeError("Array should be single pol") - - mean_fits_name,rms_fits_name=opts.rmsmean_map_filename - - mylogger.userinfo(mylog, "Skipping mean and rms image computation as external images supplied") - mylogger.userinfo(mylog, " Opening mean image: %s"%mean_fits_name) - # mean = pyfits.open(mean_fits_name, mode="readonly")[0].data - mean, hdr = read_image_from_file(mean_fits_name, img, img.indir) - CheckShape(mean); mean = mean[0,0] - mylogger.userinfo(mylog, " Opening rms image: %s"%rms_fits_name) - # rms = pyfits.open(rms_fits_name, mode="readonly")[0].data - rms, hdr = read_image_from_file(rms_fits_name, img, img.indir) - CheckShape(rms); rms = rms[0,0] - - elif (opts.rms_map is not False) or (opts.mean_map not in ['zero', 'const']): - if len(data.shape) == 2: ## 2d case - mean, rms = self.calculate_maps(img, data, mean, rms, mask, map_opts, do_adapt=do_adapt, - bright_pt_coords=isl_pos, rms_box2=img.rms_box, - logname="PyBDSM."+img.log, ncores=img.opts.ncores) - elif len(data.shape) == 3: ## 3d case - if not isinstance(mask, N.ndarray): - mask = N.zeros(data.shape[0], dtype=bool) - for i in range(data.shape[0]): - ## iterate each plane - mean, rms = self.calculate_maps(img, data[i], mean[i], rms[i], mask[i], map_opts, - do_adapt=do_adapt, bright_pt_coords=isl_pos, - rms_box2=img.rms_box, logname="PyBDSM."+img.log, - ncores=img.opts.ncores) - else: - mylog.critical('Image shape not handleable' + pol_txt) - raise RuntimeError("Can't handle array of this shape" + pol_txt) - self.output_rmsbox_size(img) - if do_adapt: - mylogger.userinfo(mylog, 'Number of sources using small scale', str(len(isl_pos))) - mylog.info('Background rms and mean images computed' + pol_txt) - - ## check if variation of rms/mean maps is significant enough: - # check_rmsmap() sets img.use_rms_map - # check_meanmap() sets img.mean_map_type - if pol == 'I': - if opts.rms_map is None and img.use_rms_map is None: - if do_adapt and len(isl_pos) > 0: - # Always use 2d map if there is at least one bright - # source and adaptive scaling is desired - img.use_rms_map = True + def CheckShape(A): + if len(A.shape)!=4: + raise RuntimeError("Array shape should be len 4 (nch,npol,nx,ny)") + if A.shape[0]!=1: + raise RuntimeError("Array should be single channel") + if A.shape[1]!=1: + raise RuntimeError("Array should be single pol") + + mean_fits_name,rms_fits_name=opts.rmsmean_map_filename + + mylogger.userinfo(mylog, "Skipping mean and rms image computation as external images supplied") + mylogger.userinfo(mylog, " Opening mean image: %s"%mean_fits_name) + # mean = pyfits.open(mean_fits_name, mode="readonly")[0].data + mean, hdr = read_image_from_file(mean_fits_name, img, img.indir) + CheckShape(mean); mean = mean[0,0] + mylogger.userinfo(mylog, " Opening rms image: %s"%rms_fits_name) + # rms = pyfits.open(rms_fits_name, mode="readonly")[0].data + rms, hdr = read_image_from_file(rms_fits_name, img, img.indir) + CheckShape(rms); rms = rms[0,0] + + elif (opts.rms_map is not False) or (opts.mean_map not in ['zero', 'const']): + if len(data.shape) == 2: ## 2d case + mean, rms = self.calculate_maps(img, data, mean, rms, mask, map_opts, do_adapt=do_adapt, + bright_pt_coords=isl_pos, rms_box2=img.rms_box, + logname="PyBDSM."+img.log, ncores=img.opts.ncores) + elif len(data.shape) == 3: ## 3d case + if not isinstance(mask, N.ndarray): + mask = N.zeros(data.shape[0], dtype=bool) + for i in range(data.shape[0]): + ## iterate each plane + mean, rms = self.calculate_maps(img, data[i], mean[i], rms[i], mask[i], map_opts, + do_adapt=do_adapt, bright_pt_coords=isl_pos, + rms_box2=img.rms_box, logname="PyBDSM."+img.log, + ncores=img.opts.ncores) + else: + mylog.critical('Image shape not handleable' + pol_txt) + raise RuntimeError("Can't handle array of this shape" + pol_txt) + self.output_rmsbox_size(img) + if do_adapt: + mylogger.userinfo(mylog, 'Number of sources using small scale', str(len(isl_pos))) + mylog.info('Background rms and mean images computed' + pol_txt) + + ## check if variation of rms/mean maps is significant enough: + # check_rmsmap() sets img.use_rms_map + # check_meanmap() sets img.mean_map_type + if pol == 'I': + if opts.rms_map is None and img.use_rms_map is None: + if do_adapt and len(isl_pos) > 0: + # Always use 2d map if there is at least one bright + # source and adaptive scaling is desired + img.use_rms_map = True + else: + self.check_rmsmap(img, rms) + 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') + else: + mylogger.userinfo(mylog, 'Using 2D map for background rms') + + if opts.mean_map == 'default' and img.mean_map_type is None: + self.check_meanmap(img, rms) + elif opts.mean_map != 'default': + img.mean_map_type = opts.mean_map + if img.mean_map_type != 'map': + mylogger.userinfo(mylog, 'Using constant background mean') else: - self.check_rmsmap(img, rms) - 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') + mylogger.userinfo(mylog, 'Using 2D map for background mean') + + ## if rms map is insignificant, or rms_map==False use const value + if img.use_rms_map is False: + if opts.rms_value is None: + rms[:] = crmss[ipol] else: - mylogger.userinfo(mylog, 'Using 2D map for background rms') - - if opts.mean_map == 'default' and img.mean_map_type is None: - self.check_meanmap(img, rms) - elif opts.mean_map != 'default': - img.mean_map_type = opts.mean_map - if img.mean_map_type != 'map': - mylogger.userinfo(mylog, 'Using constant background mean') + rms[:] = opts.rms_value + mylogger.userinfo(mylog, 'Value of background rms' + pol_txt, + '%.2e Jy/beam' % rms[0][0]) + else: + rms_min = N.nanmin(rms) + rms_max = N.nanmax(rms) + mylogger.userinfo(mylog, 'Min/max values of background rms map' + pol_txt, + '(%.2e, %.2e) Jy/beam' % (rms_min, rms_max)) + + if img.mean_map_type != 'map': + if opts.mean_map == 'zero': + val = 0.0 else: - mylogger.userinfo(mylog, 'Using 2D map for background mean') - - ## if rms map is insignificant, or rms_map==False use const value - if img.use_rms_map is False: - if opts.rms_value is None: - rms[:] = crmss[ipol] + val = img.clipped_mean + mean[:] = val + mylogger.userinfo(mylog, 'Value of background mean' + pol_txt, + str(round(val,5))+' Jy/beam') else: - rms[:] = opts.rms_value - mylogger.userinfo(mylog, 'Value of background rms' + pol_txt, - '%.2e Jy/beam' % rms[0][0]) - else: - rms_min = N.nanmin(rms) - rms_max = N.nanmax(rms) - mylogger.userinfo(mylog, 'Min/max values of background rms map' + pol_txt, - '(%.2e, %.2e) Jy/beam' % (rms_min, rms_max)) - - if img.mean_map_type != 'map': - if opts.mean_map == 'zero': - val = 0.0 + mean_min = N.nanmin(mean) + mean_max = N.nanmax(mean) + mylogger.userinfo(mylog, 'Min/max values of background mean map' + pol_txt, + '(%.2e, %.2e) Jy/beam' % (mean_min, mean_max)) + + if pol == 'I': + # Apply mask to mean_map and rms_map by setting masked values to NaN + if isinstance(mask, N.ndarray): + pix_masked = N.where(mask == True) + mean[pix_masked] = N.nan + rms[pix_masked] = N.nan + + img.mean_arr = mean + img.rms_arr = rms + + if opts.savefits_rmsim or opts.output_all: + if img.waveletimage: + resdir = img.basedir + '/wavelet/background/' + else: + resdir = img.basedir + '/background/' + if not os.path.exists(resdir): os.makedirs(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: + resdir = img.basedir + '/wavelet/background/' + else: + resdir = img.basedir + '/background/' + if not os.path.exists(resdir): os.makedirs(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: + resdir = img.basedir + '/wavelet/background/' + else: + resdir = img.basedir + '/background/' + if not os.path.exists(resdir): os.makedirs(resdir) + zero_pixels = N.where(rms <= 0.0) + rms_nonzero = rms.copy() + rms_nonzero[zero_pixels] = N.NaN + func.write_image_to_file(img.use_io, img.imagename + '.norm_I.fits', (image-mean)/rms_nonzero, img, resdir) + mylog.info('%s %s' % ('Writing ', resdir+img.imagename+'.norm_I.fits')) else: - val = img.clipped_mean - mean[:] = val - mylogger.userinfo(mylog, 'Value of background mean' + pol_txt, - str(round(val,5))+' Jy/beam') - else: - mean_min = N.nanmin(mean) - mean_max = N.nanmax(mean) - mylogger.userinfo(mylog, 'Min/max values of background mean map' + pol_txt, - '(%.2e, %.2e) Jy/beam' % (mean_min, mean_max)) - - if pol == 'I': - # Apply mask to mean_map and rms_map by setting masked values to NaN - if isinstance(mask, N.ndarray): - pix_masked = N.where(mask == True) - mean[pix_masked] = N.nan - rms[pix_masked] = N.nan - - img.mean_arr = mean - img.rms_arr = rms - - if opts.savefits_rmsim or opts.output_all: - if img.waveletimage: - resdir = img.basedir + '/wavelet/background/' - else: - resdir = img.basedir + '/background/' - if not os.path.exists(resdir): os.makedirs(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: - resdir = img.basedir + '/wavelet/background/' - else: - resdir = img.basedir + '/background/' - if not os.path.exists(resdir): os.makedirs(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: - resdir = img.basedir + '/wavelet/background/' - else: - resdir = img.basedir + '/background/' - if not os.path.exists(resdir): os.makedirs(resdir) - zero_pixels = N.where(rms <= 0.0) - rms_nonzero = rms.copy() - rms_nonzero[zero_pixels] = N.NaN - func.write_image_to_file(img.use_io, img.imagename + '.norm_I.fits', (image-mean)/rms_nonzero, img, resdir) - mylog.info('%s %s' % ('Writing ', resdir+img.imagename+'.norm_I.fits')) - else: - img.__setattr__('mean_'+pol+'_arr', mean) - img.__setattr__('rms_'+pol+'_arr', rms) + img.__setattr__('mean_'+pol+'_arr', mean) + img.__setattr__('rms_'+pol+'_arr', rms) img.completed_Ops.append('rmsimage') return img @@ -435,14 +435,14 @@ def check_meanmap(self, img, mean): # For mean map, use a higher threshold than for the rms map, as radio images # should rarely, if ever, have significant variations in the mean if stdsub > 5.0*rms_expect: - img.mean_map_type = 'map' - mylogger.userinfo(mylog, 'Variation in mean image significant') + img.mean_map_type = 'map' + mylogger.userinfo(mylog, 'Variation in mean image significant') else: - if img.confused: - img.mean_map_type = 'zero' - else: - img.mean_map_type = 'const' - mylogger.userinfo(mylog, 'Variation in mean image not significant') + if img.confused: + img.mean_map_type = 'zero' + else: + img.mean_map_type = 'const' + mylogger.userinfo(mylog, 'Variation in mean image not significant') return img @@ -468,7 +468,7 @@ def calculate_maps(self, img, data, mean, rms, mask, map_opts, do_adapt, if test: rms_ok = False if (opts.rms_box_bright is None and do_adapt) or (opts.rms_box is None and not do_adapt): - # Increase box by 20% + # Increase box by 20% if do_adapt: new_width = int(img.rms_box_bright[0]*1.2) if new_width == img.rms_box_bright[0]: @@ -781,14 +781,14 @@ def rms_mean_map(self, arr, mask=False, kappa=3, box=None, ncores=None): rms_map[co] = cr if bounds.all(): - if use_extrapolation: - ind = [-BS,arr.shape[0], -BS,arr.shape[1]] - self.for_masked(mean_map, rms_map, mask, arr, ind, - kappa, [-2, -2]) - else: - ind = [-BS,arr_pad.shape[0], -BS,arr_pad.shape[1]] - self.for_masked(mean_map, rms_map, mask_pad, arr_pad, ind, - kappa, [-1, -1]) + if use_extrapolation: + ind = [-BS,arr.shape[0], -BS,arr.shape[1]] + self.for_masked(mean_map, rms_map, mask, arr, ind, + kappa, [-2, -2]) + else: + ind = [-BS,arr_pad.shape[0], -BS,arr_pad.shape[1]] + self.for_masked(mean_map, rms_map, mask_pad, arr_pad, ind, + kappa, [-1, -1]) # Step 3: correct(extrapolate) borders of the image def correct_borders(map): @@ -913,23 +913,23 @@ def for_masked(self, mean_map, rms_map, mask, arr, ind, kappa, co): bstat = func.bstat#_cbdsm.bstat a, b, c, d = ind; i, j = co 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 - else: - pix_unmasked = N.where(mask[a:b, c:d] == False) - npix_unmasked = N.size(pix_unmasked,1) - if npix_unmasked > 20: # find clipped mean/rms - m, r, cm, cr, cnt = bstat(arr[a:b, c:d], mask[a:b, c:d], kappa) + 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 - else: - if npix_unmasked > 5: # just find simple mean/rms - cm = N.mean(arr[pix_unmasked]) - cr = N.std(arr[pix_unmasked]) - mean_map[i, j], rms_map[i, j] = cm, cr - else: # too few unmasked pixels --> set mean/rms to inf - mean_map[i, j], rms_map[i, j] = N.inf, N.inf + else: + pix_unmasked = N.where(mask[a:b, c:d] == False) + npix_unmasked = N.size(pix_unmasked,1) + if npix_unmasked > 20: # find clipped mean/rms + m, r, cm, cr, cnt = bstat(arr[a:b, c:d], mask[a:b, c:d], kappa) + if cnt > 198: cm = m; cr = r + mean_map[i, j], rms_map[i, j] = cm, cr + else: + if npix_unmasked > 5: # just find simple mean/rms + cm = N.mean(arr[pix_unmasked]) + cr = N.std(arr[pix_unmasked]) + mean_map[i, j], rms_map[i, j] = cm, cr + else: # too few unmasked pixels --> set mean/rms to inf + mean_map[i, j], rms_map[i, j] = N.inf, N.inf def for_masked_mp(self, mask, arr, ind, kappa): @@ -937,21 +937,21 @@ def for_masked_mp(self, mask, arr, ind, kappa): bstat = func.bstat #_cbdsm.bstat a, b, c, d = ind 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: - pix_unmasked = N.where(mask[a:b, c:d] == False) - npix_unmasked = N.size(pix_unmasked,1) - if npix_unmasked > 20: # find clipped mean/rms - m, r, cm, cr, cnt = bstat(arr[a:b, c:d], mask[a:b, c:d], kappa) + m, r, cm, cr, cnt = bstat(arr[a:b, c:d], mask, kappa) if cnt > 198: cm = m; cr = r - else: - if npix_unmasked > 5: # just find simple mean/rms - cm = N.mean(arr[pix_unmasked]) - cr = N.std(arr[pix_unmasked]) - else: # too few unmasked pixels --> set mean/rms to inf - cm = N.inf - cr = N.inf + else: + pix_unmasked = N.where(mask[a:b, c:d] == False) + npix_unmasked = N.size(pix_unmasked,1) + if npix_unmasked > 20: # find clipped mean/rms + m, r, cm, cr, cnt = bstat(arr[a:b, c:d], mask[a:b, c:d], kappa) + if cnt > 198: cm = m; cr = r + else: + if npix_unmasked > 5: # just find simple mean/rms + cm = N.mean(arr[pix_unmasked]) + cr = N.std(arr[pix_unmasked]) + else: # too few unmasked pixels --> set mean/rms to inf + cm = N.inf + cr = N.inf return cm, cr @@ -1023,29 +1023,29 @@ def output_rmsbox_size(self, img): do_adapt = opts.adaptive_rms_box mylog = mylogger.logging.getLogger("PyBDSM."+img.log+"RMSimage") if (opts.rms_map is not False) or (opts.mean_map not in ['zero', 'const']): - if do_adapt: - if opts.rms_box_bright is None: - mylogger.userinfo(mylog, 'Derived rms_box (box size, step size)', - '(' + str(img.rms_box_bright[0]) + ', ' + - str(img.rms_box_bright[1]) + ') pixels (small scale)') - else: - mylogger.userinfo(mylog, 'Using user-specified rms_box', - '(' + str(img.rms_box_bright[0]) + ', ' + - str(img.rms_box_bright[1]) + ') pixels (small scale)') - if opts.rms_box is None: - mylogger.userinfo(mylog, 'Derived rms_box (box size, step size)', - '(' + str(img.rms_box[0]) + ', ' + - str(img.rms_box[1]) + ') pixels (large scale)') - else: - mylogger.userinfo(mylog, 'Using user-specified rms_box', - '(' + str(img.rms_box[0]) + ', ' + - str(img.rms_box[1]) + ') pixels (large scale)') - else: - if opts.rms_box is None: - mylogger.userinfo(mylog, 'Derived rms_box (box size, step size)', - '(' + str(img.rms_box[0]) + ', ' + - str(img.rms_box[1]) + ') pixels') - else: - mylogger.userinfo(mylog, 'Using user-specified rms_box', - '(' + str(img.rms_box[0]) + ', ' + - str(img.rms_box[1]) + ') pixels') + if do_adapt: + if opts.rms_box_bright is None: + mylogger.userinfo(mylog, 'Derived rms_box (box size, step size)', + '(' + str(img.rms_box_bright[0]) + ', ' + + str(img.rms_box_bright[1]) + ') pixels (small scale)') + else: + mylogger.userinfo(mylog, 'Using user-specified rms_box', + '(' + str(img.rms_box_bright[0]) + ', ' + + str(img.rms_box_bright[1]) + ') pixels (small scale)') + if opts.rms_box is None: + mylogger.userinfo(mylog, 'Derived rms_box (box size, step size)', + '(' + str(img.rms_box[0]) + ', ' + + str(img.rms_box[1]) + ') pixels (large scale)') + else: + mylogger.userinfo(mylog, 'Using user-specified rms_box', + '(' + str(img.rms_box[0]) + ', ' + + str(img.rms_box[1]) + ') pixels (large scale)') + else: + if opts.rms_box is None: + mylogger.userinfo(mylog, 'Derived rms_box (box size, step size)', + '(' + str(img.rms_box[0]) + ', ' + + str(img.rms_box[1]) + ') pixels') + else: + mylogger.userinfo(mylog, 'Using user-specified rms_box', + '(' + str(img.rms_box[0]) + ', ' + + str(img.rms_box[1]) + ') pixels') diff --git a/bdsf/shapefit.py b/bdsf/shapefit.py index 3e932b8a..6aebe085 100644 --- a/bdsf/shapefit.py +++ b/bdsf/shapefit.py @@ -106,58 +106,58 @@ def process_island(self, isl, img, opts=None): def get_shapelet_params(self, image, mask, basis, beam_pix, fixed, ori, mode, beta=None, cen=None, nmax=None): - """ This takes as input an image, its mask (false=valid), basis="cartesian"/"polar", - fixed=(i,j,k) where i,j,k =0/1 to calculate or take as fixed for (beta, centre, nmax), - beam_pix has the beam in (pix_fwhm, pix_fwhm, deg), - beta (the scale), cen (centre of basis expansion), nmax (max order). The output - is an updated set of values of (beta, centre, nmax). If fixed is 1 and the value is not - specified as an argument, then fixed is taken as 0.""" - from math import sqrt, log, floor - from . import functions as func - import numpy as N - - if fixed[0]==1 and beta is None: fixed[0]=0 - if fixed[1]==1 and cen is None: fixed[1]=0 - if fixed[2]==1 and nmax is None: fixed[2]=0 - - if fixed[0]*fixed[1]==0: - (m1, m2, m3)=func.moment(image, mask) - - if fixed[0]==0: - try: + """ This takes as input an image, its mask (false=valid), basis="cartesian"/"polar", + fixed=(i,j,k) where i,j,k =0/1 to calculate or take as fixed for (beta, centre, nmax), + beam_pix has the beam in (pix_fwhm, pix_fwhm, deg), + beta (the scale), cen (centre of basis expansion), nmax (max order). The output + is an updated set of values of (beta, centre, nmax). If fixed is 1 and the value is not + specified as an argument, then fixed is taken as 0.""" + from math import sqrt, log, floor + from . import functions as func + import numpy as N + + if fixed[0]==1 and beta is None: fixed[0]=0 + if fixed[1]==1 and cen is None: fixed[1]=0 + if fixed[2]==1 and nmax is None: fixed[2]=0 + + if fixed[0]*fixed[1]==0: + (m1, m2, m3)=func.moment(image, mask) + + if fixed[0]==0: + try: beta = sqrt(m3[0]*m3[1])*2.0 - except ValueError: + except ValueError: beta = 0.5 - if beta == 0.0: + if beta == 0.0: beta = 0.5 - if fixed[1]==0: - cen=m2 - if fixed[2]==0: - (n, m)=image.shape - nmax=int(round(sqrt(1.0*n*n+m*m)/beam_pix[1]))-1 - nmax=min(max(nmax*2+2,10),10) # totally ad hoc - npix = N.product(image.shape)-N.sum(mask) - if nmax*nmax >= n*m : nmax = int(floor(sqrt(npix-1))) # -1 is for when n*m is a perfect square - if mode == 'fit': # make sure npara <= npix - nmax_max = int(round(0.5*(-3+sqrt(1+8*npix)))) - nmax=min(nmax, nmax_max) - - betarange=[0.5,sqrt(beta*max(n,m))] # min, max - - if fixed[1]==0: - cen=shape_findcen(image, mask, basis, beta, nmax, beam_pix) # + check_cen_shapelet - #print 'First Centre = ',cen,N.array(cen)+ori - - from time import time - t1 = time() - if fixed[0]==0: - beta, err=shape_varybeta(image, mask, basis, beta, cen, nmax, betarange, plot=False) - t2 = time() - #print 'TIME ',t2-t1, '\n' - #print 'Final Beta = ',beta, err - - if fixed[1]==0 and fixed[0]==0: - cen=shape_findcen(image, mask, basis, beta, nmax, beam_pix) # + check_cen_shapelet - #print 'Final Cen = ',N.array(cen)+ori - - return beta, cen, nmax + if fixed[1]==0: + cen=m2 + if fixed[2]==0: + (n, m)=image.shape + nmax=int(round(sqrt(1.0*n*n+m*m)/beam_pix[1]))-1 + nmax=min(max(nmax*2+2,10),10) # totally ad hoc + npix = N.product(image.shape)-N.sum(mask) + if nmax*nmax >= n*m : nmax = int(floor(sqrt(npix-1))) # -1 is for when n*m is a perfect square + if mode == 'fit': # make sure npara <= npix + nmax_max = int(round(0.5*(-3+sqrt(1+8*npix)))) + nmax=min(nmax, nmax_max) + + betarange=[0.5,sqrt(beta*max(n,m))] # min, max + + if fixed[1]==0: + cen=shape_findcen(image, mask, basis, beta, nmax, beam_pix) # + check_cen_shapelet + #print 'First Centre = ',cen,N.array(cen)+ori + + from time import time + t1 = time() + if fixed[0]==0: + beta, err=shape_varybeta(image, mask, basis, beta, cen, nmax, betarange, plot=False) + t2 = time() + #print 'TIME ',t2-t1, '\n' + #print 'Final Beta = ',beta, err + + if fixed[1]==0 and fixed[0]==0: + cen=shape_findcen(image, mask, basis, beta, nmax, beam_pix) # + check_cen_shapelet + #print 'Final Cen = ',N.array(cen)+ori + + return beta, cen, nmax diff --git a/bdsf/shapelets.py b/bdsf/shapelets.py index 88ef2278..1c6b05af 100644 --- a/bdsf/shapelets.py +++ b/bdsf/shapelets.py @@ -38,12 +38,12 @@ def decompose_shapelets(image, mask, basis, beta, centre, nmax, mode): cf[coord] = N.sum(image*B*m) if mode == 'fit': - npix = N.product(image.shape)-N.sum(mask) - npara = (nmax+1)*(nmax+2)*0.5 - cfnew = fit_shapeletbasis(image, mask, cf, Bset) - recon1 = reconstruct_shapelets(image.shape, mask, basis, beta, centre, nmax, cf) - recon2 = reconstruct_shapelets(image.shape, mask, basis, beta, centre, nmax, cfnew) - if N.std(recon2) < 1.2*N.std(recon1): cf = cfnew + npix = N.product(image.shape)-N.sum(mask) + npara = (nmax+1)*(nmax+2)*0.5 + cfnew = fit_shapeletbasis(image, mask, cf, Bset) + recon1 = reconstruct_shapelets(image.shape, mask, basis, beta, centre, nmax, cf) + recon2 = reconstruct_shapelets(image.shape, mask, basis, beta, centre, nmax, cfnew) + if N.std(recon2) < 1.2*N.std(recon1): cf = cfnew return cf @@ -185,15 +185,15 @@ def shape_findcen(image, mask, basis, beta, nmax, beam_pix): # + check_cen_shape smask2=N.array([r == 0 for r in xft2]) cen=[0.]*2 if sum(smask1) 15: # determine jmax - # Check if largest island size is - # smaller than 1/3 of image size. If so, use it to determine jmax. - min_size = min(resid.shape) - max_isl_shape = (0, 0) - for isl in img.islands: - if isl.image.shape[0] * isl.image.shape[1] > max_isl_shape[0] * max_isl_shape[1]: - max_isl_shape = isl.image.shape - if max_isl_shape != (0, 0) and min(max_isl_shape) < min(resid.shape) / 3.0: - min_size = min(max_isl_shape) * 4.0 - else: + if img.nisl == 0: + mylog.warning("No islands found. Skipping wavelet decomposition.") + img.completed_Ops.append('wavelet_atrous') + return + + mylog.info("Decomposing gaussian residual image into a-trous wavelets") + bdir = img.basedir + '/wavelet/' + if img.opts.output_all: + if not os.path.isdir(bdir): os.makedirs(bdir) + if not os.path.isdir(bdir + '/residual/'): os.makedirs(bdir + '/residual/') + if not os.path.isdir(bdir + '/model/'): os.makedirs(bdir + '/model/') + dobdsm = img.opts.atrous_bdsm_do + filter = {'tr':{'size':3, 'vec':[1. / 4, 1. / 2, 1. / 4], 'name':'Triangle'}, + 'b3':{'size':5, 'vec':[1. / 16, 1. / 4, 3. / 8, 1. / 4, 1. / 16], 'name':'B3 spline'}} + + if dobdsm: wchain, wopts = self.setpara_bdsm(img) + + n, m = img.ch0_arr.shape + + # Calculate residual image that results from normal (non-wavelet) Gaussian fitting + Op_make_residimage()(img) + resid = img.resid_gaus_arr + + lpf = img.opts.atrous_lpf + if lpf not in ['b3', 'tr']: lpf = 'b3' + jmax = img.opts.atrous_jmax + l = len(filter[lpf]['vec']) # 1st 3 is arbit and 2nd 3 is whats expected for a-trous + if jmax < 1 or jmax > 15: # determine jmax + # Check if largest island size is + # smaller than 1/3 of image size. If so, use it to determine jmax. min_size = min(resid.shape) - jmax = int(floor(log((min_size / 3.0 * 3.0 - l) / (l - 1) + 1) / log(2.0) + 1.0)) + 1 - if min_size * 0.55 <= (l + (l - 1) * (2 ** (jmax) - 1)): jmax = jmax - 1 - img.wavelet_lpf = lpf - img.wavelet_jmax = jmax - mylog.info("Using " + filter[lpf]['name'] + ' filter with J_max = ' + str(jmax)) - - img.atrous_islands = [] - img.atrous_gaussians = [] - img.atrous_sources = [] - img.atrous_opts = [] - img.resid_wavelets_arr = cp(img.resid_gaus_arr) - - im_old = img.resid_wavelets_arr - total_flux = 0.0 - ntot_wvgaus = 0 - stop_wav = False - pix_masked = N.where(N.isnan(resid) == True) - jmin = 1 - if img.opts.ncores is None: - numcores = 1 - else: - numcores = img.opts.ncores - for j in range(jmin, jmax + 1): # extra +1 is so we can do bdsm on cJ as well - mylogger.userinfo(mylog, "\nWavelet scale #" + str(j)) - im_new = self.atrous(im_old, filter[lpf]['vec'], lpf, j, numcores=numcores, use_scipy_fft=img.opts.use_scipy_fft) - im_new[pix_masked] = N.nan # since fftconvolve wont work with blanked pixels - if img.opts.atrous_sum: - w = im_new + max_isl_shape = (0, 0) + for isl in img.islands: + if isl.image.shape[0] * isl.image.shape[1] > max_isl_shape[0] * max_isl_shape[1]: + max_isl_shape = isl.image.shape + if max_isl_shape != (0, 0) and min(max_isl_shape) < min(resid.shape) / 3.0: + min_size = min(max_isl_shape) * 4.0 + else: + min_size = min(resid.shape) + jmax = int(floor(log((min_size / 3.0 * 3.0 - l) / (l - 1) + 1) / log(2.0) + 1.0)) + 1 + if min_size * 0.55 <= (l + (l - 1) * (2 ** (jmax) - 1)): jmax = jmax - 1 + img.wavelet_lpf = lpf + img.wavelet_jmax = jmax + mylog.info("Using " + filter[lpf]['name'] + ' filter with J_max = ' + str(jmax)) + + img.atrous_islands = [] + img.atrous_gaussians = [] + img.atrous_sources = [] + img.atrous_opts = [] + img.resid_wavelets_arr = cp(img.resid_gaus_arr) + + im_old = img.resid_wavelets_arr + total_flux = 0.0 + ntot_wvgaus = 0 + stop_wav = False + pix_masked = N.where(N.isnan(resid) == True) + jmin = 1 + if img.opts.ncores is None: + numcores = 1 else: - w = im_old - im_new - im_old = im_new - suffix = 'w' + repr(j) - filename = img.imagename + '.atrous.' + suffix + '.fits' - if img.opts.output_all: - 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. - if dobdsm: - wopts['filename'] = filename - wopts['basedir'] = bdir - box = img.rms_box[0] - y1 = (l + (l - 1) * (2 ** (j - 1) - 1)) - bs = max(5 * y1, box) # changed from 10 to 5 - if bs > min(n, m) / 2: - wopts['rms_map'] = False - wopts['mean_map'] = 'const' - wopts['rms_box'] = None - else: - wopts['rms_box'] = (bs, bs/3) - if hasattr(img, '_adapt_rms_isl_pos'): - bs_bright = max(5 * y1, img.rms_box_bright[0]) - if bs_bright < bs/1.5: - wopts['adaptive_rms_box'] = True - wopts['rms_box_bright'] = (bs_bright, bs_bright/3) + numcores = img.opts.ncores + for j in range(jmin, jmax + 1): # extra +1 is so we can do bdsm on cJ as well + mylogger.userinfo(mylog, "\nWavelet scale #" + str(j)) + im_new = self.atrous(im_old, filter[lpf]['vec'], lpf, j, numcores=numcores, use_scipy_fft=img.opts.use_scipy_fft) + im_new[pix_masked] = N.nan # since fftconvolve wont work with blanked pixels + if img.opts.atrous_sum: + w = im_new + else: + w = im_old - im_new + im_old = im_new + suffix = 'w' + repr(j) + filename = img.imagename + '.atrous.' + suffix + '.fits' + if img.opts.output_all: + 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. + if dobdsm: + wopts['filename'] = filename + wopts['basedir'] = bdir + box = img.rms_box[0] + y1 = (l + (l - 1) * (2 ** (j - 1) - 1)) + bs = max(5 * y1, box) # changed from 10 to 5 + if bs > min(n, m) / 2: + wopts['rms_map'] = False + wopts['mean_map'] = 'const' + wopts['rms_box'] = None else: - wopts['adaptive_rms_box'] = False - if j <= 3: - wopts['ini_gausfit'] = 'default' - else: - wopts['ini_gausfit'] = 'nobeam' - wid = (l + (l - 1) * (2 ** (j - 1) - 1))# / 3.0 - b1, b2 = img.pixel_beam()[0:2] - b1 = b1 * fwsig - b2 = b2 * fwsig - cdelt = img.wcs_obj.acdelt[:2] - - wimg = Image(wopts) - wimg.beam = (sqrt(wid * wid + b1 * b1) * cdelt[0] * 2.0, sqrt(wid * wid + b2 * b2) * cdelt[1] * 2.0, 0.0) - wimg.orig_beam = img.beam - wimg.pixel_beam = img.pixel_beam - wimg.pixel_beamarea = img.pixel_beamarea - wimg.log = 'Wavelet.' - wimg.basedir = img.basedir - wimg.extraparams['bbsprefix'] = suffix - wimg.extraparams['bbsname'] = img.imagename + '.wavelet' - wimg.extraparams['bbsappend'] = True - wimg.bbspatchnum = img.bbspatchnum - wimg.waveletimage = True - wimg.j = j - wimg.indir = img.indir - if hasattr(img, '_adapt_rms_isl_pos'): - wimg._adapt_rms_isl_pos = img._adapt_rms_isl_pos - - - self.init_image_simple(wimg, img, w, '.atrous.' + suffix) - for op in wchain: - op(wimg) - gc.collect() - if isinstance(op, Op_islands) and img.opts.atrous_orig_isl: - if wimg.nisl > 0: - - # Find islands that do not share any pixels with - # islands in original ch0 image. - good_isl = [] - - # Make original rank image boolean; rank counts from 0, with -1 being - # outside any island - orig_rankim_bool = N.array(img.pyrank + 1, dtype = bool) - - # Multiply rank images - old_islands = orig_rankim_bool * (wimg.pyrank + 1) - 1 - - # Exclude islands that don't overlap with a ch0 island. - valid_ids = set(old_islands.flatten()) - for idx, wvisl in enumerate(wimg.islands): - if idx in valid_ids: - wvisl.valid = True - good_isl.append(wvisl) + wopts['rms_box'] = (bs, bs/3) + if hasattr(img, '_adapt_rms_isl_pos'): + bs_bright = max(5 * y1, img.rms_box_bright[0]) + if bs_bright < bs/1.5: + wopts['adaptive_rms_box'] = True + wopts['rms_box_bright'] = (bs_bright, bs_bright/3) else: - wvisl.valid = False - - wimg.islands = good_isl - wimg.nisl = len(good_isl) - mylogger.userinfo(mylog, "Number of islands found", '%i' % - wimg.nisl) - - # Renumber islands: - for wvindx, wvisl in enumerate(wimg.islands): - wvisl.island_id = wvindx - - if isinstance(op, Op_gausfit): - # If opts.atrous_orig_isl then exclude Gaussians outside of - # the original ch0 islands - nwvgaus = 0 - if img.opts.atrous_orig_isl: - gaul = wimg.gaussians - tot_flux = 0.0 - - if img.ngaus == 0: - gaus_id = -1 - else: - gaus_id = img.gaussians[-1].gaus_num - for g in gaul: - if not hasattr(g, 'valid'): - g.valid = False - if not g.valid: - try: - isl_id = img.pyrank[int(g.centre_pix[0] + 1), int(g.centre_pix[1] + 1)] - except IndexError: - isl_id = -1 - if isl_id >= 0: - isl = img.islands[isl_id] - gcenter = (int(g.centre_pix[0] - isl.origin[0]), - int(g.centre_pix[1] - isl.origin[1])) - if not isl.mask_active[gcenter]: - gaus_id += 1 - gcp = Gaussian(img, g.parameters[:], isl.island_id, gaus_id) - gcp.gaus_num = gaus_id - gcp.wisland_id = g.island_id - gcp.jlevel = j - g.valid = True - isl.gaul.append(gcp) - isl.ngaus += 1 - img.gaussians.append(gcp) - nwvgaus += 1 - tot_flux += gcp.total_flux - else: - g.valid = False - g.jlevel = 0 - else: - g.valid = False - g.jlevel = 0 - vg = [] - for g in wimg.gaussians: - if g.valid: - vg.append(g) - wimg.gaussians = vg - mylogger.userinfo(mylog, "Number of valid wavelet Gaussians", str(nwvgaus)) - else: - # Keep all Gaussians and merge islands that overlap - tot_flux = check_islands_for_overlap(img, wimg) - - # Now renumber the islands and adjust the rank image before going to next wavelet image - renumber_islands(img) - - total_flux += tot_flux - if img.opts.interactive and has_pl: - dc = '\033[34;1m' - nc = '\033[0m' - print(dc + '--> Displaying islands and rms image...' + nc) - if max(wimg.ch0_arr.shape) > 4096: - print(dc + '--> Image is large. Showing islands only.' + nc) - wimg.show_fit(rms_image=False, mean_image=False, ch0_image=False, - ch0_islands=True, gresid_image=False, sresid_image=False, - gmodel_image=False, smodel_image=False, pyramid_srcs=False) - else: - wimg.show_fit() - prompt = dc + "Press enter to continue or 'q' stop fitting wavelet images : " + nc - answ = raw_input_no_history(prompt) - while answ != '': - if answ == 'q': - img.wavelet_jmax = j - stop_wav = True - break - answ = raw_input_no_history(prompt) - if len(wimg.gaussians) > 0: - img.resid_wavelets_arr = self.subtract_wvgaus(img.opts, img.resid_wavelets_arr, wimg.gaussians, wimg.islands) - if img.opts.atrous_sum: - im_old = self.subtract_wvgaus(img.opts, im_old, wimg.gaussians, wimg.islands) - 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 - pyrank[tuple(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 - - img.ngaus += ntot_wvgaus - 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('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('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('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') + wopts['adaptive_rms_box'] = False + if j <= 3: + wopts['ini_gausfit'] = 'default' + else: + wopts['ini_gausfit'] = 'nobeam' + wid = (l + (l - 1) * (2 ** (j - 1) - 1))# / 3.0 + b1, b2 = img.pixel_beam()[0:2] + b1 = b1 * fwsig + b2 = b2 * fwsig + cdelt = img.wcs_obj.acdelt[:2] + + wimg = Image(wopts) + wimg.beam = (sqrt(wid * wid + b1 * b1) * cdelt[0] * 2.0, sqrt(wid * wid + b2 * b2) * cdelt[1] * 2.0, 0.0) + wimg.orig_beam = img.beam + wimg.pixel_beam = img.pixel_beam + wimg.pixel_beamarea = img.pixel_beamarea + wimg.log = 'Wavelet.' + wimg.basedir = img.basedir + wimg.extraparams['bbsprefix'] = suffix + wimg.extraparams['bbsname'] = img.imagename + '.wavelet' + wimg.extraparams['bbsappend'] = True + wimg.bbspatchnum = img.bbspatchnum + wimg.waveletimage = True + wimg.j = j + wimg.indir = img.indir + if hasattr(img, '_adapt_rms_isl_pos'): + wimg._adapt_rms_isl_pos = img._adapt_rms_isl_pos + + + self.init_image_simple(wimg, img, w, '.atrous.' + suffix) + for op in wchain: + op(wimg) + gc.collect() + if isinstance(op, Op_islands) and img.opts.atrous_orig_isl: + if wimg.nisl > 0: + + # Find islands that do not share any pixels with + # islands in original ch0 image. + good_isl = [] + + # Make original rank image boolean; rank counts from 0, with -1 being + # outside any island + orig_rankim_bool = N.array(img.pyrank + 1, dtype = bool) + + # Multiply rank images + old_islands = orig_rankim_bool * (wimg.pyrank + 1) - 1 + + # Exclude islands that don't overlap with a ch0 island. + valid_ids = set(old_islands.flatten()) + for idx, wvisl in enumerate(wimg.islands): + if idx in valid_ids: + wvisl.valid = True + good_isl.append(wvisl) + else: + wvisl.valid = False + + wimg.islands = good_isl + wimg.nisl = len(good_isl) + mylogger.userinfo(mylog, "Number of islands found", '%i' % + wimg.nisl) + + # Renumber islands: + for wvindx, wvisl in enumerate(wimg.islands): + wvisl.island_id = wvindx + + if isinstance(op, Op_gausfit): + # If opts.atrous_orig_isl then exclude Gaussians outside of + # the original ch0 islands + nwvgaus = 0 + if img.opts.atrous_orig_isl: + gaul = wimg.gaussians + tot_flux = 0.0 + + if img.ngaus == 0: + gaus_id = -1 + else: + gaus_id = img.gaussians[-1].gaus_num + for g in gaul: + if not hasattr(g, 'valid'): + g.valid = False + if not g.valid: + try: + isl_id = img.pyrank[int(g.centre_pix[0] + 1), int(g.centre_pix[1] + 1)] + except IndexError: + isl_id = -1 + if isl_id >= 0: + isl = img.islands[isl_id] + gcenter = (int(g.centre_pix[0] - isl.origin[0]), + int(g.centre_pix[1] - isl.origin[1])) + if not isl.mask_active[gcenter]: + gaus_id += 1 + gcp = Gaussian(img, g.parameters[:], isl.island_id, gaus_id) + gcp.gaus_num = gaus_id + gcp.wisland_id = g.island_id + gcp.jlevel = j + g.valid = True + isl.gaul.append(gcp) + isl.ngaus += 1 + img.gaussians.append(gcp) + nwvgaus += 1 + tot_flux += gcp.total_flux + else: + g.valid = False + g.jlevel = 0 + else: + g.valid = False + g.jlevel = 0 + vg = [] + for g in wimg.gaussians: + if g.valid: + vg.append(g) + wimg.gaussians = vg + mylogger.userinfo(mylog, "Number of valid wavelet Gaussians", str(nwvgaus)) + else: + # Keep all Gaussians and merge islands that overlap + tot_flux = check_islands_for_overlap(img, wimg) + + # Now renumber the islands and adjust the rank image before going to next wavelet image + renumber_islands(img) + + total_flux += tot_flux + if img.opts.interactive and has_pl: + dc = '\033[34;1m' + nc = '\033[0m' + print(dc + '--> Displaying islands and rms image...' + nc) + if max(wimg.ch0_arr.shape) > 4096: + print(dc + '--> Image is large. Showing islands only.' + nc) + wimg.show_fit(rms_image=False, mean_image=False, ch0_image=False, + ch0_islands=True, gresid_image=False, sresid_image=False, + gmodel_image=False, smodel_image=False, pyramid_srcs=False) + else: + wimg.show_fit() + prompt = dc + "Press enter to continue or 'q' stop fitting wavelet images : " + nc + answ = raw_input_no_history(prompt) + while answ != '': + if answ == 'q': + img.wavelet_jmax = j + stop_wav = True + break + answ = raw_input_no_history(prompt) + if len(wimg.gaussians) > 0: + img.resid_wavelets_arr = self.subtract_wvgaus(img.opts, img.resid_wavelets_arr, wimg.gaussians, wimg.islands) + if img.opts.atrous_sum: + im_old = self.subtract_wvgaus(img.opts, im_old, wimg.gaussians, wimg.islands) + 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 + pyrank[tuple(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 + + img.ngaus += ntot_wvgaus + 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('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('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('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') ####################################################################################################### @@ -321,8 +321,8 @@ def atrous(self, image, filtvec, lpf, j, numcores=1, use_scipy_fft=True): ff = filtvec[:] for i in range(1, len(filtvec)): - ii = 1 + (2 ** (j - 1)) * (i - 1) - ff[ii:ii] = [0] * (2 ** (j - 1) - 1) + ii = 1 + (2 ** (j - 1)) * (i - 1) + ff[ii:ii] = [0] * (2 ** (j - 1) - 1) kern = N.outer(ff, ff) unmasked = N.nan_to_num(image) if use_scipy_fft: @@ -376,10 +376,10 @@ def setpara_bdsm(self, img): ops = [] for op in chain: - if isinstance(op, type): - ops.append(op()) - else: - ops.append(op) + if isinstance(op, type): + ops.append(op()) + else: + ops.append(op) return ops, opts @@ -421,18 +421,18 @@ def subtract_wvgaus(self, opts, residim, gaussians, islands): thresh = opts.fittedimage_clip for g in gaussians: - if g.valid: - C1, C2 = g.centre_pix - if hasattr(g, 'wisland_id'): - isl = islands[g.wisland_id] - else: - isl = islands[g.island_id] - b = opp.find_bbox(dummy, thresh * isl.rms, g) - bbox = N.s_[max(0, int(C1 - b)):min(shape[0], int(C1 + b + 1)), - max(0, int(C2 - b)):min(shape[1], int(C2 + b + 1))] - x_ax, y_ax = N.mgrid[bbox] - ffimg = func.gaussian_fcn(g, x_ax, y_ax) - residim[bbox] = residim[bbox] - ffimg + if g.valid: + C1, C2 = g.centre_pix + if hasattr(g, 'wisland_id'): + isl = islands[g.wisland_id] + else: + isl = islands[g.island_id] + b = opp.find_bbox(dummy, thresh * isl.rms, g) + bbox = N.s_[max(0, int(C1 - b)):min(shape[0], int(C1 + b + 1)), + max(0, int(C2 - b)):min(shape[1], int(C2 + b + 1))] + x_ax, y_ax = N.mgrid[bbox] + ffimg = func.gaussian_fcn(g, x_ax, y_ax) + residim[bbox] = residim[bbox] - ffimg return residim @@ -446,26 +446,26 @@ def morphfilter_pyramid(self, img, bdir): lpyr = [] img.npyrsrc = -1 if len(ind) > 0 : - for i in ind: - isls = img.atrous_islands[i] - for isl in isls: - if i != ind[0]: - dumr = [] - for pyrsrc in lpyr: - belongs = pyrsrc.belongs(img, isl) - if belongs: dumr.append(pyrsrc.pyr_id) - #if len(dumr) > 1: - # raise RuntimeError("Source in lower wavelet level belongs to more than one higher level.") - if len(dumr) == 1: - dumr = dumr[0] - pyrsrc = lpyr[dumr] - pyrsrc.add_level(img, i, isl) - else: - pyrsrc = Pyramid_source(img, isl, i) - lpyr.append(pyrsrc) - else: - pyrsrc = Pyramid_source(img, isl, i) - lpyr.append(pyrsrc) + for i in ind: + isls = img.atrous_islands[i] + for isl in isls: + if i != ind[0]: + dumr = [] + for pyrsrc in lpyr: + belongs = pyrsrc.belongs(img, isl) + if belongs: dumr.append(pyrsrc.pyr_id) + #if len(dumr) > 1: + # raise RuntimeError("Source in lower wavelet level belongs to more than one higher level.") + if len(dumr) == 1: + dumr = dumr[0] + pyrsrc = lpyr[dumr] + pyrsrc.add_level(img, i, isl) + else: + pyrsrc = Pyramid_source(img, isl, i) + lpyr.append(pyrsrc) + else: + pyrsrc = Pyramid_source(img, isl, i) + lpyr.append(pyrsrc) img.pyrsrcs = lpyr if img.opts.plot_pyramid and has_pl: @@ -475,45 +475,45 @@ def morphfilter_pyramid(self, img, bdir): colours = ['r', 'g', 'b', 'c', 'm', 'y', 'k'] sh = img.ch0_arr.shape for pyr in img.pyrsrcs: - for iisl, isl in enumerate(pyr.islands): - jj = pyr.jlevels[iisl] - col = colours[pyr.pyr_id % 7] - pl.subplot(a, b, jj) - ind = N.where(~isl.mask_active) - pl.plot(ind[0] + isl.origin[0], ind[1] + isl.origin[1], '.', color = col) - pl.axis([0.0, sh[0], 0.0, sh[1]]) - pl.title('J = ' + str(jj)) + for iisl, isl in enumerate(pyr.islands): + jj = pyr.jlevels[iisl] + col = colours[pyr.pyr_id % 7] + pl.subplot(a, b, jj) + ind = N.where(~isl.mask_active) + pl.plot(ind[0] + isl.origin[0], ind[1] + isl.origin[1], '.', color = col) + pl.axis([0.0, sh[0], 0.0, sh[1]]) + pl.title('J = ' + str(jj)) pl.savefig(bdir + img.imagename + '.pybdsf.atrous.pyramidsrc.png') ####################################################################################################### class Pyramid_source(object): - """ Pyramid_source is a source constructed out of multiple wavelet transform images. """ - - def __init__(self, img, island, level0): - img.npyrsrc = img.npyrsrc + 1 - self.pyr_id = img.npyrsrc - self.islands = [island] - self.jlevels = [level0] - - def belongs(self, img, isl): - from . import functions as func - # get centroid of island (as integer) - mom = func.momanalmask_gaus(isl.image, isl.mask_active, 0, 1.0, False) - cen = N.array(mom[1:3]) + isl.origin - belong = False - # check if lies within any island of self - for i, pyrisl in enumerate(self.islands): - if N.sum([pyrisl.bbox[j].start <= cen[j] < pyrisl.bbox[j].stop for j in range(2)]) == 2: + """ Pyramid_source is a source constructed out of multiple wavelet transform images. """ + + def __init__(self, img, island, level0): + img.npyrsrc = img.npyrsrc + 1 + self.pyr_id = img.npyrsrc + self.islands = [island] + self.jlevels = [level0] + + def belongs(self, img, isl): + from . import functions as func + # get centroid of island (as integer) + mom = func.momanalmask_gaus(isl.image, isl.mask_active, 0, 1.0, False) + cen = N.array(mom[1:3]) + isl.origin + belong = False + # check if lies within any island of self + for i, pyrisl in enumerate(self.islands): + if N.sum([pyrisl.bbox[j].start <= cen[j] < pyrisl.bbox[j].stop for j in range(2)]) == 2: pix = tuple([cen[j] - pyrisl.origin[j] for j in range(2)]) if not pyrisl.mask_active[pix]: - belong = True + belong = True - return belong + return belong - def add_level(self, img, level, isl): - self.islands.append(isl) - self.jlevels.append(level + 1) + def add_level(self, img, level, isl): + self.islands.append(isl) + self.jlevels.append(level + 1) Image.pyrsrcs = List(tInstance(Pyramid_source), doc = "List of Pyramidal sources") @@ -529,7 +529,7 @@ def fftconvolve(in1, in2, mode="full", pad_to_power_of_two=True, numcores=1): size = s1 + s2 - 1 if pad_to_power_of_two: - # Use 2**n-sized FFT; it might improve performance + # Use 2**n-sized FFT; it might improve performance fsize = 2 ** N.ceil(N.log2(size)) else: # Padding to a power of two might degrade performance, too diff --git a/natgrid/Test/test_natgrid.py b/natgrid/Test/test_natgrid.py index 3f7b6f87..af021f71 100755 --- a/natgrid/Test/test_natgrid.py +++ b/natgrid/Test/test_natgrid.py @@ -1,26 +1,26 @@ # Adapted for numpy/ma/cdms2 by convertcdms.py """Documentation for module natgridtest: an automatic test for natgrid, an interface to the ngmath NATGRID - - TESTING - - Typing - + + TESTING + + Typing + cdat natgridtest.py - - generates some testing of the natgridmodule using analytical functions as fields. It also writes a - hard copy of the documentation to the file natgridmodule.doc and a copy of the information describing + + generates some testing of the natgridmodule using analytical functions as fields. It also writes a + hard copy of the documentation to the file natgridmodule.doc and a copy of the information describing the nature of the tests to test.asc. For the single and the double precision interpolations from randomly spaced data to a rectangular grid on a sphere, the numerical results are written to netCDF files if there is access to the module cdms. - + DOCUMENTATION - + Without conducting the tests, documentation written to the file natgridmodule.doc can be produced after - importing the natgridtest module by typing - - natgridtest.document() - + importing the natgridtest module by typing + + natgridtest.document() + """ import sys, numpy, math, random, nat, natgridmodule @@ -34,11 +34,11 @@ def document(): #------------------------------------------------------------------------------- # - # purpose: 'document' writes documentation for the user to a file + # purpose: 'document' writes documentation for the user to a file # # usage: import natgridtest # natgridtest.document() - # + # # passed : nothing # # returned: nothing @@ -49,9 +49,9 @@ def document(): std = sys.stdout # save sys.stout to allow reassigning later sys.stdout = open( 'natgridmodule.doc', 'w') - print '**********************************************************************************************\n' + print '**********************************************************************************************\n' print '*************************** Overview of the CDAT interface to natgrid ************************\n' - print '**********************************************************************************************\n' + print '**********************************************************************************************\n' print nat.__doc__ print print @@ -69,15 +69,15 @@ def sendOutput(msg, value = None, screen = 'no'): #------------------------------------------------------------------------------ # # purpose: send a message and optionally a value a file and if screen is not 'no' - # send the same thing to the screen - # + # send the same thing to the screen + # # usage: sendOutput(msg, value = number, screen = 'yes') - # + # # passed : msg - the string to write to the output media # value - a number - # screen - a string set to something different from 'no' if the output also + # screen - a string set to something different from 'no' if the output also # goes to the screen - # + # # returned: None # #------------------------------------------------------------------------------ @@ -95,10 +95,10 @@ def sendOutput(msg, value = None, screen = 'no'): # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # ++++++++++++++++++++++++++++++++++ Autotest Calls ++++++++++++++++++++++++++++++++++ # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - + def runtests(): #----------------------------------------------------------------------------- - # + # # purpose: call test cases # #----------------------------------------------------------------------------- @@ -108,22 +108,22 @@ def runtests(): testError = 0 for n in range(1,8): - err = choose(n) - if err != 0: - #print 'test number with error :',n,err - testError = testError + 1 + err = choose(n) + if err != 0: + #print 'test number with error :',n,err + testError = testError + 1 return testError def choose(case): #------------------------------------------------------------------------------- - # - # purpose: check out natgrid # - # case 1: a simple 2D interpolation using y32 -- single precision + # purpose: check out natgrid + # + # case 1: a simple 2D interpolation using y32 -- single precision #--------------------------------------------------------------------------------- err = 0 - if case == 1: + if case == 1: sendOutput('\n******* natural neighbor linear interpolation -- single precision *****\n') # array dimensions @@ -131,7 +131,7 @@ def choose(case): nxo = 21 nyo = 21 - # input arrays and data + # input arrays and data xiList = [0.00, 1.00, 0.00, 1.00, 0.40, 0.75] yiList = [0.00, 0.00, 1.00, 1.00, 0.20, 0.65] @@ -141,17 +141,17 @@ def choose(case): yi = numpy.array(yiList, numpy.float32) dataIn = numpy.array(dataInList, numpy.float32) - # output array + # output array xo = uniformGrid(nxo, 1.0, 0.0) yo = uniformGrid(nyo, 1.0, 0.0) - r = nat.Natgrid(xi, yi, xo, yo) + r = nat.Natgrid(xi, yi, xo, yo) - dataOut = r.rgrd(dataIn) + dataOut = r.rgrd(dataIn) sendOutput('*** writing single precision linear interpolation test case to the netCDF file SingleLinearRegrid.nc') - write1D_4DField('SingleLinearRegrid', dataOut, xo, yo) + write1D_4DField('SingleLinearRegrid', dataOut, xo, yo) dataCheck = storedAnswers('linearRegrid') dataCheck = numpy.reshape(dataCheck, (nxo,nyo)) @@ -166,7 +166,7 @@ def choose(case): return err - elif case == 2: + elif case == 2: sendOutput('\n******* natural neighbor linear interpolation -- double precision *****\n') # array dimensions @@ -174,7 +174,7 @@ def choose(case): nxo = 21 nyo = 21 - # input arrays and data + # input arrays and data xiList = [0.00, 1.00, 0.00, 1.00, 0.40, 0.75] yiList = [0.00, 0.00, 1.00, 1.00, 0.20, 0.65] @@ -184,7 +184,7 @@ def choose(case): yi = numpy.array(yiList, numpy.float64) dataIn = numpy.array(dataInList, numpy.float64) - # output array + # output array xo = uniformGrid(nxo, 1.0, 0.0) yo = uniformGrid(nyo, 1.0, 0.0) @@ -192,16 +192,16 @@ def choose(case): xo = xo.astype(numpy.float64) yo = yo.astype(numpy.float64) - r = nat.Natgrid(xi, yi, xo, yo) + r = nat.Natgrid(xi, yi, xo, yo) - dataOut = r.rgrd(dataIn) + dataOut = r.rgrd(dataIn) xo = xo.astype(numpy.float32) # convert back to single precision yo = yo.astype(numpy.float32) dataOut = dataOut.astype(numpy.float32) sendOutput('*** writing double precision linear interpolation test case to the netCDF file DoubleLinearRegrid.nc') - write1D_4DField('DoubleLinearRegrid', dataOut, xo, yo) + write1D_4DField('DoubleLinearRegrid', dataOut, xo, yo) dataCheck = storedAnswers('linearRegrid') dataCheck = numpy.reshape(dataCheck, (nxo,nyo)) @@ -216,7 +216,7 @@ def choose(case): return err - elif case == 3: + elif case == 3: sendOutput('\n******* natural neighbor nonlinear interpolation -- single precision *****\n') # array dimensions @@ -224,7 +224,7 @@ def choose(case): nxo = 21 nyo = 21 - # input arrays and data + # input arrays and data xiList = [0.00, 1.00, 0.00, 1.00, 0.40, 0.75] yiList = [0.00, 0.00, 1.00, 1.00, 0.20, 0.65] @@ -234,18 +234,18 @@ def choose(case): yi = numpy.array(yiList, numpy.float32) dataIn = numpy.array(dataInList, numpy.float32) - # output array + # output array xo = uniformGrid(nxo, 1.0, 0.0) yo = uniformGrid(nyo, 1.0, 0.0) - r = nat.Natgrid(xi, yi, xo, yo) + r = nat.Natgrid(xi, yi, xo, yo) r.igr = 1 # choose nonlinear interpolation - dataOut = r.rgrd(dataIn) + dataOut = r.rgrd(dataIn) sendOutput('*** writing single precision nonlinear interpolation test case to the netCDF file SingleNonlinearRegrid.nc') - write1D_4DField('SingleNonlinearRegrid', dataOut, xo, yo) + write1D_4DField('SingleNonlinearRegrid', dataOut, xo, yo) dataCheck = storedAnswers('nonlinearRegrid') dataCheck = numpy.reshape(dataCheck, (nxo,nyo)) @@ -260,7 +260,7 @@ def choose(case): return err - elif case == 4: + elif case == 4: sendOutput('\n******* natural neighbor nonlinear interpolation -- double precision *****\n') # array dimensions @@ -268,7 +268,7 @@ def choose(case): nxo = 21 nyo = 21 - # input arrays and data + # input arrays and data xiList = [0.00, 1.00, 0.00, 1.00, 0.40, 0.75] yiList = [0.00, 0.00, 1.00, 1.00, 0.20, 0.65] @@ -278,7 +278,7 @@ def choose(case): yi = numpy.array(yiList, numpy.float64) dataIn = numpy.array(dataInList, numpy.float64) - # output array + # output array xo = uniformGrid(nxo, 1.0, 0.0) yo = uniformGrid(nyo, 1.0, 0.0) @@ -286,17 +286,17 @@ def choose(case): xo = xo.astype(numpy.float64) yo = yo.astype(numpy.float64) - r = nat.Natgrid(xi, yi, xo, yo) + r = nat.Natgrid(xi, yi, xo, yo) r.igr = 1 # choose nonlinear interpolation - dataOut = r.rgrd(dataIn) + dataOut = r.rgrd(dataIn) xo = xo.astype(numpy.float32) # convert back to single precision yo = yo.astype(numpy.float32) dataOut = dataOut.astype(numpy.float32) sendOutput('*** writing double precision nonlinear interpolation test case to the netCDF file DoubleNonlinearRegrid.nc') - write1D_4DField('DoubleNonlinearRegrid', dataOut, xo, yo) + write1D_4DField('DoubleNonlinearRegrid', dataOut, xo, yo) dataCheck = storedAnswers('nonlinearRegrid') dataCheck = numpy.reshape(dataCheck, (nxo,nyo)) @@ -312,7 +312,7 @@ def choose(case): return err - elif case == 5: + elif case == 5: sendOutput('\n******* interpolation and computation of aspects and slopes -- single precision *******\n') # array dimensions @@ -320,7 +320,7 @@ def choose(case): nxo = 21 nyo = 21 - # input array and data + # input array and data xisort, xi = randomGrid(ni, 1.2, -0.2) # xisort has random numbers monotonically increasing yisort, yi = randomGrid(ni, 1.2, -0.2) @@ -329,18 +329,18 @@ def choose(case): for i in range(ni): dataIn[i] = (xi[i] - 0.25)*(xi[i] - 0.25) + (yi[i] - 0.50)*(yi[i] - 0.50) - # output array + # output array xo = uniformGrid(nxo, 1.0, 0.0) yo = uniformGrid(nyo, 1.0, 0.0) - r = nat.Natgrid(xi, yi, xo, yo) + r = nat.Natgrid(xi, yi, xo, yo) - dataOut, aspect, slope = r.rgrd(dataIn, aspectSlope = 'yes') + dataOut, aspect, slope = r.rgrd(dataIn, aspectSlope = 'yes') sendOutput('*** writing single precision linear interpolation test case to the netCDF file AspectSlopeRegrid.nc') - write1D_4DField('AspectSlopeRegrid', dataOut, xo, yo) + write1D_4DField('AspectSlopeRegrid', dataOut, xo, yo) # Calculate the exact answer dataCheck = numpy.zeros((nxo, nyo), numpy.float32) @@ -349,7 +349,7 @@ def choose(case): dataCheck[i,j] = (xo[i] - 0.25)*(xo[i] - 0.25) + (yo[j] - 0.50)*(yo[j] - 0.50) sendOutput('*** writing exact answer to single precision interpolation test case to the netCDF file AspectSlopeExact.nc') - write1D_4DField('AspectSlopeExact', dataOut, xo, yo) + write1D_4DField('AspectSlopeExact', dataOut, xo, yo) error = rmserror(dataOut, dataCheck) # find the rms error sendOutput('\n******* compare results\n') @@ -368,15 +368,15 @@ def choose(case): sendOutput('*** writing the cosine of the aspect to xaspect.nc') sendOutput('*** writing the sine of the aspect to yaspect.nc') - write1D_4DField('xaspect', u, xo, yo) - write1D_4DField('yaspect', v, xo, yo) + write1D_4DField('xaspect', u, xo, yo) + write1D_4DField('yaspect', v, xo, yo) if error > .01: err = 1 return err - elif case == 6: + elif case == 6: sendOutput('\n******* single point mode -- single precision *****\n') # array dimensions @@ -384,7 +384,7 @@ def choose(case): nxo = 21 nyo = 21 - # input arrays and data + # input arrays and data xisort, xi = randomGrid(ni, 1.2, -0.2) # xisort has random numbers monotonically increasing yisort, yi = randomGrid(ni, 1.2, -0.2) @@ -393,19 +393,19 @@ def choose(case): for i in range(ni): dataIn[i] = (xi[i] - 0.25)*(xi[i] - 0.25) + (yi[i] - 0.50)*(yi[i] - 0.50) - # output array + # output array xo = uniformGrid(nxo, 1.0, 0.0) yo = uniformGrid(nyo, 1.0, 0.0) xn, yn = grid2Dto1D(xo, yo) - r = nat.Natgrid(xi, yi, xn, yn, listOutput = 'yes') + r = nat.Natgrid(xi, yi, xn, yn, listOutput = 'yes') r.igr = 1 # choose nonlinear interpolation - zn = r.rgrd(dataIn) + zn = r.rgrd(dataIn) xo, yo, dataOut = c1Dto2D(nxo, nyo, xn, yn, zn) sendOutput('*** writing single precision single point mode test case to the netCDF file SinglePointMode.nc') - write1D_4DField('SinglePointMode', dataOut, xo, yo) + write1D_4DField('SinglePointMode', dataOut, xo, yo) dataCheck = numpy.zeros((nxo,nyo), numpy.float32) for i in range(nxo): @@ -413,7 +413,7 @@ def choose(case): dataCheck[i,j] = (xo[i] - 0.25)*(xo[i] - 0.25) + (yo[j] - 0.50)*(yo[j] - 0.50) sendOutput('*** writing exact answer to single precision single point mode test case to the netCDF file SinglePointExact.nc') - write1D_4DField('SinglePointExact', dataOut, xo, yo) + write1D_4DField('SinglePointExact', dataOut, xo, yo) error = rmserror(dataOut, dataCheck) # find the rms error sendOutput('\n******* compare results\n') @@ -426,10 +426,10 @@ def choose(case): return err - elif case == 7: + elif case == 7: sendOutput('\n******* nonlinear interpolation of y32 with a wrap -- single precision *****\n') - # input arrays and data + # input arrays and data lati,latiSort,loni,loniSort = storedGrids() @@ -437,10 +437,10 @@ def choose(case): newOrder = (1,0) y32 = numpy.transpose(y32, newOrder) - lonLinear, latLinear, y32Linear = c2Dto1D(loni, lati, y32) # change to the linear list format + lonLinear, latLinear, y32Linear = c2Dto1D(loni, lati, y32) # change to the linear list format - # output array + # output array nlato = 71 nlono = 144 @@ -448,17 +448,17 @@ def choose(case): lono = uniformGrid(nlono, 357.5, 0.0) # start at 0. - r = nat.Natgrid(latLinear, lonLinear, lato, lono) + r = nat.Natgrid(latLinear, lonLinear, lato, lono) #r.igr = 1 # choose nonlinear interpolation - dataOut = r.rgrd(y32Linear, wrap = 'yes') + dataOut = r.rgrd(y32Linear, wrap = 'yes') dataCheck = YData(lono, lato) # longitude varies the fastest sendOutput('*** writing exact answer to single precision y32 interpolatiion test case to the netCDF file y32Exact.nc') write1D_4DField('y32Exact', dataCheck, lato, lono) # lono varies the fastest. Shape is(nlati, nloni) sendOutput('*** writing single precision y32 interpolation test case to the netCDF file y32Regrid.nc') - write1D_4DField('y32Regrid', dataOut, lato, lono) + write1D_4DField('y32Regrid', dataOut, lato, lono) error = rmserror(dataOut, dataCheck) # find the rms error sendOutput('\n******* compare results\n') @@ -480,11 +480,11 @@ def choose(case): print 'making instance for case 8' - r = nat.Natgrid(latLinear, lonLinear, lato, lono) + r = nat.Natgrid(latLinear, lonLinear, lato, lono) print 'call rgrd method for case 8' - dataOut = r.rgrd(y32Linear, wrap = 'yes') + dataOut = r.rgrd(y32Linear, wrap = 'yes') print 'returning from call rgrd method for case 8' @@ -502,7 +502,7 @@ def choose(case): def randomGrid(number, vmax, vmin): #---------------------------------------------------------------------------------------- - # + # # purpose: to construct a grid coordinate which is random but monotonically increasing # # usage: vsort, vraw = randomGrid(number, vmax, vmin) @@ -531,12 +531,12 @@ def randomGrid(number, vmax, vmin): def storedAnswers(choice): #---------------------------------------------------------------------------------------- - # + # # purpose: to store the answers to selected test cases # # usage: data = storedAnswers(choice) # - # passed : choice -- a string idetifying the desired data + # passed : choice -- a string idetifying the desired data # # returned: data # @@ -669,7 +669,7 @@ def storedAnswers(choice): def uniformGrid(number, vend, vstart): #---------------------------------------------------------------------------- - # + # # purpose: to construct a grid coordinate which is uniform # # usage: v = uniformGrid(number, vend, vstart) @@ -690,15 +690,15 @@ def uniformGrid(number, vend, vstart): def storedGrids(): """ #------------------------------------------------------------------- - # - # purpose: to construct a grid coordinate which is random # - # passed : nothing + # purpose: to construct a grid coordinate which is random + # + # passed : nothing # # returned: lati -- a 60 element latitude grid from -90. to +90. degrees - # latiSort -- lati sorted to be montonically decreasing + # latiSort -- lati sorted to be montonically decreasing # loni -- a 120 element longitude grid from 0. to 360. degrees - # loniSort -- loni sorted to be montonically increasing + # loniSort -- loni sorted to be montonically increasing # #------------------------------------------------------------------------""" latiList = [ @@ -765,15 +765,15 @@ def storedGrids(): def grid2Dto1D(x, y): """ #------------------------------------------------------------------- - # + # # purpose: to construct a linear grid from a rectangular one # - # passed : x[i] and y[j] + # passed : x[i] and y[j] # # returned: xn[n], yn[n] # #------------------------------------------------------------------------""" - + numberx = len(x) numbery = len(y) size =numberx*numbery @@ -790,15 +790,15 @@ def grid2Dto1D(x, y): def c1Dto2D(numberx, numbery, xn, yn, zn): """ #------------------------------------------------------------------- - # + # # purpose: to construct 2D z[i,j] 1D zn[n] format # # passed: xn[n], yn[n], zn[n] # - # returned : x[i], y[j] and z[i,j] + # returned : x[i], y[j] and z[i,j] # #------------------------------------------------------------------------""" - + x = numpy.zeros(numberx, numpy.float32) y = numpy.zeros(numbery, numpy.float32) @@ -813,13 +813,13 @@ def c1Dto2D(numberx, numbery, xn, yn, zn): return (x, y, z) def c2Dto1D(x, y, z): #--------------------------------------------------------------------------------------------------- - # + # # purpose: to construct 1D zn[n] from 2D z[i,j] format # # usage: xn, yn, zn = c2Dto1D(x, y, z) # - # passed: x - the array which describes the rectangular grid associated with the first z index - # y - the array which describes the rectangular grid associated with the second z index + # passed: x - the array which describes the rectangular grid associated with the first z index + # y - the array which describes the rectangular grid associated with the second z index # z - the 2D data associated with the x, y grid # # returned: xn - a list form of the x array @@ -827,7 +827,7 @@ def c2Dto1D(x, y, z): # zn - a list form of the data array (this array has the same length as xn and yn # #--------------------------------------------------------------------------------------------------- - + numberx = len(x) numbery = len(y) size =numberx*numbery @@ -844,46 +844,46 @@ def c2Dto1D(x, y, z): return (xn, yn, zn) -def write1D_4DField(varname, dataField, x, y = None, z = None, t = None): +def write1D_4DField(varname, dataField, x, y = None, z = None, t = None): #------------------------------------------------------------------------------ - # + # # purpose: write an output field which may be 1D, 2D, 3D or 4D to a NetCDF file # - # usage: write1D_4DField(varname, dataField, x, y, z = None, t = None) for a 2D write + # usage: write1D_4DField(varname, dataField, x, y, z = None, t = None) for a 2D write # # passed : varname - name of the variable and the file id # x,y,z,t - dimension vectors # dataField - the data # - # returned: None + # returned: None # #------------------------------------------------------------------------------- import cdms2 - fileObj = cdms2.createDataset(varname + '.nc') + fileObj = cdms2.createDataset(varname + '.nc') # construct the axis tuple x = x.astype(numpy.float64) - x_axis = fileObj.createAxis('x', x) + x_axis = fileObj.createAxis('x', x) axisList = [x_axis] if y is not None: y = y.astype(numpy.float64) - y_axis = fileObj.createAxis('y', y) + y_axis = fileObj.createAxis('y', y) axisList.append(y_axis) if z is not None: z = z.astype(numpy.float64) - z_axis = fileObj.createAxis('z', z) + z_axis = fileObj.createAxis('z', z) axisList.append(z_axis) if t is not None: t = t.astype(numpy.float64) - t_axis = fileObj.createAxis('t', t) + t_axis = fileObj.createAxis('t', t) axisList.append(t_axis) - + if len(axisList) == 1: axisTuple = (x_axis,) else: @@ -895,7 +895,7 @@ def write1D_4DField(varname, dataField, x, y = None, z = None, t = None): var[:] = dataField # copy in the data - fileObj.close() + fileObj.close() return None @@ -903,15 +903,15 @@ def write1D_4DField(varname, dataField, x, y = None, z = None, t = None): def YData(lonvals, latvals, data_name = 'Y32'): #---------------------------------------------------------------------------- - # - # purpose: construct Y33, Y32, Y31 or Y30 data + # + # purpose: construct Y33, Y32, Y31 or Y30 data # # usage: data = YData(lonvals, latvals, data_name = 'Y32'): # - # passed : lonvals -- longitude vactor - # latvals -- latitude vactor + # passed : lonvals -- longitude vactor + # latvals -- latitude vactor # - # returned: data + # returned: data #----------------------------------------------------------------------------- if data_name[:3] == 'Y33': @@ -931,22 +931,22 @@ def YData(lonvals, latvals, data_name = 'Y32'): def Y33(lonvals, latvals): #------------------------------------------------------------------------------ - # - # purpose: construct Y33 data + # + # purpose: construct Y33 data # # usage: y33 = Y33(lonvals, latvals) # - # passed : lonvals -- longitude vactor - # latvals -- latitude vactor + # passed : lonvals -- longitude vactor + # latvals -- latitude vactor # - # returned: data + # returned: data #------------------------------------------------------------------------------ nlon = len(lonvals) nlat = len(latvals) - phi = (math.pi/180.)*lonvals + phi = (math.pi/180.)*lonvals theta = (math.pi/180.)*latvals - + y33 = numpy.zeros( (nlat,nlon), numpy.float32) # memory fac = -(1./4.)*math.sqrt( (35./(4.*math.pi)) ) @@ -961,22 +961,22 @@ def Y33(lonvals, latvals): def Y32(lonvals, latvals): #------------------------------------------------------------------------------- - # - # purpose: construct Y32 data + # + # purpose: construct Y32 data # # usage: y32 = Y32(lonvals, latvals) # - # passed : lonvals -- longitude vactor - # latvals -- latitude vactor + # passed : lonvals -- longitude vactor + # latvals -- latitude vactor # - # returned: data + # returned: data #------------------------------------------------------------------------------- nlon = len(lonvals) nlat = len(latvals) - phi = (math.pi/180.)*lonvals + phi = (math.pi/180.)*lonvals theta = (math.pi/180.)*latvals - + y32 = numpy.zeros( (nlat,nlon), numpy.float32) # memory fac = (1./4.)*math.sqrt( (105./(4.*math.pi)) ) @@ -991,22 +991,22 @@ def Y32(lonvals, latvals): def Y31(lonvals, latvals): #-------------------------------------------------------------------------------- - # - # purpose: construct Y31 data + # + # purpose: construct Y31 data # # usage: y31 = Y31(lonvals, latvals) # - # passed : lonvals -- longitude vactor - # latvals -- latitude vactor + # passed : lonvals -- longitude vactor + # latvals -- latitude vactor # - # returned: data + # returned: data #-------------------------------------------------------------------------------- nlon = len(lonvals) nlat = len(latvals) - phi = (math.pi/180.)*lonvals + phi = (math.pi/180.)*lonvals theta = (math.pi/180.)*latvals - + y31 = numpy.zeros( (nlat,nlon), numpy.float32) # memory fac = -(1./4.)*math.sqrt( (21./(4.*math.pi)) ) @@ -1021,30 +1021,30 @@ def Y31(lonvals, latvals): def Y30(lonvals, latvals): #---------------------------------------------------------------------------------- - # - # purpose: construct Y30 data + # + # purpose: construct Y30 data # # usage: y30 = Y30(lonvals, latvals) # - # passed : lonvals -- longitude vactor - # latvals -- latitude vactor + # passed : lonvals -- longitude vactor + # latvals -- latitude vactor # - # returned: data + # returned: data #----------------------------------------------------------------------------------- nlon = len(lonvals) nlat = len(latvals) - phi = (math.pi/180.)*lonvals + phi = (math.pi/180.)*lonvals theta = (math.pi/180.)*latvals lonvals = makelon(nlon) - phi = lonvals - phi = (math.pi/180.)*lonvals + phi = lonvals + phi = (math.pi/180.)*lonvals latvals, colatvals = makelat(nlat, grid_type) latvals, colatvals = makelat(nlat) theta = (math.pi/180.)*colatvals - + y30 = numpy.zeros( (nlat,nlon), numpy.float32) # memory fac = math.sqrt( (7./(4.*math.pi)) ) @@ -1062,9 +1062,9 @@ def rmserror(data1, data2): #--------------------------------------------------------------------------------- # # purpose: compute the rms error for two data sets having the same shape - # - # passed : the two data sets - # + # + # passed : the two data sets + # # returned: rms error # #--------------------------------------------------------------------------------- @@ -1089,7 +1089,7 @@ def rmserror(data1, data2): output = open('test.asc', 'w') # global file name print 'Running the test computations' - testError = runtests() + testError = runtests() write = document() sendOutput(' ') diff --git a/natgrid/setup.py b/natgrid/setup.py index d282c77f..bfec26a5 100644 --- a/natgrid/setup.py +++ b/natgrid/setup.py @@ -10,7 +10,7 @@ url = "http://cdat.sf.net", packages = [''], package_dir = {'': 'Lib'}, - include_dirs = ['Include',], + include_dirs = ['Include',], ext_modules = [Extension('natgridmodule', sources), - ] + ] ) diff --git a/test/colourcorrection.py b/test/colourcorrection.py index e2566775..45595e20 100644 --- a/test/colourcorrection.py +++ b/test/colourcorrection.py @@ -1,11 +1,11 @@ """ This is for pybdsm for calculating spectral index. We assume a linear spectral index in log(freq) and then each channel has a flux which is bit wrong because of the colour -correction problem within that band. +correction problem within that band. Now we average n such channels. There will be another error made, partly because of the -colour correction now for channels (including discretisation) and the colour correction -of the earlier 2nd order colour correction. +colour correction now for channels (including discretisation) and the colour correction +of the earlier 2nd order colour correction. This is to see how much they differ. Refer notebook for forumlae. @@ -36,7 +36,7 @@ f_naive = N.mean(f_arr) f1 = N.power(f_arr, alpha) f2 = N.power(f_arr, alpha-2.0) - + f1 = 1.0/n*N.sum(f1) f2 = 1.0/n*N.sum(f2)*bw*bw*alpha*(alpha-1.0)/24.0 @@ -47,8 +47,8 @@ f_diff2[ifreq] = f_eff1 - f_eff2 fig = pl.figure(pl1.number) - adjustprops = dict(wspace=0.5, hspace=0.5) - fig.subplots_adjust(**adjustprops) + adjustprops = dict(wspace=0.5, hspace=0.5) + fig.subplots_adjust(**adjustprops) ax = pl.subplot(2,3,k) pl.plot(freq/1e6, f_diff1/1e3) pl.title('n='+str(n)+'; bw='+str(bw/1e6)+' MHz') @@ -57,18 +57,18 @@ pl.setp(ax.get_xticklabels(), rotation='vertical', fontsize=12) fig = pl.figure(pl2.number) - adjustprops = dict(wspace=0.5, hspace=0.5) - fig.subplots_adjust(**adjustprops) + adjustprops = dict(wspace=0.5, hspace=0.5) + fig.subplots_adjust(**adjustprops) ax = pl.subplot(2,3,k) pl.plot(freq/1e6, f_diff2) pl.title('n='+str(n)+'; bw='+str(bw/1e6)+' MHz') pl.xlabel('Freq(MHz)') pl.ylabel('Diff due to 2nd order (Hz)') pl.setp(ax.get_xticklabels(), rotation='vertical', fontsize=12) - + fig = pl.figure(pl3.number) - adjustprops = dict(wspace=0.9, hspace=0.5) - fig.subplots_adjust(**adjustprops) + adjustprops = dict(wspace=0.9, hspace=0.5) + fig.subplots_adjust(**adjustprops) ax = pl.subplot(2,3,k) f2 = f_naive+5e6 y = f_diff1*alpha/f_naive/math.log(f_naive/(f2)) @@ -84,5 +84,3 @@ pl.savefig('colourcorr_order1-2.png') pl.figure(pl3.number) pl.savefig('colourcorr_delta_spin.png') - - diff --git a/test/do_stuff.py b/test/do_stuff.py index 2d8ee063..7fef2e0d 100644 --- a/test/do_stuff.py +++ b/test/do_stuff.py @@ -11,15 +11,15 @@ def do_ws(isls, crms): thr = crms for isl in isls: - - image = isl.image*~isl.mask_active - op1, markers1 = func.watershed(image, thr=thr*3.) - - pl.figure() - pl.suptitle('Island '+str(isl.island_id)) - pl.subplot(2,2,1); pl.imshow(N.transpose(image), origin='lower', interpolation='nearest', vmin=-7*thr, vmax=15*thr); pl.title('Image') - pl.subplot(2,2,2); pl.imshow(N.transpose(op1*~isl.mask_active), origin='lower', interpolation='nearest'); pl.title('watershed1') - pl.subplot(2,2,3); pl.imshow(N.transpose(markers1*~isl.mask_active), origin='lower', interpolation='nearest'); pl.title('markers1') + + image = isl.image*~isl.mask_active + op1, markers1 = func.watershed(image, thr=thr*3.) + + pl.figure() + pl.suptitle('Island '+str(isl.island_id)) + pl.subplot(2,2,1); pl.imshow(N.transpose(image), origin='lower', interpolation='nearest', vmin=-7*thr, vmax=15*thr); pl.title('Image') + pl.subplot(2,2,2); pl.imshow(N.transpose(op1*~isl.mask_active), origin='lower', interpolation='nearest'); pl.title('watershed1') + pl.subplot(2,2,3); pl.imshow(N.transpose(markers1*~isl.mask_active), origin='lower', interpolation='nearest'); pl.title('markers1') def open_isl(isls, crms): import pylab as pl @@ -32,22 +32,18 @@ def open_isl(isls, crms): ft3 = N.ones((3,3), int) ft5 = N.ones((5,5), int) for isl in isls: - ma = ~isl.mask_active - open1 = nd.binary_opening(ma, ft1) - open2 = nd.binary_opening(ma, ft2) - open3 = nd.binary_opening(ma, ft3) - open5 = nd.binary_opening(ma, ft5) - - pl.figure() - pl.suptitle('Island '+str(isl.island_id)) - pl.subplot(2,2,1); pl.imshow(N.transpose(isl.image), origin='lower', interpolation='nearest'); pl.title('Image') - pl.subplot(2,2,2); pl.imshow(N.transpose(ma), origin='lower', interpolation='nearest'); pl.title('mask') - pl.subplot(2,2,3); pl.imshow(N.transpose(open3), origin='lower', interpolation='nearest'); pl.title('open 3x3') - pl.subplot(2,2,4); pl.imshow(N.transpose(open5), origin='lower', interpolation='nearest'); pl.title('open 5x5') - #pl.subplot(2,2,3); pl.imshow(N.transpose(open1), origin='lower', interpolation='nearest'); pl.title('open diag') - #pl.subplot(2,2,4); pl.imshow(N.transpose(open2), origin='lower', interpolation='nearest'); pl.title('open str') - pl.savefig('cyga_p_w12_bigisl_'+str(isl.island_id)+'_open.png') - - - - + ma = ~isl.mask_active + open1 = nd.binary_opening(ma, ft1) + open2 = nd.binary_opening(ma, ft2) + open3 = nd.binary_opening(ma, ft3) + open5 = nd.binary_opening(ma, ft5) + + pl.figure() + pl.suptitle('Island '+str(isl.island_id)) + pl.subplot(2,2,1); pl.imshow(N.transpose(isl.image), origin='lower', interpolation='nearest'); pl.title('Image') + pl.subplot(2,2,2); pl.imshow(N.transpose(ma), origin='lower', interpolation='nearest'); pl.title('mask') + pl.subplot(2,2,3); pl.imshow(N.transpose(open3), origin='lower', interpolation='nearest'); pl.title('open 3x3') + pl.subplot(2,2,4); pl.imshow(N.transpose(open5), origin='lower', interpolation='nearest'); pl.title('open 5x5') + #pl.subplot(2,2,3); pl.imshow(N.transpose(open1), origin='lower', interpolation='nearest'); pl.title('open diag') + #pl.subplot(2,2,4); pl.imshow(N.transpose(open2), origin='lower', interpolation='nearest'); pl.title('open str') + pl.savefig('cyga_p_w12_bigisl_'+str(isl.island_id)+'_open.png') diff --git a/test/tbdsf_process_image.py b/test/tbdsf_process_image.py index f55d3fbc..7941a70f 100644 --- a/test/tbdsf_process_image.py +++ b/test/tbdsf_process_image.py @@ -16,4 +16,3 @@ sys.exit(0) else: sys.exit(1) - diff --git a/test/test.py b/test/test.py index dcbda272..7b1a3a9a 100644 --- a/test/test.py +++ b/test/test.py @@ -27,7 +27,7 @@ def call_pybdsm(version, parameters): img = bdsm.execute(bdsm.fits_chain, parameters) return img, bdsm - + #img, bdsm = call_pybdsm('test', {'fits_name': "subim.fits", 'beam' : (0.0015, 0.0015, 0.0), 'thresh':"hard", 'atrous_do' : False}) @@ -93,4 +93,3 @@ def call_pybdsm(version, parameters): # 'polarisation_do': True, 'atrous_do' : True}) #img = bdsm.execute(bdsm.fits_chain,{'fits_name': "try.fits", 'beam': (.056, .028, 160.0), 'thresh':"hard", 'thresh_pix':20.}) - diff --git a/test/test_watershed.py b/test/test_watershed.py index aec5e52b..4536cc9f 100644 --- a/test/test_watershed.py +++ b/test/test_watershed.py @@ -10,62 +10,62 @@ from copy import deepcopy as cp for isl in img.islands: - #isl = img.islands[153] - if isl.ngaus > 1: - thr = isl.mean + img.opts.thresh_pix*isl.rms - im = isl.image; mask = isl.mask_active; av = img.clipped_mean; im1 = cp(im) - ind = N.array(N.where(~mask)).transpose() - ind = [tuple(coord) for coord in ind if im[tuple(coord)] > thr] - n, m = isl.shape; iniposn = []; inipeak = [] - for c in ind: - goodlist = [im[i,j] for i in range(c[0]-1,c[0]+2) for j in range(c[1]-1,c[1]+2) \ - if i>=0 and i=0 and j goodlist) == len(goodlist) - if peak: - iniposn.append(c); inipeak.append(im[c]) - nmulsrc = len(iniposn) - if nmulsrc > 1: - markers = N.zeros(im.shape, int) - markers[0,0] = 1 - for ipk in range(nmulsrc): - pk = inipeak[ipk]; x, y = iniposn[ipk] - markers[int(round(x)), int(round(y))] = ipk+2 - im2 = N.zeros(im.shape, int) - im1 = im1 - im1.min() - im1 = im1/im1.max()*255 - im1 = 255-im1 - nd.watershed_ift(N.array(im1, N.uint8), markers, output = im2) - fcn = MGFunction(im, isl.mask_active, 1) - fit = lmder_fit - gg1 = gg() - for ipk in range(nmulsrc): - ind = ipk+2 - mom = func.momanalmask_gaus(im, im2, ind, 1.0, True) - indd = N.where(im2==ind) - mom[3] = 3.0; mom[4]=3.0 - g = [float(N.max(im[indd])), int(round(mom[1])), int(round(mom[2])), mom[3]/fwsig, mom[4]/fwsig, mom[5]] - gg1.add_gaussian(fcn, g, dof = isl.size_active) - print g - fit(fcn, final=0, verbose=True) - print fcn.parameters - import pylab as pl - pl.figure() - pl.subplot(2,2,1);pl.imshow(N.transpose(im), interpolation='nearest', origin='lower'); pl.title(str(isl.island_id)) - pl.subplot(2,2,2);pl.imshow(N.transpose(im1), interpolation='nearest', origin='lower'); pl.title(str(isl.island_id)) - pl.subplot(2,2,3);pl.imshow(N.transpose(im2), interpolation='nearest', origin='lower'); pl.title(str(isl.island_id)) - for g in fcn.parameters: - A, x1, x2, s1, s2, th = g - s1, s2 = map(abs, [s1, s2]) - if s1 < s2: # s1 etc are sigma - ss1=s2; ss2=s1; th1 = divmod(th+90.0, 180)[1] - else: - ss1=s1; ss2=s2; th1 = divmod(th, 180)[1] - c = [A, x1, x2, ss1, ss2, th1] - x, y = N.indices(isl.shape) - x2, y2 = func.drawellipse(c) - #x2 = x2 + isl.origin[0]; y2 = y2 + isl.origin[1] - pl.subplot(2,2,4); pl.plot(x2, y2, '-r') - pl.imshow(N.transpose(im), origin='lower', interpolation='nearest') + #isl = img.islands[153] + if isl.ngaus > 1: + thr = isl.mean + img.opts.thresh_pix*isl.rms + im = isl.image; mask = isl.mask_active; av = img.clipped_mean; im1 = cp(im) + ind = N.array(N.where(~mask)).transpose() + ind = [tuple(coord) for coord in ind if im[tuple(coord)] > thr] + n, m = isl.shape; iniposn = []; inipeak = [] + for c in ind: + goodlist = [im[i,j] for i in range(c[0]-1,c[0]+2) for j in range(c[1]-1,c[1]+2) \ + if i>=0 and i=0 and j goodlist) == len(goodlist) + if peak: + iniposn.append(c); inipeak.append(im[c]) + nmulsrc = len(iniposn) + if nmulsrc > 1: + markers = N.zeros(im.shape, int) + markers[0,0] = 1 + for ipk in range(nmulsrc): + pk = inipeak[ipk]; x, y = iniposn[ipk] + markers[int(round(x)), int(round(y))] = ipk+2 + im2 = N.zeros(im.shape, int) + im1 = im1 - im1.min() + im1 = im1/im1.max()*255 + im1 = 255-im1 + nd.watershed_ift(N.array(im1, N.uint8), markers, output = im2) + fcn = MGFunction(im, isl.mask_active, 1) + fit = lmder_fit + gg1 = gg() + for ipk in range(nmulsrc): + ind = ipk+2 + mom = func.momanalmask_gaus(im, im2, ind, 1.0, True) + indd = N.where(im2==ind) + mom[3] = 3.0; mom[4]=3.0 + g = [float(N.max(im[indd])), int(round(mom[1])), int(round(mom[2])), mom[3]/fwsig, mom[4]/fwsig, mom[5]] + gg1.add_gaussian(fcn, g, dof = isl.size_active) + print g + fit(fcn, final=0, verbose=True) + print fcn.parameters + import pylab as pl + pl.figure() + pl.subplot(2,2,1);pl.imshow(N.transpose(im), interpolation='nearest', origin='lower'); pl.title(str(isl.island_id)) + pl.subplot(2,2,2);pl.imshow(N.transpose(im1), interpolation='nearest', origin='lower'); pl.title(str(isl.island_id)) + pl.subplot(2,2,3);pl.imshow(N.transpose(im2), interpolation='nearest', origin='lower'); pl.title(str(isl.island_id)) + for g in fcn.parameters: + A, x1, x2, s1, s2, th = g + s1, s2 = map(abs, [s1, s2]) + if s1 < s2: # s1 etc are sigma + ss1=s2; ss2=s1; th1 = divmod(th+90.0, 180)[1] + else: + ss1=s1; ss2=s2; th1 = divmod(th, 180)[1] + c = [A, x1, x2, ss1, ss2, th1] + x, y = N.indices(isl.shape) + x2, y2 = func.drawellipse(c) + #x2 = x2 + isl.origin[0]; y2 = y2 + isl.origin[1] + pl.subplot(2,2,4); pl.plot(x2, y2, '-r') + pl.imshow(N.transpose(im), origin='lower', interpolation='nearest') @@ -88,10 +88,10 @@ brat = [2.0, 2.0, 2.0, 2.0] markers[10,10] = 1 for ig in range(len(peaks)): - g = peaks[ig]*N.exp(-0.5*((x-posns[ig][0])*(x-posns[ig][0])+(y-posns[ig][1])*(y-posns[ig][1])) \ - /(bmaj[ig]*bmaj[ig]/brat[ig])) - image = image + g - markers[int(round(posns[ig][0])), int(round(posns[ig][1]))] = ig+2 + g = peaks[ig]*N.exp(-0.5*((x-posns[ig][0])*(x-posns[ig][0])+(y-posns[ig][1])*(y-posns[ig][1])) \ + /(bmaj[ig]*bmaj[ig]/brat[ig])) + image = image + g + markers[int(round(posns[ig][0])), int(round(posns[ig][1]))] = ig+2 image1 = image - image.min() image1 = image1/image1.max()*255 @@ -101,6 +101,3 @@ pl.figure();pl.imshow(N.transpose(image1), interpolation='nearest', origin='lower'); pl.title('input1'); pl.colorbar() pl.figure();pl.imshow(N.transpose(op1), interpolation='nearest', origin='lower'); pl.title('output1'); pl.colorbar() pl.figure();pl.imshow(N.transpose(markers), interpolation='nearest', origin='lower'); pl.title('markers'); pl.colorbar() - - -