imfit_cli.py

Go to the documentation of this file.

#
# This file was generated using xslt from its XML file
#
# Copyright 2008, Associated Universities Inc., Washington DC
#
import sys
import os
#from casac import *
import casac
import string
import time
import inspect
import gc
import numpy
from odict import odict
from taskmanager import tm
from task_imfit import imfit
class imfit_cli_:
    __name__ = "imfit"
    __async__ = {}
    rkey = None
    i_am_a_casapy_task = None
    # The existence of the i_am_a_casapy_task attribute allows help()
    # (and other) to treat casapy tasks as a special case.

    def __init__(self) :
       self.__bases__ = (imfit_cli_,)
       self.__doc__ = self.__call__.__doc__

       self.parameters={'imagename':None, 'box':None, 'region':None, 'chans':None, 'stokes':None, 'mask':None, 'includepix':None, 'excludepix':None, 'residual':None, 'model':None, 'estimates':None, 'logfile':None, 'append':None, 'newestimates':None, 'complist':None, 'overwrite':None, 'dooff':None, 'offset':None, 'fixoffset':None, 'stretch':None,  'async':None}


    def result(self, key=None):
            #### here we will scan the task-ids in __async__
            #### and add any that have completed...
            if key is not None and self.__async__.has_key(key) and self.__async__[key] is not None:
               ret = tm.retrieve(self.__async__[key])
               if ret['state'] == "done" :
                  self.__async__[key] = None
               elif ret['state'] == 'crashed' :
                  self.__async__[key] = None
               return ret
            return None


    def __call__(self, imagename=None, box=None, region=None, chans=None, stokes=None, mask=None, includepix=None, excludepix=None, residual=None, model=None, estimates=None, logfile=None, append=None, newestimates=None, complist=None, overwrite=None, dooff=None, offset=None, fixoffset=None, stretch=None,  async=None):

        """Fit one or more elliptical Gaussian components on an image region(s)
PARAMETER SUMMARY
imagename        Name of the input image
box              One or more box regions to use for fitting,
                 eg "100, 120, 200, 220, 300, 300, 400, 400"
                 to use two boxes. If both box and region
                 parameters are specified, box is used.
region           Region of interest. See help par.region for specification options.
stokes           Stokes parameter to fit. If blank, first polarization plane is used.
mask             Mask to use. See help par.mask. Default is none.
includepix       Range of pixel values to include for fitting. Array of two numeric
                 values assumed to have same units as image pixel values. Only one
                 of includepix or excludepix can be specified.
excludepix       Range of pixel values to exclude for fitting. Array of two numeric
                 values assumed to have same units as image pixel values. Only one
                 of includepix or excludepix can be specified.
residual         Name of the residual image to write.
model            Name of the model image to write.
estimates        Name of file containing initial estimates of component parameters
                 (see below for formatting details).
logfile          Name of file to write fit results.
append           If logfile exists, append to it (True) or overwrite it (False).
newestimates     File to write fit results which can be used as initial estimates
                 for next run.
complist         Name of output component list table.
overwrite        Overwrite component list table if it exists?
dooff            Simultaneously fit a zero-level offset?
offset           Initial estimate for the zero-level offset. Only used if dooff is True.
fixoffset        Hold zero-level offset constant during fit? Only used if dooff is True.
stretch          Stretch the input mask if necessary and possible. Only used if a mask is specified.
                 See help par.stretch.

OVERVIEW
imfit is used to fit one or more gaussians to sources in an image as
well as an optional zero-level offset. Fitting is limited to a single polarization
but can be performed over several contiguous spectral channels.
If the image has a clean beam, the report will contain both the convolved
and the deconvolved fit results. It is a simple wrapper around the ia.fitcomponents()
tool method.

When dooff is False, the method returns a dictionary with two keys, 'converged' and 'results'.
The value of 'converged' is a boolean array which indicates if the fit converged
on a channel by channel basis.
The value of 'results' is a dictionary representing a component list reflecting the
fit results. In the case of a deconvolved image, the sizes and position angles are those
of the source convolved with the restoring beam. For convenience, the 'results' dictionary can
be read into a compoenent list tool (default tool is named cl) using the fromrecord() method
for easier inspection using tool methods, eg

cl.fromrecord(res['results'])

If dooff is True, in addtion to the specified number of
gaussians, a zero-level offset will also be fit. The initial estimate for this
offset is specified using the offset parameter. Units are assumed to be the
same as the image brightness units.The zero level offset can be held constant during
the fit by specifying fixoffset=True. In the case of dooff=True, the returned
dictionary contains two additional keys, 'zerooff' and 'zeroofferr', which are both
arrays containing the the fitted zero level offset value and its error, respectively,
for each channel. Both sets of values are in image brightness units.
In cases where the fit did not converge, these values are set to NaN.

The region can either be specified by a box(es) or a region (there is no
garauntee which will be used if both are specified, so be certain only one is).
If specified using the box parameter, multiple boxes can be given using the format
box="blcx1, blcy1, trcx1, trcy1, blcx2, blcy2, trcx2, trcy2, ... , blcxN, blcyN, trcxN, trcyN"
where N is the number of boxes. In this case, the union of the specified boxes will be used.

Ranges of pixel values can be included or excluded from the fit via includepix and excludepix.

If specified, the task will write the residual and/or model images for
successful fits.

If an estimates file is not specified, the task will attempt to estimate
initial parameters and fit a single Gaussian. If a multiple Gaussian fit
is desired, the user must specify initial estimates via a text file
(see below for details).

The user has the option of writing the result of the fit to a log file,
and has the option of either appending to or overwriting an existing file.

The user has the option of writing the (convolved) parameters of a successful
fit to a file which can be fed back to imfit as the estimates file for a
subsequent run.

FITTING OVER MULTIPLE CHANNELS

For fitting over multiple channels, the result
of the previous successful fit is used as the estimate for the next channel. The number
of gaussians fit cannot be varied on a channel by channel basis. Thus the variation
of source structure should be reasonably smooth in frequency to produce reliable fit
results.

MASK SPECIFICATION

Mask specification can be done using an LEL expression. For example

mask = '"myimage">5' will use only pixels with values greater than 5.

INCLUDING AND EXCLUDING PIXELS

Pixels can be included or excluded from the fit based on their values
using includepix and excludepix. Note that specifying both is not permitted
and will result in  an error. If specified, both take an array of two numeric
values.

ESTIMATES

Initial estimates of fit parameters may be specified via an estimates
text file. Each line of this file should contain a set of parameters for
a single gaussian. Optionally, some of these parameters can be fixed during
the fit. The format of each line is

peak intensity, peak x-pixel value, peak y-pixel value, major axis, minor axis, position angle, fixed

The fixed parameter is optional. The peak intensity is assumed to be in the
same units as the image pixel values (eg Jy/beam). The peak coordinates are specified
by pixel values. The major and minor axes and the position angle are the convolved
parameters if the image has been convolved with a clean beam and are specified as quantities.
The fixed parameter is optional and is a string. It may contain any combination of the
following characters 'f' (peak intensity), 'x' (peak x position), 'y' (peak y position),
'a' (major axis), 'b' (minor axis), 'p' (position angle).

In addition, lines in the file starting with a "#" are considered comments.

An example of such a file is:

# peak intensity must be in map units
120, 150, 110, 23.5arcsec, 18.9arcsec, 120deg  
90, 60, 200, 46arcsec, 23arcsec, 140deg, fxp

This is a file which specifies that two gaussians are to be simultaneously fit,
and for the second gaussian the specified peak intensity, x position, and position angle
are to be held fixed during the fit.

ERROR ESTIMATES

The reported estimated errors are the maxima from two methods. In the first method, errors are calculated
according to the description given in AIPS++ Note 224: AIPS++ Least Squares Background available at
http://www.astron.nl/casacore/trunk/casacore/doc/notes/224.html . In the second method (Condon et al
1998AJ.115.1693; Richards et al. 1999MNRAS.306..954), the following algorithms are used to compute the
various errors.

R = rms of residual image
P = fitted peak intensity
S = P/R (signal to noise)

longitude error = (FWHM beam width in longitude direction)/(2*S)
latitude error = (FWHM beam width in latitude direction)/(2*S)
FWHM of major axis error = (FWHM beam along major axis)/S
FWHM of minor axis error = (FWHM beam along minor axis)/S
position angle error is calculated via standard error propogation of the FWHM major and minor axes sizes and errors
flux error is computed via standard error propogation of the peak intensity and its error and the FWHM major
and minor axes sizes and their errors.

In the case where the image has no beam, the pixel width is used as the diameter of the resolution element.

In the vast majority of (and perhaps all) cases, the errors reported come from the second method.

The deconvolved size and position angle errors are computed by taking the maximum of the absolute values of the
differences of the best fit deconvolved value of the given parameter and the deconvolved size of the eight
possible combinations of (FWHM major axis +/- major axis error), (FWHM minor axis +/- minor axis error),
and (position andle +/- position angle error).

CAUTIONARY NOTES ON ERRORS
imfit makes the following implicit assumptions:

1. The noise associated with the pixels it is given to fit is truly random and gaussian
2. There is no offset (constant (unless one is fit), gradient, or more complex function) in pixel values.
3. The things that are trying to be fit are bona-fide gaussians.

If any of these assumptions is not true, the errors reported will in general be lower limits. These
assumptions do not hold if for example, there is structure in the noise of the image (eg the residual map produced
by deconvolution has structure), the area being fit lies in a negative bowl, in extended emission, or on an image
artifact (eg ripple), or the source(s) is nongaussian (eg a disk or more complex function convolved with a gaussian;
just because something looks gaussian to the human image analysis system does not necessarily mean it is).
In addition, in the case of an image produced from interferometric observations, the errors given are appropriate
for good uv plane coverage. They may be a factor of two or more higher for sparse arrays or snapshots.
In short, imfit is not a silver bullet for determining parameters in flawed images and reported errors
should be interpreted cautiously in the vast majority of all "real life" cases. The closer the source being fit is to
fulfilling these assumptions the more realistic reported errors will be.

EXAMPLE: 

Here is how one might fit two gaussians to multiple channels of a cube using the fit
from the previous channel as the initial estimate for the next. It also illustrates
how one can specify a region in the associated continuum image as the region to use
as the fit for the channel.

default imfit
imagename = "co_cube.im"
# specify region using region from continuum
region = "continuum.im:source.rgn"
chans = "2~20"
# only use pixels with positive values in the fit
excludepix = [-1e10,0]
# estimates file contains initial parameters for two Gaussians in channel 2
estimates = "initial_estimates.txt"
logfile = "co_fit.log"
# append results to the log file for all the channels
append = "True"
imfit()


        """
        if not hasattr(self, "__globals__") or self.__globals__ == None :
           self.__globals__=sys._getframe(len(inspect.stack())-1).f_globals
        #casac = self.__globals__['casac']
        casalog = self.__globals__['casalog']
        #casalog = casac.casac.logsink()
        self.__globals__['__last_task'] = 'imfit'
        self.__globals__['taskname'] = 'imfit'
        ###
        self.__globals__['update_params'](func=self.__globals__['taskname'],printtext=False,ipython_globals=self.__globals__)
        ###
        ###
        #Handle globals or user over-ride of arguments
        #
        function_signature_defaults=dict(zip(self.__call__.func_code.co_varnames,self.__call__.func_defaults))
        useLocalDefaults = False

        for item in function_signature_defaults.iteritems():
                key,val = item
                keyVal = eval(key)
                if (keyVal == None):
                        #user hasn't set it - use global/default
                        pass
                else:
                        #user has set it - use over-ride
                        if (key != 'self') :
                           useLocalDefaults = True

        myparams = {}
        if useLocalDefaults :
           for item in function_signature_defaults.iteritems():
               key,val = item
               keyVal = eval(key)
               exec('myparams[key] = keyVal')
               self.parameters[key] = keyVal
               if (keyVal == None):
                   exec('myparams[key] = '+ key + ' = self.itsdefault(key)')
                   keyVal = eval(key)
                   if(type(keyVal) == dict) :
                      if len(keyVal) > 0 :
                         exec('myparams[key] = ' + key + ' = keyVal[len(keyVal)-1][\'value\']')
                      else :
                         exec('myparams[key] = ' + key + ' = {}')

        else :
            async = self.parameters['async']
            myparams['imagename'] = imagename = self.parameters['imagename']
            myparams['box'] = box = self.parameters['box']
            myparams['region'] = region = self.parameters['region']
            myparams['chans'] = chans = self.parameters['chans']
            myparams['stokes'] = stokes = self.parameters['stokes']
            myparams['mask'] = mask = self.parameters['mask']
            myparams['includepix'] = includepix = self.parameters['includepix']
            myparams['excludepix'] = excludepix = self.parameters['excludepix']
            myparams['residual'] = residual = self.parameters['residual']
            myparams['model'] = model = self.parameters['model']
            myparams['estimates'] = estimates = self.parameters['estimates']
            myparams['logfile'] = logfile = self.parameters['logfile']
            myparams['append'] = append = self.parameters['append']
            myparams['newestimates'] = newestimates = self.parameters['newestimates']
            myparams['complist'] = complist = self.parameters['complist']
            myparams['overwrite'] = overwrite = self.parameters['overwrite']
            myparams['dooff'] = dooff = self.parameters['dooff']
            myparams['offset'] = offset = self.parameters['offset']
            myparams['fixoffset'] = fixoffset = self.parameters['fixoffset']
            myparams['stretch'] = stretch = self.parameters['stretch']

        if type(includepix)==int: includepix=[includepix]
        if type(excludepix)==int: excludepix=[excludepix]

        result = None

