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

288 statements  

« prev     ^ index     » next       coverage.py v7.6.4, created at 2024-11-01 07:19 +0000

1########################################################################3 

2# task_immath.py 

3# 

4# Copyright (C) 2008, 2009 

5# Associated Universities, Inc. Washington DC, USA. 

6# 

7# This script is free software; you can redistribute it and/or modify it 

8# under the terms of the GNU Library General Public License as published by 

9# the Free Software Foundation; either version 2 of the License, or (at your 

10# option) any later version. 

11# 

12# This library is distributed in the hope that it will be useful, but WITHOUT 

13# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 

14# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public 

15# License for more details. 

16# 

17# You should have received a copy of the GNU Library General Public License 

18# along with this library; if not, write to the Free Software Foundation, 

19# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. 

20# 

21# Correspondence concerning AIPS++ should be adressed as follows: 

22# Internet email: casa-feedback@nrao.edu. 

23# Postal address: AIPS++ Project Office 

24# National Radio Astronomy Observatory 

25# 520 Edgemont Road 

26# Charlottesville, VA 22903-2475 USA 

27# 

28# <summary> 

29# CASA task for smoothing an image, by doing Forier-based convolution 

30# on a CASA image file. 

31# </summary> 

32# 

33# <reviewed reviwer="" date="" tests="" demos=""> 

34# </reviewed 

35# 

36# <author> 

37# Shannon Jaeger, University of Calgary (image math)  

38# Takeshi Nakazato, National Radio Astronomy Obaservatory (polarization) 

39# </author> 

40# 

41# <prerequisite> 

42# </prerequisite> 

43# 

44# <etymology> 

45# immath stands for image mathematics 

46# </etymology> 

47# 

48# <synopsis> 

49# This task evaluates mathematical expressions involving existing 

50# image files. The results of the calculations are stored in the  

51# designated output file. Options are available to specify mathematical  

52# expression directly or pre-defined expression for calculation of  

53# spectral index image, and polarization intensity and position angle  

54# images are available. The image file names imbedded in the expression or 

55# specified in the imagename parameter for the pre-defined calculations may 

56# be CASA images or FITS images. 

57# 

58# 

59# NOTE: Index values for axes start at 0 for the box and chans 

60# parameters, but at 1 when used with the indexin function 

61# in expression. Use the imhead task to see the range of 

62# values for each axes. 

63#  

64# 

65# Keyword arguments: 

66# outfile -- The file where the results of the image calculations  

67# are stored. Overwriting an existing outfile is not permitted. 

68# Default: none; Example: outfile='results.im' 

69# mode -- mode for mathematical operation 

70# Default: evalexpr 

71# Options: 'evalexpr' : evalulate a mathematical expression defined in 'expr'  

72# 'spix' : spectalindex image  

73# 'pola' : polarization position angle image  

74# 'poli' : polarization intesity image  

75# mode expandable parameters 

76# expr -- (for mode='evalexpr') A mathematical expression, with image file names. 

77# Image file names MUST be enclosed in double quotes (&quot;) 

78# Default: none  

79# Examples: 

80# Make an image that is image1.im - image2.im 

81# expr=' (&quot;image1.im&quot; - &quot;image2.im&quot; )' 

82# Clip an image below a value (0.5 in this case) 

83# expr = ' iif(&quot;image1.im&quot;>=0.5, &quot;image1.im&quot;, 0.0) ' 

84# Note: iif (a, b, c) a is the boolian expression 

85# b is the value if true 

86# c is the value if false 

87# Take the rms value of two images 

88# expr = ' sqrt(&quot;image1.im&quot; * &quot;image1.im&quot; + &quot;image2.im&quot; * &quot;image2.im&quot;) ' 

89# Note: No exponentiaion available? 

90# Build an image pixel by pixel from the minimum of (image2.im, 2*image1.im) 

91# expr='min(&quot;image2.im&quot;,2*max(&quot;image1.im&quot;))' 

92# imagename -- (for mode='spix','pola','poli') input image names  

93# Default: none; 

94# Examples: mode='spix'; imagename=['image1.im','image2.im'] will calculate  

95# an image of log(S1/S2)/log(f1/f2), where S1 and S2 are fluxes and  

