|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 ################################################ 00002 # Task to make masks. 00003 # reorganized after diesuccsion on Oct1, 2012 00004 # with Jeurgen and Urvashi 00005 # 00006 # modified by TT 00007 # based on the original code, 00008 # v1.0: 2012.03.20, U.Rau 00009 # 00010 ################################################ 00011 # Notes (to self) - TT 00012 # 1. expanding one mask to another 00013 # e.g.) expanding a continuum mask (both image mask/boolean mask) 00014 # channel mask 00015 # 2. part of copy mode func.: merging of different types of masks 00016 # e.g.) inpimage and inpmask are lists of the mask to be merged 00017 # output mask is written to either outimage or outmask as embedded 00018 # T/F mask 00019 # 3. copying mask to another or create a new one 00020 # regrid if necessary (i.e. if the coords are different) 00021 # ---------------------------------------------- 00022 # basic rules: 00023 # for mask image (1/0 mask): as it is 00024 # for internal mask : parent_imagename:mask_name 00025 # 00026 # For input, 00027 # inpimage is the casa image 00028 # - mode='list': list internal masks of inpimage 00029 # - other mode: used as a template for output 00030 # if region files are specified -> make mask specifeid with the regions on to inpimage 00031 # output is '' => modified inpimage unless overwrite=F else exception 00032 # 00033 # if inpmask='': use inpimage as input mask 00034 # if inpmask='mask0' or other embedded mask name of the inpimage, 00035 # use that T/F mask 00036 # 00037 # =expand= 00038 # case1: on the same image (outimage=''), expand mask image from 00039 # prev. run etc. No regriding. Use nearest chan mask 00040 # image to expand. 00041 # 1.a: inpimage is clean image mask (1s and 0s) 00042 # i) outimage != inpimage and outmask='' => new expanded mask image to outimage 00043 # ii) outimage != inpimage and outmask!='' => convert expanded mask image to T/F mask to store inside outimage 00044 # iii) outimage ==inpimage and outmask='' => update input mask image by expanding it 00045 # iv) outimage ==inpimage and outmask!=''=> update input image with the expanded T/F mask 00046 # 1.b: if inpmask!='', do T/F mask to 1/0 image mask conversion, then do as in 1.a 00047 00048 # case2: outimage is in diffirent coords. (need to regrid) 00049 # 00050 ################# 00051 # tests 00052 # 1. for input: mask image 00053 # 2. for input: mask/or regular image with internal mask 00054 # 3. for input: mask image; for output: mask image with different spectral grid 00055 # 4. for input: mask/regular image with internal mask; for output: image with 00056 # internal mask with different spectral grid 00057 ################### 00058 import numpy as np 00059 import shutil 00060 from taskinit import * 00061 from recipes.pixelmask2cleanmask import pixelmask2cleanmask 00062 (csys,) = gentools(['cs']) 00063 00064 def makemask(mode,inpimage, inpmask, output, overwrite, inpfreqs, outfreqs): 00065 """ 00066 make /manipulate masks 00067 00068 """ 00069 (ia,rg,im,) = gentools(['ia','rg','im']) 00070 casalog.origin('makemask') 00071 #print "params(mode,inpimage,inpmask,output,overwrite)=",mode,inpimage,inpmask,output,overwrite 00072 00073 try: 00074 # temp files 00075 tmp_maskimage='__tmp_makemaskimage' 00076 tmp_outmaskimage='__tmp_outmakemaskimage' 00077 tmp_regridim='__tmp_regridim' 00078 00079 # do parameter check first 00080 # check names of inpimage, inpmask check for existance 00081 # inpimage == output (exact match) then check overwrite 00082 # => T overwrite inpimage 00083 # => F exception 00084 00085 # check inpimage 00086 if (['list','copy','expand'].count(mode)==1): 00087 if inpimage=='': raise Exception, "inpimage is empty" 00088 if not os.path.isdir(inpimage): 00089 raise Exception, "inpimage=%s does not exist" % inpimage 00090 00091 # === list mode === 00092 if mode == 'list': 00093 inpOK=checkinput(inpimage) 00094 if inpOK: 00095 00096 ia.open(inpimage) 00097 inmasklist=ia.maskhandler('get') 00098 # now ia.maskhandler returns ['T'] if no internal mask is there... 00099 if inmasklist.count('T')!=0: 00100 inmasklist.remove('T') 00101 if len(inmasklist) ==0: 00102 casalog.post('No internal (T/F) masks were found in %s' % (inpimage),'INFO') 00103 else: 00104 defaultmaskname=ia.maskhandler('default')[0] 00105 printinmasks='' 00106 for mname in inmasklist: 00107 if mname==defaultmaskname: 00108 printinmasks+='\''+mname+'\''+'(default)' 00109 else: 00110 printinmasks+='\''+mname+'\'' 00111 if mname != inmasklist[-1]: 00112 printinmasks+=', ' 00113 00114 casalog.post('Internal (T/F) masks in %s: %s' % (inpimage, printinmasks),'INFO') 00115 ia.close() 00116 00117 elif mode == 'setdefaultmask': 00118 inpOK=checkinput(inpmask) 00119 if inpOK: 00120 (parentimage,bmask)=extractmaskname(inpmask) 00121 if bmask=='': 00122 raise Exception, "Missing an internal mask name" 00123 ia.open(parentimage) 00124 defaultmaskname=ia.maskhandler('default')[0] 00125 inmasklist=ia.maskhandler('get') 00126 if defaultmaskname==bmask: 00127 casalog.post('No change. %s is already a default internal mask' % bmask, 'INFO') 00128 else: 00129 ia.maskhandler('set',bmask) 00130 casalog.post('Set %s as a default internal mask' % bmask, 'INFO') 00131 if len(inmasklist)>1: 00132 casalog.post('Current internal masks are %s' % str(inmasklist), 'INFO') 00133 ia.close() 00134 00135 elif mode == 'delete': 00136 inpOK=checkinput(inpmask) 00137 if inpOK: 00138 (parentimage,bmask)=extractmaskname(inpmask) 00139 if bmask=='': 00140 raise Exception, "Missing an internal mask name" 00141 ia.open(parentimage) 00142 casalog.post('Deleting the internal mask, %s ' % bmask, 'INFO') 00143 defaultmaskname=ia.maskhandler('default')[0] 00144 ia.maskhandler('delete',bmask) 00145 inmasklist=ia.maskhandler('get') 00146 if inmasklist.count('T')!=0: 00147 inmasklist.remove('T') 00148 if len(inmasklist) !=0 and defaultmaskname==bmask: 00149 ia.maskhandler('set',inmasklist[0]) 00150 casalog.post('Set %s as a default internal mask' % inmasklist[0], 'INFO') 00151 if len(inmasklist)>1: 00152 casalog.post('Current internal masks are %s' % str(inmasklist), 'INFO') 00153 00154 ia.close() 00155 00156 #elif mode != 'merge': 00157 else: 00158 import commands 00159 # copy can have multiple input masks, expand has only one. 00160 # check inpimage, inpmask, output, overwrite 00161 # 00162 storeinmask = False # used to check if output to a new internal mask 00163 inpOK=checkinput(inpimage) 00164 if inpOK: 00165 (immask,inmask)=extractmaskname(inpimage) 00166 00167 # seperate text files(region files), images(with full ':'), and direct region 00168 # input mask(s) 00169 if inpmask=='': 00170 raise Exception, "Input errror. The inpmask parameter is not specified." 00171 if type(inpmask)!=list: 00172 inpmask=[inpmask] 00173 00174 # check if contains region file or region specification 00175 rgfiles=[] 00176 imgfiles=[] 00177 rglist=[] 00178 bmasks=[] 00179 for masklet in inpmask: 00180 # is text file? 00181 if type(masklet)==str: # text file or image 00182 if os.path.exists(masklet): 00183 if (commands.getoutput('file '+masklet).count('directory')): 00184 if os.path.exists(masklet+'/table.f1'): 00185 #casalog.post("%s is not in a recognized format for inpmask, ignored." % masklet, 'WARN') 00186 raise Exception, "%s is not in a recognized format for inpmask" % masklet 00187 else: 00188 # probably image file (with no mask extension) 00189 imgfiles.append(masklet) 00190 elif (commands.getoutput('file '+masklet).count('text')): 00191 rgfiles.append(masklet) 00192 else: 00193 #casalog.post("%s does not recognized format for inpmask, ignored." % masklet, 'WARN') 00194 raise Exception, "%s is not in a recognized format for inpmask" % masklet 00195 else: 00196 if masklet.count('[') and masklet.count(']'): # rough check on region specification 00197 rglist.append(masklet) 00198 else: 00199 (parentim, mask)=extractmaskname(masklet) 00200 if mask!='': 00201 bmasks.append(masklet) 00202 else: 00203 raise Exception, "%s is not an existing file/image or a region format" % masklet 00204 00205 # expand allows only a string for inpmask 00206 if mode=='expand': 00207 if type(inpmask)==list: 00208 inpmask=inpmask[0] 00209 # check for overwrite condition 00210 if output=='': 00211 if overwrite: 00212 output=inpimage 00213 else: 00214 raise Exception, "output is not specified. If you want to overwrite inpimage, please set overwrite=True" 00215 00216 if inpimage==output: 00217 if overwrite: 00218 tmp_outmaskimage=tmp_maskimage 00219 else: 00220 raise Exception, "output=inpimage. If you want to overwrite inpimage, please set overwrite=True" 00221 00222 outparentim=output 00223 outbmask='' 00224 if os.path.isdir(output): 00225 if not overwrite: 00226 raise Exception, "output=%s exists. If you want to overwrite it, please set overwrite=True" % output 00227 else: 00228 (outparentim, outbmask)=extractmaskname(output) 00229 #print "extractmaskname::: outparentim=",outparentim, "outbmask=",outbmask 00230 if outbmask!='': 00231 (parentimexist,maskexist)=checkinmask(outparentim,outbmask) 00232 if parentimexist and maskexist and not overwrite: 00233 raise Exception, "output=%s exists. If you want to overwrite it, please set overwrite=True" % output 00234 if parentimexist and not maskexist: 00235 storeinmask=True 00236 else: 00237 outparentim=output 00238 00239 00240 #print "param checks before branching out for mode==========" 00241 #print "output=",output, " is exist?=",os.path.isdir(output) 00242 #print "outparentim=",outparentim, " is exist?=",os.path.isdir(outparentim) 00243 00244 # the following code is somewhat duplicated among the modes but keep separated from each mode 00245 # for now.... copy now handle merge as well 00246 # === old copy mode === NOW combined to 'merge mode' 00247 # #if mode=='copy': 00248 # #print "Copy mode" 00249 # needregrid=True 00250 # #if outimage=='': 00251 # #overwrite 00252 # # outimage=inpimage 00253 # 00254 # if not os.path.isdir(outimage): 00255 # needregrid=False 00256 # 00257 # if inpmask!='': 00258 # # need to extract the mask and put in tmp_maskimage 00259 # pixelmask2cleanmask(imagename=inpimage, maskname=inpmask, maskimage=tmp_maskimage, usemasked=True) 00260 # else: 00261 # shutil.copytree(inpimage, tmp_maskimage) 00262 # if needregrid: 00263 # casalog.post("Regridding...",'DEBUG1') 00264 # regridmask(tmp_maskimage,outimage,tmp_outmaskimage) 00265 # # regrid may produce <1.0 pixels for mask so be sure to its all in 1.0 00266 # #ia.open(tmp_outmaskimage) 00267 # #ia.calc('iif (%s>0.0 && %s<1.0,1,%s)'%(tmp_outmaskimage,tmp_outmaskimage,tmp_outmaskimage)) 00268 # #ia.close() 00269 # #print "Copying regrid output=",tmp_outmaskimage 00270 # else: 00271 # shutil.copytree(tmp_maskimage,tmp_outmaskimage) 00272 # if outmask!='': 00273 # #convert the image mask to T/F mask 00274 # if not os.path.isdir(outimage): 00275 # shutil.copytree(inpimage,outimage) 00276 # # 00277 # ia.open(outimage) 00278 # casalog.post("convert the output image mask to T/F mask") 00279 # ia.calcmask(mask='%s<0.5' % tmp_outmaskimage,name=outmask,asdefault=True) 00280 # ia.done() 00281 # else: 00282 # # if regridded - tmp_outmaskimage is created by regridmask 00283 # # if not, tmp_outmaskimage=outimage 00284 # ia.open(tmp_outmaskimage) 00285 # ia.rename(outimage,overwrite=True) 00286 # ia.done() 00287 00288 00289 # === expand mode === 00290 if mode=='expand': 00291 (rg,) = gentools(['rg']) 00292 needtoregrid=False 00293 bychanindx=False 00294 00295 try: 00296 #print "expand mode..." 00297 # do not allow list in this mode (for inpimage and inpmask) - maybe this is redundant now 00298 if type(inpmask)==list: 00299 raise TypeError, 'A list for inpmask is not allowed for mode=expand' 00300 00301 # input image info, actually this will be output coordinates 00302 #print "inpimage=",inpimage," is exist?=",os.path.isdir(inpimage) 00303 ia.open(inpimage) 00304 inshp = ia.shape() 00305 incsys = ia.coordsys() 00306 ia.close() 00307 # prepare working input image (tmp_maskiamge) 00308 if inpmask!='': # inpmask is either image mask or T/F mask now 00309 # need to extract the mask and put in tmp_maskimage 00310 (parentimage,bmask)=extractmaskname(inpmask) 00311 if bmask!='': 00312 pixelmask2cleanmask(imagename=parentimage, maskname=bmask, maskimage=tmp_maskimage, usemasked=True) 00313 #ia.open(tmp_maskimage) 00314 else: 00315 #print "parentimage=",parentimage, " exist?=",os.path.isdir(parentimage) 00316 # copy of inpimage in tmp_maskimage 00317 ia.fromimage(outfile=tmp_maskimage, infile=parentimage) 00318 else: 00319 raise Exception, "inpmask must be specified" 00320 if ia.isopen(): ia.close() 00321 #setting up the output image (copy from inpimage or template) 00322 if not os.path.isdir(outparentim): 00323 #shutil.copytree(inpimage,tmp_outmaskimage) 00324 ia.fromshape(outfile=tmp_outmaskimage,shape=inshp, csys=incsys.torecord()) 00325 ia.close() 00326 needtoregrid=False 00327 else: 00328 shutil.copytree(outparentim,tmp_outmaskimage) 00329 00330 # if inpfreq/outfreq are channel indices (int) then 00331 # regrid in x y coords only and extract specified channel mask 00332 # to specified output channels. (no regriding in spectral axis) 00333 # if inpfreqs/outfreqs are velocity or freqs, 00334 # it assumes it is expressed in the range with minval~maxval 00335 # create subimage of the input mask with the range, 00336 # do regrid with the subimage to output. 00337 00338 # decide to regrid or not 00339 # 1. the case all channels are selected for input and output, simply regrid 00340 # 2. if inpfreqs and outfreqs are integers (= channel indices), regrid only in 00341 # first and second axes (e.g. ra,dec) and no regridding along spectral axis 00342 # 3. if inpfreqs and outfreqs are ranges in freq or vel, make subimage and regrid 00343 00344 if inshp[3]!=1 and ((inpfreqs==[] and outfreqs==[]) \ 00345 or (inpfreqs=='' and outfreqs=='')): 00346 # unless inpimage is continuum, skip chan selection part and regrid 00347 needtoregrid=True 00348 # detach input(tmp) image and open output tmp image 00349 # ia.open(tmp_outmaskimage) 00350 else: 00351 # if ia.isopen(): 00352 # if ia.name(strippath=True)!=tmp_maskimage: 00353 # ia.close() 00354 # ia.open(tmp_maskimage) 00355 # else: 00356 # ia.open(tmp_maskimage) 00357 00358 if inshp[3]!=1: casalog.post("inpmask is continuum..","INFO") 00359 # selection by channel indices (number) 00360 # if both inpfreqs and outfreqs are int skip regridding 00361 # if outfreqs is vel or freq ranges, try regridding 00362 if inpfreqs==[[]] or inpfreqs==[]: 00363 # select all channels for input 00364 inpfreqs=range(inshp[3]) 00365 00366 # check inpfreqs and outfreqs types 00367 selmode='bychan' 00368 if type(inpfreqs)==list: 00369 if type(inpfreqs[0])==int: 00370 if type(outfreqs)==list and (len(outfreqs)==0 or type(outfreqs[0])==int): 00371 selmode='bychan' 00372 elif type(outfreqs)==str: 00373 if inpfreqs[0]==0: #contintuum -allow index-type specification 00374 selmode='byvf' 00375 else: 00376 raise TypeError, "Mixed types in infreqs and outfreqs are not allowed" 00377 else: 00378 raise TypeError, "Non-integer in inpfreq is not supported" 00379 elif type(inpfreqs)==str: 00380 if type(outfreqs)!=str: 00381 raise TypeError, "Mixed types in infreqs and outfreqs" 00382 selmode='byvf' 00383 else: 00384 raise TypeError, "Wrong type for infreqs" 00385 00386 # inpfreqs and outfreqs are list of int 00387 # match literally without regridding. 00388 if selmode=='bychan': 00389 casalog.post("selection of input and output ranges by channel") 00390 00391 if ia.isopen(): 00392 ia.close() 00393 if outfreqs==[] or outfreqs==[[]]: 00394 outchans=[] 00395 else: 00396 outchans=outfreqs 00397 expandchanmask(tmp_maskimage,inpfreqs,tmp_outmaskimage,outchans) 00398 #print "done expandchanmask tmp_outmaskimge=",tmp_outmaskimage, " is exist?", os.path.isdir(tmp_outmaskimage) 00399 ia.open(tmp_outmaskimage) 00400 00401 elif selmode=='byvf': # outfreqs are quantities (freq or vel) 00402 casalog.post("selection of input/output ranges by frequencies/velocities") 00403 00404 inpfreqlist = translatefreqrange(inpfreqs,incsys) 00405 # close input image 00406 if ia.isopen(): 00407 ia.close() 00408 00409 #regrid to output image coordinates 00410 if len(inpfreqlist)==1: # continuum 00411 #do not regrid, use input image 00412 shutil.copytree(tmp_maskimage,tmp_regridim) 00413 else: 00414 regridmask(tmp_maskimage,inpimage,tmp_regridim,chanrange=inpfreqlist) 00415 # find edge masks (nonzero planes) 00416 ia.open(tmp_regridim) 00417 sh=ia.shape() 00418 chanlist=range(sh[3]) 00419 indlo=0 00420 indhi=0 00421 for i in chanlist: 00422 sl1=[0,0,0,i] 00423 sl2=[sh[0]-1,sh[1]-1,sh[2]-1,i] 00424 psum = ia.getchunk(sl1,sl2).sum() 00425 pmsum = ia.getchunk(sl1,sl2,getmask=True).sum() 00426 if pmsum!=0 and psum>0.0: 00427 indlo=i 00428 break 00429 chanlist.reverse() 00430 for i in chanlist: 00431 sl1=[0,0,0,i] 00432 sl2=[sh[0]-1,sh[1]-1,sh[2]-1,i] 00433 psum = ia.getchunk(sl1,sl2).sum() 00434 if psum>0.0: 00435 indhi=i 00436 break 00437 if indhi < indlo: 00438 raise IOError, "Incorrectly finding edges of input masks! Probably some logic error in the code!!!" 00439 else: 00440 casalog.post("Determined non-zero channel range to be "+str(indlo)+"~"+str(indhi), 'DEBUG1') 00441 00442 # find channel indices for given outfreqs 00443 ia.open(tmp_outmaskimage) 00444 ocsys=ia.coordsys() 00445 oshp=ia.shape() 00446 outfreqlist = translatefreqrange(outfreqs,ocsys) 00447 rtn=ocsys.findcoordinate('spectral') 00448 px=rtn[1][0] 00449 wc=rtn[2][0] 00450 world=ocsys.referencevalue() 00451 # assume chanrange are in freqs 00452 world['numeric'][wc]=qa.convert(qa.quantity(outfreqlist[0]),'Hz')['value'] 00453 p1 = ocsys.topixel(world)['numeric'][px] 00454 world['numeric'][wc]=qa.convert(qa.quantity(outfreqlist[1]),'Hz')['value'] 00455 p2 = ocsys.topixel(world)['numeric'][px] 00456 casalog.post("translated channel indices:"+qa.tos(outfreqlist[0])+"->p1="+str(p1)+\ 00457 " "+qa.tos(outfreqlist[0])+"-> p2="+str(p2)) 00458 if len(inpfreqs)==1: 00459 inpfreqchans=inpfreqs 00460 else: 00461 inpfreqchans=[indlo,indhi] 00462 outfreqchans=range(int(round(p1)),int(round(p2))+1) 00463 #print "inpfreqchans=",inpfreqchans 00464 #print "outfreqchans=",outfreqchans 00465 00466 expandchanmask(tmp_regridim,inpfreqchans,tmp_outmaskimage,outfreqchans) 00467 00468 # usechanims={} # list of input mask to be use for each outpfreq 00469 # for i in outfreqchans: 00470 # nearestch = findnearest(inpfreqchans,i) 00471 # usechanims[i]=nearestch 00472 # # put masks from inp image channel by channel 00473 # for j in outfreqs: 00474 # pix = refchanchunk[usechanims[j]-refchanst] 00475 # #ia.putchunk(pixels=pix,blc=[inshp[0]-1,inshp[1]-1,0,j]) 00476 # ia.putchunk(pixels=pix.transpose(),blc=[0,0,0,j]) 00477 needtoregrid=False 00478 ia.open(tmp_outmaskimage) 00479 # 00480 00481 if needtoregrid: 00482 # closing current output image 00483 if ia.isopen(): 00484 ia.close() 00485 #print "tmp_maskimage=",tmp_maskimage, " is exist?", os.path.isdir(tmp_maskimage) 00486 ia.open(tmp_maskimage) 00487 #os.system('cp -r %s beforeregrid.im' % tmp_maskimage) 00488 if os.path.isdir(tmp_outmaskimage): 00489 #print "Removing %s" % tmp_outmaskimage 00490 shutil.rmtree(tmp_outmaskimage) 00491 regridmask(tmp_maskimage,outparentim,tmp_outmaskimage) 00492 ia.remove() 00493 #print "closing after regrid" 00494 ia.open(tmp_outmaskimage) # reopen output tmp image 00495 00496 if outbmask!='': 00497 #convert the image mask to T/F mask 00498 casalog.post("Convert the image mask to T/F mask",'INFO') 00499 #ia.calcmask(mask='%s<0.5' % tmp_outmaskimage,name=outmask,asdefault=True) 00500 ia.calcmask(mask='%s=0.0' % tmp_outmaskimage,name=outbmask,asdefault=True) 00501 00502 if storeinmask: 00503 ia.open(outparentim) 00504 ia.maskhandler('copy',[tmp_outmaskimage+':'+outbmask, outbmask]) 00505 ia.maskhandler('set',outbmask) 00506 ia.done() 00507 casalog.post("Output the mask to %s in %s" % (outbmask,outparentim) ,"INFO") 00508 else: 00509 ow = False 00510 if inpimage==output: 00511 casalog.post("Updating "+output+"with new mask","INFO") 00512 else: 00513 if os.path.isdir(outparentim): 00514 casalog.post(outparentim+" exists, overwriting","INFO") 00515 ow=True 00516 else: 00517 casalog.post("Output the mask to "+outparentim ,"INFO") 00518 ia.rename(outparentim,ow) 00519 ia.done() 00520 00521 except Exception, instance: 00522 print "*** Error ***", instance 00523 if ia.isopen(): 00524 ia.close() 00525 ia.done() 00526 raise Exception, instance 00527 finally: 00528 if os.path.exists(tmp_maskimage): 00529 shutil.rmtree(tmp_maskimage) 00530 if os.path.exists(tmp_regridim): 00531 shutil.rmtree(tmp_regridim) 00532 00533 # === now called copy mode: copy/merge mode === 00534 # copy is a just special case of merge mode 00535 # CHANGE: 00536 # all input masks should be specified in inpmask 00537 # type of inpmask accepted: 1/0 mask, T/F mask, region file, and region expression in a string 00538 # already stored internally in seperate lists 00539 # rgfiles - region files 00540 # imgfiles - 1/0 image masks 00541 # rglist - region expression in strings 00542 # bmasks - T/F internal masks 00543 # 00544 # avaialble parameters: inpimage (string) , inpmask (list/string), output(string) 00545 # input inpimage as a template or image used for defining regions when it is given in inpmask 00546 # inpmask as list of all the masks to be merged (image masks, T/F internal masks, regions) 00547 00548 #was: if mode=='merge': 00549 if mode=='copy': 00550 #print "new copy mode..." 00551 sum_tmp_outfile='__tmp_outputmask' 00552 tmp_inmask='__tmp_frominmask' 00553 usedimfiles=[] 00554 usedbmasks=[] 00555 usedrgfiles=[] 00556 usedrglist=[] 00557 try: 00558 # check outparentim - image part of output and set as a template image 00559 if not os.path.isdir(outparentim) or (outparentim==inpimage): 00560 # figure out which input mask to be used as template 00561 # if inpimage is defined use the first one else try the first one 00562 # inpmask 00563 #if output=='': 00564 # if type(inpimage)==list: 00565 # raise Exception, "inputimage must be a string" 00566 # elif type(inpimage)==str: 00567 # outimage=inpimage 00568 # casalog.post("No outimage is specified. Will overwrite input image: "+outimage,'INFO') 00569 00570 #if type(inpimage)==list and len(inpimage)!=0: 00571 # tmp_template=inpimage[0] 00572 #elif inpimage!='' and inpimage!=[]: 00573 # tmp_template=inpimage # string 00574 #tmp_template=inpimage # string 00575 #else: 00576 # if type(inpmask)==list and len(inpmask)!=0: 00577 # fsep=inpmask[0].rfind(':') 00578 # if fsep ==-1: 00579 # raise IOError, "Cannot resolve inpmask name, check the format" 00580 # else: 00581 # tmp_template=inpmask[0][:inpmask[0].rfind(':')] 00582 # elif inpmask!='' and inpmask!=[]: 00583 # # this is equivalent to 'copy' the inpmask 00584 # tmp_template=inpmask #string 00585 00586 # create an empty image with the coords from inpimage 00587 makeEmptyimage(inpimage,sum_tmp_outfile) 00588 #print "making an empty image from inpimage to sum_tmp_outfile" 00589 else: 00590 #use output image - does not do zeroeing out, so output image is only modified 00591 shutil.copytree(outparentim,sum_tmp_outfile) 00592 #print "Using outpaarentim=",outparentim, " as sum_tmp_outfile=",sum_tmp_outfile 00593 00594 #if type(inpimage)==str: 00595 # inpimage=[inpimage] 00596 #if type(inpmask)==str: 00597 # inpmask=[inpmask] 00598 00599 if len(imgfiles)>0: 00600 # summing all the images 00601 casalog.post('Summing all mask images in inpimage and normalized to 1 for mask','INFO') 00602 for img in imgfiles: 00603 tmpregrid='__tmp_regrid.'+img 00604 # regrid to output image coords 00605 regridmask(img,sum_tmp_outfile,tmpregrid) 00606 #print "addimagemask... for ",tmpregrid 00607 addimagemask(sum_tmp_outfile,tmpregrid) 00608 usedimfiles.append(img) 00609 shutil.rmtree(tmpregrid) 00610 # get boolean masks 00611 # will work in image(1/0) masks 00612 00613 if len(bmasks)>0: 00614 casalog.post('Summing all T/F mask in inpmask and normalized to 1 for mask','INFO') 00615 for msk in bmasks: 00616 (imname,mskname) = extractmaskname(msk) 00617 #if msk.find(':')<0: 00618 # # assume default mask 00619 # msk=msk+':mask0' 00620 #imname=msk[:msk.rfind(':')] 00621 ia.open(imname) 00622 inmasks=ia.maskhandler('get') 00623 ia.close() 00624 #print "inmasks=",str(inmasks) 00625 if not inmasks.count(mskname): 00626 raise TypeError, mskname+" does not exist in "+imname+" -available masks:"+str(inmasks) 00627 # move T/F mask to image mask 00628 pixelmask2cleanmask(imname, mskname, tmp_inmask, True) 00629 regridmask(tmp_inmask,sum_tmp_outfile,'__tmp_fromTFmask') 00630 addimagemask(sum_tmp_outfile,'__tmp_fromTFmask') 00631 usedbmasks.append(msk) 00632 shutil.rmtree('__tmp_fromTFmask') 00633 # if overwriting to inpimage, delete the boolean mask 00634 if outparentim==inpimage and inpimage==imname: 00635 ia.open(imname) 00636 ia.maskhandler('delete',[mskname]) 00637 ia.close() 00638 ia.open(imname) 00639 ia.close() 00640 00641 00642 if len(rgfiles)>0 or len(rglist)>0: 00643 tmp_allrgmaskim='__tmp_fromAllRgn' 00644 tmp_rgmaskim='__tmp_fromRgn' 00645 # create an empty image with input image coords. 00646 #print "Using %s as a template for regions" % inpimage 00647 ia.open(inpimage) 00648 tshp=ia.shape() 00649 tcsys=ia.coordsys() 00650 ia.fromshape(tmp_allrgmaskim,shape=tshp, csys=tcsys.torecord(),overwrite=True) 00651 ia.done() 00652 if os.path.isdir(tmp_rgmaskim): 00653 shutil.rmtree(tmp_rgmaskim) 00654 shutil.copytree(tmp_allrgmaskim,tmp_rgmaskim) 00655 00656 if len(rgfiles)>0: 00657 #print "Has rgfiles..." 00658 nrgn=0 00659 for rgn in rgfiles: 00660 firstline=True 00661 subnrgn=0 00662 with open(rgn) as f: 00663 for line in f: 00664 if firstline: 00665 if line.count('#CRTF')==0: 00666 raise Exception, "Input text file does not seems to be in a correct format \ 00667 (must contains #CRTF)" 00668 firstline=False 00669 else: 00670 try: 00671 if line.count('^#'): 00672 pass 00673 else: 00674 if len(line)!=0: 00675 # reset temp mask image 00676 ia.open(tmp_rgmaskim) 00677 ia.set(pixels=0.0) 00678 ia.close() 00679 #print "tshp=",tshp 00680 #print "tcsys.torecord=",tcsys.torecord() 00681 inrgn=rg.fromtext(line, tshp, tcsys.torecord()) 00682 #print "inrgn=",inrgn 00683 im.regiontoimagemask(tmp_rgmaskim,region=inrgn) 00684 addimagemask(tmp_allrgmaskim,tmp_rgmaskim) 00685 #shutil.rmtree(tmp_rgmaskim) 00686 subnrgn +=1 00687 nrgn +=1 00688 except: 00689 break 00690 if subnrgn>0: 00691 usedrgfiles.append(rgn) 00692 casalog.post("Converted %s regions from %s region files to image mask" % (nrgn,len(rgfiles)),"INFO") 00693 00694 if len(rglist)>0: 00695 #print "Has rglist..." 00696 nrgn=0 00697 for rgn in rglist: 00698 # reset temp mask image 00699 ia.open(tmp_rgmaskim) 00700 ia.set(pixels=0.0) 00701 ia.close() 00702 inrgn=rg.fromtext(rgn, tshp, tcsys.torecord()) 00703 im.regiontoimagemask(tmp_rgmaskim,region=inrgn) 00704 addimagemask(tmp_allrgmaskim,tmp_rgmaskim) 00705 #shutil.rmtree(tmp_rgmaskim) 00706 usedrglist.append(rgn) 00707 nrgn+=1 00708 casalog.post("Converted %s regions to image mask" % (nrgn),"INFO") 00709 00710 00711 # regrid if necessary 00712 regridded_mask='__tmp_regrid_allrgnmask' 00713 regridmask(tmp_allrgmaskim, sum_tmp_outfile,regridded_mask) 00714 addimagemask(sum_tmp_outfile,regridded_mask) 00715 #shutil.rmtree('__tmp_regridded_allrgnmask') 00716 casalog.post("Added mask based on regions to output mask","INFO") 00717 #cleanup 00718 for tmpfile in [tmp_allrgmaskim,tmp_rgmaskim,regridded_mask]: 00719 if os.path.isdir(tmpfile): 00720 shutil.rmtree(tmpfile) 00721 00722 00723 if outbmask!='': 00724 casalog.post('Putting mask in T/F','INFO') 00725 ia.open(sum_tmp_outfile) 00726 #ia.calcmask(mask='%s<0.5' % sum_tmp_outfile,name=outmask,asdefault=True) 00727 # mask only pixel != 0.0 00728 ia.calcmask(mask='%s==0.0' % sum_tmp_outfile,name=outbmask,asdefault=True) 00729 ia.done() 00730 # if outfile exists initially outfile is copied to sum_tmp_outfile 00731 # if outfile does not exist initially sum_tmp_outfile is a copy of inpimage 00732 # so rename it with overwrite=T all the cases 00733 #print "open sum_tmp_outfile=",sum_tmp_outfile 00734 if storeinmask: 00735 ia.open(outparentim) 00736 ia.maskhandler('copy',[sum_tmp_outfile+':'+outbmask, outbmask]) 00737 ia.maskhandler('set',outbmask) 00738 ia.done() 00739 outputmsg="to create an output mask: %s in %s" % (outbmask,outparentim) 00740 else: 00741 ia.open(sum_tmp_outfile) 00742 ia.rename(outparentim,overwrite=True) 00743 ia.done() 00744 outputmsg="to create an output mask: %s " % outparentim 00745 00746 casalog.post("Merged masks from:","INFO") 00747 if len(usedimfiles)>0: 00748 casalog.post("mask image(s): "+str(usedimfiles),"INFO") 00749 if len(usedbmasks)>0: 00750 casalog.post("internal mask(s): "+str(usedbmasks),"INFO") 00751 if len(usedrgfiles)>0: 00752 casalog.post("region txt file(s): "+str(usedrgfiles),"INFO") 00753 if len(usedrglist)>0: 00754 casalog.post("region(s) from direct input: "+str(usedrglist),"INFO") 00755 casalog.post(outputmsg,"INFO") 00756 00757 00758 except Exception, instance: 00759 print "*** Error ***", instance 00760 raise Exception, instance 00761 finally: 00762 if os.path.exists(sum_tmp_outfile): 00763 shutil.rmtree(sum_tmp_outfile) 00764 if os.path.exists(tmp_inmask): 00765 shutil.rmtree(tmp_inmask) 00766 if type(inpimage)==list: 00767 for im in inpimage: 00768 if os.path.isdir('__tmp_regrid.'+im): 00769 shutil.rmtree('__tmp_regrid.'+im) 00770 00771 00772 00773 # === draw mode === 00774 # disabled - drawmaskinimage (working with viewer) a bit flaky 00775 # when run in succession. 00776 # if mode=='draw': 00777 # #call drawmaskinimage 00778 # from recipes.drawmaskinimage import drawmaskinimage 00779 # drawmaskinimage(inpimage,outmask) 00780 00781 00782 except Exception, instance: 00783 print '*** Error ****', instance 00784 raise Exception, instance 00785 00786 finally: 00787 # final clean up 00788 if os.path.isdir(tmp_maskimage): 00789 shutil.rmtree(tmp_maskimage) 00790 if os.path.isdir(tmp_outmaskimage): 00791 shutil.rmtree(tmp_outmaskimage) 00792 if os.path.isdir(tmp_regridim): 00793 shutil.rmtree(tmp_regridim) 00794 00795 def findnearest(arr, val): 00796 import numpy as np 00797 if type(arr)==list: 00798 arr = np.array(arr) 00799 indx = np.abs(arr - val).argmin() 00800 return arr[indx] 00801 00802 def regridmask(inputmask,template,outputmask,axes=[3,0,1],method='linear',chanrange=None): 00803 ''' 00804 Regrid input mask (image) to output mask using a template. 00805 Currently the template must be a CASA image. 00806 The default interpolation method is set to 'linear' (since 'nearest' 00807 sometime fails). 00808 ''' 00809 #print "Regrid.." 00810 #print "inputmask=",inputmask," template=",template," outputmask=",outputmask 00811 if not os.path.isdir(template): 00812 raise IOError, "template image %s does not exist" % template 00813 00814 (ia,) = gentools(['ia']) 00815 ia.open(template) 00816 ocsys = ia.coordsys() 00817 oshp = ia.shape() 00818 ia.done() 00819 ia.open(inputmask) 00820 # check axis order, if necessary re-interprete input axes correctly 00821 # assumed order of axes 00822 reforder=['Right Ascension', 'Declination', 'Stokes', 'Frequency'] 00823 axisorder=ia.summary(list=False)['axisnames'].tolist() 00824 tmp_axes=[] 00825 for axi in axes: 00826 tmp_axes.append(axisorder.index(reforder[axi])) 00827 axes=tmp_axes 00828 if type(chanrange)==list and len(chanrange)==2: 00829 incsys=ia.coordsys() 00830 spaxis=incsys.findcoordinate('spectral')[2] 00831 # create subimage based on the inpfreqs range 00832 inblc=chanrange[0] 00833 intrc=chanrange[1] 00834 casalog.post("Regridmask: spaxis=%s, inblc=%s, intrc=%s" % (spaxis,inblc,intrc), 'DEBUG1') 00835 rgn = rg.wbox(blc=inblc,trc=intrc,pixelaxes=spaxis.tolist(),csys=incsys.torecord()) 00836 else: 00837 rgn={} 00838 # for continuum case 00839 if oshp[tmp_axes[0]]==1: 00840 axes=[0,1] 00841 ir=ia.regrid(outfile=outputmask,shape=oshp,csys=ocsys.torecord(),axes=axes,region=rgn,method=method) 00842 ia.done() 00843 # to ensure to create 1/0 mask image 00844 #ir.calc('iif (%s>0.0 && %s<1.0,1,%s)'%(outputmask,outputmask,outputmask)) 00845 # treat everything not = 0.0 to be mask 00846 ir.calc('iif (abs(%s)>0.0,1,%s)'%(outputmask,outputmask)) 00847 ir.done() 00848 00849 def addimagemask(sumimage, imagetoadd, threshold=0.0): 00850 """ 00851 add image masks (assumed the images are already in the same coordinates) 00852 """ 00853 (ia,) = gentools(['ia']) 00854 #print "addimagemask: sumimage=",sumimage," imagetoadd=",imagetoadd 00855 ia.open(sumimage) 00856 ia.calc('iif ('+imagetoadd+'>'+str(threshold)+',('+sumimage+'+'+imagetoadd+')/('+sumimage+'+'+imagetoadd+'),'+sumimage+')') 00857 ia.close() 00858 00859 def expandchanmask(inimage,inchans,outimage,outchans): 00860 """ 00861 expand masks in channel direction,and insert then 00862 to output image with the same coordinates (post-regridded) 00863 only differ by channels 00864 """ 00865 # input image 00866 ia.open(inimage) 00867 inshp=ia.shape() 00868 refchanst=inchans[0] 00869 refchanen=inchans[-1] 00870 slst = [0,0,0,refchanst] 00871 slen = [inshp[0]-1,inshp[1]-1,0,refchanen] 00872 casalog.post("getting chunk at blc="+str(slst)+" trc="+str(slen),'DEBUG1') 00873 refchanchunk=ia.getchunk(blc=slst,trc=slen) 00874 refchanchunk=refchanchunk.transpose() 00875 ia.close() 00876 #print "refchanchunk:shape=",refchanchunk.shape 00877 00878 ia.open(outimage) 00879 # need find nearest inchan 00880 # store by chan indices (no regrid) 00881 outshp=ia.shape() 00882 if outchans==[]: 00883 #select all channels 00884 outchans=range(outshp[3]) 00885 usechanims={} # list of input mask to be use for each outpfreq 00886 for i in outchans: 00887 nearestch = findnearest(inchans,i) 00888 usechanims[i]=nearestch 00889 casalog.post("Mapping of channels: usechanims="+str(usechanims),'DEBUG1') 00890 for j in outchans: 00891 pix = refchanchunk[usechanims[j]-refchanst] 00892 #print "pix=",pix 00893 #print "pix.shape=",pix.shape 00894 #print "inshp=",inshp, ' j=',j 00895 #ia.putchunk(pixels=pix,blc=[inshp[0]-1,inshp[1]-1,0,j]) 00896 ia.putchunk(pixels=pix.transpose(),blc=[0,0,0,j]) 00897 #print "DONE putchunk for j=", j 00898 ia.done() 00899 00900 def translatefreqrange(freqrange,csys): 00901 """ 00902 convert the range in list 00903 mainly for frequeny and velocity range determination 00904 """ 00905 if type(freqrange)==list and type(freqrange[0])==int: 00906 #do nothing 00907 return freqrange 00908 elif type(freqrange)==str: 00909 freqlist=freqrange.split('~') 00910 for i in range(len(freqlist)): 00911 if freqlist[i].find('m/s') > -1: 00912 fq = qa.quantity(freqlist[i]) 00913 vf=csys.velocitytofrequency(value=fq['value'],velunit=fq['unit']) 00914 freqlist[i]=str(vf[0])+'Hz' 00915 return freqlist 00916 else: 00917 raise TypeError, "Cannot understand frequency range" 00918 00919 def checkinput(inpname): 00920 """ 00921 do existance check on image and internal mask 00922 """ 00923 (parentimage,tfmaskname)=extractmaskname(inpname) 00924 (parentimexist,tfmaskexist)=checkinmask(parentimage,tfmaskname) 00925 if parentimexist: 00926 if tfmaskname=='': 00927 return True # only the image 00928 else: 00929 if not tfmaskexist: 00930 ia.open(parentimage) 00931 inmasklist=ia.maskhandler('get') 00932 raise Exception, "Cannot find the internal mask, %s. Candidate mask(s) are %s" % (tfmaskname, str(inmasklist)) 00933 else: 00934 return True # image mask and internal mask 00935 else: 00936 raise Exception, "Cannot find the image=%s" % parentimage 00937 00938 00939 def checkinmask(parentimage,tfmaskname): 00940 """ 00941 check existance of the internal mask 00942 """ 00943 if os.path.isdir(parentimage): 00944 if tfmaskname!='': 00945 ia.open(parentimage) 00946 inmasks=ia.maskhandler('get') 00947 ia.close() 00948 if not any(tfmaskname in msk for msk in inmasks): 00949 return (True, False) 00950 else: 00951 return (True, True) # image mask and internal mask 00952 else: 00953 return (True, False) 00954 else: 00955 return (False,False) 00956 00957 def extractmaskname(maskname): 00958 """ 00959 split out imagename and maskname from a maskname string 00960 returns (parentimage, internalmask) 00961 """ 00962 # the image file name may contains ':' some cases 00963 # take last one in split list as an internal mask name 00964 00965 indx = maskname.find(':') 00966 for i in range(len(maskname)): 00967 if indx>-1: 00968 indx += maskname[indx+1:].find(':') 00969 indx +=1 00970 else: 00971 break 00972 if indx != -1: 00973 parentimage=maskname[:indx] 00974 maskn=maskname[indx+1:] 00975 return (parentimage,maskn) 00976 else: 00977 parentimage=maskname 00978 return (parentimage, '') 00979 00980 def makeEmptyimage(template,outimage): 00981 """ 00982 make an empty image with the coords 00983 from template 00984 """ 00985 (ia,) = gentools(['ia']) 00986 ia.open(template) 00987 inshp=ia.shape() 00988 incsys=ia.coordsys() 00989 ia.fromshape(outimage,shape=inshp,csys=incsys.torecord()) 00990 ia.done()
1.8.0
