Skip to content

Toolbox overview

Jump to bottom
villekf edited this page Jun 5, 2020 · 20 revisions

This page explains a general overview on the structure of the toolbox as well as details on how it technically works.

Toolbox overview

The general structure of OMEGA can be divided into three different layers. The top layer is the MATLAB/Octave user-interface that contains the scripts and functions necessary to call the lower layers. The middle layer is the MATLAB © MEX-interface that calls and computes the C++ code and then sends it back to the top layer. The bottom layer, which is not always used, contains the OpenCL kernels that compute the OpenCL code and then send the output data (reconstructed images) to the middle layer. The bottom layer is only used if OpenCL code is used. The middle layer can also be ignored, but this is not recommended and not enabled by default as the pure MATLAB/Octave implementation is extremely slow.

For the user, only the top layer is exposed. This is achieved in the form of scripts that are referred to as main-files. These include gate_main.m, main_PET.m, Inveon_PET_main.m, Biograph_mCT_main.m, Biograph_Vision_main.m, custom_prior_test_main.m, and forward_backward_projections_example.m. It is from these main-files that the actual functions are called. This is achieved by storing all the user selected parameters to a MATLAB/Octave struct called options, which is then input to the functions (e.g. when loading data or performing image reconstruction).

In the main-files most parameters are set either with numerical values (e.g. blocks_per_ring, span, tube_radius) or as a true/false selection, where true means that the feature is included and false that it is omitted (e.g. randoms_correction, scatter_smoothing, osem). The main-files are divided into "sections" with specific labels (MACHINE PROPERTIES, IMAGE PROPERTIES, SINOGRAM PROPERTIES, etc.) that control different aspects. Many of these sections are completely optional, e.g. CORRECTIONS section can be omitted if the user does not wish to use any corrections. The only compulsory ones are MACHINE PROPERTIES and either SINOGRAM or RAW DATA PROPERTIES. For reconstruction it is also advisable to inspect the RECONSTRUCTION PROPERTIES section, but the default values should always output a working OSEM estimate. By default all non-compulsory options are set as false (with the exception of the Inveon main-file which has several corrections enabled by default).

There are also several functions that work very independently without the need for the main-files. These include file import and export functions, visualization functions and many reconstruction algorithm functions. For help on many of these functions, you should use help function_name or alternatively doc function_name. E.g. help saveImage.

Technical aspects

This section explains how, exactly, does each component of OMEGA software work.

Precomputation phase

The corresponding m-file is lor_pixel_count_prepass.m.

Purpose

As mentioned in other help pages, the goal of the precomputation phase is to determine the number of voxels that each line of response traverses (interacts with). Line of responses that do not intersect the FOV at all will have 0 interacted voxels and as such can be safely removed. This step is also necessary for the parallel version of implementation 1 as it requires the preallocation of the system matrix in memory. This also means that, unlike with matrix free implementations, the orthogonal and volume of intersection ray tracers require their own precomputation passes since the number of voxels that are interacted with is much higher.

Execution

Dedicated codes are available for the precomputation phase. Furthermore, since the OpenCL implementations use single (32-bit float) precision numbers and the C++ versions double (64-bit double) precision numbers, different versions are available for both. This is because of floating-point rounding effects that can cause LORs that are on the very edge of a voxel to be in different voxels depending on whether the single or double precision implementations are used. While this usually affects only a very small portion of all the LORs, even a single wrong one can cause a crash due to out-of-bounds variables. Implementation 2 also has its own precomputation code, although the OpenCL kernel itself is identical to implementation 3. This is simply because of the different device/platform executions.