#
#    The following is work around to avoid a bug with current python translation
#
        mytmp = {}

        mytmp['imagename'] = imagename
        mytmp['box'] = box
        mytmp['region'] = region
        mytmp['chans'] = chans
        mytmp['stokes'] = stokes
        mytmp['mask'] = mask
        mytmp['includepix'] = includepix
        mytmp['excludepix'] = excludepix
        mytmp['residual'] = residual
        mytmp['model'] = model
        mytmp['estimates'] = estimates
        mytmp['logfile'] = logfile
        mytmp['append'] = append
        mytmp['newestimates'] = newestimates
        mytmp['complist'] = complist
        mytmp['overwrite'] = overwrite
        mytmp['dooff'] = dooff
        mytmp['offset'] = offset
        mytmp['fixoffset'] = fixoffset
        mytmp['stretch'] = stretch
        pathname='file:///'+os.environ.get('CASAPATH').split()[0]+'/share/xml/'
        trec = casac.casac.utils().torecord(pathname+'imfit.xml')

        casalog.origin('imfit')
        try :
          #if not trec.has_key('imfit') or not casac.casac.utils().verify(mytmp, trec['imfit']) :
            #return False

          casac.casac.utils().verify(mytmp, trec['imfit'], True)
          scriptstr=['']
          saveinputs = self.__globals__['saveinputs']
          saveinputs('imfit', 'imfit.last', myparams, self.__globals__,scriptstr=scriptstr)
          if async :
            count = 0
            keybase =  time.strftime("%y%m%d.%H%M%S")
            key = keybase + "_" + str(count)
            while self.__async__.has_key(key) :
               count += 1
               key = keybase + "_" + str(count)
            result = tm.execute('imfit', imagename, box, region, chans, stokes, mask, includepix, excludepix, residual, model, estimates, logfile, append, newestimates, complist, overwrite, dooff, offset, fixoffset, stretch)
            print "Use: "
            print "      tm.retrieve(return_value) # to retrieve the status"
            print 
            self.rkey = key
            self.__async__[key] = result
          else :
              tname = 'imfit'
              spaces = ' '*(18-len(tname))
              casalog.post('\n##########################################'+
                           '\n##### Begin Task: ' + tname + spaces + ' #####')
              casalog.post(scriptstr[1][1:]+'\n', 'INFO')
              result = imfit(imagename, box, region, chans, stokes, mask, includepix, excludepix, residual, model, estimates, logfile, append, newestimates, complist, overwrite, dooff, offset, fixoffset, stretch)
              casalog.post('##### End Task: ' + tname + '  ' + spaces + ' #####'+
                           '\n##########################################')

        except Exception, instance:
          if(self.__globals__.has_key('__rethrow_casa_exceptions') and self.__globals__['__rethrow_casa_exceptions']) :
             raise
          else :
             #print '**** Error **** ',instance
             tname = 'imfit'
             casalog.post('An error occurred running task '+tname+'.', 'ERROR')
             pass

        gc.collect()
        return result
