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
« 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 (")
78# Default: none
79# Examples:
80# Make an image that is image1.im - image2.im
81# expr=' ("image1.im" - "image2.im" )'
82# Clip an image below a value (0.5 in this case)
83# expr = ' iif("image1.im">=0.5, "image1.im", 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("image1.im" * "image1.im" + "image2.im" * "image2.im") '
89# Note: No exponentiaion available?
90# Build an image pixel by pixel from the minimum of (image2.im, 2*image1.im)
91# expr='min("image2.im",2*max("image1.im"))'
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
165import os
166import re
167import shutil
168import sys
170from casatools import image, imagepol, regionmanager, coordsys, quanta
171from casatasks import casalog
172from .ialib import write_image_history
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)
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()
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
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()
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]
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()
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
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)
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)
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 )
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] + ')'
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
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
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
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
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 )
552def _immath_parse( expr='' ):
553 retValue=[]
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="'"
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;
572 retValue.append( current[start:end] )
573 current=current[end+1:]
575 return retValue
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())
587 tmpvars.sort()
588 tmpvars.reverse()
589 for varname in tmpvars:
590 expr = expr.replace(varname, '"' + tmpfiles[varname] + '"')
591 return(expr)
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]