Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_makemask.py: 83%

1################################################ 

2# Task to make masks. 

3# reorganized after diesuccsion on Oct1, 2012 

4# with Jeurgen and Urvashi 

5# 

6# modified by TT 

7# based on the original code,  

8# v1.0: 2012.03.20, U.Rau 

9# 

10################################################ 

11# Notes (to self) - TT  

12# 1. expanding one mask to another  

13# e.g.) expanding a continuum mask (both image mask/boolean mask) 

14# channel mask  

15# 2. part of copy mode func.: merging of different types of masks  

16# e.g.) inpimage and inpmask are lists of the mask to be merged 

17# output mask is written to either outimage or outmask as embedded 

18# T/F mask  

19# 3. copying mask to another or create a new one 

20# regrid if necessary (i.e. if the coords are different)  

21# ---------------------------------------------- 

22# basic rules: 

23# for mask image (1/0 mask): as it is 

24# for internal mask : parent_imagename:mask_name 

25# 

26# For input, 

27# inpimage is the casa image 

28# - mode='list': list internal masks of inpimage 

29# - other mode: used as a template for output 

30# if region files are specified -> make mask specifeid with the regions on to inpimage 

31# output is '' => modified inpimage unless overwrite=F else exception 

32# 

33# if inpmask='': use inpimage as input mask 

34# if inpmask='mask0' or other embedded mask name of the inpimage,  

35# use that T/F mask 

36#  

37# =expand= 

38# case1: on the same image (outimage=''), expand mask image from  

39# prev. run etc. No regriding. Use nearest chan mask 

40# image to expand.  

41# 1.a: inpimage is clean image mask (1s and 0s) 

42# i) outimage != inpimage and outmask='' => new expanded mask image to outimage 

43# ii) outimage != inpimage and outmask!='' => convert expanded mask image to T/F mask to store inside outimage 

44# iii) outimage ==inpimage and outmask='' => update input mask image by expanding it  

45# iv) outimage ==inpimage and outmask!=''=> update input image with the expanded T/F mask 

46# 1.b: if inpmask!='', do T/F mask to 1/0 image mask conversion, then do as in 1.a  

47 

48# case2: outimage is in diffirent coords. (need to regrid) 

49# 

50################# 

51# tests  

52# 1. for input: mask image  

53# 2. for input: mask/or regular image with internal mask 

54# 3. for input: mask image; for output: mask image with different spectral grid 

55# 4. for input: mask/regular image with internal mask; for output: image with  

56# internal mask with different spectral grid 

57################### 

58import os 

59import shutil 

60import numpy as np 

61 

62from casatools import image, regionmanager, imager, table, quanta 

63from casatasks import casalog 

64from .imtools import pixelmask2cleanmask 

65 

66import subprocess 

67subprocess_getoutput = subprocess.getoutput 

68 

69_ia = image() 

70_rg = regionmanager() 

71_qa = quanta() 

72 

73pid = str(os.getpid()) 

74debug = False 

75#debug = True  

76 

77def makemask(mode,inpimage, inpmask, output, overwrite, inpfreqs, outfreqs): 

78 """ 

79 make /manipulate masks 

80 

81 """ 

82 _ia = image() 

83 _rg = regionmanager() 

84 _im = imager() 

85 casalog.origin('makemask') 

86 #casalog.post("params(mode,inpimage,inpmask,output,overwrite)=",mode,inpimage,inpmask,output,overwrite) 

87 

88 try: 

89 # temp files 

90 tmp_maskimage='__tmp_makemaskimage_'+pid 

91 tmp_outmaskimage='__tmp_outmakemaskimage_'+pid 

92 tmp_regridim='__tmp_regridim_'+pid 

93 

94 #intial cleanup to make sure nothing left from previous runs 

95 tmplist = [tmp_maskimage,tmp_outmaskimage,tmp_regridim] 

96 cleanuptempfiles(tmplist) 

97 

98 # do parameter check first 

99 # check names of inpimage, inpmask check for existance 

100 # inpimage == output (exact match) then check overwrite 

101 # => T overwrite inpimage 

102 # => F exception 

103 

104 # check inpimage  

105 if (['list','copy','expand'].count(mode)==1): 

106 if inpimage=='': raise ValueError("inpimage is empty") 

107 if not os.path.isdir(inpimage): 

108 raise RuntimeError("inpimage=%s does not exist" % inpimage) 

109 

110 # === list mode === 

111 if mode == 'list': 

112 inpOK=checkinput(inpimage) 

113 if inpOK: 

114 

115 if _ia.isopen(): _ia.close() 

116 _ia.open(inpimage) 

117 inmasklist=_ia.maskhandler('get') 

118 # now ia.maskhandler returns ['T'] if no internal mask is there... 

119 if inmasklist.count('T')!=0: 

120 inmasklist.remove('T') 

121 if len(inmasklist) ==0: 

122 casalog.post('No internal (T/F) masks were found in %s' % (inpimage),'INFO') 

123 else: 