#
#
#
    def paramgui(self, useGlobals=True, ipython_globals=None):
        """
        Opens a parameter GUI for this task.  If useGlobals is true, then any relevant global parameter settings are used.
        """
        import paramgui
        if not hasattr(self, "__globals__") or self.__globals__ == None :
           self.__globals__=sys._getframe(len(inspect.stack())-1).f_globals

        if useGlobals:
            if ipython_globals == None:
                myf=self.__globals__
            else:
                myf=ipython_globals

            paramgui.setGlobals(myf)
        else:
            paramgui.setGlobals({})

        paramgui.runTask('imfit', myf['_ip'])
        paramgui.setGlobals({})

#
#
#
    def defaults(self, param=None, ipython_globals=None, paramvalue=None, subparam=None):
        if not hasattr(self, "__globals__") or self.__globals__ == None :
           self.__globals__=sys._getframe(len(inspect.stack())-1).f_globals
        if ipython_globals == None:
            myf=self.__globals__
        else:
            myf=ipython_globals

        a = odict()
        a['imagename']  = ''
        a['box']  = ''
        a['region']  = ''
        a['chans']  = ''
        a['stokes']  = ''
        a['mask']  = ''
        a['includepix']  = []
        a['excludepix']  = []
        a['residual']  = ''
        a['model']  = ''
        a['estimates']  = ''
        a['logfile']  = ''
        a['newestimates']  = ''
        a['complist']  = ''
        a['dooff']  = False

        a['async']=False
        a['mask'] = {
                    0:odict([{'notvalue':''}, {'stretch':False}])}
        a['logfile'] = {
                    0:odict([{'notvalue':''}, {'append':True}])}
        a['complist'] = {
                    0:odict([{'notvalue':''}, {'overwrite':False}])}
        a['dooff'] = {
                    0:odict([{'notvalue':False}, {'offset':0.0}, {'fixoffset':False}])}

