SDCLN: Task to deconvolve image with CLEAN & Steer/Dewdney/Ito
INPUTS
Dirty image
INNAME Image name (name)
INCLASS Image name (class)
INSEQ 0.0 9999.0 Image name (seq. #)
INDISK 0.0 9.0 Image disk drive #
BLC 0.0 Used to specify plane.
Beam image
IN2NAME Image name (name)
IN2CLASS Image name (class)
IN2SEQ 0.0 9999.0 Image name (seq. #)
IN2DISK 0.0 9.0 Image disk drive #
Clean image
OUTNAME Image name (name)
OUTCLASS Image name (class)
OUTSEQ -1.0 9999.0 Image name (seq. #)
OUTDISK 0.0 9.0 Image disk drive #
INVERS -1.0 46655.0 CC file ver. no.
GAIN * 0.0 2.0 CLEAN loop gain
FLUX * 0.0 Min. residual in window (Jy)
PHAT 0.0 100.0 Prussian hat size.
NCOUNT * -1.0 # pixels > CUTOFF*peak to
switch from CLEAN to SDI
0 => no SDI method
CUTOFF 0.0 1.0 Cutoff for inclusion in SDI
STFACTOR * 0.0 1.0 SDI loop gain
NITER * Maximum # of CLEAN components
BITER Begin at BITER component
BMAJ * -999.9 FWHM(asec) maj. axis CLEAN
restoring beam.
BMIN * -999.9 FWHM(asec) min. axis CLEAN
restoring beam.
BPA * -360.0 360.0 CLEAN beam position angle
NBOXES 0.0 50.0 Number of boxes for CLEAN
CLBOX 0.0 4096.0 Four coordinates for each box
FACTOR * -5.0 5.0 Speedup factor see HELP SDCLN
MINPATCH * 0.0 1024.0 Min. beam halfwidth in AP.
MAXPIXEL * 0.0 500000.0 Maximum pixels searched in
each major cycle.
DOTV * -1.0 1.0 TV display of residuals?
BADDISK -1.0 1000.0 Disks to avoid for scratch.
HELP SECTION
SDCLN
Type: Task
Use: SDCLN does an AP-based CLEAN of the type devised by B. G.
Clark. CLEAN components are found by CLEANing the
brightest points in the residual map with a small area of
the dirty beam. More exact CLEANing is then done by
transforming the CLEAN components, multiplying by the
transform of the dirty beam, then transforming back to the
map plane and subtracting the result from the dirty map.
SDCLN may be restarted if further CLEANing is desired. The
dirty beam may be modified by adding a delta function of
amplitude PHAT centered on its peak. TJC claims that this
encourages smoothness in the clean map. SDCLN differs from
APCLN in that it offers the option to switch from the
standard CLEAN method of finding components to a Steer -
Dewdney - Ito (SDI) method. In the latter, all pixels
having residual flux > CUTOFF times the current peak
residual are taken as components. This method should avoid
the "ripples" produced by the standard CLEAN on extended
emission.
A single plane in a cube may be processed per run; the
resulting CLEAN (or residual) image will be put in the
corresponding plane of the output image. Other planes in
the output image will not be affected, i.e. are likely to
contain garbage. When several runs of SDCLN are made to
CLEAN a cube be sure to completely specify the output
image. A clean components table will be created with a
sequence number corresponding to the plane number. At the
end of each execution of SDCLN, the output components table
is compressed so that a pixel will appear only once in the
table. This avoids the necessity of running CCMRG.
SDCLN now partially supports assymmetric beams, see the
EXPLAIN section for more details. SDCLN supports TELL
including some change of parameters: see SHOW SDCLN.
Adverbs:
INNAME......The dirty map image name. Standard defaults.
INCLASS.....The dirty map image class. Standard defaults.
INSEQ.......The dirty map image sequence number. 0 => high
INDISK......The dirty map image disk drive no. 0 => any
BLC.........The third and higher dimensions are used to
specify the plane in the dirty and residual
images to be processed. NOTE: process in
increasing plane number.
IN2NAME.....The dirty beam image name. blank => actual
INNAME, otherwise standard defaults.
IN2CLASS....The dirty beam image class. Standard defaults
except blank => 'RBM001' if INCLASS = 'RIM001'
'LBM001' if INCLASS = 'LIM001'
'IBM001' if INCLASS = anything else
IN2SEQ......The dirty beam image sequence number. 0 => actual
INSEQIN2DISK.....The dirty beam image disk drive no. 0 => any
OUTNAME.....The clean map image name. blank => Standard
defaults based on INNAME.
OUTCLASS....The clean map image class. Standard behavior with
default = 'QCL001' if INCLASS = 'QIM001'
'UCL001' if INCLASS = 'UIM001'
'VCL001' if INCLASS = 'VIM001'
'RCL001' if INCLASS = 'RIM001'
'LCL001' if INCLASS = 'LIM001'
'ICL001' if INCLASS = anything else
OUTSEQ......The clean map image seq. no., 0=> highest unique
If >0; image will be created if new,
overwritten if image name exists.
OUTDISK.....Clean Map disk drive no., 0=> highest with space
INVERS......The CC file version number. 0=> highest
current one. -1 => create new one. If none
exist, one will be created in any case.
GAIN........The CLEAN loop gain. 0 => 0.10
Can be reset with TELL.
FLUX........Stop CLEAN when abs(resid. image max in CLEAN
window) < FLUX (Jy). Can be reset with TELL.
PHAT........The offset added to the pixel at the beam peak.
Should be about 0.1 to 0.5
NCOUNT......The CLEAN method is used until the number of
pixels included in a SDI iteration exceeds
NCOUNT. 0 => do not use the SDI method.
Can be reset with TELL.
CUTOFF......All pixels > CUTOFF times the current peak
residual are considered components. If the 2nd
highest point in the beam is B2, CUTOFF must be
> B2. 0 => 1.01 * B2.
STFACTOR....The loop gain (fraction of each pixel's flux taken
as a component flux) during the SDI phase of the
deconvolution. It must be <= the default.
0 => lower of 0.9 * (1-CUTOFF) and GAIN.
Can be reset to a lower value with TELL.
NITER.......CLEAN iteration limit. 0 => 200
If <0 then the CLEAN will stop at the first
negative component or ABS (NITER) whichever is
encountered first. Can be reset with TELL.
BITER.......Begin new CLEAN at iteration BITER.
If >0, you must completely specify the
OUTNAME, OUTCLASS and OUTSEQ of the clean map.
BMAJ........The FWHM (asec) major axis of the restoring beam.
If 0; value obtained from dirty image.
If <0; output will contain the residual image.
Can be reset with TELL.
BMIN........The FWHM (asec) minor axis of the restoring beam.
Can be reset with TELL.
BPA.........The position angle in the unrotated image of BMAJ.
Can be reset with TELL.
NBOXES......Number (<=50) of rectangular search boxes. 0 => 1
CLBOX.......A 4x50 array with the BLC and TRC of each box.
0 => use inner quarter of input map (first box)
0 => ignore box (boxes 2 - 50)
FACTOR......FACTOR>0 causes deeper CLEAN in each major cycle,
speeding CLEAN, maybe "eating" extended structure.
FACTOR=0 => the normal Clark CLEAN. FACTOR=-0.3 is
good for deep CLEANs of extended structure.
Can be reset with TELL.
MINPATCH....Minimum half width of the portion of the beam
which is used in the AP minor CLEAN. < 3 => 127.
Use 51 for deep CLEANs of extended sources.
Can be reset with TELL. Upper limit 1024.
MAXPIXEL....The maximum number of pixels that are searched for
components inside the ``AP'' in each major cycle.
< 1000 => 20050. This number affects the cpu
usage significantly. Too many causes the task to
search over many points it will never use. Too
few causes the task to do many more small major
cycles, also at great expense. Use this with
caution, but big wins are possible using larger
sizes on very large cleans. Can be changed by
TELL.
DOTV........Display residuals on TV channel 1. >= 0 => yes.
If true, you may stop the CLEANing with TV button
D after each residual map is displayed.
Can be reset with TELL.
BADDISK.....This array contains the numbers of disks on which
it is desired that scratch files not be located.
BADDISK has no effect on input and output maps.
EXPLAIN SECTION
SDCLN: Task which cleans an image using the BGC-AP algorithm
RELATED PROGRAMS: UVMAP,APCLN,IMAGR,CALIB,UVSUB,CCMOD,PRTCC
PURPOSE
SDCLN deconvolves a dirty beam from a dirty map image using
the CLEAN algorithm [Hogbom 1974] as modified to take advantage
of the Array Processor [Clark 1980]. At some stage, specified
by the user, SDCLN will switch to a Steer-Dewdney-Ito (SDI)
[1984] method for finding the model components.
CLEAN iteratively constructs discrete approximants F(n) to
a solution F of the convolution equation:
B^F = D (1)
where D denotes the discrete representation of the dirty map and
B of the dirty beam, the symbol ^ here denoting convolution. The
initial approximant is F(0)=0 everywhere. At the n'th iteration,
CLEAN searches for the extremum of the residual map R determined
at the (n-1)'th iteration:
R(n-1) = D - B^F(n-1) (2)
A delta-function "CLEAN component", centered at this extremum,
and of amplitude g (the loop GAIN) times its value, is added to
F(n-1) to yield F(n). The search over R is restricted to an
area A called the "CLEAN window". A is specified as a number
NBOXES of rectangular sub-areas, or CLBOXes.
Iterations continue until either the number of iterations
n reaches a preset limit N (=ABS(NITER)), or the absolute value
of the extremum in the CLEAN window of the residual map
decreases to a preset value FLUX. If NITER is negative then the
CLEAN will stop at the first negative component. This option
should be used with caution and never be used on the last CLEAN
of an image. CLEAN uses negative components to correct when
too much flux has been subtracted from a pixel and does not
necessarily mean that the noise has been reached.
To diminish any spurious high spatial frequency features in
the solution, F(N) is normally convolved with a "hypothetical"
Gaussian "Clean Beam" H to construct a final "Clean Map" C:
C = H^F(N) + R(N) (3)
The clean beam H may be specified by the user through the
parameters BMAJ, BMIN, BPA, or it may be defaulted to an
elliptical Gaussian fitted to the central region of the dirty
beam B. SDCLN writes the array of "Clean Components" F(N) to
the CC extension of the clean map image file.
The Clark algorithm speeds up the deconvolution process by
splitting it into "major" and "minor" iteration cycles. At the
beginning of the m'th major cycle, it loads into the AP a
RESTRICTED residual map R'(m) containing only the LARGEST
(positive and negative) values in the current residual map R(m).
It then performs a "minor" cycle of iterations wherein new CLEAN
components are sought with (a) the restricted residual map R'(m)
replacing the full residual map R and (b) the dirty beam B being
approximated by its values inside a small area (the "beam
PATCH") with zeroes outside.
A minor cycle is terminated at iteration n' when the peak
in the restricted residual map R'(n') falls to a given multiple
[Clark 1980] of the largest value that was ignored in R(m) when
R'(m) was passed to the the AP. At the end of the cycle of
minor iterations, the current clean component list F(n') is
Fourier transformed, multiplied by the transform of the beam B,
transformed back into the map plane, and then subtracted from
the dirty map D, thereby performing step (2) EXACTLY with the
components list F(n') obtained at the end of the minor cycle.
Errors introduced in the minor cycle through use of the
restricted beam patch are corrected to some extent at this step.
This ends the m'th major cycle, the (m+1)th beginning when the
new restricted residual map R'(m+1) is loaded into the AP.
CLEANing ends (with the transform steps used at the end of a
major cycle) when either the total number of minor iterations
reaches NITER, or the residual value being CLEANed at a minor
iteration reaches FLUX.
The SDI method of finding components differs from the
CLEAN method described above. When there are many pixels in
the residual map of intensity nearly equal to that of the
maximum residual, it no longer makes sense to select just one
pixel at a time. The SDI method selects instead all pixels
having residual intensity greater than some fraction of the peak
residual. For this method to work correctly on a point source
located at a pixel, this fraction must be greater than the 2nd
highest point in the beam and SDCLN forces the fraction (set by
the user as CUTOFF) to be in the range Bp2-0.98. Once a pixel
has been selected by this method, it should in most cases be
included in every major cycle. For this reason, the loop gain
must be < (1-CUTOFF). The actual loop gain at any major cycle
is smaller than this to account for the convolution of the other
selected pixels onto the pixel having the maximum value. In
SDCLN, each selected pixel is regarded as a clean component and
written to the clean component file. This selection is the
equivalent of the minor cycle iterations described above. SDCLN
then removes the components found in the minor cycles by the
same method used for CLEAN. Major cycles continue until the
number of components found is greater than or equal abs(NITER).
SDI components are restored to the output image by the same
convolution used for the CLEAN components. Because the same
pixels will appear many times in the output components file,
this file is sorted, merged and resorted to combine all
components at the same pixels. This much smaller CC file will
speed later processing by CALIB, UVSUB, and all other tasks which
must use the CC model data.
An alternative Maximum entropy deconvolution program is
available in AIPS as task VM but it is experimental and
slow.
COMMENTS
Don't Run SDCLN ...
===================
If there is so little signal in your map that no sidelobes
of any source in it exceed the thermal noise, then no sidelobe
deconvolution is necessary, and running SDCLN is a waste of your
time and of CPU cycles.
General - You can help SDCLN when you run UVMAP
===============================================
Other things being equal, the accuracy of the deconvolution
process is greatest when the shape of the dirty beam is well
sampled. When mapping complicated fields, it is often necessary
to compromise between cell size and field of view; if you are
going to CLEAN a map image, you should set up your parameters in
UVMAP so that there will be at least three or four cells across
the main lobe of the dirty beam.
It is also important to make the map size large enough that
no strong sources whose sidelobes will affect your map have been
aliased by the FFT (preventing SDCLN from handling them
correctly). Consider making a strongly tapered map of a wide
field around your source at low resolution to diagnose confusion
before running SDCLN on a high resolution map (especially when
processing snapshot data from the lower VLA frequencies).
It is helpful to regard CLEAN as an attempt to interpolate
missing samples in the (u,v) plane. The accuracy of the
interpolation is greatest where the original sampling is dense
or where the visibility function varies slowly. The accuracy is
least where you ask CLEAN to EXTRApolate into poorly sampled or
unsampled regions of the (u,v) plane where the visibility
function changes rapidly.
One such region is the center of the (u,v) plane in any VLA
map made from data where all of the fringe visibilities were
less than the integrated flux density of the source. You can
help CLEAN to guess what may have happened in the center of the
(u,v) plane (and thus to guess what the more extended structure
on your map should look like) by including a zero-spacing flux
density when you make your map with UVMAP. This gives CLEAN
a datum to "aim at" in the center of the (u,v) plane. Extended
structure can often be reconstructed well by deep CLEANing when
the zero-spacing flux density is between 100 percent and 125 percent of the
average visibility amplitude at the shortest spacings. If your
data do not meet this criterion, there may be no RELIABLE way
for you to reconstruct the more extended structure. (Some cases
with higher ratios of zero-spacing flux density to maximum
visibility amplitude can be successfully CLEANed, but
success is difficult to predict). If you see an increase in the
visibility amplitudes on the few shortest baselines in your
data, but not to near the integrated flux density, you may
get better maps of the FINE structure from SDCLN by excluding
these innermost baselines when you run UVMAP.
Another unsampled region lurks in the outer (u,v) plane in
many VLA maps of sources at declinations south of +50, if the
source has complicated fine structure. To see why, consult the
plots of (u,v) coverage for the VLA in Section 4 of the "Green
Book" [Hjellming 1982]. At lower declinations, some sectors of
the outer (u,v) plane are left poorly sampled, or unsampled,
even by "full synthesis" mapping. (There are missing sectors in
the outer (u,v) plane in ANY snapshot map). If the visibility
function of your source has much structure in the unsampled
sectors, CLEAN may work poorly on a high resolution map unless
it gets good "clues" about the source structure from the
well-sampled domain. If the clues are weak, badly extrapolated
visibilities in the unsampled regions can cause high frequency
ripples on the CLEAN map. In such cases, CLEAN may give maps
with better dynamic range if you are not too resolution-greedy,
and restrict your data to the well-sampled "core" of the (u,v)
plane when you run UVMAP.
Before applying CLEAN, examine your (u,v) coverage and
think whether you will be asking the algorithm to guess what
happened in such unsampled regions.
Frailties, Foibles and Follies
==============================
There are excellent discussions of CLEAN's built-in
idiosyncracies by Schwarz (1978, 1979), by Clark (1982) and by
Cornwell (1982).
Another way of looking at CLEAN is to think of it as
attempting to answer the question "What is the distribution of
amplitudes at the CLEAN component positions [F(N)] which best
fits the visibility data, if we define the sky to be blank
everywhere else ?" The algorithm can then be thought of as
a "search" for places where F should be non-zero, and an
adjustment of the intensities in F(N) to obtain the "best"
agreement with the data.
The reconvolution of F(N) with the hypothetical "clean
beam" H produces a "clean map" C whose transform is no longer a
"best fit" to the data (due to differences between the
transforms of H and of the dirty beam B). The merit of the
reconvolution is that it produces maps whose noise properties
are pleasing to the eye. It may also be used to "cover up"
instabilities in CLEAN stemming from poor extrapolation into the
unsampled regions of the (u,v) plane, by making H significantly
wider than the main lobe of B.
Note also that step (3) of the standard CLEAN combines this
reconvolution with the residual map, which contains faint sky
features convolved with the DIRTY beam B. If there is
significant signal remaining in the residual map, the effective
resolution of the Clean Map C varies with brightness. You must
therefore be particularly careful when comparing Clean maps made
at different frequencies or in different polarizations; you
should CLEAN all such maps sufficiently deeply that the
integrated flux density in the CLEAN components F(N) is much
greater than that in the residual map R(N).
A recurrent question about CLEAN concerns the uniqueness of
the Clean Map. In the absence of noise, CLEAN could adjust the
amplitudes of the components in F(N) to minimize the rms
difference between the observed visibility function and the
transform of F(N) [Schwarz 1978, 1979]. If the number of
degrees of freedom in F(N) is less than the number in the data,
CLEAN can (and in many practical cases does) converge on a
solution that is sensibly independent of your input parameters.
Noise and approximations in the algorithms complicate this
[Cornwell 1982], but realise that the solution CANNOT be unique
if the number of positions at which you end up with CLEAN
components exceeds the number of independent (u,v) data points.
Be suspicious if your Clean Map contains structures which
resemble those of the dirty beam. This may mean either that you
have not CLEANed deeply enough, or that SDCLN has had
difficulty in some unsampled sector of the (u,v) plane in your
case. This test is particularly important in the case of
snapshot maps, for which the sidelobes of the dirty beam have a
pronounced star (snowflake) structure.
GAIN and NITER
The depth to which SDCLN carries out its deconvolution is
approximately measured by the product NITER*GAIN. The value of
NITER and the execution time needed to reach a given CLEANing
depth are minimized by setting GAIN = 1.0, but setting
GAIN > 0.5 is recommended only when removing the sidelobes of
a single bright unresolved component from surrounding fainter
structure.
When CLEANing diffuse emission, GAIN = 0.1 (the
default) will be much better, for the following reason. The
search step of the algorithm begins its work at the highest peak
(which in an extended source may be a random noise "spike").
After one iteration, the spike is replaced in R by a negative
beamshape, so the next highest peaks are more likely to be found
where the spike would have had its biggest negative sidelobes
[see the diagram on p.11 of Clark (1982)]. If GAIN is high,
subsequent iterations tend to populate F(n) around the negative
sidelobe regions of the highest peaks. This "feedback" can be
turned off by making GAIN small enough that the negative
sidelobes of the first peaks examined in an extended structure
are lost in the noise, i.e. GAIN * (worst negative sidelobe
level) < signal-to-noise on the extended structure. In practice
setting GAIN << 0.1 makes SDCLN unacceptably slow (NITER too
large for a given CLEANing depth) so a compromise is needed.
GAINs in the range 0.1 to 0.25 are most commonly used.
If the source has some very bright compact features
embedded in weaker diffuse emission, it is illuminating to
examine the Clean Map when only the brightest structure has been
CLEANed, to check whether subsequent CLEANing of weaker diffuse
emission causes it to "go lumpy" via the sidelobe feedback
effect. This could be done with GAIN = 0.3-0.5, using either
NITER or the CLBOX selection to ensure that the search does not
stray into the extended emission. Then SDCLN can be restarted
with lower GAIN, higher NITER and wider CLBOXes to tackle the
diffuse structure. If the weak emission "goes lumpy" you may
need to rerun SDCLN with different combinations of GAIN and
NITER to find the most effective one for your particular case.
Ultimately you will stop SDCLN when the new CLEAN component
intensities approach the noise level on your map. On a map of
Stokes I, the appropriate stopping point will be indicated by
comparable numbers of positive and negative components appearing
in the CC list. On maps of Stokes Q and U, which can and will
be legitimately negative, you need to know the expected
sensitivity level of your observation to judge how deep to go.
It is NEVER worth increasing NITER and restarting SDCLN
once many negative CLEAN components appear in the CC list of an
I map. When this occurs, you are using SDCLN to shuffle the
noise in the residual map, which is about as sensible as
rearranging the deck chairs on the Titanic after it hit the
iceberg.
FLUX
This provides an alternative to NITER for terminating the
iterations in a given run of SDCLN. In practice, most users
prefer to control SDCLN by limiting each execution with NITER.
FLUX should then be set to your expected rms noise level times
the dynamic range of the dirty beam (peak/worst sidelobe), to
ensure that you do not inadvertently waste time iterating to a
high value of NITER while in fact CLEANing emission whose
sidelobes are lost in the noise. The FLUX criterion is applied
only to the residuals inside the CLEAN window.
PHAT
When dealing with two-dimensional extended structures,
CLEAN can produce artifacts in the form of low-level high
frequency stripes running through the brighter structure. These
stripes derive from poor interpolations into unsampled or poorly
sampled regions of the (u,v) plane. [When dealing with quasi
one-dimensional sources (jets), the artifacts resemble knots
(which may not be so readily recognized as spurious)].
SDCLN can invoke a modification of CLEAN that is intended to
bias it towards generating smoother solutions to the
deconvolution problem while preserving the requirement that the
transform of the CLEAN components list fits the data. The
mechanism for introducing this bias is the addition to the dirty
beam of a delta-function (or "spike") of small amplitude
(PHAT) while searching for the CLEAN components. [The beam
used for the deconvolution thereby resembles the helmet worn by
German military officers in World War I, hence the name
"Prussian Helmet Clean"]. The theory underlying the algorithm
is given by Cornwell (1982, 1983), where it is described as the
Smoothness Stabilized CLEAN (SSC).
PHAT has given good results in some cases (e.g. Sag A at
20cm - see Cornwell 1983), reducing the amplitude of fringelike
artifacts and producing more believable images.
If you think SDCLN has left CLEAN artifacts in what should
be a smooth source, you have nothing to lose, and possibly
something to gain, by using PHAT for comparison purposes. Tim
Cornwell recommends that lower loop gains than normal, and
deeper CLEANing than normal, be used when using PHAT. Use
GAIN = 0.05 where you would have used 0.1 in the regular SDCLN,
for example. He also notes that his original recommendation for
PHAT = **2/2*