96# f1 and f2 are frequencies 

97# mode='pola'; imagename=['imageQ.im','imageU.im'] will calculate  

98# an image of polarization angle distribution, where imageQ.im and  

99# imageU.im are Stokes Q and U images, respectively. Calculate 0.5*arctan(U/Q). 

100# mode='poli'; imagename=['imageQ.im','imageU.im','imageV.im'] will calculate 

101# total polarization intensity image, where imageQ.im, imageU.im, imageV.im 

102# are Stokes Q, U, and V images, respectively. 

103# sigma - (for mode='poli') standard deviation of noise of Stokes images with unit such as 

104# Jy/beam to correct for bias  

105# Default: '0.0Jy/beam' (= no debiasing) 

106# mask -- Name of mask applied to each image in the calculation 

107# Default '' means no mask; Example: mask='orion.mask'.  

108# region -- File path to an ImageRegion file. 

109# An ImageRegion file can be created with the CASA 

110# viewer's region manager. Typically ImageRegion files 

111# will have the suffix '.rgn'. If a region file is given 

112# then the box, chans, and stokes selections whill be  

113# ignored. 

114# Default: none 

115# Example: region='myimage.im.rgn' 

116# box -- A box region on the directional plane 

117# Only pixel values acceptable at this time. 

118# Default: none (whole 2-D plane); Example: box='10,10,50,50' 

119# chans -- channel numbers, velocity, and/or frequency 

120# Only channel numbers acceptable at this time. 

121# Default: none (all); Example: chans='3~20'  

122# stokes -- Stokes parameters to image, may or may not be separated 

123# by commas but best if you use commas. 

124# Default: none (all); Example: stokes='IQUV'; 

125# Options: 'I','Q','U','V','RR','RL','LR','LL','XX','YX','XY','YY', ...  

126# 

127# Available functions in the <i>expr</i> and <i>mask</i> paramters: 

128# pi(), e(), sin(), sinh(), asinh(), cos(), cosh(), tan(), tanh(), 

129# atan(), exp(), log(), log10(), pow(), sqrt(), complex(), conj() 

130# real(), imag(), abs(), arg(), phase(), aplitude(), min(), max() 

131# round(), isgn(), floor(), ceil(), rebin(), spectralindex(), pa(),  

132# iif(), indexin(), replace(), ... 

133# 

134# For a full description of the allowed syntax see the  

135# Lattice Expression Language (LEL) documentation on the at: 

136# http://aips2.nrao.edu/docs/notes/223/223.html 

137# 

138# NOTE: where indexing and axis numbering are used in the above 

139# functions they are 1-based, ie. numbering starts at 1. 

140# 

141# </synopsis>  

142# 

143# <example> 

144# <srcblock> 

145# </srcblock 

146# 

147# </example> 

148# 

149# <motivation> 

150# To provide a user-friendly task interface to imagecalc and ??? 

151# as well as an more user-friendling math syntax then what is 

152# provided by the CASA Lattice Exprssion Language. 

153# </motivation> 

154# 

155# <todo> 

156# Crystal wanted different masks for different inputs 

157# but unlikely that its really needed. 

158# 

159# Add an "overwrite" output file parameter 

160# 

161# Add polygon and circle region selection  

162# </todo> 

163########################################################################3 

164 

165import os 

166import re 

167import shutil 

168import sys 

169 

170from casatools import image, imagepol, regionmanager, coordsys, quanta 

171from casatasks import casalog 

172from .ialib import write_image_history 

173 

174def immath( 

175 imagename, mode, outfile, expr, varnames, sigma, 

176 polithresh, mask, region, box, chans, stokes, stretch, 

177 imagemd, prec 

178): 

179 # Tell CASA who will be reporting 

180 casalog.origin('immath') 

181 tmpFilePrefix='_immath_tmp' + str(os.getpid()) + '_' 

182 try: 

183 myia = image() 

184 myia.dohistory(False) 

185 outia = None 

186 _immath_initial_cleanup(tmpFilePrefix) 

187 outfile = _immath_check_outfile(outfile) 

188 # Find the list of filenames in the expression 

189 # also do a quick check to see if all of the files 

