AIPS++ Newsletter
November 2000
AIPS++ Tools
|
Main Newsletter Index
Tool Articles:
Polarimetric Analysis of Images
Mosaicwizard
Vis Plane Calibration Components
GBT Observing Interface
|
Polarimetric Analysis of Images
Neil Killeen - ATNF
The Image tool is the basic AIPS++ tool with which you manipulate images.
Working with it is the Imagepol tool; its job is dedicated polarimetric
analysis of images.
Creating Imagepol tools
The imagepol tool is constructed from either an AIPS++ disk image
file, or from an Image tool. This image must have a Stokes coordinate
or an error will occur. The idea is that you will generally
provide an Image holding RA, DEC, Frequency and Stokes coordinates.
Here are some examples:
- include 'imagepol.g'
- p1 := imagepol('pks1333.iquv') # 1
- im1 := image('pks1333.iquv') # 2
- p2 := imagepol(im1) # 3
- p2.summary() # 4
The first example makes the Imagepol tool directly from a disk file
which has Stokes I, Q, U and V. The second line makes an Image tool
from the same file and then we construct the Imagepol tool from that
Image tool. We use the summary function
of the Imagepol tool to summarise the image (it's just the Image tool
function summary that is really invoked).
What can I do with it?
So far we have seen how to construct Imagepol tools, and
the use of the summary function. Let us assume from now
on that we constructed an Imagepol tool from an image
with RA,DEC, Frequency (many channels), and Stokes (IQUV say)
- include 'imagepol.g'
- p := imagepol('pks1333.iquv')
Recovering specific Stokes
You can recover a particular Stokes Image :
- i := p.stokesi() # Make Image tools
- v := p.stokesv()
- q := p.stokes('q')
-
- i.statistics() # Run Image tool functions
- v.view()
- i2 := i.subimage('pks1333.i') # Copy to disk file
-
- imagedones() # Destroy all Image tools when fed up
You can use either the dedicated functions such as stokesi, stokesq
or use the stokes function which takes an argument saying which
Stokes you want.
Note that the variables i, q and v are in fact Image tools.
These Image tools are actually `virtual reference images' (see the
August 2000 Newsletter);
they reference bits of the image from which the
Imagepol tool was originally constructed. You can use all of the
Image tool functions on them as usual (a couple of examples are shown).
Generating Secondary Images
You can generate all of the usual secondary polarimetric analysis
secondary images such as polarized intensity, linearly polarized intensity
position angle and so on. There are dedicated functions and
a generic one pol to which you
supply an argument saying which quantity you want (useful for
scripts). Here are a few examples.
- lpi := p.linpolint() # linearly polarized intensity
- lpa := p.linpolposang() # linearly polarized position angle
- tpi := p.totpolint)() # total polarized intensity
- flp := p.fraclinpol() # fractional linear polarization
- ftp := p.fractotpol() # fractional total polarization
- lpi2 := p.pol('lpi') # Specify which polarized quantity
-
- lpi.statistics() # Run Image tool function
Again all of the returned variables are virtual Image tools upon which
you can run all the Image tool functions.
Generating Error Images
You can also recover the error images assuming simplistic propagation
of Gaussian errors. In AIPS++ we do not (yet ?) associate an
error value with each image pixel value and propagate errors
automatically (very hard problem...). The errors are stored
in their own image as needed. The standard deviation of the thermal
noise is either worked out for you (with a clipped-mean algorithm)
from the data, or you can supply it if you know it better.
Some of the functions return a scalar (as the error is constant
over the image), some return an Image tool.
- sig := p.sigma(clip=5) # A scalar; best guess at thermal noise
- sigi := p.sigmastokes('i') # A scalar; standard deviation of noise
- sigi := p.sigmastokesi();
- sigv := p.sigmastokesv(); # A scalar
-
- siglpi := p.sigmalinpolint() # A scalar; linearly pol intensity
- siglpa := p.sigmalinpolposang() # An Image tool;
# linearly polarized position angle
- sigflp := p.sigmafraclinpol() # An Image tool; fractional linear pol
The philosophy behind the handling of errors is to defer it as long as
possible. Many of you are probably familiar with AIPS and Miriad where
we blank output images based on a variety of statistical tests. For
example, you might blank the linear polarization position angle image
when its error is greater than some value. Doing that means you have to
keep on regenerating the output image every time you want to try a
different blanking value or method.
So to avoid this we have tried to defer this step until you actually use
the images. For example
- lpa := p.linpolposang() # linearly polarized position angle
- siglpa := p.sigmalinpolposang() # and its error image
-
- lpa.statistics(mask='$siglpa<10') # Blank when error > 10 degrees
- lpa.view(mask='$siglpa<20') # Blank when error > 20 degrees
You can see we have used the mask argument of the
statistics and
view functions to do this on-the-fly. This
is very convenient. Note that we have used the $ syntax (e.g.
mask='$siglpa < 20') because the error Image tools are virtual (there is
no disk file associated with them); the only way to get at their values
is via the tool.
You can of course always copy virtual images to disk if you wish:
- siglpa2 := siglpa.subimage('lpa')
- siglpa2.done()
Finding the Rotation Measure
Now we come to the vexatious question of the Rotation Measure.
The Imagepol tool offers you two algorithms. The first is a
'traditional' approach of fitting position angle as a function of
frequency. The algorithm used is that of Leahy et al (Astronomy &
Astrophysics, 156, 234).
Recall that our Imagepol tool was constructed from an image holding (probably)
RA, DEC, Frequency, and Stokes coordinates.
- ok := p.rm(rm='rm', rmerr='rmerr', rmmax=800, maxpaerr=10)
For each spatial pixel, a spectrum (frequency) of position angles is
computed (from Q and U stored in the image) and fit. It is very
important to specify the rmmax for this algorithm. All RM algorithms of
this nature struggle with the n-pi ambiguity and this argument allows
you to constrain the process of handling it in some vaguely useful way.
We have also specified here a maximum allowed error in the position
angle; otherwise that frequency is rejected for that spatial pixel. We
output two images, the RM image and its error image. Note that these
are disk image files, not virtual images.
- ok := p.rm(rm='rm', rmerr='rmerr', rmmax=800, maxpaerr=10)
The Imagepol tool also offers a new and novel Fourier-based algorithm
(
Killeen, Fluke, Zhao and Ekers, 2000, submitted).
This algorithm generates an
output image which is the polarized intensity as a function of Rotation
Measure. It does not suffer from ambiguity and it can recover multiple
RM screens unresolved by the beam (as long as they are not along the
same line-of-sight).
- p := imagepoltestimage(outfile='iquv.im', rm="1e5 5e5", nx=8,
+ ny=8, nf=256, f0=1.4e9, bw=8e6)
- p.frm(amp='amp', pa='pa') # Fourier transform;
- amp := image('amp') # Look at polarized intensity
- amp.view() # And reorder to put RM along X-axis
In this example, we have used the other Imagepol constructor
imagepoltestimage to generate a test image with the specified Rotation
Measure values. The generated test image is 4-dimensional (RA/DEC; 8 by 8
pixels), Stokes (4 pixels; IQUV) and Frequency (256 pixels). The source
is just a constant I (if you don't add noise all spatial pixels will be
identical) and V. Q and U vary with frequency according to the given Rotation
Measures.
The Fourier algorithm Fourier transforms the Complex polarization
P = Q + iU image as a function of frequency to a function of Rotation Measure.
In the example we save images of the polarized intensity and position angle.
Figure T1: [Click on figure for higher resolution.]
The figure shows the result of viewing the polarized intensity
where we have reordered the display to put Rotation Measure along
the first axis. You can see the two RM components recovered;
the one to the left is the larger one.
What's missing
Presently there is no specialized display for polarimetric quantities.
For example, the traditional vector overlay of linearly polarization.
This will be developed for release 1.5. One can also imagine doing
interesting things with Complex representations of polarimetric
quantities. The Image tool does not support Complex images yet
(although of course we can create them and write them to disk with C++
code) and this path is also not yet well persued.
Functions for depolarization ratio and errors will also be in place for
release 1.5
|
Main Newsletter Index
Tool Articles:
Polarimetric Analysis of Images
Mosaicwizard
Vis Plane Calibration Components
GBT Observing Interface
|
The Mosaicwizard in AIPS++
Mark Holdaway - NRAO/Tucson
Glish, the scripting language and command interpreter for AIPS++,
allows the astronomer to have a great deal of control over a series of
complex processing steps, running tools and automating decision-making
based on the tools' output. But glish also permits astronomical
programmers to bundle up their understanding of complex data
processing into simple scripts or wizards which are very easy for
non-experts to use. Mosaicwizard and imagerwizard, included in release 1.4 of
AIPS++, are good examples of the simplicity which can be achieved using
glish in this manner.
To start the mosaicwizard, start up AIPS++ and type at the command line:
- include 'mosaicwizard.g';
- mosaicwizard();
Two new windows will appear, one for the mosaicwizard GUI, and one for
the scripter, which displays commands emitted by the mosaicwizard.
The mosaicwizard GUI guides you through the mosaicing process, asking
for any required information and making reasonable defaults whenever
possible. Each step of the process has a title (such as "Select an
AIPS++ MeasurementSet or UVFITS file"), along with a few sentences
that explain in more depth about what is required and what will be
happening. Below the explanatory message there will usually be one or
more widgets to accept your input. After you have given your input,
the green "Next" button moves you
along to the next step. The status line near the bottom of the GUI
tells you brief messages about what is going on, such as problems with
any input you've given the wizard, or whether an image file has been
created or not.
The first step of the mosaicwizard is to select a MeasurementSet or
UVFITS file to process. If you don't have one, just leave the input
field , and by default AIPS++ will create a MeasurementSet of a
seven-pointing VLA D array observation of Cas A at 8 GHz. If you
start from a FITS file, as in the case of the Cas A default, handling
MeasurementSet arcana will add some extra processing time to a few of
the initial steps, so you probably want to keep the MeasurementSet and
use it directly in subsequent work with the mosaicwizard. If you
already have MeasurementSets in the current directory, click on the
wrench icon to the right of the input field, and a file catalog will
show you all available MeasurementSets. You can then select a MeasurementSet
and click on "send & dismiss" to send the selected MeasurementSet to the
mosaicwizard GUI.
The next steps include reading in an optional image that can be used
as an initial starting model, selecting the spectral window to image,
and selecting the fields you wish to image. after completing this
input, you are ready to start thinking about how we want to make our
mosaic image.
Because mosaics often have several pointings, cover a large angle on
the sky and have many pixels in their final images, and they often
reproduce extended structure which takes a long time to deconvolve,
mosaics are famous for being slow. However, if you make a small
mosaic at low resolution, they aren't slow at all. By default, the
mosaicwizard will first make a low resolution image with pixels 4
times larger than full resolution, using only the inner 0.25 of the
(u,v) plane. The low resolution mosaic is followed by an intermediate
resolution image with pixels 2 times larger than full resolution,
using only the inner 0.5 of the (u,v) plane. Finally, a full
resolution image is made using all of the (u,v) data. In most cases,
the low or intermediate resolution image can be used as a starting
model for the next higher resolution image, hence saving on
deconvolution time: the high resolution image doesn't need to
deconvolve all the flux from scratch, it just needs to incrementally
fill in the fine details that were missing from the low resolution
model. The first input on the deconvolution page, "(u,v) scaling
parameter", controls the maximum baseline used, as a
fraction of the maximum baseline present in the MeasurementSet.
In addition to the (u,v) scaling, you can specify the visibility
weighting (by the way, when the (u,v) scaling is less than 1.0, a
Gaussian taper is also applied in the (u,v) plane so that there is no
sharp cutoff in the Fourier plane coverage), the number of Clean or
Mem iterations, the deconvolution algorithm (choose between the
mfclark, mfhogbom, and mfmultiscale Clean algorithms, or the
mfemptiness or mfentropy MEM algorithms), and the presence of a
progress display. The multiscale and MEM algorithms require
additional input on a supplemental page.
Note that nowhere does the user input the image cell size or the image
size. The cell size is set by the reciprocal of the maximum baseline
length used for each resolutions' image (divided by 3 for super
Nyquist sampling). Initially, the image size is set so that it covers
all selected pointings plus half a primary beam's worth of slack on
all sides. In later stages, or at higher resolution, the user is
invited to generate a mask region which will determine the image size.
While the low resolution deconvolution is progressing, you may want to
focus your attention on the mosaicwizard scripter window showing the
most important mosaicing commands, the AIPS++ Log Message window, and
the PGPlotter window (only if display progress was invoked) which
shows the peak image residual and deconvolved flux as a function of
iteration number. Because AIPS++ deconvolves the entire mosaic image
with a single approximate PSF, the multifield (ie, mf) algorithms are
implemented with a major/minor cycle framework: deconvolution can only
proceed so deeply until errors will be made due to the differences
between the approximate PSF and the actual PSF for each field. At the
end of each cycle, a User Choice window will ask you if you want to
proceed with the deconvolution. (If you don't choose after 30
seconds, the default is to continue.) At the start of the next major
cycle, the residuals will likely be a bit higher due to errors made
from cleaning with an inexact PSF, but these errors are corrected in
just a few iterations.
Figure T2: [Click on figure for higher resolution.]
The progress display in imager shows algorithm-specific indicators
of the deconvolution progress. For Multi-Scale Clean, the three panels
are the peak positive residual, the peak negative residual, and the
total flux as a function of iteration number. The different scales' data are
color coded. Each major cycle can clearly be seen by its early correction
of errors made in the previous cycle.
When deconvolution ends (by manually stopping at a major cycle, or by
achieving the requested number of iterations), the image is displayed
with the viewer, along with instructions for how to create an optional
rectangular or polygonal region. This mask region will determine the
angular size of the next higher resolution image and will also restrict
modeled emission to reside within that region.
Figure T3: [Click on figure for higher resolution.]
When deconvolution is complete, the restored low resolution image is
displayed by the viewer. You have the option of defining a rectangular or
polygonal mask for the next stage of deconvolution at high resolution.
When you click on "Next" to leave the viewer, the mosaicwizard leads
you back to set the deconvolution parameters once again. The (u,v)
scaling has automatically increased by a factor of 2. After you adjust
the deconvolution parameters to your liking, you are off making the intermediate
resolution image, then inspecting that image with the viewer, perhaps
prescribing a new mask. And then on to the full resolution mosaic, and
you are finished.
There is nothing special about increasing the resolution in steps of
2. You can jump the (u,v) scaling however you want. You could, of
course, make your initial image at full resolution, but I prefer to
make my mistakes at low resolution where they aren't quite so
glaringly obvious and they don't take so much time.
We are planning a number of enhancements to the mosaicwizard, such as the
addition of total power data, dealing with problematic point sources, and
handling spectral line data. Look for a more capable version of the mosaicwizard
in the release 1.5 of AIPS++.
For more detailed information about the
mosaicwizard, see the
reference manual
under Synthesis:imager:mosaicwizard. For more detailed information about
multifield imaging in AIPS++, see the multifield chapter of
getting results.
|
Main Newsletter Index
Tool Articles:
Polarimetric Analysis of Images
Mosaicwizard
Vis Plane Calibration Components
GBT Observing Interface
|
Synthesis III: Visibility plane calibration components
A. J. Kemball - NRAO/Socorro
The previous article in this
series provided a description of the Measurement Equation (ME), which
is the calibration formalism used in AIPS++. As described in that
article, the ME allows antenna-based uv- and image-plane calibration
terms, as well as additive and multiplicative baseline-based
corrections. This article considers the uv-plane calibration
components in more detail for the case of synthesis reduction, and
their identification with well-known physical effects.
The antenna-based uv-plane calibration components are represented in
the ME by [2,2] Jones matrices, Jvisi, in a
selected polarization representation, using either a circular or
linear polarization decomposition. The previous article described how
these calibration components are applied in the ME formalism. Further
information regarding this formalism can be found in
Cornwell (AIPS++ Note 183) and
Noordam (AIPS++ Note 182).
We enumerate the different types of Jvisi components here:
- Gi: This diagonal term represents the composite complex
electronic gain for all components located after the feed in the
signal path.
- Bi: This term denotes the bandpass response, and allows a
dependence on frequency in the complex gain terms.
- Di: This term represents the instrumental polarization
response of the feed.
- Ti: This term denotes atmospheric corrections.
The terms listed here can be solved for; there are also two terms of
this class in the ME which are either known a priori or pre-computed,
namely: Ci, the polarization configuration matrix, which performs
polarization conversion operations, and the parallactic angle
correction Pi. Note that a subset of the solvable, and non-solvable
terms may also may a dependence on angular position, and have
image-plane counterparts in the ME to reflect this. The image-plane
effects will be the subject of a future article in this series.
The letter designation [G, B, D, T] is used to identify and
select calibration component types in the user interface to the
imaging and calibration software in AIPS++. The two AIPS++ tools
most directly concerned with calibration and imaging are
calibrater and imager. The use of these tools is fully
described in the AIPS++ cookbook, particularly in the chapters on
synthesis calibration and imaging. These chapters detail all the steps
required to calibrate and image uv-data, and the associated
calibrater and imager tool functions required.
An important point to note is that the tools calibrater and
imager represent the uv- and image-plane sides of the Measurement
Equation respectively. Therefore imager offers functions
concerned with imaging, image-plane calibration components (such as
the primary beam correction) and deconvolution, while calibrater
is concerned with the application and solution for uv-plane
calibration components, as enumerated above. The tools communicate via
the MODEL_DATA and CORRECTED_DATA columns in the Measurement Set.
The imager tool fills the MODEL_DATA column by integrating
across the ME from right to left, starting from the source model or
image. For different ways in which to do this see imager.setjy()
and imager.ft(). The calibrater tool fills the
CORRECTED_DATA column by integrating across the ME from left to
right, and applying the known calibration components to the observed
data in this process. The associated tool functions for the latter
process are calibrater.setapply() and calibrater.correct().
The solvers can form chi squared at any point in the ME at which a
calibration component is being solved for; the solvers will be the
subject of the next article in this series.
|
Main Newsletter Index
Tool Articles:
Polarimetric Analysis of Images
Mosaicwizard
Vis Plane Calibration Components
GBT Observing Interface
|
GBT's Observer Interface Written in Glish
Rick Fisher - NRAO/Green Bank
Several years ago I threw together a prototype graphical interface for
the GBT specifically tailored to observers' control of the telescope.
This prototype was written in glish/Tk plus a bit of C++, mainly
because I was familiar with glish from early aips++ experience. Glish
does not benefit from the enormous user community of a language and
toolbox like Tcl/Tk, but this does not appear to have been a big
disadvantage. Much of the GUI design involved detailed creation of
widget combinations at a fairly low level to make optimum use of
screen space and quite a bit of widget interconnection to make the GUI
respond to interdependent parameters in a way that might be expected
by an observer.
The GBT observer interface prototype seemed to fit the requirements of
observer's well enough to warrant writing the real interface, now
called "GBT Observe," in glish/Tk after a redesign of the code
structure. The observer interface now consists of a simulator for
each GBT hardware device, written entirely in glish/Tk, a parameter
interface to the GBT's monitor and control software, a GUI for
interactive observing, and a basic telescope control scripting syntax
for programmed observing. The interactive and programmed modes of
observing are integrated so that one can use one or the other as
appropriate at any given moment. The observing script syntax parser
has been written by Darrell Schiebel in C++, and a few utilities and
an interface to the JPL solar system ephemeris have also been written
in C++. Nearly all of the other code is in glish simply because code
development is easier in the interpreted language, and computation
efficiency has seldom been an issue.
Telescope observing procedures, from simple tracking and on-offs to
raster maps, are written in glish and built into the observer
interface. These should suffice for the great majority of GBT
observing, but an observer is free to write his or her own procedures
using the templates provided. The scripting syntax for programmed
observing intentionally excludes arithmetic operators and most flow
control to keep the syntax simple. The more powerful language
features needed for procedure writing are already available in glish,
which most observers will know from their use of aips++.
GBT Observe is still under construction, mostly adding interfaces
to the many hardware modules in the system. The core interface to the
antenna and a couple of back-ends has been in use for about a year on
the GBT mockup in the Jansky lab. A few key interfaces need to be
completed in time for GBT commissioning, and detailed control of some
of the hardware not normally set up directly by observers will be
added as time permits.
A few web documents describing the basic layout from the observer's
point of view and the scripting language syntax may be found at
More complete observer documentation is being adapted for the main
GBT web pages.
A bit of programmer's documentation on the layout of the code
directories and the details of adding a hardware device interface to
GBT Observe is posted at
Since a GBT simulator is part of the observer's interface it
should be possible to export the entire interface to a user's home
computer for testing observing scripts and user-written procedures
before using them on the telescope. We still need to work out the
details of exporting the associated C++ glish clients. One way might
be to make the interface code part of the aips++ distribution.
Mark Holdaway
|
|