124 defaultmaskname=_ia.maskhandler('default')[0] 

125 printinmasks='' 

126 for mname in inmasklist: 

127 if mname==defaultmaskname: 

128 printinmasks+='\''+mname+'\''+'(default)' 

129 else: 

130 printinmasks+='\''+mname+'\'' 

131 if mname != inmasklist[-1]: 

132 printinmasks+=', ' 

133 

134 casalog.post('Internal (T/F) masks in %s: %s' % (inpimage, printinmasks),'INFO') 

135 _ia.close() 

136 

137 # === setdefaultmask mode === 

138 elif mode == 'setdefaultmask': 

139 inpOK=checkinput(inpmask) 

140 if inpOK: 

141 (parentimage,bmask)=extractmaskname(inpmask) 

142 if bmask=='': 

143 raise RuntimeError("Missing an internal mask name") 

144 if _ia.isopen(): _ia.close() 

145 _ia.open(parentimage) 

146 defaultmaskname=_ia.maskhandler('default')[0] 

147 inmasklist=_ia.maskhandler('get') 

148 if defaultmaskname==bmask: 

149 casalog.post('No change. %s is already a default internal mask' % bmask, 'INFO') 

150 else: 

151 _ia.maskhandler('set',bmask) 

152 casalog.post('Set %s as a default internal mask' % bmask, 'INFO') 

153 if len(inmasklist)>1: 

154 casalog.post('Current internal masks are %s' % str(inmasklist), 'INFO') 

155 _ia.close() 

156 

157 # === delete mode === 

158 elif mode == 'delete': 

159 inpOK=checkinput(inpmask) 

160 if inpOK: 

161 (parentimage,bmask)=extractmaskname(inpmask) 

162 if bmask=='': 

163 raise RuntimeError("Missing an internal mask name") 

164 _ia.open(parentimage) 

165 casalog.post('Deleting the internal mask, %s ' % bmask, 'INFO') 

166 defaultmaskname=_ia.maskhandler('default')[0] 

167 _ia.maskhandler('delete',bmask) 

168 inmasklist=_ia.maskhandler('get') 

169 if inmasklist.count('T')!=0: 

170 inmasklist.remove('T') 

171 if len(inmasklist) !=0 and defaultmaskname==bmask: 

172 _ia.maskhandler('set',inmasklist[0]) 

173 casalog.post('Set %s as a default internal mask' % inmasklist[0], 'INFO') 

174 if len(inmasklist)>1: 

175 casalog.post('Current internal masks are %s' % str(inmasklist), 'INFO') 

176 

177 _ia.close() 

178 

179 else: 

180 # PREPROCESS STAGE for mode='copy' and 'expand' 

181 #DEBUG 

182 #casalog.post("mode=",mode) 

183 # copy can have multiple input masks, expand has only one. 

184 # check inpimage, inpmask, output, overwrite 

185 #  

186 storeinmask = False # used to check if output to a new internal mask 

187 isNewOutfile = False 

188 

189 inpOK=checkinput(inpimage) 

190 if inpOK: 

191 (immask,inmask)=extractmaskname(inpimage) 

192 

193 # seperate text files(region files), images(with full ':'), and direct region  

194 # input mask(s) 

195 if inpmask=='': 

196 raise ValueError("Input errror. The inpmask parameter is not specified.") 

197 if type(inpmask)!=list: 

198 inpmask=[inpmask] 

199 

200 # check if inpmask contains region file or region specification 

201 rgfiles=[] 

202 imgfiles=[] 

203 rglist=[] 

204 bmasks=[] 

205 for masklet in inpmask: 

206 # is text file? 

207 if type(masklet)==str: # text file or image 

208 if os.path.exists(masklet): 

209 # region file or image 

210 if (subprocess_getoutput('file '+masklet).count('directory')): 

211 if os.path.exists(masklet+'/table.f1'): 

212 #casalog.post("%s is not in a recognized format for inpmask, ignored." % masklet, 'WARN')  

213 raise ValueError("%s is not in a recognized format for inpmask" % masklet) 

214 else: 

215 # probably image file (with no mask extension) 

216 imgfiles.append(masklet) 

217 # text file (CRTF) 

218 elif (subprocess_getoutput('file '+masklet).count('text')): 

219 rgfiles.append(masklet) 

220 else: 

221 #casalog.post("%s does not recognized format for inpmask, ignored." % masklet, 'WARN') 

222 raise ValueError("%s is not in a recognized format for inpmask" % masklet) 

223 else: 

224 # direct string specification 

225 if masklet.count('[') and masklet.count(']'): # rough check on region specification  

226 rglist.append(masklet) 

227 # extract internal mask from the input image 

228 else: 

229 (parentim, mask)=extractmaskname(masklet) 

230 if mask!='': 

231 bmasks.append(masklet) 

232 else: 

233 raise ValueError("%s is not an existing file/image or a region format" % masklet) 

234 

235 # expand allows only a string for inpmask 

236 if mode=='expand': 

237 if type(inpmask)==list: 