190 # exist 

191 tmpfilenames = '' 

192 filenames = imagename 

193 if mode=='evalexpr': 

194 tmpfilenames = _immath_parse(expr) 

195 if isinstance(filenames, str): 

196 filenames = [filenames] 

197 varnames = _immath_varnames(varnames, filenames, tmpfilenames) 

198 filenames = _immath_filenames(filenames, tmpfilenames, varnames, mode) 

199 expr = expr.replace(' ', '') 

200 if mode == 'spix': 

201 expr = _immath_dospix(len(filenames), varnames) 

202 if mode == 'pola': 

203 _immath_new_pola( 

204 filenames, outfile, tmpFilePrefix, mask, region, 

205 box, chans, stokes, stretch, polithresh, myia 

206 ) 

207 return True 

208 elif mode == 'poli' or mode == 'lpoli' or mode == 'tpoli': 

209 _immath_new_poli( 

210 filenames, outfile, tmpFilePrefix, mask, region, 

211 box, chans, stokes, stretch, sigma, myia, mode 

212 ) 

213 return True 

214 elif mode == 'evalexpr' or mode == 'spix': 

215 if box or chans or stokes or region or mask: 

216 (subImages, file_map) = _immath_createsubimages( 

217 box, chans, stokes, region, mask, 

218 stretch, filenames, myia, tmpFilePrefix 

219 ) 

220 if imagemd: 

221 casalog.post( 

222 "Specifying region, box, chan, or stokes will " 

223 + "create smaller sub-images. The image " 

224 + "metadata specified in imagemd will have to " 

225 + "conform to the output, not the input image " 

226 + "dimensions. Please check your output image " 

227 + "for accurate header definition.", 'WARN' 

228 ) 

229 (expr, varnames, subImages) = _immath_updateexpr( 

230 expr, varnames, subImages, filenames, file_map 

231 ) 

232 outia = _immath_compute( 

233 imagename, expr, outfile, imagemd, myia, prec 

234 ) 

235 else: 

236 # If the user didn't give any region or mask information 

237 # then just evaluated the expression with the filenames in it. 

238 outia = _immath_dofull( 

239 imagename, imagemd, outfile, mode, expr, 

240 varnames, filenames, myia, prec 

241 ) 

242 else: 

243 raise(Exception, "Unsupported mode " + str(mode)) 

244 try: 

245 vars = locals() 

246 param_names = immath.__code__.co_varnames[:immath.__code__.co_argcount] 

247 param_vals = [vars[p] for p in param_names] 

248 write_image_history( 

249 outia, sys._getframe().f_code.co_name, 

250 param_names, param_vals, casalog 

251 ) 

252 except Exception as instance: 

253 casalog.post("*** Error \'%s\' updating HISTORY" % (instance), 'WARN') 

254 return True 

255 except Exception as error: 

256 if mode == 'evalexpr': 

257 casalog.post("Unable to process expression " + expr, 'SEVERE') 

258 else: 

259 casalog.post("Error running immath", 'SEVERE') 

260 casalog.post("Exception caught was: " + str(error), 'SEVERE') 

261 raise 

262 finally: 

263 if myia: 

264 myia.done() 

265 if outia: 

266 outia.done() 

267 _immath_cleanup(tmpFilePrefix) 

268 

269def _immath_concat_stokes(filenames, target, myia): 

270 myia.open(filenames[0]) 

271 stokes_axis = myia.coordsys().findaxisbyname("stokes") 

272 myia.done() 

273 casalog.post("Concatenating images along stokes axis") 

274 myia = myia.imageconcat( 

275 outfile=target, infiles=filenames, axis=stokes_axis 

276 ) 

277 myia.done() 

278 

279def _immath_getregion(region, box, chans, stokes, mode, myia, target): 

280 myreg = region 

281 if (type(region) != type({})): 

282 if stokes: 

283 casalog.post( 

284 "Ignoring stokes parameters selection for mode='" 

285 + mode + "'." 

286 ,'WARN' 

287 ) 

288 stokes='' 

289 myia.open(target) 

290 myrg = regionmanager() 

