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.

HSI_CLEAN has serval input parameters. the default values can be seen with #IDL> o->print, /clean, /control #we store this structure for later use: ;IDL> cp_flare1 = o->get( /clean, /control )

IDL> help, hsi_clean_parameters(),/str

once HSI_CLEAN has run, the clean parameters can be seen with

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

we store this structure for later use:

IDL> cp_flare1 = hsi_clean_parameters()

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

We can set these parameters with

IDL> o -> set,image_algorithm='clean'
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. if we are planning to reuse the same clean parameters later again, we can save the changes in the structure cp_flare1:

IDL> cp_flare1.niter=1
IDL> cp_flare1.frac=0.3
IDL> cp_flare1.show_maps=1

OK, now let's run HSI_CLEAN.

1. ITERATION

run HSI_CLEAN

IDL> cleaned_map = o -> getdata(image_alg='clean',_extra=cp_flare1)

The input parameters are given with the _extra keyword. if the input parameters are set with o->set,/xxxxx, the _extra keyword is not needed. the input parameters can also be given with cleaned_map = o -> getdata(image_alg='clean',clean_iter=33) etc.

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.

  • 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> cp_flare1.more_iter=1
IDL> cp_flare1.niter=50
IDL> cleaned_map = o -> getdata(_extra=cp_flare1)

JANUARY 8, 2002: bug alert! _extra keyword not working. andre is (luckely) working on it. for the moment, please set the keyword with o->set

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> *cp_flare1.clean_box=cb
IDL> cp_flare1.more_iter=0
IDL> cleaned_map = o -> getdata(image_alg='clean',_extra=cp_flare1)

BUG ALERT (see above): use o->set, .....
IDL> o -> set,clean_box=cb
IDL> o -> set,clean_more_iter=0
IDL> cleaned_map = o -> getdata()

BUG ALERT january 2002: mark_box not working .....

Or 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 slected by coursor:

IDL> cp_flare1.mark_box=1
IDL> cp_flare1.more_iter=0
IDL> cleaned_map = o -> getdata(image_alg='clean',_extra=cp_flare1)


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

BUG ALERT january 2002: mark_box not working .....

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