238 inpmask=inpmask[0] 

239 # check for overwrite condition 

240 if output=='': 

241 if overwrite: 

242 output=inpimage 

243 else: 

244 raise ValueError("output is not specified. If you want to overwrite inpimage, please set overwrite=True") 

245 

246 if inpimage==output: 

247 #if overwrite: 

248 # tmp_outmaskimage=tmp_maskimage 

249 #else: 

250 if not overwrite: 

251 raise ValueError("output=inpimage. If you want to overwrite inpimage, please set overwrite=True") 

252 

253 outparentim=output 

254 outbmask='' 

255 if os.path.isdir(output): 

256 if not overwrite: 

257 raise RuntimeError("output=%s exists. If you want to overwrite it, please set overwrite=True" % output) 

258 else: 

259 (outparentim, outbmask)=extractmaskname(output) 

260 if outbmask!='': 

261 (parentimexist,maskexist)=checkinmask(outparentim,outbmask) 

262 if parentimexist and maskexist: 

263 if not overwrite: 

264 raise ValueError("output=%s exists. If you want to overwrite it, please set overwrite=True" % output) 

265 else: 

266 casalog.post("Will overwrite the existing internal mask, %s in %s" % (outbmask,outparentim)) 

267 storeinmask=True 

268 

269 #if parentimexist and not maskexist: 

270 else: 

271 storeinmask = True 

272 if not parentimexist: isNewOutfile=True 

273 else: 

274 outparentim=output 

275 

276 #casalog.post("param checks before branching out for mode==========")) 

277 #casalog.post("storeinmask = ",storeinmask) 

278 #casalog.post("output=",output, " is exist?=",os.path.isdir(output)) 

279 #casalog.post("outparentim=",outparentim, " is exist?=",os.path.isdir(outparentim)) 

280 

281 # MAIN PROCESS for 'copy' or 'expand' mode 

282 # the following code is somewhat duplicated among the modes but keep separated from each mode 

283 # for now.... copy now handle merge as well 

284 # === old copy mode === NOW combined to 'merge mode' 

285# #if mode=='copy': 

286# #casalog.post("Copy mode") 

287# needregrid=True 

288# #if outimage=='': 

289# #overwrite 

290# # outimage=inpimage 

291# 

292# if not os.path.isdir(outimage): 

293# needregrid=False 

294# 

295# if inpmask!='': 

296# # need to extract the mask and put in tmp_maskimage 

297# pixelmask2cleanmask(imagename=inpimage, maskname=inpmask, maskimage=tmp_maskimage, usemasked=True)  

298# else: 

299# shutil.copytree(inpimage, tmp_maskimage) 

300# if needregrid: 

301# casalog.post("Regridding...",'DEBUG1') 

302# regridmask(tmp_maskimage,outimage,tmp_outmaskimage) 

303# # regrid may produce <1.0 pixels for mask so be sure to its all in 1.0 

304# #_ia.open(tmp_outmaskimage)  

305# #_ia.calc('iif (%s>0.0 && %s<1.0,1,%s)'%(tmp_outmaskimage,tmp_outmaskimage,tmp_outmaskimage)) 

306# #_ia.close() 

307# #casalog.post("Copying regrid output=",tmp_outmaskimage) 

308# else: 

309# shutil.copytree(tmp_maskimage,tmp_outmaskimage) 

310# if outmask!='': 

311# #convert the image mask to T/F mask 

312# if not os.path.isdir(outimage): 

313# shutil.copytree(inpimage,outimage) 

314# # 

315# _ia.open(outimage) 

316# casalog.post("convert the output image mask to T/F mask") 

317# _ia.calcmask(mask='%s<0.5' % tmp_outmaskimage,name=outmask,asdefault=True) 

318# _ia.done() 

319# else: 

320# # if regridded - tmp_outmaskimage is created by regridmask 

321# # if not, tmp_outmaskimage=outimage 

322# _ia.open(tmp_outmaskimage) 

323# _ia.rename(outimage,overwrite=True) 

324# _ia.done() 

325 

326 

327 # === expand mode === 

328 if mode=='expand': 

329 _rg = regionmanager() 

330 needtoregrid=False 

331 bychanindx=False 

332 

333 try: 

334 # These coordsys objects will need to be closed, if created 

335 incsys = None 

336 inmaskcsys = None 

337 ocsys = None 

338 

339 #casalog.post("expand mode main processing blocks...") 

340 # do not allow list in this mode (for inpimage and inpmask) - maybe this is redundant now 

341 if type(inpmask)==list: 

342 raise TypeError('A list for inpmask is not allowed for mode=expand') 

343 

344 # input image info, actually this will be output coordinates 

345 _ia.open(inpimage) 

346 inshp = _ia.shape() 

347 incsys = _ia.coordsys() 

348 _ia.close() 

349 #casalog.post("inpimage=",inpimage," is exist?=",os.path.isdir(inpimage)) 

350 #casalog.post(" inshp for inpimage=",inshp) 

351 

