PET Tutorial
This tutorial applies only for PET data. For CT help, see Using CT data.
Two tutorials are presented, a simple one using gate_main_simple.m and a more advanced one utilizing gate_main.m. For non-GATE data (main_PET.m) the flow is similar except there are no GATE specific features. When using main_PET.m or Inveon_PET_main.m, Biograph_mCT_main.m, Biograph_Vision_main.m without GATE data, the user will be prompted for the input data.
Both files (gate_main_simple.m and gate_main.m) are immediately ready to use, once the GATE data has been obtained. Only the path to search for the data needs to be changed (line 308/310 in gate_main_simple.m and line 704/706 in gate_main.m). However, if you use different device than the one in the example, you need to modify at least device/sinogram specific settings accordingly.
Ready-made simulated GATE data can be found from (both ASCII and ROOT format): https://doi.org/10.5281/zenodo.3526859
For reconstruction purposes, ready-made example raw data and sinograms can be found from: https://doi.org/10.5281/zenodo.3522199
Below is an explanation on what the different parameters actually do in gate_main_simple.m. Once you have set all the necessary parameters (when using the demo data from above, simply modify the input data path) simply run the file. The reconstructed images are found in f_osem matrix (gate_main_simple.m only).
The most important part of OMEGA is the Scanner properties section. This has already been filled for the example case, but here are detailed explanations as to where the values are derived from.
Blocks per ring
In the GATE macro camera.mac this is found from line 64. I.e. it is the ring/R-sector repeater. The R-sector in this case contains all the blocks/buckets in the axial direction and is repeated radially along the ring. Inside these R-sectors are the modules that are the individual blocks/buckets.
For non-GATE data, this should be the number of blocks/buckets the scanner has transaxially. This multiplied with the below crystals per block should equal the total number of detectors on one ring (transaxially).
Linear multiply
This is the axial repetition factor, from line 59 in camera.mac. Each module (block/bucket) is repeated this many times in the axial direction. The blocks have small gaps between them (0.2 mm).
For non-GATE data this should be the number of blocks/buckets axially. This multiplied with the below crystals per block should equal the number of crystal rings (e.g. the number of crystals along a line axially).
Transaxial multiplier
Basically same as above, but for transaxial case. In the example case it is simply 1. However, in certain cases there can be transaxial repeaters and in such case you should set that value here.
Crystals per block, transaxial
This is simply the number of crystals/detectors in the module/block/bucket transaxially. Line 53 in camera.mac.
If you have a block with 20x10 crystals, set crystals per block to 20.
Crystals per block, axial
This is the axial crystals per block. In the example case it is the same as the transaxial one.
If you have a block with 20x10 crystals, set axial crystals per block to 10.
Crystal pitch
Both the transaxial (x- and y-direction) as well as the axial (z-direction) crystal pitches need to be input. These are the distances between the centers of adjacent crystals in the same block. Line 55 in camera.mac. There are small gaps between the crystals.
Ring diameter
This is the total diameter of the scanner bore, i.e. the distance between perpendicular detectors. Usually this should be the Rmin, see line 6 in camera.mac.
Transaxial FOV
This is the length of a one side of the transaxial FOV (needs to be rectangular). Both x- and y-directions can be selected individually, though square FOVs are recommended. This parameter can be selected freely, but the current value was selected such that the size is of similar size as that in the Inveon PET scanner (which uses the vendor specified FOV size).
Axial FOV
Same as above, but for axial (z) direction. Has been selected to be slightly smaller than the actual scanner length.
Machine name
Simply the name for your scanner. Used only to name the saved mat-files.
Before loading GATE data, you should specify first if you want to load some GATE specific data. These include trues, scatter, randoms and the source image. For more detailed information on extracting GATE trues, scatter and randoms see Extracting GATE scatter, randoms and trues data.
By default, the main-file saves trues, scatter, randoms and the source images. These can be, respectively, disabled by setting the values on lines 80, 94, 121 and 132 to false. Optionally, the trues or scatter can be reconstructed instead of prompts. For trues this is achieved by setting line 87 to true, and for scatter line 113.
Three GATE input data are supported: ASCII, LMF and ROOT. Of these ASCII is used by default (line 148). The coincidence mask is correct if no modifications have been made. The numbers come straight from the ASCII mask order as specified by GATE (see GATE documentation).
LMF data can be enabled on line 168. For LMF data you need to specify the number of bits dedicated for each geometric element (R-sectors, modules, submodules, crystals, layers). This is output by GATE at the start of the simulation, see the image below:
Furthermore, you need the coincidence window length as well the "clock time step" that can be found from the output cch-files.
ROOT data is enabled from line 216 and doesn’t need any specific information. ROOT data is, however, unstable on MATLAB R2018b and earlier and is not supported on Windows. See Known issues and limitations for more information.
The voxel counts for each slice and the slice count can be specified on lines 236-242. Slice count is usually the number of crystal rings x 2 - 1.
Span factor (line 269) determines the axial compression, that is how many sinograms are combined in the axial direction. Currently 3 is the smallest span factor supported (if you wish to use span of 1, use raw data). Higher span values compress the sinogram data more and result in faster reconstruction, but can have a negative impact on the image quality.
Ring difference is related to the span value and can be at most the number of crystal rings - 1. Oblique sinograms are created from the ring distance specified by the ring difference. I.e. with maximum ring difference the coincidences between the first and last crystal ring are included. For more information on span and ring difference see Data Acquisition in PET Imaging.
The number of angular positions (views) is the first dimension of a single sinogram slice (usually depicted with an s). Primarily you should use the same value as the device you are modelling, but if that is not available you can use the function ndist_max to determine optimal values. This function outputs the orthogonal distance of every line of response as well as displays the view counts for the case of fully encompassing the FOV (i.e. a circle with the FOV square inside) or for a case where the sinogram FOV is fully inside the FOV (i.e. a circle that just fits inside the FOV square). It is recommended the value is set somewhere between these two values.
The number of angles (tangential positions) determine the second dimension of a single sinogram slice (usually depicted with a φ). This should be the number of detectors on a ring divided by two. Smaller values are supported, but they need to be values that are obtained from the number of detectors on a ring divided by a value divisible by two. Smaller values than the number of detectors per ring divided by two will be considered as sinogram mashing (e.g. number of detectors per ring divided by four).
The name (line 301) is used for naming purposes only. I.e. sinograms are saved for specific examination and scanner.
Only sinograms can be computed by setting line 316 to true and running the file. This loads the GATE data and then forms the sinograms, but does not continue to image reconstruction.
Likewise, only reconstruction can be enabled on line 323. Running with this true, will skip the data load and sinogram formation steps. Previously created measurement data will be automatically loaded if such exists (name and machine name match).
Status messages can be turned off by setting line 328 variable to false.
By default, the reconstruced images contain the total number of detected counts. This can be converted to Bq/mL with the scaleImages function (line 441), but you have to provide the total time of the measurement. options.dx/dy/dz are automatically computed.
If you want to export the reconstructed image, you can use saveImage function. For example, you can export the reconstructed image as NIfTI format with the following command saveImage(f_osem, 'nifti', [], options).
This is based on gate_main.m. The process is identical to that of the simple version, except that the non-simple one gives more possibilities for user adjustments. Below is an explanation on what the different parameters present in all the other main-files except for the simple case. Once you have set all the necessary parameters (again for the demo data, simply modify the input data path) simply run the file. The reconstructed images are found in pz cell matrix, where each cell represents each algorithm/prior. See visualize_pet.m for the order of the algorithms. The last cell element contains various variables that were set for the reconstruction (such as number of iterations, subsets, regularization parameters, etc.).
Almost identical to the simple version, but the number of pseudo rings/detectors can be adjusted. Aside from the pseudo rings/detectors, all other variables not present in the simple version are calculated automatically.
For an example on how to use pseudo detectors and rings, see Biograph_mCT_main.m.
Pseudo rings
These can be adjusted in options.pseudot. If your scanner has pseudo rings, input the number of pseudo rings here. If no pseudo rings are present (as is with most scanners), use 0 or empty array [].
Pseudo detectors
Pseudo detectors can be added to options.det_w_pseudo. Normally, without pseudo detectors, this is computed as options.det_w_pseudo = options.blocks_per_ring*(options.cryst_per_block); which is the number of detectors on a one crystal ring. Usually pseudo detectors are a one additional detector on each block, which means that this becomes options.det_w_pseudo = options.blocks_per_ring*(options.cryst_per_block + 1); If no pseudo detectors are present, options.det_w_pseudo should be the same as options.det_per_ring.
Ring difference
If you use raw data, you can specify the maximum ring difference here. Default includes all the rings, but decreasing the value removes axial coincidences. For example, if ring difference is number of rings - 2, then the coincidences between the first ring and the last two rings are not included.
Increased sampling
Raw data sampling can be increased by increasing options.sampling_raw. E.g. options.sampling_raw = 2 doubles the length of the both dimensions of the raw data matrix. This is achieved through interpolation that can be selected with options.sampling_interpolation_method_raw. All interpolation methods supported by interp2 are available. Unlike with sinogram data, the sampling is increased in both dimensions.
This can be used to remove aliasing artifacts caused by too low sampling rate. It is not recommended to use this, if no aliasing artifacts are present and even then it is better to test PSF reconstruction and/or orthogonal distance-based ray tracer first.
Segment table
Table of sinogram segments. This is automatically computed, but can be filled manually as well.
Total number of sinograms
This should be the total number of sinograms, i.e. the sum of the segment table.
Number of sinograms used in reconstruction
options.NSinos can be used to utilize N FIRST sinograms in the image reconstruction. E.g. if you want to use only the direct plane sinograms, then this should equal the number of crystal rings * 2 - 1.
Angular cut-off
options.ndist_side is used to determine from which side is one angular position removed. E.g. if you have 128 views, as with the Inveon PET, that means that you have an even number of views. However, since the smallest orthogonal distance between a LOR and the center of the FOV is usually zero, that means that there should be 64 views on both sides + the center one. This would result in 129 views. Usually, however, one of the views is removed from either side.
Increased sampling
Sinogram sampling can be increased by increasing options.sampling. E.g. options.sampling = 2 doubles the length of the first dimension of the sinogram. This is achieved through interpolation that can be selected with options.sampling_interpolation_method. All interpolation methods supported by interp1 are available. The sampling is only increased in the first (views) dimension, not on the second (angles).
This can be used to remove aliasing artifacts caused by too low sampling rate. It is not recommended to use this, if no aliasing artifacts are present and even then it is better to test arc correction and/or PSF reconstruction and/or orthogonal distance-based ray tracer first.
Gap filling
If pseudo detectors are present on your scanner, these will not obviously receive any measurements. Due to this, the location of these pseudo detectors cause gaps in the formed sinogram. These gaps can be, however, filled if the fill_sinogram_gaps is set to true.
Two different gap filling methods can be selected. The first is fillmissing (built-in function in MATLAB), that uses 1D interpolation. Since it uses 1D interpolation, the gap filling is performed separately for each transaxial dimension. You can also specify the type of interpolation used, default is 'linear'. See the help of fillmissing for details on the interpolation methods.
The second interpolation method is inpaint_nans which is available from MathWorks file exchange. This is a 2D interpolation method. As with fillmissing it also supports different interpolation methods that you can specify. See the help for inpaint_nans for more information.
To apply different corrections to a sinogram previously saved, use options.SinM = form_sinograms(options, true) (e.g. apply additional true input after options).
When using sinogram data that has been previously saved, and options.corrections_during_reconstruction = false, applying any correction will automatically use the corrected sinogram, even if the sinogram did not actually contain that correction. E.g. if you have previously constructed a sinogram with ONLY normalization correction, both the raw uncorrected sinogram as well as the normalization corrected sinograms have been saved. Using this same data (same machine name and name) and setting all corrections to false will cause the uncorrected sinogram to be used during the reconstruction (if the sinogram is not re-created, i.e. options.only_reconstructions = true). Assuming that the sinogram is NOT re-created, setting any correction to true (e.g. randoms correction) will cause the corrected sinogram to be used although the specified correction may not have been included (in this case it would be the normalized sinogram).
The above does NOT apply when options.corrections_during_reconstruction = true as the uncorrected sinogram is always used and corrections are not applied to measurement data. The above also does not apply to raw data. For raw data the corrections are performed in reconstructions_main if options.corrections_during_reconstruction = false.
Attenuation correction is always applied during the reconstruction and as such does not affect the measurement data in any case.
Arc correction is always applied before the reconstruction, but is never saved. Thus arc corrected sinogram cannot be obtained automatically (but can be obtained manually, but saving the options.SinM variable after arc_correction has been run) and is always computed regardless of other choices as long options.arc_correction = true and options.precompute_lor = false.
Randoms correction
Randoms correction in OMEGA is performed by using the delayed coincidence window data. For GATE, Inveon and Biograph mCT/Vision data the delayed coincidences are automatically (and only) collected when randoms correction is enabled. This delayed coincidence data is then used for the randoms corrections either in the sinogram/raw data matrix or during reconstruction (ordinary Poisson). For other devices, the user will be prompted for the location of the randoms correction (mat) file if randoms correction is enabled and then used similarly as in the GATE/Inveon/Biograph case.
If randoms correction is selected, the randoms can be optionally smoothed (with fixed 7x7 moving mean smoothing) and/or perform variance reduction on the randoms data. These are performed before the corrections are applied to the measurement data/before reconstruction.
Scatter correction
Scatter correction is not inherently provided in OMEGA. However, since you can extract scatter data from GATE simulations, you can use that data for scatter correction. If scatter correction is selected, the user will be prompted for the scatter data (scn or mat file). Data in other format could be converted with loadMeasurementData (data is saved to options.SinM) or imported with importData.
Like randoms, scatter data can also be similarly smoothed and undergoe variance reduction. Scatter data can also be normalized if normalization correction is enabled and normalize_scatter is selected. Furthermore, scatter data can also be included as a part of the system matrix (as a diagonal matrix with the diagonal containing the scatter data). This can be enabled by setting options.subtract_scatter = false.
Attenuation correction
Attenuation is achieved by using attenuation images scaled for 511 keV. CT and 122 keV attenuation (e.g. germanium phantom) can be converted to 511 keV attenuation images with attenuationCT_to_511.m and attenuation122_to_511.m, respectively. To use attenuation correction, the attenuation_correction needs to be set to true as well as provide the name (and if not present on MATLAB path also the full path) of the mat-file containing the attenuation images to attenuation_datafile.
An experimental CT DICOM image to PET attenuation map is available with the function create_atten_matrix_CT.m. For converting of other formats it is recommended to try to use importData.
Attenuation, unlike all other corrections, is always performed during the reconstruction. I.e. when computing the system matrix.
Normalization correction
The normalization correction is component based and contains the following components: axial block profile factors, axial geometric factors, intrinsic detector efficiencies, transaxial block profile, and transaxial geometric factors. These are grouped into four different groups, with axial block profile factors and axial geometric factors being one, detector efficiencies one, transaxial block profile one and transaxial geometric factors one. Each or any of these can be deselected. The corresponding function normalization_coefficients.m can also output any of the specific normalization coefficients.
First output is the normalization coefficient matrix, second the normalized input data (sinogram or raw data file), third the axial geometric coefficients, fourth axial block profiles, fifth the transaxial block profile, sixth the detector efficiencies and seventh the transaxial geometric factors.
The process for sinogram data or raw data is slightly different. Sinogram data supports only fan-sum algorithm, while raw data also supports single-plane Casey (detector efficiency).
Transaxial geometric factors can be unreliable if small (non-FOV covering) sources are used.
Arc correction
Arc correction is currently an experimental feature. It is not recommended except with Inveon data. Arc correction also currently only works with sinogram data.
Arc correction is simply enabled by setting options.arc_correction = true. Increased sampling is supported.
Arc correction automatically uses parallel computing toolbox (parfor) if it is available.
Dynamic imaging in OMEGA means a time-series of images. For input data that means a time-series of sinograms/raw data contained in cells. For reconstructed images, the results are a time-series of the reconstructed images.
For dynamic imaging the start time, end time, total time and the number of time steps can be specified. The amount of time per time-step is currently fixed; variable rate is not supported at the moment.
No inherent TAC or Patlak plotting is available (currently).
For details on TOF, see Using TOF data.
Setting the precompute option to true, enables precomputation when the gate_main.m, Inveon_PET_main.m, Biograph_mCT_main.m, Biograph_Vision_main.m or main_PET.m is run. The precomputation phase includes the computation of the number of voxels that each LOR traverses. This allows for more efficient computation as LORs which are not inside the FOV can be ignored. This is especially true when using raw data. With sinogram data you may not get any benefit from precomputation, depending on the number of views (precomputed data might actually be slower, especially on Octave). This is highly recommended when using implementation 1 as without precomputation it is NOT parallelized. With precomputation the memory for the system matrix can be preallocated and, as such, the matrix creation can be parallelized. To actually use the precomputed data, you need to set options.precompute_lor = true.
You can load ONLY GATE/list-mode data by selecting options.only_sinos = true. In this case, running the main-file will load only the selected data type and will not perform any reconstructions. Normalization coefficients, however, are computed if they were selected (options.compute_normalization). Not applicable to main_PET.m.
You can skip ONLY data load by setting options.no_data_load = true. All other steps, e.g. normalization coefficient calculations, are performed. Not applicable to main_PET.m.
If implementation 1 is used, the entire system matrix (no subsets) can be computed by setting precompute_obs_matrix. This will most likely use significant amount of memory, however, as the matrix is output without any reduction to its size. This is the only way to compute MLEM with implementation 1.
You can perform ONLY image reconstruction step by setting options.only_reconstructions = true. This, however, requires an earlier measurement data load by OMEGA. The sinogram/raw data file specified by options.name and and the type specified by options.use_machine are automatically loaded as long as the files exist. For example, you could do a first pass with options.only_sinos = true and load GATE/list-mode data. Then you can perform any number of reconstructions with the same data by setting options.only_reconstructions = true and skipping the data load phase completely. This simply uses the previously loaded data. You cannot, however, make modifications to the sinogram (and most likely scanner) properties without doing the data load again. Not applicable to main_PET.m.
You can disable most of the status messages by setting options.verbose = false. They are enabled by default in all main-files.
Implementations
Four different implementations can be selected.
Implementation 1
Implementation 1 computes the system matrix as a sparse (MATLAB) matrix. This matrix is then used in MATLAB to compute the selected algorithm(s). This matrix can be either the full system matrix (not recommended) or a subset. Full system matrix is obtained by setting options.precompute_obs_matrix = true in misc properties.
If options.precompute_lor = false (misc properties), the system matrix is computed such that there is only partial memory preallocation and as such there is no parallellization. If options.precompute_lor = true, then the matrix can be preallocated after the precomputation phase and thus allowing parallellization and more efficient memory usage. If you use implementation 1, it is highly recommended to do the precomputation phase (options.precompute_lor = true).
This implementation uses the most amount of memory and is also, most likely, the slowest.
Implementation 2
This is a matrix-free reconstruction method. Implementation 2 uses both OpenCL and ArrayFire and as such you need to have both OpenCL and ArrayFire libraries installed and on your path. Implementation 2 supports all algorithms and priors that implementation 1 supports.
As with implementation 1, you can run 2 with both precomputation on (options.precompute_lor = true) or off (options.precompute_lor = false). The difference in implementation 2, however, is much smaller. Precomputation only leads to somewhat faster reconstruction, especially with raw data, but is otherwise identical in performance.
Unlike implementation 1, 2 has some additional properties you can set. First is the device used (options.use_device). Default is device 0 and, depending on the OpenCL supported devices, might also be the only device available. This is often also a GPU. You can query the available device numbers with ArrayFire_OpenCL_device_info(). Any devices shown with the aforementioned function can be used, though devices with less than 2GB of memory are not recommended.
Second is the usage of 64-bit atomics (options.use_64bit_atomics). This is on by default and is recommended when using GPUs. If you use CPUs, it is recommended to turn this to false as 64-bit atomics are not supported on CPU and as such cause slight delay in the kernel creation process.
Lastly is the ability to use CUDA (options.use_CUDA), that uses CUDA code instead of OpenCL code. This is recommended only for improved Siddon as other projectors tend to be slower than on OpenCL. CUDA support is also experimental and should not use 64-bit atomics.
Implementation 3
This is a matrix-free reconstruction method. Implementation 3 is a pure OpenCL method, meaning that ArrayFire libraries are not required. Implementation 3 also supports multi-device (heterogeneous or multi-GPU) computing. Only OSEM and MLEM are supported (though all projectors are supported). Precomputation works as with implementation 2.
Similarily to implementation 2, you can set the device used with the same parameter (options.use_device), however, unlike implementation 2 you do not select a single device, but rather a platform. Platform contains all the supported devices from the same vendor. You can view available platforms with OpenCL_device_info(). Some computing devices (especially CPUs) can be in multiple platforms. Selecting a platform will, by default, use all devices available on that platform. E.g. if you both a GPU and a CPU on the same platform, then both will be used. If you have two GPUs from the same vendor, both will be used, etc. Multi-device computing from different vendors are not supported (e.g. you can’t use both an AMD and a Nvidia GPU at the same time). In multi-GPU/device case, devices with less than 2GB memory are ignored (not used).
The amount of data distributed between the CPU and GPU in heterogeneous computing can be adjusted with options.cpu_to_gpu_factor. E.g. if options.cpu_to_gpu_factor = 2.5 then 2.5 times more data is given to the GPU. Alternatively, if options.cpu_to_gpu_factor = 0, then in multi-device platform ONLY the GPU with the highest amount of memory is used.
Implementation 4
This is a matrix-free reconstruction method. Implementation 4 is a pure CPU implementation using OpenMP for parallellization. It behaves similarly to implementations 2 and 3, except that OpenCL is not required and double precision (64-bit) values are used. All algorithms except MRAMLA and MBSREM are supported. Though only one subset-based algorithm and one MLEM method can be used at the same time (with one prior). E.g. you can use MLEM and OSL-OSEM with NLM, but you can’t use MLEM with OSEM and OSL-OSEM, or OSL-OSEM with MRP and NLM.
Precomputation works just as with implementations 2 and 3.
There are no additional parameters for implementation 4. All cores are always used and sometimes all threads as well.
Projectors
Three different projectors can be selected, the improved Siddon’s ray tracer, the orthogonal distance-based ray tracer or the volume of intersection based ray tracer. All projectors can also utilize PSF reconstruction.
First is selected by setting options.projector_type = 1, which is also the default. Second with options.projector_type = 2 and third with options.projector_type = 3. Implementation 1, with no precomputation, supports also the regular Siddon’s algorithm (options.projector_type = 0), but it is not recommended.
Improved Siddon can be used with any number of rays if no precomputation is used (options.precompute_lor = false). These can be set with options.n_rays_transaxial and options.n_rays_axial. More than 1 ray is not supported with implementation 1.
Orthogonal distance-based ray tracer can be used in either 2.5D or 3D mode. 2.5D mode is enabled by setting options.tube_width_z = 0 and options.tube_width_xy to greater than zero. 3D mode is enabled if options.tube_width_z is nonzero. In 3D-mode options.tube_width_xy is ignored.
Volume of intersection based ray tracer allows you to specify the radius of the tube of response (cylinder) with options.tube_radius and the relative size of the spherical voxels. This projector approximates the voxels as spheres and with options.voxel_radius you can specify their relative size. If options.voxel_radius = 1 then the spherical voxels are just large enough to contain the entire regular voxel inside them. Lower values cause some of the original voxel to go outside the sphere while larger values cause more space between the original voxel and the sphere boundary.
PSF (point spread function) reconstruction can be enabled for any projector. options.use_psf controls whether PSF is included or not. The full width at half maximum (FWHM) of the PSF can be controlled with options.FWHM for each of the three dimensions separately. Furthermore, an optional deblurring step can be enabled for PSF reconstruction with options.use_deblurring. The total number of deblurring iterations is the same as the number of subsets times the number of iterations. The deblurring phase is performed only after the reconstruction is completed and is performed for all iterations. The function gaussianKernel is used to form the Gaussian PSF.
Reconstruction settings
Subset type can be selected with options.subset_type.
Initial value can be set to options.x0. Dimensions should match with the values in image properties.
options.epps is a small value that prevents division by zero. Shouldn’t be need to adjust.
Misc settings
Setting options.use_Shuffle = true uses the file exchange function shuffle, when options.subset_type = 3 (subsets selected randomly). Using this is optional, but recommended (reduces memory usage). This function needs to be downloaded from the file exchange: https://se.mathworks.com/matlabcentral/fileexchange/27076-shuffle
Fast sparse can be enabled with options.use_fsparse = true. This is also optional, but speeds up sparse matrix generation (no longer applicable in MATLAB R2020a and up). Download from https://github.com/stefanengblom/stenglib
When computing any median root prior based priors, the prior is normally computed as (x - median(x)) / median(x). Setting options.med_no_norm = true changes this such that the prior is computed as x - median(x). This can lead to improved image quality, but is NOT the way it is presented in the literature.
Reconstruction algorithms
Select the desired algorithms by setting them to true. Deselect them by setting them to false.
Algorithm properties
Adjust the various algorithm properties here.
Here you can input custom detector coordinates. These enable the use of non-cylindrical devices or continuous detectors (e.g. CPET in GATE). Transaxial coordinates are specified as options.x and options.y while axial are options.z (or alternatively options.z_det). Simply save your custom coordinates to the specified x, y and z variables in the options struct. For sinogram data, the x and y sizes should be [sinogram_elements 2] and z [number_of_sinograms 2]. For raw data [detectors_per_ring 1] and [number_of_rings 1], respectively.
It is also possible to use pure event-by-event list-mode reconstruction with custom coordinates. In this case each event should have its own detector coordinates and as such each coordinate should be of size [total_events 2]. Currently only implementation 4 supports event-by-event reconstruction and it is also ONLY supported with custom detectors.
When using main_PET.m you will be prompted for the input data when you either run the script, run the section or run the loadMeasurementData function.
The input data can be either a mat-file (raw data supports ONLY mat files), NIfTI, Analyze, Interfile, DICOM or MetaImage file. The file type is detected from the file extension. .hdr, .img. .nii and .gz are considered either NIfTI or Analyze format images, .i33 and .h33 Interfile, .dcm DICOM, .mhd and .mha MetaImage and all others as raw data.
DICOM support requires image processing toolbox. NIfTI either image processing toolbox or "Tools for NIfTI and ANALYZE image" toolbox from MathWorks file exchange https://se.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image
Sinogram data is loaded from options.SinM variable and raw data from options.coincidences.
Always check the correct orientation of the input data before proceeding with reconstruction. E.g. for sinogram data the input data should have the dimensions of [options.Ndist options.Nang options.TotSinos].