### This function sets the default values but also will return the list of
### parameters or the default value of a given parameter
        if(param == None):
                myf['__set_default_parameters'](a)
        elif(param == 'paramkeys'):
                return a.keys()
        else:
            if(paramvalue==None and subparam==None):
               if(a.has_key(param)):
                  return a[param]
               else:
                  return self.itsdefault(param)
            else:
               retval=a[param]
               if(type(a[param])==dict):
                  for k in range(len(a[param])):
                     valornotval='value'
                     if(a[param][k].has_key('notvalue')):
                        valornotval='notvalue'
                     if((a[param][k][valornotval])==paramvalue):
                        retval=a[param][k].copy()
                        retval.pop(valornotval)
                        if(subparam != None):
                           if(retval.has_key(subparam)):
                              retval=retval[subparam]
                           else:
                              retval=self.itsdefault(subparam)
                     else:
                        retval=self.itsdefault(subparam)
               return retval


#
#
    def check_params(self, param=None, value=None, ipython_globals=None):
      if ipython_globals == None:
          myf=self.__globals__
      else:
          myf=ipython_globals
#      print 'param:', param, 'value:', value
      try :
         if str(type(value)) != "<type 'instance'>" :
            value0 = value
            value = myf['cu'].expandparam(param, value)
            matchtype = False
            if(type(value) == numpy.ndarray):
               if(type(value) == type(value0)):
                  myf[param] = value.tolist()
               else:
                  #print 'value:', value, 'value0:', value0
                  #print 'type(value):', type(value), 'type(value0):', type(value0)
                  myf[param] = value0
                  if type(value0) != list :
                     matchtype = True
            else :
               myf[param] = value
            value = myf['cu'].verifyparam({param:value})
            if matchtype:
               value = False
      except Exception, instance:
         #ignore the exception and just return it unchecked
         myf[param] = value
      return value
