|
Other RHESSI |
 |
 |
 |
 |
 |
 |
HESSI Forward-Fitting Tutorial
This tutorial is intended for HESSI users who want to use the forward-fitting 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). We start with the standard example of simulated data that is generated in the HESSI software. We explain various steps how the forward-fitting tool can be used to reconstruct and improve an image. The resulting numbers serve as a rough guide, but may vary in different runs, because the data simulator contains a random generator and Poisson noise.First we initiate an object reference:
o = hsi_image()
The minimum call to reconstruct an image with the forward-fitting algorithm and to display the image isim = o -> getdata(image_algorithm='forwardfit',/plot)
The resulting image shows a gaussian source at location [600,200] arcsecs offset from the Sun, reconstructed with detectors 4-8 and in the energy range 12-25 keV. Let us check the image control parameters, which we extract into the structure p,p = o -> get()
and check parameters that are relevant for optimization of our image processing: print,p.det_index_mask
print,p.energy_band
print,p.time_range
print,p.xyoffset
print,p.image_dim
print,p.rmap_dim
print,p.pixel_size
print,p.sim_pixel_size
print,p.sim_xyoffset
print,p.sim_model_ptr(0)
print,p.sim_model_ptr(1)
The simulation model parameter SIM_MODEL_PTR
is a structure that contains the amplitude, the x and y gaussian
widths, the x and y position, and the tilt angle. In our standard example we see that the
image contains two gaussian sources with equal amplitude (1.0), elliptical shapes (x-width=0.5,
y-width=1.0), one source at position [x,y]=[-11.5",-11.5"], and the second source at position
[x,y]=[+0.5",+0.5"], relative to the xyoffset=[600",200"].
Given the information about the simulated image we know that 2 gaussian components should be used for the forward-fitting model. In real data, where the true source structure is not known, we would use either the chi-square information from a preliminary forward-fitting run or another quick imaging method (e.g. back-projection) to guess a suitable number of source components. A perfect model would yield a chi-square in the range of chi2=0.95...1.05. If the resulting chi2 is higher, it is justified to increase the number of gaussian components. In the run above we find, e.g. a value of chi2=1.085, which suggests to increase the number of source components from 1 to 2. So, we change the default value n_gaussians=1 explicitely to n_gaussians=2 in the keyword:
im = o -> getdata(image_algorithm='forwardfit',n_gaussians=2,/plot)
In an arbitrary test run we find (and you should find similar values) the following results for 2 Gaussian components (see the last two lines on screen display, each line contains CHI, CHI-decrement, CEOFF (amplitude, width, x-position, y-position) for one Gaussian component: CHI=1.074, gaussian component #1: amplitude=1.0, width=4.28", x,y=[2.31", 4.15"], gaussian component #2: amplitude=0.43, width=1.09, x,y=[-8.95,-10.04]. Comparing these solutions with the simulation parameters we find agreement in the location of the source components and widhts within ~2", which is satisfactory, given the pixel-size of 4" used in the map grid.The accuracy can be improved for a better signal-to-noise ratio. Because in the standard default setting only the detectors #4-8 are used, we might increase the photon statistics somewhat be including all detectors #1-9 in the forward-fitting. The algorithm tests the ratio between the expected modulation depth and the Poisson noise statistics in the time profiles of each detector and decides by itself which noisy detectors to suppress. So, we enable all detectors:
o->set,det_index_mask=bytarr(9)+1B
im = o -> getdata(image_algorithm='forwardfit',n_gaussians=2,/plot)
The display shows that detectors #2-9 have been found with significant modulation in
the time profiles. The repeated run yields the following results:
CHI=0.970, gaussian component #1: amplitude=1.0, width=2.60",
x,y=[1.61", 2.72"], gaussian component #2: amplitude=0.32, width=1.17, x,y=[-9.89,-9.92].
Another improvement could be achieved by including photons over a broader energy range, assuming that the centroid of the source location does not change significantly.
o->set,energy_band=[6.,300.]
o->set,det_index_mask=bytarr(9)+1B
im = o -> getdata(image_algorithm='forwardfit',n_gaussians=2,/plot)
In order to obtain some quantitative feeling about the accuracy and fidelity
of the image reconstruction, the steps above can be repeated with different
random representations (just initiate a new object o=hsi_image(), and the subsequent
call GETDATA will generate a new image with different random noise).
If the source area is not too big, the size of the reconstructed image can be reduced, either by reducing the image dimension and/or the pixel size, e.g.
o->set,image_dim=[32,32]
o->set,pixel_size=[2.,2.]
o->set,energy_band=[6.,300.]
o->set,det_index_mask=bytarr(9)+1B
im = o -> getdata(image_algorithm='forwardfit',n_gaussians=2,/plot)
For this variant we find the most satisfactory agreement with the simulated data,
within a fraction of an arcsecond:
CHI=1.06, gaussian component #1: amplitude=1.0, width=0.97",
x,y=[-11.38,-10.72], gaussian component #2: amplitude=0.95, width=1.01, x,y=[0.71,1.80].
Note that the computation times for all these runs are of order of 1 minute on workstations, and up to factor of 3 faster on modern laptops.
Further insights into the fidelity of image reconstruction can be obtained by reconstructing the same image with other algorithms, e.g.