352 # prepare working input mask image (tmp_maskimage) 

353 if debug: casalog.post("prepare working input image (tmp_maskimage)...") 

354 if inpmask!='': # inpmask is either image mask or T/F mask now 

355 # need to extract the mask and put in tmp_maskimage 

356 # Note: changed usemasked=F, so that True (unmasked) part to be used. CAS-  

357 # ==> tmp_maskimage is an input mask image 

358 (parentimage,bmask)=extractmaskname(inpmask) 

359 if bmask!='': 

360 pixelmask2cleanmask(imagename=parentimage, maskname=bmask, maskimage=tmp_maskimage, usemasked=False) 

361 #_ia.open(tmp_maskimage) 

362 else: 

363 if debug: 

364 casalog.post("parentimage=" + parentimage + " exist?=" + os.path.isdir(parentimage)) 

365 casalog.post("tmp_maskimage=" + tmp_maskimage + " exist?=" + os.path.isdir(tmp_maskimage)) 

366 # copy of inpimage in tmp_maskimage 

367 _ia.fromimage(outfile=tmp_maskimage, infile=parentimage) 

368 else: 

369 raise ValueError("inpmask must be specified") 

370 if _ia.isopen(): _ia.close() 

371 #setting up the output image (copy from inpimage or template) 

372 if not os.path.isdir(outparentim): 

373 #shutil.copytree(inpimage,tmp_outmaskimage) 

374 _ia.fromshape(outfile=tmp_outmaskimage,shape=inshp, csys=incsys.torecord()) 

375 _ia.close() 

376 needtoregrid=False 

377 else: 

378 # if inpimage == output, tmp_outmaskimage is already created...  

379 if not os.path.isdir(tmp_outmaskimage): 

380 shutil.copytree(outparentim,tmp_outmaskimage) 

381 if debug: casalog.post("done setting up the out image=" + tmp_outmaskimage) 

382 # if inpfreq/outfreq are channel indices (int) then 

383 # regrid in x y coords only and extract specified channel mask 

384 # to specified output channels. (no regriding in spectral axis) 

385 # if inpfreqs/outfreqs are velocity or freqs,  

386 # it assumes it is expressed in the range with minval~maxval 

387 # create subimage of the input mask with the range, 

388 # do regrid with the subimage to output. 

389 

390 # decide to regrid or not 

391 # 1. the case all channels are selected for input and output, simply regrid 

392 # 2. if inpfreqs and outfreqs are integers (= channel indices), regrid only in 

393 # first and second axes (e.g. ra,dec) and no regridding along spectral axis 

394 # 3. if inpfreqs and outfreqs are ranges in freq or vel, make subimage and regrid 

395 _ia.open(tmp_maskimage) 

396 inmaskshp = _ia.shape() 

397 inmaskcsys = _ia.coordsys() 

398 _ia.close() 

399 regridmethod = 'linear' 

400 if 'spectral2' in inmaskcsys.torecord(): 

401 inspecdelt = inmaskcsys.torecord()['spectral2']['wcs']['cdelt'] 

402 _ia.open(tmp_outmaskimage) 

403 ocsys=_ia.coordsys() 

404 oshp=_ia.shape() 

405 _ia.close() 

406 outspecdelt = ocsys.torecord()['spectral2']['wcs']['cdelt'] 

407 if outspecdelt < inspecdelt: 

408 regridmethod='nearest' 

409 

410 if inmaskshp[3]!=1 and ((inpfreqs==[] and outfreqs==[]) \ 

411 or (inpfreqs=='' and outfreqs=='')): 

412 # unless inpimage is continuum, skip chan selection part and regrid  

413 needtoregrid=True 

414 # detach input(tmp) image and open output tmp image 

415 # _ia.open(tmp_outmaskimage) 

416 else: 

417 # if _ia.isopen(): 

418 # if _ia.name(strippath=True)!=tmp_maskimage: 

419 # _ia.close() 

420 # _ia.open(tmp_maskimage) 

421 # else: 

422 # _ia.open(tmp_maskimage) 

423 

424 #if inshp[3]!=1: casalog.post("inpmask is continuum..","INFO") 

425 if inmaskshp[3]==1: casalog.post("inpmask is continuum..","INFO") 

426 # selection by channel indices (number)  

427 # if both inpfreqs and outfreqs are int skip regridding 

428 # if outfreqs is vel or freq ranges, try regridding  

429 if inpfreqs==[[]] or inpfreqs==[]: 

430 # select all channels for input 

431 inpfreqs=list(range(inmaskshp[3])) 

432 

433 # check inpfreqs and outfreqs types 

434 # index based 

435 selmode='bychan' 

436 if type(inpfreqs)==list: 

437 if type(inpfreqs[0])==int: 

438 if type(outfreqs)==list and (len(outfreqs)==0 or type(outfreqs[0])==int): 

439 selmode='bychan' 

440 elif type(outfreqs)==str: 

441 #if inpfreqs[0]==0: #contintuum -allow index-type specification 