The codes themselves are stripped versions of the actual projectors. In all cases, the ray tracing itself is always performed, no parameters outside of the ones necessary for the ray tracing are computed. For C++ version, there is a different file/function for this (improved_Siddon_algorithm_discard.cpp), but for OpenCL versions all the code are in "master" kernel file (multidevice_kernel.cl) where the correct lines are compiled during runtime by taking advantage of preprocessor directives (#ifdef, #else, etc.). For precomputation phase, the preprocessor directive that is defined is FIND_LORS. No precomputation code is available for CUDA.

Data load

The corresponding m-files are load_data.m, load_data_mCT.m and load_data_Vision.m.

Execution

Data load is divided into three different sections: GATE data, list-mode data and sinogram data.

GATE data

For GATE data, the data import is separate for LMF, ASCII and ROOT.

In LMF, the data import is done in a C++ file (gate_lmf_matlab.cpp) and each binary packet is read sequentially. In LMF, unlike other methods, the coincidences are formed manually from the singles. This is done by first checking if the time difference between the consecutive singles is within the coincidence window length. If yes, then the singles are assumed to be coincidences and all the specified values (e.g. detector number, whether the coincidence is true, random or scattered event, source coordinates) are saved. If not, then the first single is discarded and the latter single becomes the first single and compared with the next single, etc. For time steps, there is first a check on whether the current time is larger than the time of the next time step. If yes, then the current event index is stored. This index indicates where the time was large enough for the next time step to begin and is then used to separate the measurement data into the time bins. Of all the three data types, LMF is the most limited one as it does not support delayed coincidences or scattered events except for Compton scatter in phantom.

ASCII data import is different from both LMF and ROOT in sense that it is done purely in MATLAB/Octave. It uses either readmatrix (if available) or importdata to import the ASCII text into a matrix containing all the rows and columns. If any of the columns on any of the rows is corrupted/missing (replaced by NaN value), the code will detect these, discard them and inform the line number where this occurred. All the chosen variables are then stored in individual vectors/matrices. Time steps are handled similarly to LMF.

ROOT data import is handled in a C++ MEX-file. Unlike LMF, ROOT has three different MEX-versions. One uses the "traditional" C MEX-interface (GATE_root_matlab_C.cpp) and is intended for MATLAB version 2018b and earlier, the second uses the new C++ MEX-interface (GATE_root_matlab.cpp) and is for MATLAB 2019a and newer and lastly there is a dedicated version for Octave as well (GATE_root_matlab_oct.cpp). All the ROOT import functions first open the trees (coincidences and delays if selected) and then the desired branches. There are also checks in place that first guarantee that a specific branch is available, e.g. the scatter data. If not a message is displayed, but the data load will still continue, unless it is the detector information that is missing as it is vital for the data load.

List-mode data

List-mode data, more specifically Siemens Inveon, Biograph mCT and Biograph Vision list-mode data, is loaded in a separate MEX-file. For Inveon, the source code is available (inveon_list2matlab.cpp), but for mCT and Vision only the MEX-files themselves are distributed (i.e. a closed source release). For Inveon, the code loops through all the bit-packets, determines whether they are prompt, delay or time tags and then extracts the corresponding information. Static and dynamic cases are handled a bit differently; in the former case the counts are stored in one detectors x detectors sized matrix, while in the latter they are stored event-by-event basis. In dynamic case, however, the events are stored in a same type of (sparse) matrix as in static case for each time step with the use of accumarray function. Time steps are handled as with LMF data, where the index is stored where the time exceeded the previous time step.

List-mode data is saved with _listmode in the end of the filename. For 32-bit list-mode data (mCT only) _listmode_sinogram is added to the end of the filename.

Sinogram data

This applies only to Biograph mCT and Vision. For Inveon, the sinogram data load occurs in form_sinograms.m. Uncompressed mCT and Vision sinograms (.s or .ptd files) can be loaded, but TOF data is currently lost, i.e. all the TOF bins are summed together. Corrections can be applied normally. The data itself is saved with _machine_sinogram in the end of the filename.

Saving data

GATE data and list-mode data go through the same procedures when saving data. All steps are repeated for the selected number of time steps, where first the sinogram is created (if sinogram data is used) and then the raw data is stored. For GATE data, trues, randoms and scatter are stored as well if selected.

Forming sinograms

The corresponding m-file is form_sinograms.m.

Execution

Sinograms can be formed from saved raw data, during data load (no need to load the raw data separately) and also by simply modifying the corrections applied to the sinogram (e.g. no actual new sinogram is created). When sinograms are formed, a raw uncorrected sinogram is always created regardless of the corrections applied. This is saved as raw_SinM.

The first step of sinogram creation is the formation of an "initial Michelogram". This is an intermediate step between the raw data format and the Michelogram/sinogram format. The raw data is divided into vectors that contain the future Michelogram bins. This is performed in initial_michelogram.m.

Next step is the formation of the Michelograms by selecting the data points that are within the predetermined orthogonal distance from the center of the field-of-view. These are saved as unsigned 16-bit integers and performed for all the selected data types (trues, prompts, delays, etc.).

After this, the next step performs the axial compression, though using span of 1 (no axial compression) is also possible. However, span of 1 is only supported with prompts.

Last step are the corrections if any have been selected and when options.corrections_during_reconstruction = false. Corrections are handled in the following order: Randoms (variance reduction, then smoothing) -→ Scatter without normalization (variance reduction, then smoothing) -→ normalization correction -→ Scatter when using normalized scatter (variance reduction, then smoothing) -→ global correction factor -→ Sinogram gap filling. If any of the corrections are set as false, then that step is omitted. Only prompts go through corrections. Scatter can be applied only with normalization separately applied to it or without separate normalization.

All the separate sinograms are saved in a same mat-file with the sinogram dimensions in the name. Included are also structs that contain whether certain corrections were applied (appliedCorrections) and what corrections were applied to scatter or randoms (ScatterProp, RandomsProp). In appliedCorrections normalization is stored as a boolean variable (false means no normalization), randoms and scatter as char (empty array means no corrections, otherwise they can be e.g. "randoms correction with smoothing"), gap filling as boolean, mashing factor as an integer and lastly the user specified global correction factor. The prop-structs contain booleans indicating whether variance reduction and/or smoothing was applied.

Randoms correction is applied as randoms subtraction from the delayed coincidences data. Scatter correction can be applied either as a subtraction by setting options.subtract_scatter = true, or alternatively by multiplication. In the latter case the scatter data is multiplied with the sinogram. For dynamic case, the delayed coincidences are stored as in the static case, but then divided by the amount of time steps.

When the function is used to modify the applied corrections (e.g. form_sinograms(options, true)), the sinogram creation step is skipped and the uncorrected sinogram is loaded. By default, form_sinograms assumes that the sinogram needs to be created, i.e. the boolean value after options needs to be true in order to perform only corrections. Any sinogram, no matter where created, can be corrected like this. However, the data needs to saved as raw_SinM in a mat-file with the same name as the current machine properties (e.g. [options.machine_name '_' options.name '_sinograms_combined_static_' num2str(options.Ndist) 'x' num2str(options.Nang) 'x' num2str(options.TotSinos) '_span' num2str(options.span) '.mat'] for static data and [options.machine_name '_' options.name '_sinograms_combined_' num2str(options.partitions) 'timepoints_for_total_of_ ' num2str(options.tot_time) 's_' num2str(options.Ndist) 'x' num2str(options.Nang) 'x' num2str(options.TotSinos) '_span' num2str(options.span) '.mat'] for dynamic).

In the bottom of the m-file, there is a separate section for loading Inveon Acquisition Workplace created sinograms. These sinograms automatically have randoms corrections applied. All other corrections can be applied just as with raw data.

Clone this wiki locally