|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 from taskinit import * 00002 from clean import clean 00003 from cleanhelper import * 00004 from immath import immath # Not needed? 00005 from sys import stdout # Not needed? 00006 from math import ceil 00007 from math import floor 00008 from math import sqrt 00009 import numpy 00010 import shutil 00011 00012 ##################################################################### 00013 # iterative cleaning while automatically selecting new clean regions 00014 ##################################################################### 00015 00016 # Some places in code commented by #TMP# indicate temporary fixes that 00017 # should be removed/modified once certain features or bugs in CASA have 00018 # been implemented or fixed. 00019 00020 def autoclean(vis, imagename, field, spw, selectdata, timerange, uvrange, 00021 antenna, scan, mode, nchan, start, width, interpolation, doconcat, 00022 psfmode, imagermode, cyclefactor, cyclespeedup, imsize, cell, 00023 phasecenter, restfreq, stokes, weighting, robust, noise, npixels, 00024 interactive, mask, modelimage, uvtaper, outertaper, innertaper, 00025 niter, npercycle, npercycle_speedup, gain, pbcor, minpb, 00026 clean_threshold, Nrms, eps_maxres, useabsresid, allow_maxres_inc, 00027 island_rms, diag, peak_rms, gain_threshold, Npeak, shape, 00028 boxstretch, irregsize): 00029 00030 casalog.origin('autoclean') 00031 00032 orig_npercycle = npercycle 00033 00034 # set up initial mask values if input by user 00035 inputmask = None 00036 if(mask): 00037 casalog.post('Using input mask from user.') 00038 imCln = imtool() 00039 imset = cleanhelper(imCln, vis) 00040 imset.defineimages(imsize=imsize, cell=cell, stokes=stokes, 00041 mode=mode, spw=spw, nchan=1, 00042 start=start, width=width, 00043 restfreq=restfreq, field=field, 00044 phasecenter=phasecenter) 00045 imset.datselweightfilter(field=field, spw=spw, timerange=timerange, 00046 uvrange=uvrange, antenna=antenna, scan=scan, 00047 wgttype=weighting, robust=robust, noise=noise, 00048 npixels=npixels, 00049 #mosweight=mosweight, 00050 innertaper=innertaper, outertaper=outertaper, 00051 calready=False) 00052 inputmask = '__temporary.mask' 00053 imset.makemaskimage(outputmask=inputmask, maskobject=mask) 00054 imCln.close() 00055 ia.open(inputmask) 00056 maskVals = ia.getchunk() 00057 ia.close() 00058 ia.removefile(inputmask) 00059 00060 # begin loop over all channels 00061 for ichan in xrange(nchan): 00062 00063 if imagermode!='mfs': 00064 casalog.post('***** Beginning channel %d of %d' % (ichan, nchan)) 00065 00066 if (mode=='mfs') or (nchan==1): 00067 thisImage = imagename 00068 thisStart = start 00069 else: 00070 thisImage = imagename + '.channel.' + str(ichan) 00071 if mode=='channel': 00072 thisStart = start + ichan*width 00073 else: 00074 thisStart = qa.add(start, qa.mul(width, qa.quantity(ichan))) 00075 00076 # create dirty image from which to choose initial clean regions 00077 clean(vis=vis, imagename=thisImage, field=field, spw=spw, 00078 selectdata=selectdata, timerange=timerange, uvrange=uvrange, 00079 antenna=antenna, scan=scan, mode=mode, nchan=1, start=thisStart, 00080 width=width, interpolation=interpolation, psfmode=psfmode, 00081 imagermode=imagermode, cyclefactor=cyclefactor, 00082 cyclespeedup=cyclespeedup, imsize=imsize, cell=cell, 00083 phasecenter=phasecenter, restfreq=restfreq, stokes=stokes, 00084 weighting=weighting, robust=robust, noise=noise, npixels=npixels, 00085 modelimage=modelimage, uvtaper=uvtaper, outertaper=outertaper, 00086 innertaper=innertaper, niter=0, gain=gain, minpb=minpb) 00087 casalog.origin('autoclean') 00088 00089 # make the mask image 00090 ia_tool = ia.newimagefromimage(infile=thisImage+'.image', 00091 outfile=thisImage+'.mask', 00092 overwrite=True) 00093 if(mask): 00094 ia_tool.putchunk(pixels=maskVals) 00095 else: 00096 ia_tool.set(pixels=0.0) 00097 ia_tool.close() 00098 00099 # keep track of changes in maximum residual over iterations 00100 maxres_increased = 0 00101 post_max = 5000 00102 # iterations variable stores number of clean internal iterations 00103 iterations = 0 00104 # new_npercycle: number of minor cycles this iteration (can change) 00105 new_npercycle = orig_npercycle 00106 while iterations < niter: 00107 cleanmore, nregions = \ 00108 autowindow(imagename=thisImage, Npeak=Npeak, 00109 shape=shape, irregsize=irregsize, 00110 boxstretch=boxstretch, Nrms=Nrms, 00111 island_rms=island_rms, peak_rms=peak_rms, 00112 gain_threshold=gain_threshold, 00113 clean_threshold=clean_threshold, diag=diag, 00114 useabsresid=useabsresid, ichan=ichan) 00115 00116 # if first iteration, make sure CLEAN is cleaning 00117 if not(cleanmore): 00118 if not(iterations): 00119 casalog.post('Input parameters do not induce cleaning of ' 00120 'channel %d.' % ichan) 00121 break 00122 00123 # did autowindow make new clean regions? 00124 if(nregions): 00125 new_npercycle = orig_npercycle 00126 else: 00127 casalog.post('Will continue cleaning with current clean regions.') 00128 new_npercycle = int(round(new_npercycle * npercycle_speedup)) 00129 00130 # maximum absolute residual before next clean iteration 00131 if (eps_maxres) or (allow_maxres_inc >= 0): 00132 pre_max = post_max 00133 00134 # niter-iterations is the *total* number of iterations left, and might 00135 # be less than npercycle. Use the minimum of (npercycle, niter-iterations) 00136 casalog.post('Starting clean at iteration %d' % iterations) 00137 00138 #TMP# CAS-1623: an improvement will be made to CLEAN to allow it to 00139 # start from an existing residual image. The call to clean may need 00140 # to be changed to reflect this. 00141 clean(vis=vis, imagename=thisImage, field=field, spw=spw, 00142 selectdata=selectdata, timerange=timerange, uvrange=uvrange, 00143 antenna=antenna, scan=scan, mode=mode, nchan=1, width=width, 00144 start=thisStart, interpolation=interpolation, psfmode=psfmode, 00145 imagermode=imagermode, imsize=imsize, 00146 cyclefactor=cyclefactor, cyclespeedup=cyclespeedup, 00147 cell=cell, phasecenter=phasecenter, restfreq=restfreq, 00148 stokes=stokes, weighting=weighting, robust=robust, 00149 noise=noise, npixels=npixels, modelimage=modelimage, 00150 uvtaper=uvtaper, outertaper=outertaper, innertaper=innertaper, 00151 niter=min(new_npercycle,niter-iterations), 00152 mask=thisImage+'.mask', gain=gain, minpb=minpb, npercycle=niter, 00153 interactive=bool(nregions and interactive)) 00154 casalog.origin('autoclean') 00155 iterations += new_npercycle 00156 00157 # stop if we've reached maximum number of clean minor cycles 00158 if iterations >= niter: 00159 casalog.post('Reached maximum number of CLEAN cycles (%d)' % niter) 00160 break 00161 00162 # check maximum residual value 00163 if (eps_maxres) or (allow_maxres_inc >= 0): 00164 residualImage = thisImage+'.residual' 00165 ia.open(residualImage) 00166 stats = ia.statistics(list=False) 00167 ia.close() 00168 if(useabsresid): 00169 post_max = max([ stats['max'][0], abs(stats['min'][0]) ]) 00170 else: 00171 post_max = stats['max'][0] 00172 # check fractional change of maximum residual 00173 if(eps_maxres): 00174 if abs(post_max-pre_max)/pre_max < eps_maxres: 00175 casalog.post('Maximum residual has changed by less ' 00176 'than %.3f; stopping' % eps_maxres) 00177 break 00178 if post_max >= pre_max: 00179 maxres_increased += 1 00180 casalog.post('Maximum residual has increased.')#, 'INFO1') 00181 if maxres_increased > allow_maxres_inc: 00182 casalog.post('Maximum residual has increased %d times; ' 00183 'stopping.' % maxres_increased) 00184 break 00185 00186 # finished with all channels: concatenate if there are multiple channels 00187 if (mode!='mfs') and (nchan!=1) and (doconcat): 00188 concat_images(imagename, '.image', nchan) 00189 concat_images(imagename, '.mask', nchan) 00190 concat_images(imagename, '.flux', nchan) 00191 concat_images(imagename, '.psf', nchan) 00192 concat_images(imagename, '.model', nchan) 00193 concat_images(imagename, '.residual', nchan) 00194 concat_regions(imagename, '.rgn', nchan) 00195 # The individual .channel. tables are no longer needed 00196 shutil.rmtree(imagename + '.channel.*') 00197 00198 if(pbcor): 00199 # user wants primary beam corrected .image file and .residual file 00200 ia.open(imagename+'.image') 00201 ia.calc(imagename+'.image/'+imagename+'.flux') 00202 ia.open(imagename+'.residual') 00203 ia.calc(imagename+'.residual/'+imagename+'.flux') 00204 ia.close() 00205 00206 00207 # Autowindow selects the clean region based on peaks in each island above threshold 00208 def autowindow(imagename='', island_rms=0, gain_threshold=0, peak_rms=0, Nrms=0, 00209 boxstretch=0, clean_threshold=0, Npeak=None, shape=0, 00210 irregsize=100, diag=False, useabsresid=False, ichan=0): 00211 00212 maskImage = imagename+'.mask' 00213 residualImage = imagename+'.residual' 00214 00215 # what is the rms of the residual image outside the current clean regions? 00216 # option: could modify this to get rms by iterative sigma clipping 00217 ia.open(residualImage) 00218 rms = ia.statistics(mask=maskImage+'==0',list=False)['rms'][0] 00219 max_residual = ia.statistics(list=False)['max'][0] 00220 if(useabsresid): 00221 min_residual = ia.statistics(list=False)['min'][0] 00222 if abs(min_residual) > max_residual: 00223 max_residual = abs(min_residual) 00224 ia.close() 00225 00226 # threshold values for selecting clean regions 00227 threshold = max(peak_rms*rms, max_residual*gain_threshold) 00228 casalog.post('peak_threshold = %f' % threshold, 'INFO1') 00229 casalog.post('island_threshold = %f' % (island_rms*rms), 'INFO1') 00230 casalog.post('max_residual (%.5f) is %.1f times rms (%.5f)' % 00231 (max_residual, max_residual/rms, rms)) 00232 # If no units for clean_threshold, assume mJy for consistency with CLEAN. 00233 # But convert to Jy, because that's what units the images are in. 00234 clean_threshold = qa.convert(qa.quantity(clean_threshold,'mJy'),'Jy') 00235 clean_threshold_Jy = qa.getvalue(clean_threshold) 00236 00237 # if max_res. is below thresholds, no need to continue 00238 if (max_residual < clean_threshold_Jy) or (max_residual < rms*Nrms): 00239 if max_residual < clean_threshold_Jy: 00240 casalog.post('Max_residual is less than clean_threshold.') 00241 if max_residual < rms*Nrms: 00242 casalog.post('Max_residual is less than '+str(Nrms)+' times rms.') 00243 return 0, 0 00244 00245 # if max_res. below peak threshold, need more cleaning with current regions 00246 if max_residual < peak_rms*rms: 00247 casalog.post('Max_residual is less than box threshold for peaks.') 00248 return 1, 0 00249 00250 # select the new clean regions; put them in the mask 00251 casalog.post('Selecting clean regions.') 00252 Nregions = get_islands(imagename, Npeak, island_rms*rms, threshold, shape, 00253 boxstretch, irregsize, diag, ichan) 00254 00255 if not(Nregions): 00256 casalog.post('No new peaks passed selection.') 00257 return 1, 0 00258 00259 return 1, Nregions 00260 00261 00262 # Starting with peak in image, find islands: contiguous pixels above threshold. 00263 # Tests peak: bright enough to cleanbox? Then adds chosen region shape to mask. 00264 # Continues with next peak pixel until Npeak peaks have been found. 00265 # Won't always add Npeak new clean regions if current clean regions still 00266 # contain peaks. But isolated bright pixels (<2.5*peak_threshold) are ignored. 00267 def get_islands(imagename='', Npeak=3, island_threshold=0, peak_threshold=0, 00268 shape=0, boxstretch=0, irregsize=100, diag=False, ichan=0): 00269 00270 maskImage = imagename+'.mask' 00271 residualImage = imagename+'.residual' 00272 00273 # types for the recarrays that will store pixel and island info. 00274 island_dtype = [('box', '4i4'), ('peak_flux', 'f4')] 00275 pix_dtype = [('x','i4'),('y','i4'),('tmp_mask','f4'),('cln_mask','f4')] 00276 00277 # Get mask for all current clean regions which we'll call cln_mask 00278 ia.open(maskImage) 00279 cln_mask = ia.getregion().squeeze() 00280 ia.close() 00281 # Find all pixels above the threshold; make temporary mask: tmp_mask 00282 ia.open(residualImage) 00283 tmp_mask = ia.getregion(mask=residualImage+'>'+str(island_threshold), 00284 getmask=True).squeeze() 00285 # pixel values 00286 pixelVals = ia.getregion().squeeze() 00287 ia.close() 00288 # store pixel positions, mask values, clean region status in recarray 00289 grid = numpy.indices(tmp_mask.shape) 00290 xyMask = numpy.rec.fromarrays([grid[0], grid[1], tmp_mask, cln_mask], 00291 dtype=pix_dtype) 00292 nx, ny = tmp_mask.shape 00293 00294 # keep going until we've found Npeak islands 00295 # or there are no more pixels above the island_threshold 00296 # or the peak is less than the peak_threshold 00297 Nregions = 0 00298 Nkept = 0 00299 while Nregions < Npeak: 00300 00301 if not(xyMask['tmp_mask'].max()): 00302 # no more pixels above island threshold: we're done 00303 break 00304 00305 # find the next peak and its location 00306 xok, yok = numpy.where(xyMask['tmp_mask']) 00307 pixok = pixelVals[xok,yok] 00308 peak = pixok.max() 00309 00310 if peak < peak_threshold: 00311 # peak is below peak threshold for clean boxing: we're done 00312 break 00313 00314 # store location of peak 00315 peakind = pixok == peak 00316 xpeak = xok[peakind][0] 00317 ypeak = yok[peakind][0] 00318 x = xpeak 00319 y = ypeak 00320 00321 # since we've already checked this pixel, set its tmp mask to 0 00322 xyMask[x,y]['tmp_mask'] = 0 00323 listPix = [xyMask[x,y]] 00324 00325 # find all above-threshold contiguous pixels of this island 00326 # python lets us loop over items in a list while adding to the list! 00327 for pixel in listPix: 00328 x = pixel['x'] 00329 y = pixel['y'] 00330 # search the pixels surrounding the pixel of interest 00331 xLook1 = max(0,x-1) # in case we're at the image edge... 00332 xLook2 = min(x+2,nx-1) # | 00333 yLook1 = max(0,y-1) # | 00334 yLook2 = min(y+2,ny-1) # V 00335 # add new above-threshold pixels to the list for this island 00336 if(diag): 00337 contig_pix = xyMask[xLook1:xLook2, yLook1:yLook2] 00338 listPix += [pix for pix in contig_pix.ravel() if(pix['tmp_mask'])] 00339 # since we've already added these pixels, set their tmp mask to 0 00340 contig_pix['tmp_mask'] = 0 00341 else: 00342 contig_pix = [] 00343 contig_pix += xyMask[xLook1:xLook2, y] 00344 contig_pix += xyMask[x, yLook1:yLook2] 00345 listPix += [pix for pix in contig_pix if(pix['tmp_mask'])] 00346 # since we've already added these pixels, set their tmp mask to 0 00347 xyMask[xLook1:xLook2, y]['tmp_mask'] = 0 00348 xyMask[x, yLook1:yLook2]['tmp_mask'] = 0 00349 00350 if peak < peak_threshold: 00351 # reject islands with peak < peak_threshold 00352 continue 00353 00354 # found all pixels in this island; get bounding box 00355 islandPix = numpy.rec.fromrecords(listPix, dtype=pix_dtype) 00356 xmin = islandPix['x'].min() 00357 xmax = islandPix['x'].max() 00358 ymin = islandPix['y'].min() 00359 ymax = islandPix['y'].max() 00360 00361 if (xmin==xmax) or (ymin==ymax): 00362 # reject islands that are only 1 pixel wide, unless very bright 00363 if peak < 2.5 * peak_threshold: 00364 continue 00365 00366 Nregions += 1 # This island should be in a clean region. 00367 00368 # Is the peak already in a clean region? 00369 if(islandPix[0]['cln_mask']): 00370 casalog.post('Peak of %f ' % peak + 'at pixel ' + 00371 str([xpeak,ypeak]) + ' is already in mask.', 'INFO1') 00372 continue 00373 00374 # This is a new island for clean boxing. Prepare to mask! 00375 Nkept += 1 00376 newIsland = numpy.array(([xmin, ymin, xmax, ymax], peak), 00377 dtype=island_dtype) 00378 if (irregsize) and (min([xmax-xmin, ymax-ymin]) >= irregsize): 00379 # user wants clean region to be outline of island 00380 # for large islands 00381 mask_island(imagename, newIsland, islandPix) 00382 else: 00383 # region will be circle or box 00384 mask_region(imagename, newIsland, shape, boxstretch, ichan) 00385 return Nkept 00386 00387 00388 # If user prefers, large islands (size >= irregsize) will be masked 'as-is' 00389 # instead of surrounding them by a box-shaped or circular clean region 00390 def mask_island(imagename='', island=None, pixels=None): 00391 ia.open(imagename+'.mask') 00392 csys = ia.coordsys() # Not needed? 00393 mask = ia.getregion() 00394 mask[pixels['x'], pixels['y']] = 1 00395 ia.putchunk(mask) 00396 ia.close() 00397 casalog.post('Adding irregular region for peak of %f ' % 00398 island['peak_flux'] + 'inside box ' + str(island['box'])) 00399 #TMP# for now, can only write boxes out to a region file. So for a 00400 # temporary fix, will continue putting irregular region in mask, 00401 # but will also output a box covering the island to the .rgn file. 00402 box_to_regionfile(imagename, island['box']) 00403 00404 00405 # Chooses appropriate shape (box or circle) for clean region 00406 def mask_region(imagename='', island=None, shape=0, boxstretch=0, ichan=0): 00407 00408 peak_flux = island['peak_flux'] 00409 box = island['box'] 00410 00411 if shape==2: 00412 xsize = box[2] - box[0] 00413 ysize = box[3] - box[1] 00414 if abs(xsize-ysize) <= 1: 00415 shape = 0 00416 else: 00417 shape = 1 00418 00419 if(shape): 00420 add_box(imagename, box, peak_flux, boxstretch) 00421 else: 00422 add_circle(imagename, box, peak_flux, boxstretch) 00423 00424 00425 # Add a circular region to the mask image 00426 def add_circle(imagename, box, peak_flux, boxstretch=0): 00427 #TMP# At some point im.regiontoimagemask should allow circles that 00428 # specify a particular channel 00429 xsize = box[2] - box[0] 00430 ysize = box[3] - box[1] 00431 radius = max(sqrt(xsize**2+ysize**2)/2. + boxstretch, 1) 00432 xcen = numpy.average([box[0], box[2]]) 00433 ycen = numpy.average([box[1], box[3]]) 00434 circle = [radius, xcen, ycen] 00435 im.regiontoimagemask(mask=imagename+'.mask', circles=circle) 00436 casalog.post('Adding circle for peak of %f with center (%.1f,%.1f) and ' 00437 'radius %.1f' % (peak_flux, xcen, ycen, radius)) 00438 #TMP# can't write out circular regions to a region file, 00439 # so for now add a box that encompasses the circle 00440 # Since we increased size of circle to include potential "corners" 00441 # of emission outside a smaller (max(xsize,ysize)) radius, must 00442 # now make a box that completely encompasses the circle 00443 newbox = [floor(xcen - radius), floor(ycen - radius), 00444 ceil(xcen + radius), ceil(ycen + radius)] 00445 box_to_regionfile(imagename, newbox) 00446 00447 00448 # Add a box region to the mask image; save new region in region file 00449 def add_box(imagename, box, peak_flux, boxstretch=0): 00450 box[0:2] -= boxstretch 00451 box[2:4] += boxstretch 00452 # in case we used boxstretch < 0 and a one-pixel sized box: 00453 if box[0] > box[2]: 00454 box[0] += boxstretch 00455 box[2] -= boxstretch 00456 if box[1] > box[3]: 00457 box[1] += boxstretch 00458 box[3] -= boxstretch 00459 im.regiontoimagemask(mask=imagename+'.mask', boxes=box) 00460 casalog.post('Adding box for peak of %f ' % peak_flux + 'with coordinates ' 00461 + str(box)) 00462 box_to_regionfile(imagename, box) 00463 00464 00465 # concatenate multiple-channel images 00466 def concat_images(imagename='', suffix='', number=0, relax=False): 00467 imagelist = [imagename + '.channel.' + i + suffix 00468 for i in map(str, range(number))] 00469 ia.imageconcat(imagename + suffix, imagelist, overwrite=True, relax=relax) 00470 ia.close() 00471 00472 00473 # read in multiple region (.rgn) files; save to one file 00474 def concat_regions(imagename='', suffix='', number=0): 00475 regions = {} 00476 for file in [imagename + '.channel.' + i + suffix 00477 for i in map(str, range(number))]: 00478 if(os.path.exists(file)): 00479 regions[file] = rg.fromfiletorecord(file) 00480 nregion = len(regions) 00481 if(nregion): 00482 regionfile = imagename + '.rgn' 00483 if(os.path.exists(regionfile)): 00484 os.system('rm -f '+regionfile) 00485 if nregion > 1: 00486 rg.tofile(regionfile, rg.makeunion(regions)) 00487 else: 00488 rg.tofile(regionfile, regions[file]) 00489 00490 def box_to_regionfile(imagename, box): 00491 ia.open(imagename+'.image') 00492 csys = ia.coordsys() 00493 # need pixel corners, not pixel centers, for region file 00494 blccoord = [box[0]-0.5, box[1]-0.5] 00495 trccoord = [box[2]+0.5, box[3]+0.5] 00496 blc = ia.toworld(blccoord, 's')['string'] 00497 trc = ia.toworld(trccoord, 's')['string'] 00498 ia.close() 00499 newregion = rg.wbox(blc=blc, trc=trc, csys=csys.torecord()) 00500 regionfile = imagename+'.rgn' 00501 if(os.path.exists(regionfile)): 00502 oldregion = rg.fromfiletorecord(regionfile) 00503 regions = {'0':oldregion, '1':newregion} 00504 unionregion = rg.makeunion(regions) 00505 os.system('rm -f '+regionfile) 00506 else: 00507 unionregion = newregion 00508 rg.tofile(regionfile, unionregion)
1.8.0