442 if len(inpfreqs)==1: #contintuum -allow index-type specification 

443 selmode='byvf' 

444 else: 

445 raise TypeError("Mixed types in infreqs and outfreqs are not allowed") 

446 else: 

447 raise TypeError("Non-integer in inpfreq is not supported") 

448 # by velocity or frequency 

449 elif type(inpfreqs)==str: 

450 if type(outfreqs)!=str: 

451 raise TypeError("Mixed types in infreqs and outfreqs") 

452 selmode='byvf' 

453 else: 

454 raise TypeError("Wrong type for infreqs") 

455 

456 # inpfreqs and outfreqs are list of int 

457 # match literally without regridding. 

458 if selmode=='bychan': 

459 casalog.post("selection of input and output ranges by channel") 

460 

461 if _ia.isopen(): 

462 _ia.close() 

463 if outfreqs==[] or outfreqs==[[]]: 

464 outchans=[] 

465 else: 

466 outchans=outfreqs 

467 expandchanmask(tmp_maskimage,inpfreqs,tmp_outmaskimage,outchans) 

468 _ia.open(tmp_outmaskimage) 

469 

470 elif selmode=='byvf': # outfreqs are quantities (freq or vel) 

471 casalog.post("selection of input/output ranges by frequencies/velocities") 

472 

473 # do it for input mask image (not the template ) 

474 inpfreqlist = translatefreqrange(inpfreqs,inmaskcsys) 

475 #casalog.post("inpfreqlist=",inpfreqlist) 

476 # close input image 

477 if _ia.isopen(): 

478 _ia.close() 

479 

480 #regrid to output image coordinates 

481 if len(inpfreqlist)==1: # continuum 

482 #do not regrid, use input image 

483 shutil.copytree(tmp_maskimage,tmp_regridim) 

484 else: 

485 if debug: casalog.post("regrid the mask,tmp_maskimage=" + tmp_maskimage + " tmp_regridim=" + tmp_regridim) 

486 #shutil.copytree(tmp_maskimage,'before_regrid_tmp_maskimage') 

487 regridmask(tmp_maskimage,inpimage,tmp_regridim,chanrange=inpfreqlist,method=regridmethod) 

488 #regridmask(tmp_maskimage,inpimage,tmp_regridim,chanrange=inpfreqlist) 

489 # find edge masks (nonzero planes) 

490 ##shutil.copytree(tmp_regridim,'saved_'+tmp_regridim) 

491 if _ia.isopen(): _ia.close() 

492 _ia.open(tmp_regridim) 

493 sh=_ia.shape() 

494 chanlist = list(range(sh[3])) 

495 indlo=0 

496 indhi=0 

497 for i in chanlist: 

498 sl1=[0,0,0,i] 

499 sl2=[sh[0]-1,sh[1]-1,sh[2]-1,i] 

500 psum = _ia.getchunk(sl1,sl2).sum() 

501 pmsum = _ia.getchunk(sl1,sl2,getmask=True).sum() 

502 if pmsum!=0 and psum>0.0: 

503 indlo=i 

504 break 

505 chanlist.reverse() 

506 for i in chanlist: 

507 sl1=[0,0,0,i] 

508 sl2=[sh[0]-1,sh[1]-1,sh[2]-1,i] 

509 psum = _ia.getchunk(sl1,sl2).sum() 

510 if psum>0.0: 

511 indhi=i 

512 break 

513 if indhi < indlo: 

514 raise RuntimeError("Incorrectly finding edges of input masks! Probably some logic error in the code!!!") 

515 else: 

516 casalog.post("Determined non-zero channel range to be "+str(indlo)+"~"+str(indhi), 'DEBUG1') 

517 

518 # find channel indices for given outfreqs 

519 #_ia.open(tmp_outmaskimage) 

520 #ocsys=_ia.coordsys() 

521 #oshp=_ia.shape()  

522 outfreqlist = translatefreqrange(outfreqs,ocsys) 

523 rtn=ocsys.findcoordinate('spectral') 

524 px=rtn['pixel'][0] 

525 wc=rtn['world'][0] 

526 world=ocsys.referencevalue() 

527 # assume chanrange are in freqs 

528 world['numeric'][wc]=_qa.convert(_qa.quantity(outfreqlist[0]),'Hz')['value'] 

529 p1 = ocsys.topixel(world)['numeric'][px] 

530 world['numeric'][wc]=_qa.convert(_qa.quantity(outfreqlist[1]),'Hz')['value'] 

531 p2 = ocsys.topixel(world)['numeric'][px] 

532 casalog.post("translated channel indices:"+_qa.tos(outfreqlist[0])+"->p1="+str(p1)+\ 

533 " "+_qa.tos(outfreqlist[0])+"-> p2="+str(p2)) 

534 if len(inpfreqs)==1: 

535 inpfreqchans=inpfreqs 

536 elif inpfreqs.find('~'): 

537 inpfreqchans=list(range(indlo,indhi+1)) 

538 else: 

539 inpfreqchans=[indlo,indhi] 

