 News |
FAQ
|
| Search
|
Home
|
| Getting Started | Documentation | Glish | Learn More | Programming | Contact Us |
|
| Version 1.6 Build 363 |
|
Tim Cornwell and Kumar Golap
This chapter describes how to make and deconvolve images starting from a calibrated MeasurementSet. There is also some discussion of how to do self-calibration.
The major AIPS++ tools that you will need are:
In addition, you will need to use the Viewer to examine your images, and perhaps an MS tool to load and save MeasurementSets.
The User Reference Manual contains the primary documentation for all these tools and their functions.
To make an image, a reasonably well-calibrated MeasurementSet is required. You may either import calibrated data from another system or calibrate the data in AIPS++. See the Calibration chapter for more details on the latter.
If the dataset was calibrated in another system then you'll need to load it from a FITS file.
include 'ms.g'
x := fitstoms('pks1331-0.ms', 'PKS1331-09.UVFITS') # Convert from FITS
x.summary() # Summarise to logger
x.done() # Finished with tool
The tool for imaging and deconvolution is Imager. It has been designed to do everything that is required: gridding/degridding of visibility data, Fourier transforms, deconvolution using various methods, etc. A complete overview of the capabilities of Imager can be found in the Reference Manual.
To use an Imager tool one must first construct it from a MeasurementSet and then configure it (set its ``state''). The steps are:
include 'imager.g'
imgr := imager('/home/data/mydata.ms') # Construct
imgr.summary() # Summarize state of tool
imgr.done() # Destroy
The summary function is always useful to summarize the current state of your Imager tool. You may find you have forgotten what functions you ran and need reminding what you state you have set.
Here are some examples (refer to the Reference Manual for more details).
imgr.setdata(mode='channel', nchan=15, # Select 15 channels, second
spwid=[2], fieldid=3) # spectral window and third field
#
imgr.setdata(mode='channel', nchan=256, step=2) # Select every second channel to make 256
#
imgr.setdata(mode='velocity', nchan=64, mstart='20km/s', # Select data in radio velocity space
mstep='-100m/s');
#
imgr.setdata(mode='opticalvelocity', nchan=64, # Select data in optical velocity space
mstart='20km/s', mstep='-100m/s');
Note that the channel selected can be further combined in the Imager.setimage function (e.g. to grid blocks of the selected channels together).
Other functions can also further reduce the data selected. For example, Imager.filter and Imager.uvrange.
This function controls things like image size, cell size, spectral and polarization selectionsampling, phasecenter, and number of facets (for wide-field imaging). Images constructed afterwards (using e.g. the Imager.clean function) will have these parameters.
One important point to note is that the image size does NOT have to be a power of two. Just choose some reasonably composite even number close to the number that you want. Do not use a prime number since the Fourier Transform will then be very slow. The combination of functions setimage and setdata determines how spectral averaging is done in making spectral line images (see the example at the end of this document).
include 'measures.g'
imgr.setimage(nx=600, ny=600, cellx='0.5arcsec', celly='0.5arcsec',
fieldid=2, mode='channel')
#
phctr := dm.direction('J2000', '19h50m00', '50d00m00')
imgr.setimage(nx=600, ny=600, cellx='0.5arcsec', celly='0.5arcsec',
doshift=T, phasecenter=phctr, mode='mfs')
In the first example, the phase center is that of the specfied field (source). In addition, each selected channels (in function setdata) will be imaged as one channel of the output image.
In the second example, we specify the phase center explicitly and set doshift=T to indicate that the phase center must be taken from the phasecenter argument (Ed. note - the interface that specifies the phase center needs to be rationalized some).
The above steps show how to set up your Imager tool as desiresd. In addition, before imaging one may wish to set some values in the MeasurementSet itself. This is necessary for visibility weighting. The visibility imaging weights are computed prior to the actual imaging and stored in a column of the MeasurementSet, ``IMAGING_WEIGHT''.
The first time an Imager tool is attached to a MeasurementSet it initializes the imaging weights to the natural weighting scheme. Thereafter, whatever weighting scheme you last used, is the one you will get if you don't explicitly run one of the weighting functions.
The values in IMAGING_WEIGHT column are set, changed, and examined by a number of functions.
imgr.weight(type='natural') # Natural weighting imgr.weight(type='uniform') # Uniform over entire field of view imgr.weight(type='uniform', npixels=300) # Uniform over specified size imgr.weight(type='briggs', rmode='norm', robust=0.5) # An example of Briggs weighting
imgr.filter(type='gaussian', bmaj='2.5arcsec', bmin='2.5arcsec', bpa='60deg')
In Imager, we have chosen separate out the calculation of the weighting from the gridding and Fourier transformation to give more control over this important step. Note that this also gives the ability to set the weights by hand using getcol and putcol functions of the table tool. For example, to replace the weights by their square root:
t:=table('3C273XC1.ms', readonly=F)
wt:=t.getcol('IMAGING_WEIGHT')
wt[wt>0.0] := sqrt(wt[wt>0.0])
t.putcol('IMAGING_WEIGHT')
t.flush();
t.close()
Before performing such a step, one should close the imager and then construct it again afterwards.
WARNING: All the weighting schemes modify the MeasurementSet and thus are conserved even after destroying the Imager tool and may no longer be suitable for subsequent imaging. You can of course reset the weighting scheme with the weight function.
Now that your Imager tool and the MeasurementSet are set up, you can make an image.
Imagine we want to make an image over a field of view of 10 arcmin. The advise function will return the maximum pixel size required and number of pixels. If the number of facets required is larger than 1 then one needs a wide-field imaging algorithm as described in the wide-field chapter of Getting Results.
myimager.advise(fieldofview='10arcmin')
You might also like to make the point spread function.
imgr.makeimage(type='corrected', image='dirty.image') imgr.makeimage(type='psf', image='dirty.beam')
At this point, you are ready to deconvolve the point spread function (usually called the dirty beam) from the dirty image. The goal of this step is to correct the dirty image for the sidelobes in the point spread function, which arise from the limited Fourier plane coverage in your observation. You have two choices for the deconvolution: use the Deconvolver tool on the dirty image and point spread function or use the various deconvolution functions of Imager. In the latter, the dirty image and point spread function are recomputed as needed, whereas Deconvolver works with a fixed dirty image and point spread function. If you want to mosaic, or image wide- or multiple fields then you must use Imager.
The options for deconvolving are :
If you use a multi-scale algorithm, you should also set up the scale sizes via the Imager.setscales function.
The Pixon algorithm is available via an IDL library courtesy of the Pixon LLC. This means that you must have the library and IDL installed on your computer. The Pixon library is available free of charge for your personal scientific use direct from Pixon LLC. IDL is available commercially from RSI.
For all deconvolution approaches (except Pixons)a Componentmodel tool may be specified, in which case it is removed from the data prior to imaging. Thus the resulting deconvolved image does not include the flux modelled by the Componentmodel.
The deconvolution methods do not keep a list of CLEAN components in sequence. Instead the CLEAN components are immediately added to the model image. This allows interoperability of the CLEAN, MEM and NNLS algorithms as well the deconvolution functions in the Deconvolver tool. It also allows editing of the resulting images using the Image tool.
The deconvolution functions can return the residual and restored images, as well as the updated model image. In addition, there are Imager.residual and Imager.restore functions that also compute residual and restored images. The distinction is that in the former, the residual is that calculated by the deconvolution method at the end of iteration, whereas the latter returns the residual calculated by going back to the original visibility data (except for the cases of multi-field and wide-field algorithms).
imgr.clean(algorithm='clark', niter=2000, gain=0.2, model='model_2000',
image='restored_2000', residual='residual_2000')
In the above example we run a Clark CLEAN deconvolution with Imager. 2000 CLEAN components with a gain of 0.2 are requested.
If the model image does not pre-exist it will be created and filled with the CLEAN delta function components found. If it does pre-exist (it may be non-empty (perhaps the result of a prior deconvolution) and its Fourier transform is first subracted. This is also how you would continue to CLEAN deeper from a prior run of clean, mem, or nnls.
Note that the model image is updated when this function finishes. So if it was non-empty on input, you should save a copy first, if you wish to preserve that model.
In deconvolution, one must often specify the region of the dirty image over which the sky brightness is expected to be non-zero. In CLEAN, possible components can only come from this region. Designation of regions to deconvolve is done via images of the same shape as the dirty image. The value of these mask images ranges from 0 to 1 (but usually we use 0 [no components] or 1 [allow components]). There are a number of ways of making such mask images.
The most powerful approach is probably to use Imager.regionmask, since regions have a number of advantages: they can be polygonal, they can be combined in various ways, they can be defined in world coordinates, and they can be saved to a table for later retrieval and use. The Regionmanager tool allows these and other operations on regions.
A compound region holding the union of lots of interactively specified simple (boxes or polygons) regions can be made via the Viewer. Please see the documentation in Viewerdisplaypanelviewer:vdpgui to see how to do this.
NNLS requires two such masks images : one for the region where flux is to be allowed and one for the region to be used as constraints. Typically the latter region will be bigger than the former. The product of the pixels in the two regions should not exceed about 5-10 million.
# Imagine you have a Region tool called 'r1' in Glish
imgr.regionmask(mask='mask_image', region=r1) # Make mask image from region
imgr.clean(algorithm='clark', niter=2000, gain=0.2,
model='model_2000', image='restored_2000',
residual='residual_2000', mask='mask_image') # Use it !
For a quick look at an image, one may wish to use the Imagerwizard function. This is built using various tools including Imager. It makes many decisions, such as cell size, automatically, and simplifies obtaining an initial imaging.
include 'imagerwizard.g'
ok := imagerwizard('data.ms')
The wizard leads you step by step through the imaging process via interactive GUIs. It also writes a Glish script which you can then look at, edit and re-run if you wish.
A priori, survey, or model fitting information can provide the user with the basis for a component model of bright, unresolved or somewhat resolved, confusing sources. These sources can be removed from the visibilities prior to imaging mosaicing by setting them in the Imager.clean or Imager.mem functions. Primary beams are accounted for (see below).
However, these components are assumed to be small compared to the primary beam, and only the primary beam value at the center of a Gaussian component is applied to the Gaussian component (ie, a Gaussian component stays Gaussian after application of the primary beam, even if it's wings extend beyond the primary beam pattern). Hence, errors will result if highly extended components are included in the component list. If your application requires a model with extended Gaussian components, for the time being you may want to turn the extended components into an image with the Image tool's Image.modify function.
Component models and primary beams are dealt with as follows:
I wish I knew.
Delta function CLEAN algorithms, such as the Clark or Hogbom algorithms, must use a modest gain factor (traditionally 0.1) in order to image extended emission in a reasonable manner. Otherwise, emissions and sidelobes get confused and error striping can result.
Multi-scale CLEAN algorithms attempt to recognize that the emission in the sky usually comes on a variety of scale sizes; it decomposes the image into Gaussians of the specified scale sizes. This means it does not spend huge amounts of time CLEANing large extended sources pretending that they are really collections of delta functions.
Function Imager.setscales allows you to choose these scales, either with an automatic method (you say how many scales you want) or by direct specification of the scale sizes in pixels.
Multi-scale CLEAN does not suffer from the same problems as delta function CLEAN and can use a larger value of the gain factor, in the range of 0.4 and 0.8. Sometimes an oscillatory pattern will occur in the total cleaned flux with the peak residuals bouncing back and forth between positive and negative. If this occurs, try reducing the gain factor or try reducing the size of the largest scale.
The progress display for multi-scale CLEAN shows the flux and peak residuals for each cleaning scale. For any extended brightness distribution, there will be initially more flux on the large scales than on the small scales, and this will be reflected by much larger residuals on the large scales. Consequently, the flux is built up much more quickly by CLEANing on these large scales. However, in the early stages of deconvolution, we may never get anywhere close to the stated CLEANing threshold if we are dominated by extended flux. The exact manner in which the larger scales' residuals will be inflated depends upon the details of the brightness distribution.
imgr.setscales(nscales) # Specify number of scales
imgr.setimage(nx=600, ny=600, cellx='0.5arcsec',
celly='0.5arcsec', mode='mfs')
imgr.clean(algorithm='msclean', # Deconvolve
model='model',
image='restored',
niter=100, gain=0.3)
The imagermultiscale function provides a different sort of packaging of multi-scale algorithms and works on mosaicing or single field data.
Conceptually, the imagermultiscale() function offers a second order multi-scale method, or doing multi-scale CLEAN as the deconvolution kernel inside of a cascade of multi-resolution CLEANs each working on increasingly larger images with smaller pixels.
The multi-scale CLEAN algorithm addresses an essential problem in imaging: you often have emission on a wide range of scale sizes which must be imaged simultaneously. While it is vastly superior to delta-function CLEAN algorithms at reconstructing extended emission, it is not incredibly efficient as it seeks to deconvolve very large structures with very small pixels.
A simple scripting approach to multi-resolution imaging which overcomes this inefficiency is provided by the imagermultiscale function. The user provides a vector of image sizes and the corresponding pixel sizes, and the function performs a multi-scale CLEAN on the first image (with the largest pixels). When the first CLEAN stops, its model is regridded onto the second image size and is used as a starting model for the next round of Multi-Scale CLEANing. Of course, since the regridded model will have smaller pixels, it will be in error; but the large scale structure will be pretty well in place. So, each round of Multi-Scale CLEANing will concentrate on the small scale structure, though each will also be able to make larger scale adjustments as well.
If the imagermultiscale function is too restrictive for your needs, you are encouraged to pull it out of imager.g and edit your own function. For example, you might want to have some fine control over which channels to image. This function could also easily be converted to a multi-scale MEM function.
include 'imager.g'
#
ok := imagermultiscale(msname='orion.ms', imsizes=[128, 256, 512],
cellsizes=[16,8,4], scales=[0, 5, 15],
nitermult=0.15, niterpower=1.2, fields=[2],
centerfield=[2],
spwid=[1])
Wide-field imaging and mosaicing are covered in separate chapters of this document.
Sometimes, the lack of total-power data in an interferometric observation (single- or multi-field) can result in image defects on the size scale of the interferometer baselines. For example, in a maximum entropy algorithm with positivity constraint, the ``negative bowl'' will manifest itself in the model brightness distribution by a flat region of near zero emission surrounded by a spurious positive ``splash'' around the bowl; the ``splash'' artifact sometimes has small scale structure.
Currently, the only method available in Imager to combine total power measurements with interferometer data is via a feathering technique.
Since feathering is an intrinsically post-deconvolution method, it cannot help such a case. However, there are two ways to help out the deconvolution algorithms. First, convert the single dish image to the same image grid as the interferometric mosaic by using the Image.regrid function. Then either rescale or deconvolve the image into units of Janskys per pixel. This image can be used as a prior image in maximum entropy or maximum emptiness and will prevent small scale image artifacts due to the lack of total power data. However the prior image is not elevated to the level of a hard constraint, so a subsequent pass of the interferometer image and the single dish image through the Imager.feather function will ensure that the large scale structure in the final image is consistent with the single dish image.
imgr.feather(image='feathered.image', # Specify the single dish
highres='casa.vlaonly', # voltage pattern via
lowres='casa.sd', # a table
usedefaultvp=F,
vptable='sd.vptable');
Self-calibration is the process of using a model of the source to improve the calibration of the data. One makes an initial model using perhaps componentmodels or a deconvolved image. This is then used to predict the visibilities, and the antenna gains are derived by a least squares fit between the observed and modelled visibilities. Once the antenna gains thus derived have been corrected, the model may be rederived and thus improved. This process is iterated until a satisfactory improvement has been acheived.
In AIPS++ this process is done with an Imager tool and a Calibrater tool working together.
The steps are:
This process (steps 3 onwards) can be repeated until the image stops improving. Usually it is a good idea to start with phase-only solutions (controlled in Calibrater.setsolve) and then allow amplitude solutions after the phase-only solutions have ceased improving.
The types of solution allowed are:
More details may be found in the Calibrater documentation.
Diagnosis of the success of the self-calibration is ultimately judged by the improvement in image quality (measured from the off-source noise level). One may also expect that the rms errors per time interval reported by Calibrater.solve will decrease, and that the residual visibility as displayed by Imager:plotvis will diminish. In addition, it is also advisable to examine the gain solutions using Calibrater.plotcal. These should be small if the initial calibration was of good quality.
include 'imager.g'
include 'os.g'
#
imgr := imager('mydata.ms') # Construct
#
imgr.setdata(mode='channel', nchan=15, start=1, # Select data (15 channels, 1 field)
fieldid=3)
#
imgr.setimage(nx=512, ny=512, cellx='2arcsec', # Set image parameters
celly='2arcsec', mode='mfs')
#
imgr.clean(algorithm='clark', niter=10000, # CLEAN 10000 components
model='model', image='restored',
residual='residual')
#
dos.copy ('residual', 'residual_10000') # Save copies
dos.copy ('restored', 'restored_10000')
#
imgr.clean(algorithm='clark', niter=5000, # Restart and CLEAN another
model='model', image='restored', # 5000 components. Updates
residual='residual') # images
#
rest1 := image('restored_10000') # Display restored images
rest1.view()
rest2 := image('restored')
rest2.view()
#
resid1 := image('residual_10000') # Compare statistics of
resid1.statistics() # residual images
resid2 := image('residual')
resid2.statistics()
The following script demonstrate of how we can image multiple regions of the sky simultaneously (in AIPS this is implemented in task MX). This is a technique used to deal with a source which is outside the main field of view being imaged but whose sidelobes are seen in the image we want to make. It is also sometimes useful to remove some kinds of interference which tend to appear as a source at the pole.
It is assumed in this script that there is only one fieldid (i.e. one source).
imgr:=imager('mydata.ms') # Construct
#
imgr.setdata(mode='channel', nchan=25, # Select data: 25 channels,
start=4, step=1, spwid=[1,2]) # 2 spectral windows.
#
mydir:=dm.direction('J2000', '19h58m49', '40d45m29') # Direction of outlier (CygnusA)
imgr.setimage(nx=200, ny=200, # Set image parameters for outlier
cellx='19arcsec',celly='19arcsec',
doshift=T, phasecenter=mydir,
mode='mfs', spwid=[1,2])
imgr.make('cygA') # Make empty image for outlier
#
#
imgr.setimage(nx=600, ny=600, # Set image parameters
cellx='19arcsec', celly='19arcsec', # for main field direction
spwid=[1,2])
imgr.make('mainfield') # Make empty image for main field
#
imgr.uvrange(uvmin=10, uvmax=5000) # Select uv range
imgr.weight(type='briggs', rmode='norm', # Weighting
robust=0.0, npixels=300)
imgr.clean(algorithm='mfclark', # Deconvolve
niter=1000, gain=0.05,
model="mainfield cygA", # Vector of strings
image="mainfield.restored cygA.restored",
residual="mainfield.residual cygA.residual")
#
Note that the combination of Imager.setimage and Imager.make can be used as many times as is necessary to define all the outlier fields.
We had to explicitly make the empty model images before passing them to Imager.clean. Otherwise Imager would recognize only the last call to Imager.setimage and would thus fail to handle the multiple fields.
In the running of the Imager.clean function we used the Clark algorithm. In addition, we must specify algorithm='mfclark' because we are making a multi-field image.
The following is an example of making a spectral-line image. We demonstrate how Imager.setdata and Imager.setimage can be used to average data channels to make one of the spectral-line image channels.
The MeasurementSet is created from a FITS file that comes with the AIPS++ installation.
include 'imager.g'
include 'ms.g'
include 'sysinfo.g'
#
aipsroot:=sysinfo().root()
fitsfile:=spaste(aipsroot, '/data/demo/BLLAC.fits')
myms := fitstoms('mydata.ms', fitsfile) # Generate MeasurementSet
myms.summary()
myms.done()
#
imgr:=imager('mydata.ms') # Construct Imager
#
imgr.setdata(mode='channel', nchan=40, start=4, # Select data: 40 channels
step=1, spwid=[1]) # from the first spectral window
#
imgr.setimage(nx=600, ny=600, # Set imaging parameters. We
cellx='19arcsec', celly='19arcsec', # output 20 channels from the 40
doshift=F, spwid=[1], mode='channel', # by combining pairs of channels
nchan=20, step=2 )
#
imgr.weight(type='uniform') # Weighting
#
imgr.clean(algorithm='clark', niter=1000, gain=0.05, # Deconvolve
model='model',
image='restored',
residual='residual')
#
imgr.done()
Note that if the channel images were needed in the velocity domain, then argument mode in both functions Imager.setdata and Imager.setimage must be put to velocity and instead of the start and step parameters, use the mstart and mstep parameters.
In this example we assume we have bright point sources and low level broad scale structures to image. It is well known that CLEAN does not do well with large scale structures and that MEM does not do well with bright compact sources in presence of large scale structures. We can CLEAN on the bright point sources followed by MEM on the low level structures. What follows is an example of how to do that.
imgr:=imager('mydata.ms')
imgr.setimage(nx=600, ny=600, cellx='2arcsec', celly='2arcsec',
doshift=F, spwid=[1])
#
# Use clean with a mask image around the bright point sources
# to restrict clean to deconvolve these sources only. Alternatively,
# if the point sources are brightest, just CLEAN a little bit
# to remove the bulk of them
#
imgr.clean(algorithm='clark', niter=1000, gain=0.1,
model='cleanmodel',
mask='pointsource.mask',
image='clean.restored',
residual='clean.residual')
#
# Now run MEM on the data. Note that we use the model image that we
# have got from clean, so that the point sources will be subtracted
# from the visibilities and MEM will deconvolve the residual
#
imgr.mem(model='cleanmodel', niter=40, sigma='0.001Jy',
image='final.restored', residual='final.residual')
#
imgr.done()
The Deconvolver tool has functions which have very similar meanings to to those in Imager. However, no MeasurementSet is used with the Deconvolver; it is based on image plane deconvolution only. It needs a point spread function and a dirty image.
deco := deconvolver('3C273XC1.dirty', '3C273XC1.psf')
deco.setscales(scalemethod='uservector', # Scales of delta function,
uservector=[0, 3, 10]) # 3 and 10 pixels wide
#
deco.clean(algorithm='msclean', # deconvolve
model='3C273XC1.clean',
niter=100, gain=0.3)
deco.done()
Note that it is not very useful to use CLEAN components that are wider than the largest structure in the image.
The following is a script of a multi step CLEAN and self-calibration. We use here a test MeasurementSet that is in the AIPS++ system.
include 'imager.g';
include 'calibrater.g';
include 'os.g'
#
imagermaketestms(msfile='3C273XC1.ms') # Make MS
imgr:=imager(filename='3C273XC1.ms') # Construct Imager tool
imgr.summary()
imgr.setdata(mode="none" , nchan=1, start=1, step=1, # Data selection
spwid=1, fieldid=1);
imgr.setimage(nx=600, ny=600, cellx='0.6arcsec', # Set image parameters
celly='0.6arcsec', stokes="I" )
imgr.uvrange(uvmin=1.0, uvmax=500000);
#
dos.remove('model0')
dos.remove('restored0')
dos.remove('residual0')
#
imgr.clean(algorithm="hogbom" , niter=1000, # Deconvolve 1000
model="model0" , fixed=F,
image="restored0" , residual="residual0")
#
dos.remove('model0_5000')
dos.remove('restored0_5000')
dos.remove('residual0_5000')
imgr.clean(algorithm="hogbom" , niter=5000, # Deconvolve
model="model0_5000" , fixed=F,
image="restored0_5000" , residual="residual0_5000")
#
cal := calibrater(filename="3C273XC1.ms") # Construct Calibrater tool
cal.setsolve(type="G" , t=100, preavg=0.0, # Set selfcal solution parameters
phaseonly=T, refant=-1, # phaseonly
table="selfcal.gcal",append=F);
cal.setapply(type="G" , t=0.0, table="selfcal.gcal", # Set parameters for application of
select=''); # the calibration
# The next step Fourier transforms the model from the CLEAN
# above and makes model visbilities in the MeasurementSet. The
# previously created Calibrater tool fits the model visibilities with
# the observed ones and hence calculates the necessary gains to make
# the best fit. The solution is applied.
ok := imgr.selfcal(caltool=cal, model="model0", complist='');
cal.plotcal(plottype="PHASE", tablename="selfcal.gcal" , # Inspect solution
antennas=[], fields=[], polarization=1, spwids=[],
timeslot=1, multiplot=T, nx=3, ny=3)
#
dos.remove('model1')
dos.remove('restored1')
dos.remove('residual1')
imgr.clean(algorithm="hogbom" , niter=2000, # Deconvolve 2000
model="model1" , fixed=F,
image="restored1" , residual="residual1")
cal.setsolve(type="G" , t=30, # Setup for phase and amplitude
preavg=0.0, phaseonly=F, refant=-1,
table="selcal.gcal" , append=F);
#
ok := imgr.selfcal(caltool=cal, model="model1") # Solve and apply for a gain selfcal
#
dos.remove('model2')
dos.remove('restored2')
dos.remove('residual2')
imgr.clean(algorithm="hogbom" , niter=5000, # Final deconvolution, 5000
model="model2" , fixed=F,
image="restored2" , residual="residual2")
#
r := drm.quarter()
r0_1000 := image('restored0')
r0_5000 := image('restored0_5000')
r1_2000 := image('restored1')
r2_5000 := image('restored2')
#
r0_1000.view(region=r, includepix=[-.003,.003])
r0_5000.view(region=r, includepix=[-.003,.003])
r1_2000.view(region=r, includepix=[-.003,.003])
r2_5000.view(region=r, includepix=[-.003,.003])