291 myreg = myrg.frombcs( 

292 csys=myia.coordsys().torecord(), shape=myia.shape(), box=box, 

293 chans=chans, stokes=stokes, stokescontrol="a", region=region 

294 ) 

295 myia.done() 

296 myrg.done() 

297 return myreg 

298 

299def _immath_new_pola( 

300 filenames, outfile, tmpFilePrefix, mask, region, 

301 box, chans, stokes, stretch, polithresh, myia 

302): 

303 target = filenames[0] 

304 if len(filenames) > 1: 

305 target = tmpFilePrefix + "_concat_for_pola" 

306 _immath_concat_stokes(filenames, target, myia) 

307 myreg = _immath_getregion(region, box, chans, stokes, "pola", myia, target) 

308 mypo = imagepol() 

309 myqa = quanta() 

310 if (polithresh): 

311 if (mask != ""): 

312 mask = "" 

313 casalog.post("Ignoring mask parameter in favor of polithresh parameter", 'WARN') 

314 if (myqa.getunit(polithresh) != ""): 

315 initUnit = myqa.getunit(polithresh) 

316 myia.dohistory(False) 

317 myia.open(filenames[0]) 

318 bunit = myia.brightnessunit() 

319 polithresh = myqa.convert(polithresh, bunit) 

320 myia.done() 

321 if (myqa.getunit(polithresh) != bunit): 

322 raise Exception("Units of polithresh " + initUnit \ 

323 + " do not conform to input image units of " + bunit \ 

324 + " so cannot perform thresholding. Please correct units and try again.") 

325 polithresh = myqa.getvalue(polithresh)[0] 

326 lpol = tmpFilePrefix + "_lpol" 

327 mypo.open(target) 

328 myia = mypo.linpolint(debias=False, outfile=lpol, region=myreg) 

329 myia.done() 

330 mypo.done() 

331 mypo.open(target) 

332 myia = mypo.linpolposang( 

333 outfile=outfile, region=myreg, mask=mask, stretch=stretch 

334 ) 

335 mypo.done() 

336 if (polithresh): 

337 myexpr = "'" + lpol + "' >= " + str(polithresh) 

338 myia.dohistory(False) 

339 myia.calcmask(name='mask0', mask=myexpr) 

340 casalog.post( 

341 'Calculated mask based on linear polarization threshold ' 

342 + str(polithresh), 

343 'INFO' 

344 ) 

345 myia.done() 

346 

347def _immath_new_poli( 

348 filenames, outfile, tmpFilePrefix, mask, region, 

349 box, chans, stokes, stretch, sigma, myia, mode 

350): 

351 target = filenames[0] 

352 if len(filenames) > 1: 

353 target = tmpFilePrefix + "_concat_for_poli" 

354 _immath_concat_stokes(filenames, target, myia) 

355 if mode == 'tpoli': 

356 myia.open(target) 

357 csys = myia.coordsys() 

358 myia.done() 

359 stokes = csys.stokes() 

360 csys.done() 

361 for p in ["Q", "U", "V"]: 

362 if stokes.count(p) == 0: 

363 raise Exception( 

364 "Stokes " + p + " is required for mode tpoli but was not found" 

365 ) 

366 debias = False 

367 newsigma = 0 

368 myqa = quanta() 

369 if sigma: 

370 qsigma = myqa.quantity(sigma) 

371 if myqa.getvalue(qsigma)[0] > 0: 

372 debias = True 

373 sigmaunit = myqa.getunit(qsigma) 

374 try: 

375 myia.open(filenames[0]) 

376 iunit = myia.brightnessunit() 

377 myia.done() 

378 except: 

379 raise Exception('Unable to get brightness unit from image file ' + filenames[0]) 

380 if sigmaunit != iunit: 

381 newsigma = myqa.convert(qsigma,iunit) 

382 else: 

383 newsigma = sigma 

384 myreg = _immath_getregion(region, box, chans, stokes, "poli", myia, target) 

385 mypo = imagepol() 

386 mypo.open(target) 

387 # for some annoying reason, qa.getvalue() returns an array in this context 

388 numeric_sigma = myqa.getvalue(myqa.quantity(newsigma))[0] 

389 

390 if mode == 'tpoli' or mode == 'poli': 

