|
Other RHESSI |
 |
 |
 |
 |
 |
 |
HESSI Imaging Overview
Updated 12-April-2000 by cmj@ssl.berkeley.edu
Introduction
This document describes the command line interface to the HESSI imaging software that most general investigators will use. In addition to the documentation given here, additional software documentation is available on the web. In particular, several HESSI data analysis software manuals can be found at:
http://hessi.ssl.berkeley.edu/software/
though much of the information here is mainly of use to programmers writing modules for HESSI data analysis. This manual deals only with HESSI imaging software and assumes you have already installed and are using the Solar SoftWare (SSW) package. Instructions for downloading and installing SSW can be found at the above link.
HESSI Imaging Software
The main command line interface to the HESSI imaging software is the routine hessi_image.pro. At the IDL command line, one can create an image simply by typing:
IDL> hessi_image,im
With this call, the software will simulate HESSI data using a default model and then analyze the data using the HESSI imaging software. The resulting image is returned in the variable im. This image can be plotted with the traditional SSW routine plot_image.pro. The simple call would be:
IDL> plot_image,im
By selecting the keyword plot in the call to hessi_image, the program will automatically display the final image at the end of its processing, along with a summary of the information that went into creating the image. The call would be:
IDL> hessi_image,im,/plot
The above command can also be given without the return variable, im, in which case the resulting image is only displayed to the screen. Real HESSI telemetry data will reside in FITS files. An image can be constructed from such FITS files by specifying a filename in the call to hessi_image. For example, if HESSI telemetry data resides in the file telemetry.fits (in the current working directory), this data can be analyzed with the call:
IDL> hessi_image,im,filename='telemetry.fits'
Examples of simulated HESSI telemetry data files are available on the
web at http://hessi.ssl.berkeley.edu/data/archive.html,
under
the heading Test Data (not currently available). It is
important to note that the HESSI imaging software has memory. Once
invoked, repeated calls to
hessi_image will not do anything new
unless something (parameter, filename, etc.) is explicitly changed.
The current imaging parameters can be returned in an optional second return variable. For example, the call
IDL> hessi_image,im,pars,/plot
will display the final image, return the image in the variable im, and will return the imaging control parameters in the structure pars. This structure can be examined by typing:
IDL> help,/structure,pars
An example of the resulting help display is:
** Structure <2cdb08>, 53 tags, length=80576, refs=1:
IMAGE_ALGORITHM STRING 'BACK PROJECTION'
ALGORITHM_USED STRING 'HSI_MODUL_PATTERN'
VERBOSE
BYTE 1
NEED_UPDATE BYTE
0
DEBUG
BYTE 1
PLOT
BYTE 0
MODPAT_FILE STRING
'modul_pattern_file.dat'
IMAGE_DIM INT
Array[2]
PIXEL_SIZE FLOAT
Array[2]
VRATE
POINTER <NullPointer>
VIMAGE
POINTER <NullPointer>
LONGSCALE BYTE
0
XYPSF
LONG
2048
XYPSF_PTR LONG
2048
MODPAT_SET_MAX INT
0
N_MODPAT LONG
0
HARMONIC_MASK BYTE
Array[27, 5]
HARMONIC_MAX BYTE
1
SASZERO
BYTE 0
TIME_BIN_DEF LONG
Array[9]
TIME_BIN_MIN INT
1024
XYOFFSET FLOAT
Array[2]
RADIUS
FLOAT 100.000
SPIN_PERIOD FLOAT
4.00000
JITTER_AMP FLOAT
10.0000
JITTER_PERIOD FLOAT
0.200000
UT_PHASE0 DOUBLE
6.7867200e+08
DUMMY
BYTE 0
UT_REF
DOUBLE 6.7867200e+08
TIME_RANGE DOUBLE
Array[2]
IN_TIME_UNIT INT
1
ENERGY_BAND DOUBLE
Array[2]
A2D_INDEX_MASK BYTE
Array[27]
TIME_UNIT LONG
0
N_EVENT
LONG
0
SIM_A2D_INDEX_MASK
BYTE = Array[27]
SIM_PIXEL_SIZE FLOAT
1.00000
SIM_F1D_PIXEL_SIZE
FLOAT = 0.500000
SIM_XYOFFSET FLOAT
Array[2]
SIM_MODEL_PTR STRUCT ->
GAUSSIAN_SOURCE_STR Array[2]
SIM_TIME_RANGE DOUBLE Array[2]
SIM_PHOTONS_PER_COLL
FLOAT = 1000.00
SIM_COUNT_DIST_PTR
FLOAT = Array[10000]
SIM_TIME_PROFILE_PTR
STRUCT = -> <Anonymous> Array[1]
SIM_TIME_UNIT INT
1
SIM_UT_REF DOUBLE
6.7867200e+08
SIM_BACKGROUND FLOAT
0.00000
SIM_BKGD_TIME_PROFILE_PTR
STRUCT = -> <Anonymous> Array[1]
SIM_BKGD_COUNT_DIST_PTR
FLOAT = Array[10000]
SIM_ENERGY_BAND FLOAT Array[2]
SIM_MAX_SIZE FLOAT
128.000
SIM_JUST_BACKGROUND
BYTE = 0
SIM_SASZERO BYTE
0
The keywords shown above reflect two distinct jobs being handled by hessi_image: 1. image reconstruction and 2. data simulation. All of the tags which do not begin with sim_ either control some aspect of the image reconstruction or give information regarding the reconstructed image. All of the tags which do begin with sim_ control some aspect of the simulated data. Each of these structure tags can passed in as a keyword to hessi_image in order to change how the image is reconstructed. For example, suppose you wanted to use the mem sato image reconstruction method, you could simply type:
IDL> hessi_image,im,pars,/plot,image_algorithm='mem sato'
Although, in general it is quite time consuming to simply invoke mem sato with the other default imaging parameters. Substantial time savings is achieved by limiting the number of collimators used via the a2d_index_mask keyword (see below). Again, it is also important to note that the software remembers any changes made to the imaging control parameters in all subsequent usage. In the examples given above, the mem sato algortihm will continue to be used for each new image until it is explicitly changed, perhaps with the call:
IDL> hessi_image,im,image_algorithm='back projection'
HESSI Image Reconstruction Algorithms
Currently, there are 4 different image reconstruction algorithms available for use with hessi_image. These are: 1.) back projection + clean, 2.) mem sato, 3.) pixons, and 4.) forward fitting. Here, we give a brief description of each algorithm, its basic usage, some hints about appropriate usage, and a list of important keywords used with each algorithm. Additional keywords can be found in the list of Image Control Keywords given below.
1. Back Projection + Clean is the default image control algorithm with clean disabled. Back projection simply operates by multiplying the calibrated event list by the collimator modulation patterns to construct a dirty map of the image. This map can then be improved by cleaning, which is set through the clean_alg keyword. Thus, the IDL call:
IDL> hessi_image,im,filename='telemetry.fits'
will construct a dirty map of the image using the data file telemetry.fits. To clean the resulting image, you can simply type:
IDL> hessi_image,im,/clean_alg
if you have already made the first call to create the dirty map, otherwise you can do the whole thing with one IDL call by typing
IDL> hessi_image,im,filename='telemetry.fits',/clean_alg
All the basic image control keywords related to image size, position, pixel size etc. can be used with back projection. These keywords are identified by giving the usage as all in the list of image control keywords given below. The use of clean is unique to the back projection image algorithm. The clean algorithm used here is a straight forward iterative procedure that takes the dirty map, finds the brightest pixel, uses some fraction of this value and the modulation paterns to compute its side lobes and subtracts the result from the dirty map. The process is then repeated until the brightest pixel is negative, or the maximum number of iterations has been reached. The maximum number of iterations is one of a number of clean parameters that can be specified. The parameters related to clean are passed back as part of the imaging control parameters. Thus, in a call such as:
IDL> hessi_image,im,pars,/plot,/clean_alg
the following tag names on the structure pars define the clean related image parameters. These parameters with their default values are:
niter - (default = 30) max number of clean iterations which
will be performed
negative_max_test - (default = 1) a flag for stopping when the
highest absolute value comes from a negative peak.
frac - (default = 0.1) fraction to multiply the maximum by when
subtracting the PSF from the dirty map.
sigma_factor - (default = 0.206) factor by which the pitch of
the grid parameters will be multiplied by to generate the assumed width
(Gaussian sigma) of a point source.
sigma_beam - (default = pitch*sigma_factor) the assumed
width (Gaussian sigma) produced in the image by a point source.
taperpsf - (default = fltarr[9]+1.0) tapering weights
which will be applied to each collimator.
chi_sq_crit - (default = -1.0) the maximum chi square value
obtained by clean.
lambda - (default = -1.0) clean tweaking parameter #1.
mu - (default = -1.0) clean tweaking parameter #2.
nu - (default = -1.0) clean tweaking parameter #3.
sigma - (default = -1.0) clean tweaking parameter #4.
tau - (default = -1.0) clean tweaking parameter #5.
Simple back projection without clean is a fast image reconstruction technique which is good for initial examination of an image. For eaxmple, if you are looking for specific instances (say flares with at least 3 roughly equally bright sources), back projection without clean is useful for looking at a bunch of images to select out times for additional analysis. Once you have honed in on a chunk of data you are particularly interested in, you will almost certainly want to sharpen the image, at least by using the clean_alg keyword, and probably by using one of the more sophisticated image reconstruction algortihms.
2. Mem Sato is an optional image reconstruction algorithm utilizing the maximum entropy method described by Sato, Takeo, and Kazuo 1999, PASPJ, 51, 127. In general, mem sato gives very sharp images, but is time consuming to run. This is especially true is if all the collimators are selected for the imaging. The finest two collimators (0 and 1) require substantial computing time since they have the narrowest PSFs. Eliminating these collimators from the calculation will cause a slight loss of resolution, but will dramatically increase the speed of the image reconstruction. This is done through the a2d_index_mask keyword. For example, if you wish to use mem sato to construct an image of the default model using only collimators 3 - 5 (indices 2 - 4), you would type:
IDL> hessi_image,im,pars,/plot,image_algorithm='mem sato',a2d_index_mask=[0,0,1,1,1,0,0,0,0]
3. Pixon based image reconstruction tries to maximize the image resolution in locations where the information content of the data is high, but permits the resolution to be low where there is little information. Thus, the 'pixels' in a pixon code can vary in shape and size accross the image. The method used here is the fractal pixon method of Puetter and Pina 1994 (Experimental Astronomy, 3, 293: see also Pina and Puetter 1993, PASP, 105, 630 and Metcalf et al. 1996, ApJ, 466, 585). In principle, the pixon method is superior to both back projection and MEM in terms of noise suppression and recovered photometry, but is even more time consuming. It is suggested that users utilize a different technique to decide the key issues of time coverage, energy range, etc. for their images, and then use the pixon method to create final images for analysis and display. The pixon specific optional keywords accepted by hessi_image are:
sensitivity - (default = 1.0) The sensitivity in creating the
pixon map. 0.0 is most sensitive, >3.0 is very insensitive. The default
was arrived at after much testing and simulation. Don't mess with it unless
you know what you are doing. If sensitivity is set too small, you
may get spurious sources; too large and you may miss real sources.
guess - (default guess is computed by the quick look software)
The initial guess for reconstructed image. If guess is not nx by
ny or if all elements of guess are LE 0, then the guess is not used.
resolution - (default = 1) Minimum resolution in pixels.
pixon_sizes - (default = powers of sqrt(2.0)) An integer list
of pixon sizes (resolutions) to use. Slower, but more robust, would
be something like pixon_sizes=indgen(17).
snr - (default = 0, can be set to 1) Weight the residuals by
the signal-to-noise ratio.
chisqr - (default = 0, can be set to 1) Use chi square statistics
rather that poisson statistics. Sigma will be used if set. If sigma
is not set, the Poisson statistics ala /poisson are used.
poisson - (default = 0, can be set to 1) Use Poisson statistics
rather than chi^2. If this keyword is set, sigma is not used.
iterate - (a return variable) Returns with the number of iterations
used at the final resolution.
outresolution - (a return variable) Returns with the final resolution,
i.e. the minimum pixon size.
notty - (default = 0, can be set to 1) Not running on a tty.
Set this for background jobs.
quiet - (default = 0, can be set to 1) Work quietly.
btot_in - (a return variable) An estimate of the total # of
counts in the image. The default is a weighted sum of the counts
in the SC's.
residual - (a return array) The data residuals
gof - (a return variable) The goodness of fit
error - (a return array) Outputs the estimated uncertainty on
the reconstructed image
4. Forward Fitting assumes that an image can be represented by
a number of superimposed 2D gaussians, which can be either spherical, elliptical,
or curved. The algorithm fits the modulation patterns calculated
from a modelmap (specified by Gaussian parameters) to the calibrated event
list of the telemetry file. A detailed description of the algorithm
and its testing can be found in Aschwanden et al. 2000 (ApJ, submitted).
The primary advantages of the forward fitting image algorithm is its speed
and relatively robust estimation of the
source parameters. It's drawbacks are that it currently works
well only for images with 3 or fewer discreet 2D Gaussian like sources.
More complicated real sources will be shoehorned into the assumed 2D Gaussian
shapes. The user can specify the number of sources to find, or the
routine will make an initial guess based on an image computed using back
projection. The forward fitting specific image control keywords
along with their default values are:
telefile - (string variable, default not set) Observed or simulated
telemetry file (*.FITS)
savefile - (string variable, default not set) IDL saveset file
to which the time profiles of the observed calibrated event list and fitted
modulation patterns will be saved along with a number of additional parameters
useful for test displays.
simfile - (string variable, default not set) file containing
simulation parameters (see FWD_SIMULATE)
n_par - (default = 4) The number of free parameters per gaussian,
n_par=4,6,7,...
n_gaussians - (default = 1) The number of gaussian components
to fit.
min_sigma - (default = 1.02) Keyword to select detector with
significant modulation.
nitmax - (default = 10) Maximum number of iteration cycles.
If nitmax is set to zero, the input model from /SIMFILE is directly
used.
minsep - (default = 25.0) The minimum separation of sources
(in arcseconds) in gaussian source decmposition.
fixed_pixel - (default = 0, can be set to 1) Option to disable
automated pixel adjustement in FWD_MINIMUM
ftol - (default = 1.0e-4) Tolerance value for convergence of
forward-fitting algorithm.
plot - (default = 0, can be set to 1) Enables testplots
from procedures FWD_DECOMPOSE and FWD_MINIMUM. NOTE: This keyword
currently conflicts with basic plot keyword sent to hessi_image
to simply make it plot the final image at the end of the reconstruction.
verbose - (default = 0, can be set to 1) Set to print comments
on screen
HESSI Imaging Control Keywords
Here, a description is given of each of the imaging control keywords, along with their default values.
a2d_index_mask - 27 element byte array - Controls which collimator's data is used in the image reconstruction. A value of 1 indicates the collimator will be used, and a value of 0 indicates that it will not. The first 9 elements refer to the front segments of collimators 1-9, the next 9 elements of the array refer to the low energy portion of the rear segments of collimators 1-9, and the last 9 elements refer to the high energy portion of the rear segments for collimators 1-9. The energy at which the changeover occurs in the rear segments is not yet known. As an example, to contstruct an image using only the front segments of collimators 3,4,5, and 6, a2d_index_mask should be set equal to [0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]. The default is to use the front segements of all 9 collimators in the first harmonic: that is the first 9 elements are set to 1 and all others are set to 0.
algorithm_used - a returned informational string variable - gives the filename of the HESSI imaging program used to reconstruct the image.
clean - byte - I think this is redundant with clean_alg given below, and may not survive. See clean_alg for a description of what it is supposed to do.
clean_alg - byte - Controls whether clean is performed on images produced using 'back projection' and does nothing when other imaging algorithms are used. Default value is 0, corresponding to no cleaning, while 1 turns clean on.
debug - byte -
energy_band - double precision 2 element array - Controls the energy band used to construct the image. Units are in keV, with defaults of [6.,100.].
filename - string - Is used to specify the filename containing the HESSI telemetry data which will be used to construct the image.
harmonic_max - byte - maximum harmonic to be used in the constructing the image. Default is 1, maximum is 5.
harmonic_mask - a 27 X 5 byte array - an array giving the mask for the harmonics to be used. The 27 elements reference the a2d indices (see a2d_index_mask). The 5 elements reference the 5 potential harmonics that can be used: element 0 is the 1st harmonic, element 1 is the 2nd harmonic, and so on. A value of 1 in the mask will use the given harmonic for the given collimator, and a value of 0 will force the given collimator/harmonic combination to not be used in the image reconstruction. The default is to use only the 1st harmonic in collimators 1-9, so the mask has elements [0:8,0] = 1, while all other elements are set to 0.
image_algorithm - a string variable - Is used to set what imaging algorithm is to be used in the image reconstruction. The default is 'back projection'. The two currently available algorithms are 'back projection' and 'mem sato', with 'pixons' and 'forward fitting' to be available soon. When using the 'mem sato' algorithm, substantial improvements in speed result if the first 2 elements of the a2d_index_mask are set to 0.
image_dim - a 2 element integer array - Sets the size of the image to construct. Units are in pixels, and it is generally recommended to use square images in the reconstruction. The total size (in arcseconds) of the resulting image is then the size of the image given by image_dim multiplied by the pixel_size described below. The default for image_dim is [64,64].
in_time_unit - integer -
jitter_amp - float -
jitter_period - float -
lightcurve - probably will not survive as a keyword to hessi_image, but will become its own module.
longscale - byte -
modpat_file - string - filename to which the modulation patterns are written.
modpat_set_max - integer -
need_update - byte -
n_event - long -
n_modpat - long -
noclean - Can be given to turn off cleaning. Is redundant with setting clean_alg = 0.
out_imgalg_pars
pixel_size - a 2 element floating point array - Sets the size of the pixels used in the image reconstruction. Both an x- and y-size are specified (hence the two elements, x- first), allowing for the use of non-square pixels. Units are in arcseconds and the default is [1.,1.].
plot - byte - A flag to tell the program to plot the final image automatically. Default is 0 - no plotting. Set to 1 to see a plot of the final image.
radius - float - The distance of the current radial arc segment from the center of the Sun.
saszero - byte - Used in currently available simulation mode. Default is 0, which means no aspect solution is used, and assumes the spin axis is aligned with Sun center. When real data becomes available, the SAS solution will be taken from it.
spectrum - probably will not survive as a keyword to hessi_image, but will become its own module.
spin_period - float - The spin period of the satellite. The nominal (default) period is 4 seconds, but this will be updated in flight and the software will read the current value from the telemetry data files.
time_bin_def - a 9 element long integer array - the time bins used for each collimator. Units are in milliseconds.
time_bin_min - integer - the minimum allowed time bin for any collimator, in units of milliseconds. Default is 512.
time_range - double precision 2 element array - Controls the time range for the data used to construct the image. Units are in seconds from the start time of the data file. The default is [0.,4.], ideally corresponding to one full rotation of the spacecraft. In reality, the spacecraft spin period is likely to be slightly different, requiring a setting like [0.,3.943] to specify one full rotation.
time_unit - integer - Default value is 1.
ut_phase0 - double -
ut_ref - double precision number - The UT given in milliseconds from the beginning of the year. This is the start time taken from the data file on disk or is assigned a default value of 6.7867200e+008 for simulated data.
verbose - byte - Is a flag to set whether the software spews information as it is churning along. Default is 1, which is to spew information to the screen. Setting to zero will minimize the amount of information printed to the screen.
vimage - pointer - can be used to return an estimated count rate from a specified image? Again, not generally set by the user.
vrate - pointer - mainly used for internal testing. Is set only if for some reason a photon rate other than that observed is to be used? This keyword will not generally be set by the user.
xaxis - not really working at this point, but could perhaps be used to return the X-axis of the resulting image.
xyoffset - a 2 element floating point array - Sets the center of the region to be imaged. Units are in arcseconds from the telescope axis, which is nominally supposed to be pointing at the disk center of the Sun. The default pointing is [600.,200.].
xypsf - a longword integer -
xypsf_ptr - long -
yaxis - not really working at this point, but could perhaps be used to return the X-axis of the resulting image.
Simulating HESSI Data
As described above, images can be constructed from HESSI telemetry data files, or the user can use hessi_image to actually simulate HESSI data and then analyze it. A more sophisticated simulation tool is under development; however, much can be done with hessi_image. The default model consists of two point sources (one centered, one in the lower left portion of the image), assuming a photon rate of 1000 photons per second per collimator. The simulation control keywords all have the prefix sim_ to set them apart. Most models of HESSI data use Gaussian sources of various sizes and aspect ratios. The suite of HESSI software has a convenient means built in to simulate Gaussian sources. One must specify the source amplitude (normalized to unit maximum - the overall flux normalazation is handled in the number of photons specified), location with respect to the center of the simulation box, width in the X and Y directions, and the tilt angle. This information is entered into a source structure for use by the software. For example, suppose you wish to specify a single Gaussian source offset +10 arcseconds in the X-direction and +2 arcseconds in the Y-direction from the center of the image, with a X width (sigma) of 5.0 arcseconds and a Y width of 2.0 arcseconds, and a tilt angle of 45.0 degrees, the following commands should be entered:
IDL> hsi_struct_def
(note: if you have already run hessi_image, you do not need to type the above command)
IDL> gg = {gaussian_source_str}
IDL> gg.amplitude = 1.0
IDL> gg.xysigma = [5.,2.]
IDL> gg.xypos = [10.,2.]
IDL> gg.tilt_angle_deg = 45.
IDL> hessi_image,im,model=gg
Note that the keyword specifying the model does not have the leading sim_ in it. This is a bug and will hopefully be fixed soon.
You can specify multiple Gaussian sources as well. Suppose you wish to specify 2 sources, both with widths of [5.0,2.0], with one located at an offset [20.0,-3.0] and one half as strong as this at [-5.0,15.0], both with 0.0 tilt angles, you would type:
IDL> gg1 = {gaussian_source_str}
IDL> gg1.amplitude = 1.0
IDL> gg1.xysigma = [5.,2.]
IDL> gg1.xypos = [20.,-3.]
IDL> gg1.tilt_angle_deg = 0.
IDL> gg2 = {gaussian_source_str}
IDL> gg2.amplitude = 0.5
IDL> gg2.xysigma = [5.,2.]
IDL> gg2.xypos = [20.,-3.]
IDL> gg2.tilt_angle_deg = 0.
IDL> model = [gg1,gg2]
IDL> hessi_image,im,model=model
Another means of specifying point sources is to simply create a model
image loaded with 0s everywhere except for the locations
of the point sources, which should be loaded with a normalized intensity.
For example, to create a 64 X 64 image with 3 point
equal strength point sources at positions [10,20], [30,12], and [15,40]
the following syntax can be used:
IDL> model=fltarr(64,64)
IDL> model[10,20]=1.0
IDL> model[30,12]=1.0
IDL> model[15,40]=1.0
IDL> hessi_image,im,sim_model=model
The strength of the point sources again is given in normailized units, while the overall flux level is given by the keyword sim_photons. In the above example, all three point source are of equal strength. If you wanted to make the one at [30,12] equal to 1/4 of the other two, you would have typed:
IDL> model[30,12]=0.25
instead. The simulated data will be constructed assuming the default of 1000 photons per second per collimator incident on the spacecraft, which is not particularly good statistics, especially when spread over 3 sources. The overall flux level in the simulation could be increased with the call:
IDL> hessi_image,im,sim_model=model,sim_photons=5000
HESSI Simulation Keywords
The following is a list of keywords that can be specified to control the simulated data to be produced.
sim_model - typically given as an array with zeroes everywhere except at the location of point sources. The non zero values should be specified in normalized units, with a maximum of 1.0. Alternatively, one can use the keyword model (without the leading sim_) to specify single Gaussian components, controlling the amplitude, x and y standard deviations, and the tilt angle of the Gaussian.
model - Can be used to pass in a structure (labelled gaussian_source_str) which specifies a 2d Gaussian source. The Gaussian is specified by its amplitude, x and y standard deviation (in arcseconds), its position relative to the center of the image (in arcseconds), and its tilt angle. See above for examples.
I
The result will show the image of this Gaussian source and load it into the variable im.
sim_photons - floating point number - Is used to control the overall flux level of the simulated data. The units are in photons per second per collimator incident on the spacecraft, with a default value of 1000.
sim_a2d_index_mask
sim_pixel_size - Is used to set the size of the pixels when using a model image of some size (e.g. 64 X 64). The default sim_pixel_size is 2.0 arcseconds, which is different from the default imaging pixel_size.
sim_xyoffset
sim_time_range
sim_time_unit
sim_ut_ref
sim_background
sim_energy_band
sim_max_size
