Hessi CLEAN tutotrial

This tutorial is intended for HESSI users who want to use the CLEAN method in the IDL SSW environment to reconstruct HESSI images. For basic information on image processing see tutorial written by Chris Johns-Krull ( http://hessi.ssl.berkeley.edu/~cmj/hessi/doc.html). Please, study first Chris Johns-Krull's tutorial.

How HSI_CLEAN works, can be be most easily explained by running the program.

The Standard Example: two gaussian sources

First, we are making a backprojection image:

IDL> o = hsi_image()
IDL> o -> set,image_algorithm='bproj'
IDL> dirty_map = o -> getdata()
IDL> loadct , 5
IDL> plot_image , dirty_map

Adapting terms from radio astronomy, we call the backprojection map the 'DIRTY MAP'. The dirty map is the map we want to clean.

let's set the image algorithm to HSI_CLEAN:

IDL> o -> set,image_algorithm='clean'

HSI_CLEAN has serval input parameters. the default values can be seen with

IDL> help, hsi_clean_parameters(),/str

once image_algorithm has set, the clean parameters can be seen with

IDL> o->print, /clean, /control

if image_algorithm was never set to 'clean', the clean parameters are not set, i.e. there are -1

information about the clean parameters can also be found at http://hesperia.gsfc.nasa.gov/ssw/hessi/doc/hsi_params_clean.htm

to save the clean parameters for a possible later use:

IDL> cp_flare1 = o->get( /clean, /control )

In the beginning, the following three input parameters are of main interested:

NITER: gives the MAXIMAL number of iteration HSI_CLEAN will perform.

FRAC: in each iteration, only the value of FRAC of the maximum in the map is assumed to be real. In radio astronomy, this is called the GAIN.

SHOW_MAPS: the progress of the clean process is displayed, if SHOW_MAPS is set.

We can set these parameters with

IDL> o -> set,clean_niter=1
IDL> o -> set,clean_frac=0.3
IDL> o -> set,clean_show_maps=1

REMARK: print,o->get(/clean_niter) will show the value in clean_niter. however, if image_alg is not set to 'clean' the clean paramters cannot be seen. once clean was called, it will work.

we will run HSI_CLEAN first for only one iteration (so that we can see what is going on) with a relatively large gain of 0.3. OK, now let's run HSI_CLEAN.

1. ITERATION

run HSI_CLEAN

IDL> cleaned_map = o -> getdata()

alternatively to using o->set, the input parameters can also be given with cleaned_map = o -> getdata(image_alg='clean',clean_iter=33, ....) etc.

REMARK: if no keyword is changed, getdata will not run again. to force to rerun the calculation use cleaned_map = o -> getdata(/force)

Since the SHOW_MAPS keyword is set, a window opens showing the progress of clean (your simulation looks different, and the default parameters have too. sorry for that).

in the top left corner the input image is display, the DIRTY MAP.
in the each iteration, the maximum of the dirty MAP is searched first. the location of the maximum is displayed in the top right corner ('SELECTED MAXIMA') and in the bottom right corner ('CLEAN_COMPONENTS'). the clean components additionally show the strength of the maxima times the gain (gain=FRAC).
in the next step, the point spread function (PSF) is calculated at the location of the maximum.
the PSF multiplied by the gain is then subtracted from the dirty map. the resulting map is called 'RESIDUAL MAP' and is displayed on the top second from the left.
then the PSF is approximated by a gaussian shaped function (CLEAN PSF or CLEAN BEAM) supressing side lobes and the wings.
the 'clean components' is then convolved with the CLEAN PSF and displayed in the top second from right. this is the 'COMPONENT MAP'.
the sum of the component map and the residual map give the 'CLEANED MAP'.
this is the end of the first iteration.
since we set NITER to 1, only one iteration was performed and HSI_CLEAN stoped.

some more information about the displayed maps:

The dirty map, residual map (top), component map (top), and cleaned map are all shown with the same scaling, whereas the residual map and component map shown in the bottom panels are each scaled from minimum to maximum.

images are in polar coordinates

the units of the cleaned map are photons/subcollimator/pixel. so we can sum over the flare area and get photons/subcollimator.

on the right, some information about the clean process are plotted. from top to bottom:

current iteration
maximum in residual map
ratio of the maximum in the cleaned map and the standard deviation of the residual map
ratio of the maximum in the cleaned map and the maximum of the residual map
the sum of all clean components
gain, e.g. the keyword FRAC

MORE ITERATIONS

let's run some more iterations. We can set the keyword 'MORE_ITER' to tell HSI_CLEAN to continue cleaning (so it starts with iteration 2, since iteration 1 is already done). Additionally, we increase the number of iterations, NITER, to 50:

IDL> o -> set,clean_more_iter=1
IDL> o -> set,clean_niter=50
IDL> cleaned_map = o -> getdata()

HSI_CLEAN was now running for 7 iterations (could be different in the simulation you are running). it stoped because the extremum in the residual map is NEGATIVE, i.e. the residual map contains bigger negative values ('noise') than positive values (signal). So it makes sense to stop cleaning. However, we can still see that there is some remaining flux left. the flare source can be still seen in the residual map.

we can now help clean by tell it where we expect the source to be. we are doing that by defining so-called 'CLEAN BOXES' or 'CLEAN WINDOWS'.

CLEAN WINDOWS

At the moment, only rectangular clean windows can be set (March 2001). Clean windows can be set directly or by clicking with the mouse.

if one knows the coordinates (arcseconds from solar center, o->plot, makes a lot of the current images including the axis in arcsec from solar center), one can set clean boxes the following way (if we want to run clean again from iteration 1, we additionally have to set back the keyword MORE_ITER):

IDL> nbox=2
IDL> cb=fltarr(4,nbox)
IDL> cb(*,0)=[584.,183.,599.,199.]
IDL> cb(*,1)=[595.,195.,606.,207.]
IDL> o -> set,clean_box=cb
IDL> o -> set,clean_more_iter=0
IDL> cleaned_map = o -> getdata()

Or, much more convenient, one can set the keyword MARK_BOX. HSI_CLEAN then displays the dirty map and a hopefully selfexplaining menu. an unlimited number of clean boxes can be selected by coursor:

IDL> o -> set,clean_mark_box=1
IDL> o -> set,clean_more_iter=0
IDL> cleaned_map = o -> getdata()

Clean was now running for 17 iterations. the residual map does not show any remining flux within the clean boxes. we can even see that their are larger maxima outside the clean boxes indicating that we have maybe even cleaned a little bit too much.

the residual map gives us a measure for the quality of the component map. it gives us error estimates to the component map.

HOW TO ACCESS THE DIFFERENT MAPS

all results of HSI_CLEAN (residual map, clean components, etc) can be accessed with

IDL> help, o -> get(class_name='hsi_clean',/info),/str

all the maps are in polar coordinates.

the clean components for each iteration (pixel number (1dimensional), x-coordinate, y-coordinate, maximum) for each iteration are saved in

IDL> print,o->get(class_name='hsi_clean',/info,/clean_components)

PLAY AROUND

play around with the following parameters:

use more photons for the simulations by setting

o->set,sim_photon=5000.

the default value was 1000. HSI_CLEAN gives much nicer results for images with more photons.

try out to clean images taken with different detectors:

o->set,det_index_mask=[0,0,0,0,0,1,0,0,0]

if you use small grids only, change the pixel size too.

clean more carefully by using a smaller gain, e.g of 0.1 or even less.

try to clean a spatial resolved source. make a wide gauss'schen source and try to clean it. using a low gain gives best results.

try out to clean a single point source and compare with the point spread function. the PSF can be calculated with the following routine:

psf = o ->GetData( CLASS_NAME = 'HSI_PSF', /SUM, xy_pixel=xy )

where xy is the 1 dimensional pixel coordinate, e.g. center pixel in 65*67 image is 33*65+32=2177 (PSF is in polar coordinates). PSF at different location in the images are slightly different. most of the calculation time in HSI_CLEAN is used to calcuated PSFs.

create a more fancy model source and see what clean can do.

QUESTIONS/COMMENTS/ETC. to Säm Krucker, krucker@ssl.berkeley.edu, 510 643 3101