391 myia = mypo.totpolint( 

392 debias=debias, sigma=numeric_sigma, outfile=outfile, 

393 region=myreg, mask=mask, stretch=stretch 

394 ) 

395 elif mode == 'lpoli': 

396 myia = mypo.linpolint( 

397 debias=debias, sigma=numeric_sigma, outfile=outfile, 

398 region=myreg, mask=mask, stretch=stretch 

399 ) 

400 else: 

401 raise Exception("Logic Error: Unhandled mode " + mode) 

402 myia.done() 

403 mypo.done() 

404 

405def _immath_compute( 

406 imagename, expr, outfile, imagemd, myia, prec 

407): 

408 # Do the calculation 

409 res = myia.imagecalc( 

410 pixels=expr, outfile=outfile, 

411 imagemd=_immath_translate_imagemd(imagename, imagemd), prec=prec 

412 ) 

413 res.dohistory(False) 

414 # modify stokes type for polarization intensity image 

415 return res 

416 

417def _immath_updateexpr(expr, varnames, subImages, filenames, file_map): 

418 # Make sure no problems happened 

419 if len(filenames) != len(subImages): 

420 raise Exception( 

421 'Unable to create subimages for all image names given' 

422 ) 

423 # because real file names also have to be mapped to a corresponding subimage, CAS-1830 

424 for k in file_map.keys(): 

425 # we require actual image names to be in quotes when used in the expression 

426 varnames.extend(["'" + k + "'", '"' + k + '"']) 

427 subImages.extend(2 * [file_map[k]]) 

428 # Put the subimage names into the expression 

429 try: 

430 expr = _immath_expr_from_varnames(expr, varnames, subImages) 

431 except Exception as e: 

432 casalog.post( 

433 "Unable to construct pixel expression aborting immath: " + str(e), 

434 'SEVERE' 

435 ) 

436 raise 

437 return (expr, varnames, subImages) 

438 

439def _immath_createsubimages( 

440 box, chans, stokes, region, mask, 

441 stretch, filenames, myia, tmpFilePrefix 

442): 

443 subImages = [] 

444 file_map = {} 

445 i = 0 

446 for image in filenames: 

447 try: 

448 myia.open(image) 

449 myrg = regionmanager() 

450 reg = myrg.frombcs(csys=myia.coordsys().torecord(), 

451 shape=myia.shape(), box=box, chans=chans, stokes=stokes, 

452 stokescontrol="a", region=region 

453 ) 

454 tmpFile = tmpFilePrefix + str(i) 

455 subim = myia.subimage( 

456 region=reg, mask=mask, outfile=tmpFile, stretch=stretch 

457 ) 

458 subim.done() 

459 file_map[image] = tmpFile 

460 subImages.append(tmpFile) 

461 myia.done() 

462 i = i + 1 

463 except Exception as e: 

464 raise Exception( 

465 'Unable to apply region to file: ' + image 

466 ) 

467 finally: 

468 myia.done() 

469 return (subImages, file_map) 

470 

471def _immath_dofull( 

472 imagename, imagemd, outfile, mode, expr, 

473 varnames, filenames, myia, prec 

474): 

475 expr = _immath_expr_from_varnames(expr, varnames, filenames) 

476 return _immath_compute( 

477 imagename, expr, outfile, imagemd, myia, prec 

478 ) 

479 

480def _immath_dospix(nfiles, varnames): 

481 if nfiles != 2: 

482 raise Exception('Requires two images at different frequencies') 

483 return 'spectralindex(' + varnames[0] + ', ' + varnames[1] + ')' 

484 

485def _immath_filenames(filenames, tmpfilenames, varnames, mode): 

486 ignoreimagename = False 

487 if mode=='evalexpr': 

488 varnamesSet = set(varnames) 

489 count = 0 

490 for imname in tmpfilenames: 

491 # check if it is one of varnames, if not check the files in expr exist  

492 if(not varnamesSet.issuperset(imname)): 

493 if( not os.path.exists(imname)): 

494 raise Exception('Image data set not found - please verify ' + imname) 

495 else: 

496 count = count + 1 

497 if len(tmpfilenames) == count: 

498 ignoreimagename = True 

499 filenames = tmpfilenames 