540 outfreqchans=list(range(int(round(p1)),int(round(p2))+1)) 

541 #casalog.post("inpfreqchans=",inpfreqchans) 

542 #casalog.post("outfreqchans=",outfreqchans) 

543 expandchanmask(tmp_regridim,inpfreqchans,tmp_outmaskimage,outfreqchans) 

544 #shutil.copytree(tmp_regridim,'my_tmp_regrid')  

545 #shutil.copytree(tmp_outmaskimage,'my_tmp_outmaskimage')  

546 

547# usechanims={} # list of input mask to be use for each outpfreq 

548# for i in outfreqchans: 

549# nearestch = findnearest(inpfreqchans,i) 

550# usechanims[i]=nearestch  

551# # put masks from inp image channel by channel 

552# for j in outfreqs: 

553# pix = refchanchunk[usechanims[j]-refchanst] 

554# #_ia.putchunk(pixels=pix,blc=[inshp[0]-1,inshp[1]-1,0,j]) 

555# _ia.putchunk(pixels=pix.transpose(),blc=[0,0,0,j]) 

556 needtoregrid=False 

557 

558 if _ia.isopen(): _ia.close() 

559 _ia.open(tmp_outmaskimage) 

560 #  

561 

562 if needtoregrid: 

563 # closing current output image 

564 if _ia.isopen(): 

565 _ia.close() 

566 _ia.open(tmp_maskimage) 

567 #os.system('cp -r %s beforeregrid.im' % tmp_maskimage) 

568 if os.path.isdir(tmp_outmaskimage): 

569 #casalog.post("Removing %s" % tmp_outmaskimage) 

570 shutil.rmtree(tmp_outmaskimage) 

571 #regridmask(tmp_maskimage,outparentim,tmp_outmaskimage) 

572 regridmask(tmp_maskimage,inpimage,tmp_outmaskimage,method=regridmethod) 

573 _ia.remove() 

574 #casalog.post("closing after regrid") 

575 _ia.open(tmp_outmaskimage) # reopen output tmp image 

576 

577 # for debugging 

578 #os.system('cp -r '+outparentim+" expandmode-copy-"+outparentim) 

579 #os.system('cp -r '+tmp_outmaskimage+" expandmode-copy-"+tmp_outmaskimage) 

580 if outbmask!='': 

581 #convert the image mask to T/F mask 

582 casalog.post("Convert the image mask to T/F mask",'INFO') 

583 # regions will be masked if == 0.0 for a new outfile, if outfile exists  

584 # the pixel values inside specified mask is preserved and the rest is masked 

585 if os.path.isdir(outparentim): 

586 _ia.calcmask(mask='%s==1.0' % tmp_outmaskimage,name=outbmask,asdefault=True) 

587 else: 

588 _ia.calcmask(mask='%s!=0.0' % tmp_outmaskimage,name=outbmask,asdefault=True) 

589 if storeinmask: 

590 isNewFile=False 

591 if not os.path.isdir(outparentim): 

592 makeEmptyimage(inpimage,outparentim) 

593 isNewFile=True 

594 _ia.open(outparentim) 

595 if isNewFile: 

596 _ia.set(1) 

597 # if output image exist its image pixel values will not be normalized the region 

598 # outside input mask will be masked. 

599 #check  

600 curinmasks = _ia.maskhandler('get') 

601 if outbmask in curinmasks: 

602 if not overwrite: 

603 raise RuntimeError("Internal mask,"+outbmask+" exists. Please set overwrite=True.") 

604 else: 

605 _ia.maskhandler('delete',outbmask) 

606 

607 _ia.maskhandler('copy',[tmp_outmaskimage+':'+outbmask, outbmask]) 

608 _ia.maskhandler('set',outbmask) 

609 _ia.done() 

610 casalog.post("Output the mask to %s in %s" % (outbmask,outparentim) ,"INFO") 

611 else: 

612 ow = False 

613 if inpimage==output: 

614 casalog.post("Updating "+output+" with new mask","INFO") 

615 ow=True 

616 else: 

617 if os.path.isdir(outparentim): 

618 casalog.post(outparentim+" exists, overwriting","INFO") 

619 ow=True 

620 else: 

621 casalog.post("Output the mask to "+outparentim ,"INFO") 

622 _ia.rename(outparentim,ow) 

623 _ia.done() 

624 

625 except Exception as instance: 

626 casalog.post("*** Error (1) *** %s" % instance, 'ERROR') 

627 if _ia.isopen(): 

628 _ia.close() 

629 _ia.done() 

630 raise 

631 finally: 

632 if inmaskcsys: 

633 inmaskcsys.done() 

634 if incsys: 

635 incsys.done() 

636 if ocsys: 

637 ocsys.done() 

638 if os.path.exists(tmp_maskimage): 

639 shutil.rmtree(tmp_maskimage) 

640 if os.path.exists(tmp_regridim): 

641 shutil.rmtree(tmp_regridim) 

642 if os.path.exists(tmp_outmaskimage): 

