Toolbox overview

Table of Contents

• Toolbox overview
• Technical aspects
  • Precomputation phase
    • Purpose
    • Execution
  • Data load
    • Execution

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

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

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 and is intended for MATLAB version 2018b and earlier, the second uses the new C-MEX interface and is for MATLAB 2019a and newer and lastly there is dedicated version for Octave as well.