#
#
    def description(self, key='imfit', subkey=None):
        desc={'imfit': 'Fit one or more elliptical Gaussian components on an image region(s)',
               'imagename': 'Name of the input image',
               'box': 'Specify one or more box regions for the fit.',
               'region': 'Region. See help par.region for specs.',
               'chans': 'Spectral channels on which to perform fit.',
               'stokes': 'Stokes parameter to fit. If blank, first stokes plane is used.',
               'mask': 'Mask to use. See help par.mask. Default is none.',
               'includepix': 'Range of pixel values to include for fitting.',
               'excludepix': 'Range of pixel values to exclude for fitting.',
               'residual': 'Name of output residual image.',
               'model': 'Name of output model image.',
               'estimates': 'Name of file containing initial estimates of component parameters.',
               'logfile': 'Name of file to write fit results.',
               'append': 'If logfile exists, append to it if True or overwrite it if False',
               'newestimates': 'File to write fit results which can be used as initial estimates for next run.',
               'complist': 'Name of output component list table.',
               'overwrite': 'Overwrite component list table if it exists?',
               'dooff': 'Also fit a zero level offset? Default is False',
               'offset': 'Initial estimate of zero-level offset. Only used if doff is True. Default is 0.0',
               'fixoffset': 'Keep the zero level offset fixed during fit? Default is False ',
               'stretch': 'Stretch the mask if necessary and possible? See help par.stretch ',

               'async': 'If true the taskname must be started using imfit(...)'
              }

#
# Set subfields defaults if needed
#

        if(desc.has_key(key)) :
           return desc[key]

    def itsdefault(self, paramname) :
        a = {}
        a['imagename']  = ''
        a['box']  = ''
        a['region']  = ''
        a['chans']  = ''
        a['stokes']  = ''
        a['mask']  = ''
        a['includepix']  = []
        a['excludepix']  = []
        a['residual']  = ''
        a['model']  = ''
        a['estimates']  = ''
        a['logfile']  = ''
        a['append']  = True
        a['newestimates']  = ''
        a['complist']  = ''
        a['overwrite']  = False
        a['dooff']  = False
        a['offset']  = 0.0
        a['fixoffset']  = False
        a['stretch']  = False

        #a = sys._getframe(len(inspect.stack())-1).f_globals

        if self.parameters['mask']  != '':
            a['stretch'] = False

        if self.parameters['logfile']  != '':
            a['append'] = True

        if self.parameters['complist']  != '':
            a['overwrite'] = False

        if self.parameters['dooff']  != False:
            a['offset'] = 0.0
            a['fixoffset'] = False

        if a.has_key(paramname) :
              return a[paramname]
imfit_cli = imfit_cli_()