643 shutil.rmtree(tmp_outmaskimage) 

644 

645 # === main process for copy mode: also does merge of masks === 

646 # copy is a just special case of merge mode 

647 # CHANGE: 

648 # all input masks should be specified in inpmask  

649 # type of inpmask accepted: 1/0 mask, T/F mask, region file, and region expression in a string  

650 # already stored internally in seperate lists 

651 # rgfiles - region files (binary or CRTF-format text file) 

652 # imgfiles - 1/0 image masks 

653 # rglist - region expression in strings 

654 # bmasks - T/F internal masks 

655 # 

656 # avaialble parameters: inpimage (string) , inpmask (list/string), output(string) 

657 # input inpimage as a template or image used for defining regions when it is given in inpmask  

658 # inpmask as list of all the masks to be merged (image masks, T/F internal masks, regions) 

659 

660 if mode=='copy': 

661 sum_tmp_outfile='__tmp_outputmask_'+pid 

662 tmp_inmask='__tmp_frominmask_'+pid 

663 tmp_allrgmaskim='__tmp_fromAllRgn_'+pid 

664 tmp_rgmaskim='__tmp_fromRgn_'+pid 

665 # making sure to remove pre-existing temp files 

666 cleanuptempfiles([sum_tmp_outfile, tmp_inmask, tmp_allrgmaskim, tmp_rgmaskim]) 

667 usedimfiles=[] 

668 usedbmasks=[] 

669 usedrgfiles=[] 

670 usedrglist=[] 

671 # This coordsys object used in several places below will need to be closed, if created 

672 tcsys = None 

673 try: 

674 # check outparentim - image part of output and set as a template image 

675 if not (os.path.isdir(outparentim) or (outparentim==inpimage)): 

676 # Output is either a new image or the same as inpimage 

677 

678 # figure out which input mask to be used as template 

679 # if inpimage is defined use the first one else try the first one 

680 # inpmask 

681 #if output=='': 

682 # if type(inpimage)==list: 

683 # raise Exception, "inputimage must be a string" 

684 # elif type(inpimage)==str: 

685 # outimage=inpimage 

686 # casalog.post("No outimage is specified. Will overwrite input image: "+outimage,'INFO') 

687 

688 #if type(inpimage)==list and len(inpimage)!=0: 

689 # tmp_template=inpimage[0] 

690 #elif inpimage!='' and inpimage!=[]: 

691 # tmp_template=inpimage # string 

692 #tmp_template=inpimage # string 

693 #else: 

694 # if type(inpmask)==list and len(inpmask)!=0: 

695 # fsep=inpmask[0].rfind(':') 

696 # if fsep ==-1: 

697 # raise IOError, "Cannot resolve inpmask name, check the format" 

698 # else: 

699 # tmp_template=inpmask[0][:inpmask[0].rfind(':')] 

700 # elif inpmask!='' and inpmask!=[]: 

701 # # this is equivalent to 'copy' the inpmask 

702 # tmp_template=inpmask #string 

703 

704 # create an empty image with the coords from inpimage 

705 makeEmptyimage(inpimage,sum_tmp_outfile) 

706 #casalog.post("making an empty image from inpimage to sum_tmp_outfile") 

707 else: 

708 #use output image(either the same as the input image or other existing image)  

709 # - does not do zeroeing out, so output image is only modified *for outbmask!=''*  

710 shutil.copytree(outparentim,sum_tmp_outfile) 

711 # temporary clear out the internal masks from the working image 

712 if _ia.isopen(): _ia.close() 

713 _ia.open(sum_tmp_outfile) 

714 if (len(imgfiles) or len(rglist) or len(rgfiles)): 

715 # do zeroeing out working base image for masks  

716 _ia.set(0) 

717 origmasks = _ia.maskhandler('get') 

718 _ia.maskhandler('delete',origmasks) 

719 _ia.close() 

720 

721 if len(imgfiles)>0: 

722 # summing all the images 

723 casalog.post('Summing all mask images in inpmask and normalized to 1 for mask','INFO') 

724 for img in imgfiles: 

725 #tmpregrid='__tmp_regrid.'+img 

726 dirname = os.path.dirname(img) 

727 basename = os.path.basename(img) 

728 tmpregrid= dirname+'/'+'__tmp_regrid.'+basename if len(dirname) else '__tmp_regrid.'+basename 

729 tmpregrid+='_'+pid 

730 if os.path.exists(tmpregrid): 

731 shutil.rmtree(tmpregrid) 

732 # regrid to output image coords 

733 try: 

734 regridmask(img,sum_tmp_outfile,tmpregrid) 

735 addimagemask(sum_tmp_outfile,tmpregrid) 

736 usedimfiles.append(img) 

737 finally: 

738 shutil.rmtree(tmpregrid) 

739 # get boolean masks 

740 # will work in image(1/0) masks 

741 

742 if debug: 

743 casalog.post(("imgfiles=" + imgfiles)) 

744 shutil.copytree(sum_tmp_outfile,sum_tmp_outfile+"_imagemaskCombined") 