500 if not ignoreimagename: 

501 for i in range(len(filenames)): 

502 if not os.path.exists(filenames[i]): 

503 casalog.post("Image data set not found - please verify " +filenames[i], "SEVERE") 

504 raise Exception('Image data set not found - please verify '+filenames[i]) 

505 return filenames 

506 

507def _immath_varnames(varnames, filenames, tmpfilenames): 

508 # Construct the variable name list. We append to the list the 

509 # default variable names if the user hasn't supplied a full suite. 

510 if not isinstance(varnames, list): 

511 name0 = varnames 

512 varnames = [] 

513 if name0: 

514 varnames.append(name0) 

515 nfile = max(len(filenames),len(tmpfilenames)) 

516 for i in range(len(varnames), nfile): 

517 varnames.append('IM'+str(i)) 

518 casalog.post( 'Variable name list is: '+str(varnames), 'DEBUG1' ) 

519 return varnames 

520 

521def _immath_initial_cleanup(tmpFilePrefix): 

522 try: 

523 _immath_cleanup(tmpFilePrefix) 

524 except Exception as e: 

525 casalog.post( 'Unable to cleanup working directory '+os.getcwd()+'\n'+str(e), 'SEVERE' ) 

526 raise 

527 

528def _immath_check_outfile(outfile): 

529 if not outfile: 

530 outfile = 'immath_results.im' 

531 casalog.post( 

532 "The outfile parameter is empty, consequently the " 

533 + "resultant image will be saved on disk in an image named " 

534 + outfile, 'WARN' 

535 ) 

536 if (os.path.exists(outfile)): 

537 raise Exception( 

538 'Output file '+ outfile 

539 + ' exists. immath can not proceed, please ' 

540 + 'remove it or change the output file name.' 

541 ) 

542 return outfile 

543 

544def _immath_cleanup(filePrefix): 

545 # Remove any old tmp files that may be left from 

546 # a previous run of immath 

547 fileList=os.listdir( os.getcwd() ) 

548 for file in fileList: 

549 if ( file.startswith( filePrefix ) ): 

550 shutil.rmtree( file ) 

551 

552def _immath_parse( expr='' ): 

553 retValue=[] 

554 

555 # Find out if the names are surrounded by single or double quotes 

556 quote='' 

557 if ( expr.find('"') > -1 ): 

558 quote='"' 

559 if ( expr.find("'") > -1 ): 

560 quote="'" 

561 

562 current=expr; 

563 while( current.count(quote) > 1 ): 

564 start = current.index(quote)+1 

565 end = current[start:].index(quote)+start 

566 if ( retValue.count( current[start:end] ) > 0 ) : 

567 # We already have this file name so we won't add it 

568 # to the list again. This saves us work and disk space. 

569 current=current[end+1:] 

570 continue; 

571 

572 retValue.append( current[start:end] ) 

573 current=current[end+1:] 

574 

575 return retValue 

576 

577# it is important to sort the varnames in reverse order before doing 

578# the substitutions to assure the substitution set is performed correctly 

579# CAS-1678 

580def _immath_expr_from_varnames(expr, varnames, filenames): 

581 tmpfiles = {} 

582 for i in range(len(filenames)): 

583 tmpfiles[varnames[i]] = filenames[i] 

584 # python 3 requires explicit list conversion 

585 tmpvars = list(tmpfiles.keys()) 

586 

587 tmpvars.sort() 

588 tmpvars.reverse() 

589 for varname in tmpvars: 

590 expr = expr.replace(varname, '"' + tmpfiles[varname] + '"') 

591 return(expr) 

592 

593def _immath_translate_imagemd(imagename, imagemd): 

594 # may IM0 etc is a real file name 

595 if os.path.exists(imagemd): 

596 return imagemd 

597 # match IM0, IM1, ... etc 

598 m = re.match("^IM(\d+)$", imagemd) 

599 if not m: 

600 return imagemd 

601 idx = int(m.groups()[0]) 

602 if idx == 0 and type(imagename) == str: 

603 return imagename 

604 # if here, then imagename is an array of strings 

605 if idx >= len(imagename): 

606 # out of range 

607 return imagemd 

608 return imagename[idx] 

609