Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/imfit.py: 56%
27 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-10-31 17:39 +0000
1##################### generated by xml-casa (v2) from imfit.xml #####################
2##################### 78593ac558383844ca5bf3587a31cf5b ##############################
3from __future__ import absolute_import
4import numpy
5from casatools.typecheck import CasaValidator as _val_ctor
6_pc = _val_ctor( )
7from casatools.coercetype import coerce as _coerce
8from casatools.errors import create_error_string
9from .private.task_imfit import imfit as _imfit_t
10from casatasks.private.task_logging import start_log as _start_log
11from casatasks.private.task_logging import end_log as _end_log
12from casatasks.private.task_logging import except_log as _except_log
14class _imfit:
15 """
16 imfit ---- Fit one or more elliptical Gaussian components on an image region(s)
18 --------- parameter descriptions ---------------------------------------------
20 imagename Name of the input image
21 box Rectangular region(s) to select in direction plane. Default is to use the entire direction plane.
22 region Region selection. Default is to use the full image.
23 chans Channels to use. Default is to use all channels.
24 stokes Stokes planes to use. Default is to use first Stokes plane.
25 mask Mask to use. Default is none.
26 includepix Range of pixel values to include for fitting.
27 excludepix Range of pixel values to exclude for fitting.
28 residual Name of output residual image.
29 model Name of output model image.
30 estimates Name of file containing initial estimates of component parameters.
31 logfile Name of file to write fit results.
32 append If logfile exists, append to it if True or overwrite it if False
33 newestimates File to write fit results which can be used as initial estimates for next run.
34 complist Name of output component list table.
35 overwrite Overwrite component list table if it exists?
36 dooff Also fit a zero level offset? Default is False
37 offset Initial estimate of zero-level offset. Only used if doff is True. Default is 0.0
38 fixoffset Keep the zero level offset fixed during fit? Default is False
39 stretch Stretch the mask if necessary and possible?
40 rms RMS to use in calculation of uncertainties. Numeric or valid quantity (record or string). If numeric, it is given units of the input image. If quantity, units must conform to image units. If not positive, the rms of the residual image, in the region of the fit, is used.
41 noisefwhm Noise correlation beam FWHM. If numeric value, interpreted as pixel widths. If quantity (dictionary, string), it must have angular units.
42 summary File name to which to write table of fit parameters.
43 [1;42mRETURNS[1;m void
45 --------- examples -----------------------------------------------------------
48 PARAMETER SUMMARY
49 imagename Name of the input image
50 box Rectangular region(s) to select in direction plane.
51 for details. Default is to use the entire direction plane.
52 eg "100, 120, 200, 220, 300, 300, 400, 400" to use two boxes.
53 region Region selection. Default is to use
54 the full image.
55 chans Channels to use. Default is to use all
56 channels.
57 stokes Stokes planes to use. Default is to
58 use first Stokes plane.
59 mask Mask to use. Default is none.
60 includepix Range of pixel values to include for fitting. Array of two numeric
61 values assumed to have same units as image pixel values. Only one
62 of includepix or excludepix can be specified.
63 excludepix Range of pixel values to exclude for fitting. Array of two numeric
64 values assumed to have same units as image pixel values. Only one
65 of includepix or excludepix can be specified.
66 residual Name of the residual image to write.
67 model Name of the model image to write.
68 estimates Name of file containing initial estimates of component parameters
69 (see below for formatting details).
70 logfile Name of file to write fit results.
71 append If logfile exists, append to it (True) or overwrite it (False).
72 newestimates File to write fit results which can be used as initial estimates
73 for next run.
74 complist Name of output component list table.
75 overwrite Overwrite component list table if it exists?
76 dooff Simultaneously fit a zero-level offset?
77 offset Initial estimate for the zero-level offset. Only used if dooff is True.
78 fixoffset Hold zero-level offset constant during fit? Only used if dooff is True.
79 stretch Stretch the input mask if necessary and possible. Only used if a mask is specified.
81 rms RMS to use in calculation of various uncertainties, assumed to have units of the input
82 image. If not positve, the rms of the residual image is used.
83 noisefwhm Noise correlation beam FWHM. If numeric value, interpreted as pixel widths. If
84 quantity (dictionary, string), it must have angular units.
85 summary File name to which to write table of fit parameters.
87 OVERVIEW
88 This application is used to fit one or more two dimensional gaussians to sources in an image as
89 well as an optional zero-level offset. Fitting is limited to a single polarization
90 but can be performed over several contiguous spectral channels.
91 If the image has a clean beam, the report and returned dictionary will contain both the convolved
92 and the deconvolved fit results.
94 When dooff is False, the method returns a dictionary with three keys, 'converged', 'results',
95 and 'deconvolved'. The value of 'converged' is a boolean array which indicates if the fit
96 converged on a channel by channel basis. The value of 'results' is a dictionary representing
97 a component list reflecting the fit results. In the case of an image containing beam information,
98 the sizes and position angles in the 'results' dictionary are those of the source(s) convolved
99 with the restoring beam, while the same parameters in the 'deconvolved' dictionary represent the
100 source sizes deconvolved from the beam. In the case where the image does not contain a beam,
101 'deconvolved' will be absent. Both the 'results' and 'deconvolved' dictionaries can
102 be read into a component list tool (default tool is named cl) using the fromrecord() method
103 for easier inspection using tool methods, eg
105 cl.fromrecord(res['results'])
107 although this currently only works if the flux density units are conformant with Jy.
109 There are also values in each component subdictionary not used by cl.fromrecord() but meant to
110 supply additional information. There is a 'peak' subdictionary for each component that provides the
111 peak intensity of the component. It is present for both 'results' and 'deconvolved' components.
112 There is also a 'sum' subdictionary for each component indicated the simple sum of pixel values in
113 the the original image enclosed by the fitted ellipse. There is a 'channel' entry in the 'spectrum'
114 subdictionary which provides the zero-based channel number in the input image for which the solution
115 applies. In addtion, if the image has a beam(s), then there will be a 'beam' subdictionary associated
116 with each component in both the 'results' and 'deconvolved' dictionaries. This subdictionary will
117 have three keys: 'beamarcsec' will be a subdictionary giving the beam dimensions in arcsec,
118 'beampixels' will have the value of the beam area expressed in pixels, and 'beamster' will have the
119 value of the beam area epressed in steradians. Also, if the image has a beam(s), in the component level
120 dictionaries will be an 'ispoint' entry with an associated boolean value describing if the component
121 is consistent with a point source.
123 If dooff is True, in addtion to the specified number of
124 gaussians, a zero-level offset will also be fit. The initial estimate for this
125 offset is specified using the offset parameter. Units are assumed to be the
126 same as the image brightness units. The zero level offset can be held constant during
127 the fit by specifying fixoffset=True. In the case of dooff=True, the returned
128 dictionary contains two additional keys, 'zerooff' and 'zeroofferr', which are both
129 dictionaries containing 'unit' and 'value' keys. The values associated with the 'value'
130 keys are arrays containing the the fitted zero level offset value and its error, respectively,
131 for each channel. In cases where the fit did not converge, these values are set to NaN.
132 The value associated with 'unit' is just the i`mage brightness unit.
134 The region can either be specified by a box(es) or a region.
135 Ranges of pixel values can be included or excluded from the fit. If specified using
136 the box parameter, multiple boxes can be given using the format
137 box="blcx1, blcy1, trcx1, trcy1, blcx2, blcy2, trcx2, trcy2, ... , blcxN, blcyN, trcxN, trcyN"
138 where N is the number of boxes. In this case, the union of the specified boxes will be used.
140 If specified, the residual and/or model images for successful fits will be written.
142 If an estimates file is not specified, an attempt is made to estimate
143 initial parameters and fit a single Gaussian. If a multiple Gaussian fit
144 is desired, the user must specify initial estimates via a text file
145 (see below for details).
147 The user has the option of writing the result of the fit to a log file,
148 and has the option of either appending to or overwriting an existing file.
150 The user has the option of writing the (convolved) parameters of a successful
151 fit to a file which can be fed back to fitcomponents() as the estimates file for a
152 subsequent run.
154 The user has the option of writing the fit results in tabular format to a file whose
155 name is specified using the summary parameter.
157 If specified and positive, the value of rms is used to calculate the parameter uncertainties,
158 otherwise, the rms in the selected region in the relevant channel is used for these calculations.
160 The noisefwhm parameter represents the noise-correlation beam FWHM. If specified as a quantity,
161 it should have angular units. If specified as a numerical value, it is set equal to that number
162 of pixels. If specified and greater than or equal to the pixel size, it is used to calculate
163 parameter uncertainties using the correlated noise equations (see below). If it is specified but
164 less than a pixel width, the the uncorrelated noise equations (see below) are used to
165 compute the parameter uncertainties. If it is not specified and the image has a restoring beam(s),
166 the the correlated noise equations are used to compute parameter uncertainties using the
167 geometric mean of the relevant beam major and minor axes as the noise-correlation beam FWHM. If
168 noisefwhm is not specified and the image does not have a restoring beam, then the uncorrelated
169 noise equations are used to compute the parameter uncertainties.
171 SUPPORTED UNITS
173 Currently only images with brightness units conformant with Jy/beam, Jy.km/s/beam, and K are fully
174 supported for fitting. If your image has some other base brightness unit, that unit will be assumed
175 to be equivalent to Jy/pixel and results will be calculated accordingly. In particular,
176 the flux density (reported as Integrated Flux in the logger and associated with the "flux" key
177 in the returned component subdictionary(ies)) for such a case represents the sum of pixel values.
179 Note also that converting the returned results subdictionary to a component list via cl.fromrecord() currently
180 only works properly if the flux density units in the results dictionary are conformant with Jy.
181 If you need to be able to run cl.fromrecord() on the resulting dictionary you can first modify the
182 flux density units by hand to be (some prefix)Jy and then run cl.fromrecord() on that dictionary,
183 bearing in mind your unit conversion.
185 If the input image has units of K, the flux density of components will be reported in units
186 of [prefix]K*rad*rad, where prefix is an SI prefix used so that the numerical value is between
187 1 and 1000. To convert to units of K*beam, determine the area of the appropriate beam,
188 which is given by pi/(4*ln(2))*bmaj*bmin, where bmaj and bmin are the major and minor axes
189 of the beam, and convert to steradians (=rad*rad). This value is included in the beam portion
190 of the component subdictionary (key 'beamster'). Then divide the numerical value of the
191 logged flux density by the beam area in steradians. So, for example
193 begin{verbatim}
194 # run on an image with K brightness units
195 res = imfit(...)
196 # get the I flux density in K*beam of component 0
197 comp = res['results']['component0']
198 flux_density_kbeam = comp['flux']['value'][0]/comp['beam']['beamster']
199 end{verbatim}
201 FITTING OVER MULTIPLE CHANNELS
203 For fitting over multiple channels, the result of the previous successful fit is used as
204 the estimate for the next channel. The number of gaussians fit cannot be varied on a channel
205 by channel basis. Thus the variation of source structure should be reasonably smooth in
206 frequency to produce reliable fit results.
208 MASK SPECIFICATION
210 Mask specification can be done using an LEL expression. For example
212 mask = '"myimage">5' will use only pixels with values greater than 5.
214 INCLUDING AND EXCLUDING PIXELS
216 Pixels can be included or excluded from the fit based on their values
217 using these parameters. Note that specifying both is not permitted and
218 will cause an error. If specified, both take an array of two numeric
219 values.
221 ESTIMATES
223 Initial estimates of fit parameters may be specified via an estimates
224 text file. Each line of this file should contain a set of parameters for
225 a single gaussian. Optionally, some of these parameters can be fixed during
226 the fit. The format of each line is
228 peak intensity, peak x-pixel value, peak y-pixel value, major axis, minor axis, position angle, fixed
230 The fixed parameter is optional. The peak intensity is assumed to be in the
231 same units as the image pixel values (eg Jy/beam). The peak coordinates are specified
232 in pixel coordinates. The major and minor axes and the position angle are the convolved
233 parameters if the image has been convolved with a clean beam and are specified as quantities.
234 The fixed parameter is optional and is a string. It may contain any combination of the
235 following characters 'f' (peak intensity), 'x' (peak x position), 'y' (peak y position),
236 'a' (major axis), 'b' (minor axis), 'p' (position angle).
238 In addition, lines in the file starting with a # are considered comments.
240 An example of such a file is:
242 begin{verbatim}
243 # peak intensity must be in map units
244 120, 150, 110, 23.5arcsec, 18.9arcsec, 120deg
245 90, 60, 200, 46arcsec, 23arcsec, 140deg, fxp
246 end{verbatim}
249 This is a file which specifies that two gaussians are to be simultaneously fit,
250 and for the second gaussian the specified peak intensity, x position, and position angle
251 are to be held fixed during the fit.
253 ERROR ESTIMATES
255 Error estimates are based on the work of Condon 1997, PASP, 109, 166. Key assumptions made are:
256 * The given model (elliptical Gaussian, or elliptical Gaussian plus constant offset) is an
257 adequate representation of the data
258 * An accurate estimate of the pixel noise is provided or can be derived (see above). For the
259 case of correlated noise (e.g., a CLEAN map), the fit region should contain many "beams" or
260 an independent value of rms should be provided.
261 * The signal-to-noise ratio (SNR) or the Gaussian component is large. This is necessary because
262 a Taylor series is used to linearize the problem. Condon (1997) states that the fractional
263 bias in the fitted amplitude due to this assumption is of order 1/(S*S), where S is the overall
264 SNR of the Gaussian with respect to the given data set (defined more precisely below). For a 5
265 sigma "detection" of the Gaussian, this is a 4% effect.
266 * All (or practically all) of the flux in the component being fit falls within the selected region.
267 If a constant offset term is simultaneously fit and not fixed, the region of interest should be
268 even larger. The derivations of the expressions summarized in this note assume an effectively
269 infinite region.
271 Two sets of equations are used to calculate the parameter uncertainties, based on if
272 the noise is correlated or uncorrelated. The rules governing which set of equations are
273 used have been described above in the description of the noisefwhm parameter.
275 In the case of uncorrelated noise, the equations used are
277 f(A) = f(I) = f(M) = f(m) = k*s(x)/M = k*s(y)/m = (s(p)/sqrt(2))*((M*M - m*m)/(M*m))
278 = sqrt(2)/S
280 where s(z) is the uncertainty associated with parameter z, f(z) = s(z)/abs(z) is the
281 fractional uncertainty associated with parameter z, A is the peak intensity, I is the flux
282 density, M and m are the FWHM major and minor axes, p is the position angle of the
283 component, and k = sqrt(8*ln(2)). s(x) and s(y) are the direction
284 uncertainties of the component measured along the major and minor axes; the resulting
285 uncertainties measured along the principle axes of the image direction coordinate are
286 calculated by propagation of errors using the 2D rotation matrix which enacts the rotation through
287 the position angle plus 90 degrees. S is the overall signal to noise ratio of the component,
288 which, for the uncorrelated noise case is given by
290 S = (A/(k*h*r))*sqrt(pi*M*m)
292 where h is the pixel width of the direction coordinate and r is the rms noise (see the
293 discussion above for the rules governing how the value of r is determined).
295 For the correlated noise case, the same equations are used to determine the uncertainties
296 as in the uncorrelated noise case, except for the uncertainty in I (see below). However,
297 S is given by
299 S = (A/(2*r*N)) * sqrt(M*m) * (1 + ((N*N/(M*M)))**(a/2)) * (1 + ((N*N/(m*m)))**(b/2))
301 where N is the noise-correlation beam FWHM (see discussion of the noisefwhm parameter for
302 rules governing how this value is determined). "**" indicates exponentiation and a and b
303 depend on which uncertainty is being calculated. For sigma(A), a = b = 3/2. For M and x,
304 a = 5/2 and b = 1/2. For m, y, and p, a = 1/2 and b = 5/2. f(I) is calculated in the
305 correlated noise case according to
307 f(I) = sqrt( f(A)*f(A) + (N*N/(M*m))*(f(M*f(M) + f(m)*f(m))) )
309 Note well the following caveats:
310 * Fixing Gaussian component parameters will tend to cause the parameter uncertainties reported for free
311 parameters to be overestimated.
312 * Fitting a zero level offset that is not fixed will tend to cause the reported parameter
313 uncertainties to be slightly underestimated.
314 * The parameter uncertainties will be inaccurate at low SNR (a ~10% for SNR = 3).
315 * If the fitted region is not considerably larger than the largest component that is fit,
316 parameter uncertainties may be mis-estimated.
317 * An accurate rms noise measurement, r, for the region in question must be supplied.
318 Alternatively, a sufficiently large signal-free region must be present in the selected region
319 (at least about 25 noise beams in area) to auto-derive such an estimate.
320 * If the image noise is not statistically independent from pixel to pixel, a reasonably accurate noise
321 correlation scale, N, must be provided. If the noise correlation function is not approximately Gaussian,
322 the correlation length can be estimated using
324 N = sqrt(2*ln(2)/pi)* double-integral(dx dy C(x,y))/sqrt(double-integral(dx dy C(x, y) * C(x,y)))
326 where C(x,y) is the associated noise-smoothing function
327 * If fitted model components have significan spatial overlap, the parameter uncertainties are likely to
328 be mis-estimated (i.e., correlations between the parameters of separate components are not accounted
329 for).
330 * If the image being analyzed is an interferometric image with poor uv sampling, the parameter
331 uncertainties may be significantly underestimated.
333 The deconvolved size and position angle errors are computed by taking the maximum of the absolute values of the
334 differences of the best fit deconvolved value of the given parameter and the deconvolved size of the eight
335 possible combinations of (FWHM major axis +/- major axis error), (FWHM minor axis +/- minor axis error),
336 and (position andle +/- position angle error). If the source cannot be deconvolved from the beam (if the best
337 fit convolved source size cannot be deconvolved from the beam), upper limits on the deconvolved source size
338 are sometimes reported. These limits simply come from the maximum major and minor axes of the deconvolved
339 gaussians taken from trying all eight of the aforementioned combinations. In the case none of these combinations
340 produces a deconvolved size, no upper limit is reported.
342 EXAMPLE:
344 Here is how one might fit two gaussians to multiple channels of a cube using the fit
345 from the previous channel as the initial estimate for the next. It also illustrates
346 how one can specify a region in the associated continuum image as the region to use
347 as the fit for the channel.
350 begin{verbatim}
351 default imfit
352 imagename = "co_cube.im"
353 # specify region using region from continuum
354 region = "continuum.im:source.rgn"
355 chans = "2~20"
356 # only use pixels with positive values in the fit
357 excludepix = [-1e10,0]
358 # estimates file contains initial parameters for two Gaussians in channel 2
359 estimates = "initial_estimates.txt"
360 logfile = "co_fit.log"
361 # append results to the log file for all the channels
362 append = "True"
363 imfit()
366 """
368 _info_group_ = """analysis"""
369 _info_desc_ = """Fit one or more elliptical Gaussian components on an image region(s)"""
371 def __call__( self, imagename='', box='', region='', chans='', stokes='', mask='', includepix=[ ], excludepix=[ ], residual='', model='', estimates='', logfile='', append=True, newestimates='', complist='', overwrite=False, dooff=False, offset=float(0.0), fixoffset=False, stretch=False, rms=int(0), noisefwhm='', summary='' ):
372 schema = {'imagename': {'type': 'cReqPath', 'coerce': _coerce.expand_path}, 'box': {'type': 'cStr', 'coerce': _coerce.to_str}, 'region': {'type': 'cVariant', 'coerce': [_coerce.to_variant]}, 'chans': {'type': 'cVariant', 'coerce': [_coerce.to_variant]}, 'stokes': {'type': 'cStr', 'coerce': _coerce.to_str}, 'mask': {'type': 'cStr', 'coerce': _coerce.to_str}, 'includepix': {'type': 'cFloatVec', 'coerce': [_coerce.to_list,_coerce.to_floatvec]}, 'excludepix': {'type': 'cFloatVec', 'coerce': [_coerce.to_list,_coerce.to_floatvec]}, 'residual': {'type': 'cStr', 'coerce': _coerce.to_str}, 'model': {'type': 'cStr', 'coerce': _coerce.to_str}, 'estimates': {'type': 'cStr', 'coerce': _coerce.to_str}, 'logfile': {'type': 'cStr', 'coerce': _coerce.to_str}, 'append': {'type': 'cBool'}, 'newestimates': {'type': 'cStr', 'coerce': _coerce.to_str}, 'complist': {'type': 'cStr', 'coerce': _coerce.to_str}, 'overwrite': {'type': 'cBool'}, 'dooff': {'type': 'cBool'}, 'offset': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'fixoffset': {'type': 'cBool'}, 'stretch': {'type': 'cBool'}, 'rms': {'anyof': [{'type': 'cInt'}, {'type': 'cFloat', 'coerce': _coerce.to_float}, {'type': 'cDict'}, {'type': 'cStr', 'coerce': _coerce.to_str}]}, 'noisefwhm': {'anyof': [{'type': 'cInt'}, {'type': 'cFloat', 'coerce': _coerce.to_float}, {'type': 'cDict'}, {'type': 'cStr', 'coerce': _coerce.to_str}]}, 'summary': {'type': 'cStr', 'coerce': _coerce.to_str}}
373 doc = {'imagename': imagename, 'box': box, 'region': region, 'chans': chans, 'stokes': stokes, 'mask': mask, 'includepix': includepix, 'excludepix': excludepix, 'residual': residual, 'model': model, 'estimates': estimates, 'logfile': logfile, 'append': append, 'newestimates': newestimates, 'complist': complist, 'overwrite': overwrite, 'dooff': dooff, 'offset': offset, 'fixoffset': fixoffset, 'stretch': stretch, 'rms': rms, 'noisefwhm': noisefwhm, 'summary': summary}
374 assert _pc.validate(doc,schema), create_error_string(_pc.errors)
375 _logging_state_ = _start_log( 'imfit', [ 'imagename=' + repr(_pc.document['imagename']), 'box=' + repr(_pc.document['box']), 'region=' + repr(_pc.document['region']), 'chans=' + repr(_pc.document['chans']), 'stokes=' + repr(_pc.document['stokes']), 'mask=' + repr(_pc.document['mask']), 'includepix=' + repr(_pc.document['includepix']), 'excludepix=' + repr(_pc.document['excludepix']), 'residual=' + repr(_pc.document['residual']), 'model=' + repr(_pc.document['model']), 'estimates=' + repr(_pc.document['estimates']), 'logfile=' + repr(_pc.document['logfile']), 'append=' + repr(_pc.document['append']), 'newestimates=' + repr(_pc.document['newestimates']), 'complist=' + repr(_pc.document['complist']), 'overwrite=' + repr(_pc.document['overwrite']), 'dooff=' + repr(_pc.document['dooff']), 'offset=' + repr(_pc.document['offset']), 'fixoffset=' + repr(_pc.document['fixoffset']), 'stretch=' + repr(_pc.document['stretch']), 'rms=' + repr(_pc.document['rms']), 'noisefwhm=' + repr(_pc.document['noisefwhm']), 'summary=' + repr(_pc.document['summary']) ] )
376 task_result = None
377 try:
378 task_result = _imfit_t( _pc.document['imagename'], _pc.document['box'], _pc.document['region'], _pc.document['chans'], _pc.document['stokes'], _pc.document['mask'], _pc.document['includepix'], _pc.document['excludepix'], _pc.document['residual'], _pc.document['model'], _pc.document['estimates'], _pc.document['logfile'], _pc.document['append'], _pc.document['newestimates'], _pc.document['complist'], _pc.document['overwrite'], _pc.document['dooff'], _pc.document['offset'], _pc.document['fixoffset'], _pc.document['stretch'], _pc.document['rms'], _pc.document['noisefwhm'], _pc.document['summary'] )
379 except Exception as exc:
380 _except_log('imfit', exc)
381 raise
382 finally:
383 task_result = _end_log( _logging_state_, 'imfit', task_result )
384 return task_result
386imfit = _imfit( )