745 

746 if len(bmasks)>0: 

747 casalog.post('Summing all T/F mask in inpmask and normalized to 1 for mask','INFO') 

748 for msk in bmasks: 

749 (imname,mskname) = extractmaskname(msk) 

750 #if msk.find(':')<0: 

751 # # assume default mask 

752 # msk=msk+':mask0' 

753 #imname=msk[:msk.rfind(':')] 

754 if _ia.isopen(): _ia.close() 

755 _ia.open(imname) 

756 inmasks=_ia.maskhandler('get') 

757 _ia.close() 

758 if not inmasks.count(mskname): 

759 raise TypeError(mskname+" does not exist in "+imname+" -available masks:"+str(inmasks)) 

760 # move T/F mask to image mask 

761 # changed to usemasked=False as of CAS-5443  

762 

763 pixelmask2cleanmask(imname, mskname, tmp_inmask, False) 

764 cleanuptempfiles(['__tmp_fromTFmask_'+pid]) 

765 regridmask(tmp_inmask,sum_tmp_outfile,'__tmp_fromTFmask_'+pid) 

766 # for T/F mask do AND operation 

767 _ia.open(sum_tmp_outfile) 

768 if _ia.statistics()['sum'][0] != 0: 

769 multiplyimagemask(sum_tmp_outfile,'__tmp_fromTFmask_'+pid) 

770 else: 

771 addimagemask(sum_tmp_outfile,'__tmp_fromTFmask_'+pid) 

772 _ia.close() 

773 usedbmasks.append(msk) 

774 # need this temp file for the process later 

775 ###shutil.rmtree('__tmp_fromTFmask')  

776 if os.path.isdir(tmp_inmask): 

777 shutil.rmtree(tmp_inmask) 

778 # if overwriting to inpimage and if not writing to in-mask, delete the boolean mask 

779 if outparentim==inpimage and inpimage==imname: 

780 if outbmask=="": 

781 _ia.open(imname) 

782 _ia.maskhandler('delete',[mskname]) 

783 _ia.close() 

784 _ia.open(imname) 

785 _ia.close() 

786 

787 if debug: casalog.post("check rgfiles and rglist") 

788 

789 if len(rgfiles)>0 or len(rglist)>0: 

790 # create an empty image with input image coords. 

791 #casalog.post("Using %s as a template for regions" % inpimage ) 

792 if _ia.isopen(): _ia.close() 

793 _ia.open(inpimage) 

794 tshp=_ia.shape() 

795 tcsys=_ia.coordsys() 

796 _ia.fromshape(tmp_allrgmaskim,shape=tshp, csys=tcsys.torecord(),overwrite=True) 

797 _ia.done() 

798 if os.path.isdir(tmp_rgmaskim): 

799 shutil.rmtree(tmp_rgmaskim) 

800 shutil.copytree(tmp_allrgmaskim,tmp_rgmaskim) 

801 

802 if len(rgfiles)>0: 

803 # Here only CRTF format file is expected  

804 nrgn=0 

805 for rgn in rgfiles: 

806 subnrgn=0 

807 with open(rgn) as f: 

808 # check if it has '#CRTF' in the first line 

809 firstline=f.readline() 

810 if firstline.count('#CRTF')==0: 

811 raise Exception("Input text file does not seems to be in a correct format \ 

812 (must contains #CRTF in the first line)") 

813 else: 

814 try: 

815 # reset temp mask image 

816 if _ia.isopen(): _ia.close() 

817 _ia.open(tmp_rgmaskim) 

818 _ia.set(pixels=0.0) 

819 _ia.close() 

820 #casalog.post... "tshp=",tshp 

821 #casalog.post... "tcsys.torecord=",tcsys.torecord() 

822 inrgn=_rg.fromtextfile(rgn, tshp, tcsys.torecord()) 

823 #casalog.post... "inrgn=",inrgn 

824 _im.regiontoimagemask(tmp_rgmaskim,region=inrgn) 

825 addimagemask(tmp_allrgmaskim,tmp_rgmaskim) 

826 #shutil.rmtree(tmp_rgmaskim) 

827 subnrgn +=1 

828 nrgn +=1 

829 except: 

830 break 

831 if subnrgn>0: 

832 usedrgfiles.append(rgn) 

833 casalog.post("Converted %s regions from %s region files to image mask" % (nrgn,len(rgfiles)),"INFO") 

834 

835 if debug: casalog.post("processing rglist...") 

836 if len(rglist)>0: 

837 #casalog.post( "Has rglist...") 

838 nrgn=0 

839 for rgn in rglist: 

840 # reset temp mask image 

841 if _ia.isopen(): _ia.close() 

842 _ia.open(tmp_rgmaskim) 

843 _ia.set(pixels=0.0) 

844 _ia.close() 

845 inrgn=_rg.fromtext(rgn, tshp, tcsys.torecord()) 

846 _im.regiontoimagemask(tmp_rgmaskim,region=inrgn)