|
 |
|
|
SimSET provides several ‘user functions’ where the user can intercept a photon or event during its processing and modify it, thus altering how it will be processed further downstream. Normally the functions do nothing, but they are provided as convenient places for the user to alter program execution. These functions are found in the source code files ColUsr.c (for the collimator module), DetUsr.c (for the detector module), PhgUsrBin.c (for the binning module), addrandUsr.c (for the randoms processing), and their associated header (.h) files. The functions are called at certain critical points in decay/photon processing, providing a safe means of modifying or adding to SimSET. Furthermore, user functions can be ported to new versions of SimSET easily, unlike changes to the main modules. This web page lists the user functions, gives some of the reasons to use one, and gives ideas of how to use the functions.
General design and usage
Each of the user function modules contains an initialization function that gets called at the beginning of a simulation, processing functions that get called for each photon or event, and a termination function that gets called at the end of the simulation. The user may use the initialization and termination functions to allocate/deallocate ‘local global’ and global variables, initialize variables, read in user function parameters, error check input parameters, write output, report results, etc. However, the meat of the user functions is in the processing functions: it is here that the user can alter the flow of the simulation, and these functions tend to dictate what needs doing in the initialization and termination functions. The processing functions are passed each photon (or coincidence) with all its simulation history. The function can reject the photon (or coincidence) outright, or can modify some of the fields in the photon to change how it is treated downstream. Below we list each of these modules, where they are called, and indicate possible uses for each.
[top of page] [examples]
Collimator user module
The collimator user module, ColUsr.c and ColUsr.h, contains four functions: ColUsrInitialize, ColUsrPETPhotons, ColUsrSPECTPhotons, and ColUsrTerminate. Global #defines, variables and structures can be declared in ColUsr.h (if they are needed outside the collimator user module, for instance in one of the other user modules) or at the beginning of ColUsr.c (for 'local globals').
ColUsrInitialize is called by the collimator initialization function, ColInitialize in Collimator.c, after the collimator parameters have been read in. ColUsrInitialize has access to all the collimator parameters, which can be found in the data structure that is passed to the function, *ColRunTimeParamsPtr (the data type is given in ColParams.h). While this structure is passed as a pointer, the user should be circumspect about altering its contents: the contents would be better modified using the collimator parameter file. This function is intended for allocation of memory and initialization of variables to be used in ColUsrPETPhotons and ColUsrSPECTPhotons. The user can also create an extra parameter file and use this function to read its contents.
Depending on whether the simulation is a PET or SPECT simulation, ColUsrPETPhotons or ColUsrSPECTPhotons is called for each photon as it enters the collimator. They are passed pointers to the collimator parameters (ColRunTimeParamsPtr), the decay currently being processed (decayPtr) and the photon currently being tracked (photonPtr). (The data types for the latter two are in Photon.h.) The user can modify any of the decay or photon data, thus altering downstream processing, or can discontinue the tracking of the current photon by setting the function return value to false.
ColUsrTerminate is called after all photon tracking and event processing as the first step of the collimator termination function, ColTerminate, in Collimator.c. ColUsrTerminate receives a pointer to the current collimator parameters, ColRunTimeParamsPtr. The user should deallocate any memory allocated in the other ColUsr functions, and can report results or write output.
[top of page] [examples]
Detector user module
The detector user module, DetUsr.c and DetUsr.h, contains four functions: DetUsrInitialize, DetUsrPETPhotons, DetUsrSPECTPhotons, and DetUsrTerminate. Global #defines, variables and structures can be declared in DetUsr.h (if they are needed outside the detector user module, for instance in one of the other user modules) or at the beginning of DetUsr.c (for 'local globals').
DetUsrInitialize is called by the detector initialization function, DetInitialize in Detector.c, after the detector parameters have been read in. DetUsrInitialize has access to all the detector parameters, which can be found in the data structure that is passed to the function, *DetRunTimeParamsPtr (the data type is given in DetParams.h). While this structure is passed as a pointer, the user should be circumspect about altering its contents: the contents would be better modified using the detector parameter file. This function is intended for allocation of memory and initialization of variables to be used in DetUsrPETPhotons and DetUsrSPECTPhotons. The user can also create an extra parameter file and use this function to read its contents.
Depending on whether the simulation is a PET or SPECT simulation, DetUsrPETPhotons or DetUsrSPECTPhotons is called for each photon as it enters the detector. They are passed pointers to the detector parameters (DetRunTimeParamsPtr), the decay currently being processed (decayPtr) and the photon currently being tracked (photonPtr). (The data types for the latter two are in Photon.h.) The user can modify any of the decay or photon data, thus altering downstream processing, or can discontinue the tracking of the current photon by setting the function return value to false.
DetUsrTerminate is called after all photon tracking and event processing as the first step of the detector termination function, DetTerminate, in Detector.c. DetUsrTerminate receives a pointer to the current detector parameters, DetRunTimeParamsPtr. The user should deallocate any memory allocated in the other DetUsr functions, and can report results or write output.
[top of page] [examples]
Binning user module
The binning user module, PhgUsrBin.c and PhgUsrBin.h, contains six functions: PhgUsrBinInitialize, PhgUsrBinPETPhotons, PhgUsrBinSPECTPhotons, PhgUsrBinPETPhotons2, PhgUsrBinSPECTPhotons2, and PhgUsrBinTerminate. Global #defines, variables and structures can be declared in PhgUsrBin.h (if they are needed outside the binning user module, for instance in one of the other user modules) or at the beginning of PhgUsrBin.c (for 'local globals').
PhgUsrBinInitialize is called by the binning initialization function, PhgBinInitialize in PhgBin.c, after the binning parameters have been read in. PhgUsrBinInitialize has access to all the binning parameters, which can be found in a data structure that is passed to the function, *binParams (the data type is given in PhgParams.h). While this structure is passed as a pointer, the user should be circumspect about altering its contents: the contents would be better modified using the binning parameter file. PhgUsrBinInitialize also receives a pointer to the histogram data arrays, *binData (the data type is given in phg.h): in most cases these arrays will be initialized to zeroes at this point, but if the user has selected "add_to_existing_img = true" the renormalized histograms from previous runs will already be stored there. Again, while this structure is passed as a pointer, the user should be circumspect about altering the data. PhgUsrBinInitialize is intended for allocation of memory and initialization of variables to be used in PhgUsrBinPETPhotons, PhgUsrBinSPECTPhotons, PhgUsrBinPETPhotons2 and PhgUsrBinSPECTPhotons. The user can also create an extra parameter file and use this function to read its contents.
PhgUsrBin.c has two functions for processing PET photons/events and two for processing SPECT photons/events. PhgUsrBinPETPhotons or PhgUsrBinSPECTPhotons is called for each photon as it enters the binning function, before any steps are taken to accept/reject it or compute histogram indices for it. The functions are passed pointers to the binning parameters (*binParams), the histogram data arrays (*binData), the decay currently being processed (*decay). For PET pointers to the current blue and pink photons are also passed, *bluePhoton and *pinkPhoton; for SPECT, a pointer to the current photon is passed (*photon). (The data types for the photons and decays are in Photon.h.) The user can modify any of the decay or photon data, thus altering downstream processing, or can discontinue the tracking of the current photon(s) by setting the function return value to false.
The second set of functions for processing PET photons/events, PhgUsrBinPETPhotons2 and PhgUsrBinSPECTPhotons2, are called just before the output histograms are incremented for the current event. Note that these functions are not called if the event is, for some reason (e.g., energy below the cutoff), rejected within the PhgBin.c binning function. These functions are passed the same arguments as PhgUsrBinPETPhotons and PhgUsrBinSPECTPhotons, plus plus pointers to all the indices for binning (please see the source code file!), and a pointer to the index into the histograms to be updated (*imageIndex). Please note that *imageIndex must be changed if the histogram position for the event is to be changed: changing one of the component indices will not change this position as *imageIndex is used for incrementing the array. PhgUsrBinPETPhotons2 also receives pointers to the values that will be used when incrementing the output weight and weight-squared histograms, *coincidenceWt and *coincidenceSqWt.
PhgUsrBinTerminate is called after all photon tracking and event processing as the first step of the binning termination function, PhgBinTerminate, in PhgBin.c. PhgBinTerminate receives a pointer to the current binning parameters, *binParams, and to the histogram data arrays, *binData. The user should deallocate any memory allocated in the other PhgUsrBin functions, and can report results or write output.
[top of page] [examples]
Randoms user module
The random coincidence software includes a user module for the randoms processing, addrandUsr.c and addrandUsr.h, containing four functions: addrandUsrInitialize, addrandUsrModDecays1, addrandUsrModDecays2, and addrandUsrTerminate. Global #defines, variables and structures can be declared in addrandUsr.h (if they are needed outside the randoms user module, for instance in one of the other user modules) or at the beginning of addrandUsr.c (for 'local globals').
addrandUsrInitialize is called by the addrandoms.c initialization function, adrandInitialize, after the randoms processing parameters have been read in but before processing of the time-sorted history list starts. addrandUsrInitialize has access to all the parameters for randoms processing through a pointer *detParams (the data type is given in DetParams.h). Note that this structure may not be completely filled in, or may be filled in differently than it was for the detector module. However, the parameters needed for randoms processing will be filled in. While this structure is passed as a pointer, the user should be circumspect about altering its contents: the contents would be better modified using the addrandoms parameter file. addrandUsrInitialize is intended for allocation of memory and initialization of variables to be used in addrandUsrModDecays1 and addrandUsrModDecays2. The user can also create an extra parameter file and use this function to read its contents.
addrandUsr.c has two functions for processing PET photons/events, addrandUsrModDecays1 and addrandUsrModDecays2. addrandUsrModDecays1 is called before the photons/decays are time-windowed. All the decays within a time-window of each other (or actually longer: the window is extended each time another photon arrives until the time elapsed between succesive photons is greater than the time window) are passed in the timeWindowDetectionsTy structure (declared in addrandUsr.h). The photons and decays may be altered/deleted as desired, e.g. to account for deadtime, different triples handling than SimSET's (currently all triples are deleted). This function has no return value.
addrandUsrModDecays2 is called for each coincidence after the photons/decays are time-windowed. Only the decays that are actually going to be written out will be passed to this routine - many of the decays/time windows passed to addrandUsrModDecays1 will not generate a call to this function. The function returns a boolean that determines if the coincidence will be accepted. This gives the user a chance to change/reject the coincidences - whether random, scatter or true - that are getting written out.
addrandUsrTerminate is called after the entire time-sorted history file has been processed as the first step of the addrandoms termination function, adrandTerminate, in addrandoms.c. The user should deallocate any memory allocated in the other addrandUsr functions, and can report results or write output.
[top of page] [examples]
Ideas for user functions
Below we offer an example of a user function; we would also like to solicit examples from users to add to this page.
User functions are most useful in cases where a few changes to the attributes of a photon or decay can improve the realism of the simulation or add a new functionality. They can even be used to supplant the entire functionality of their module with a new model, for instance a new collimator. Below we suggest a few possible uses for each of the user function modules. Note that the ColUsr.c and DetUsr.c photon tracking functions are called only as the photon enters the collimator or detector, respectively; thus, to modify a photon as the photon exits the collimator or detector, one must use the user functions in the next module, e.g. DetUsr.c and PhgUsrBin.c respectively.
Collimator user functions: ColUsr.c
- Unsupported collimators. If the collimator to be simulated is different from those SimSET supports, the user could implement their own collimator here.
- Modifications to existing collimators. For instance, for a collimator like one of those currently offered, but which doesn't cover the full angular range, the user could reject photons that would fail to hit the collimator in ColUsrPETPhotons or ColUsrSPECTPhotons.
Detector user functions: DetUsr.c
- Pinhole collimator. One might be able to simulate a pinhole collimator by using a PET collimator simulation, creating a collimator that used the axial profile of the pinhole to define the axial composition of the PET septa. The transaxial and off-axis components of the pinhole would then be approximated in the DetUsr.c function DetUsrSPECTPhotons. Photons are passed to this function before they enter the detector. Projecting the photon position/direction back to the collimator ring, one can determine where the 'pinhole' would have to have been for the photon to pass through; this value can be stored in the photonPtr field detectorAngle. The user binning module, PhgUsrBin, would then use this value, along with the detected position of the photon in a planar detector, to determine the correct position to bin the photon. This pinhole collimator approximation would not accurately reflect collimator penetration or collimator scatter, though some tabular approximations of these effects could be used to improve the modeling. Note, as mentioned above that this modification of the collimator is implemented in the DetUsr and PhgUsrBin files, rather than in the ColUsr files.
Binning user functions: PhgUsrBin.c
- Create a new binning variable. Users can either do their own binning or alter the photon(s) to trick SimSET into binning the new variable into an otherwise unused variable. For example, a user not binning by photon energy can, nevertheless, set SimSET up to bin by energy, but then use PhgUsrBinPETPhotons or PhgUsrBinSPECTPhotons to change the photon energy to a value corresponding to the desired bin for the new binning variable.
- Change bin boundaries. SimSET uses evenly spaced bins, e.g. if the user specifies two energy bins from 110 to 160 keV, the bins will be 110-135 and 135-160. If, however, the user wished to have the bins be 110-125 and 125-160, photons with energy 125-135 could be changed to have energy 120 in PhgUsrBinPETPhotons or PhgUsrBinSPECTPhotons.
- Change detected position. Most of the SimSET detectors position photons at the energy-weighted centroid of the interactions in the crystal. This is not an accurate representation of how some detectors work. In PhgUsrBinPETPhotons or PhgUsrBinSPECTPhotons the user has access to all the interaction information (position and energy deposited for each interaction) and can alter the detected position before the photon (or coincidence) is binned. (An example of this option is given below.)
- Weight events based on their history or bin. For example, when imaging 140 keV photons, one might want to give more credence (a higher weight) to a photon detected at 145 keV than 135 keV (a scatter photon is less likely to register the former than a true photon). This would be done in PhgUsrBinPETPhotons2 or PhgUsrBinSPECTPhotons2.
- Approximating detector gaps. This is a possibility that our group used quite successfully to approximate block detectors using the cylindrical detector model before the block detector software was available. We determined in PhgUsrBinPETPhotons (or PhgUsrBinSPECTPhotons) if a photon had interacted in parts of the cylindrical detector that would have been in the gaps between the block detectors we were trying to simulate and altered the detected position and energy of the photon accordingly. While this particular approximation is no longer necessary, similar approximations might be useful for structures in the detector that are not well represented by right-rectangular boxes. Note, as mentioned above that this modification of the detector is implemented in the PhgUsrBin files, rather than in the DetUsr files.
- Split block detectors into individual crystals. An example of this is given below. As noted in the bug reports and elsewhere on the News page, the current version of block detectors does not allow blocks to be subdivided into smaller crystals. We developed a work-around for this, as described below.
[top of page] [examples]
Example: splitting blocks into crystals
(This example works only with version 2.9! We are leaving it here for now as it is the only example we have written up to this point.)
As noted in the bug reports for version 2.9 and elsewhere on the News page, the version 2.9 of block detectors did not allow blocks to be subdivided into smaller crystals.
We used version 2.9 to simulate the General Electric DSTE positron emission tomograph. To do so, we defined a block-based tomograph based on the DSTE with the blocks represented as solid blocks of crystal; the actual blocks are subdivided into a 6x8 array of crystals. Without a user function, SimSET would track a photon through the detector blocks correctly, compute an energy-weighted centroid of the interactions in the block where the most energy was deposited (so far, so good: this is a reasonable model of the detection/electronics), and then assign the center of the block as the detected position (not so good - only lines-of-response between the centers of blocks will be populated with events).
To design the user functions, we first had to decide how we wanted the output data binned. We decided to bin by crystal pair. SimSET already provides crystal-crystal binning, but it allocates an output array that has its size based on the number of active elements described in the number of crystals defined in the parameter files. As the parameter files we used don't break the blocks down into individual crystals (the problem the user function is trying to solve), this means that the array would be too small if we used this option. We address this by specifying the output as a (z1,z2,angle,distance) array for PET or (z,angle) array for SPECT. The user must set the number of bins for the z variables to the number of rings of crystals axially, and the number of distances and angles for PET, or angles for SPECT, to the number of crystals in each ring.
To solve this problem we wrote PhgUsrBin_crystal.c and PhgUsrBin_crystal.h. The difference between the header file and the default PhgUsrBin.h is 8 '#defines' giving the number of detector blocks axially and transaxially (this could have been taken from the ring parameter files, but we chose to specify it here as well), the number of crystals within a block axially and transaxially, the radius at of the block centers (the detected photon position is placed on the cylinder of this radius), the number of the first active ring with crystal (i.e., some rings may just consist of shielding, this is the first ring which contains active crystal), the number of the first block in this ring which is active, and the depth within the block at which photons' detected position should be set (by default SimSET sets this to the center of the crystal radially, but most tomographs set this to some average depth, e.g. the depth of the first interaction):
#define PUB_BLOCS_AXIAL 4 /* num blocks axially */
#define PUB_BLOCS_RAD 70 /* num blocks radially */
#define PUB_XTALS_AXIAL 6 /* num crystals/block axially */
#define PUB_XTALS_TRANS 8 /* num crystals/block transaxially */
#define PUB_RAD_BLOC_CENT 45.805 /* radius to block cent */
#define PUB_FIRST_ACTIVE_RING 3 /* first ring with active crystal */
#define PUB_FIRST_ACTIVE_BLOCK 0 /* first block in PUB_FIRST_ACTIVE_RING which is active crystal */
#define PUB_INTERACTION_DEPTH 1.0 /* depth in crystal to position events */It is possible that the only modifications a user would have to make to these files to model another tomograph would be to change these #defines. However, the PhgUsrBin_crystal files make many assumptions, e.g. that all the active blocks have the same dimensions. We describe the changes to PhgUsrBin_crystal.c below as an example of how to use the user functions and in case further modifications are necessary for other tomographs.
The function PhgUsrBinInitialize checks that the number of z bins has been set to the number of rings of crystal, and that the number of distances and angles are set to the number of crystals in a ring, calculates a number of constants giving the dimensions of the blocks and the numbers of crystals.
The function PhgUsrBinPETPhotons calculates the energy-weighted centroids of the interactions in the detection blocks for the blue and pink photons, and then sets the detected positions of the photons to the center of the nearest crystal in the axial and transaxial directions, and to the depth within the crystal specified by PUB_INTERACTION_DEPTH. It also calculates new detection crystal numbers for the two photons based on the actual number of crystals in the tomograph. (SimSET, by default, currently assigns the active block number instead as the detection crystal, as the blocks cannot be split into crystals.) This function is complicated by the fact that the positions of the detector interactions are recorded in tomograph coordinates (rather than block-specific coordinates - this is something we may change in future versions). For each photon we compute the energy-weighted centroid in tomograph coordinates, convert the centroid to block coordinates, calculate the crystal in which the centroid lies, assign that crystal's position as the detected position after converting it back to tomograph coordinates.
The function PhgUsrBinPETPhotons2 is called after PhgBinPETPhotons has calculated the distance-angle bin, the z indices, and the image index for the coincidence. It replaces the calculated z-axis bins with the detector crystal ring numbers for the blue and pink photons (in the zUp- and zDownIndex respectively), and replaces the calculated angle and distance indices with the blue and pink photon crystal numbers within the ring respectively (i.e. the active crystal number mod the number of crystals per ring). These new index values are used to compute a new pointer to the output data bin, *imageIndex. Please note that we do not try and make the output upper triangular: there is no difference between a coincidence with (zUp, zDown, angle, distance) = (zBlue, zPink, dBlue, dPink) or = (zPink, zBlue, dPink, dBlue) (it just depends on which photon was blue and which pink).
The functions PhgUsrBinSPECTPhotons and PhgUsrBinSPECTPhotons2, which are used for single photon binning, have been modified in a manner similar to that described above for PhgUsrBinPETPhotons and PhgUsrBinPETPhotons2 respectively. The main difference is that, as only one photon's position need be stored, the crystal position within the ring is stored in the angle index - the distance index is not used.
The function PhgUsrBinTerminate is unmodified.
[top of page] [examples]
- Generic contains samples of all types of detectors.
- Planar SPECT detector parameters file.
[top of page] [examples]
Last revised by: Robert Harrison Revision date: 9 April 2011
