Extracting the system matrix

This page details how you can extract only the system matrix (or a subset of it) and use it manually in your own algorithms/code.

Extracting the system matrix

Currently there are two ways to extract the system matrix. A class based method using similar methods as when computing the forward/backward projections or a function-based method that has been included since the first release. Currently the class based method is recommended. For CT data, only the class based method is supported. For more information on extracting the system matrix in CT, see bottom of page Computing the forward and/or backward projections.

Class based example

This is based on forward_backward_projections_example.m for PET and e.g. for CT. The initial part shows how to compute the forward and backward projections, but "System matrix (OSEM) example" shows how to extract the system matrix (or a subset of it) by using the same class. First you need to build the class object with e.g. forwardBackwardProject(options) for PET and forwardBackwardProjectCT(options) for CT and then you can obtain the system matrix simply with formMatrix(A) or a subset of it with formMatrix(A, subsetNumber). Note that if you are using precomputed data (options.precompute_lor = true) then this is the TRANSPOSE of the system matrix. For CT data, precomputed data is always used so the matrix is always a transpose of the system matrix.

Function based example

When extracting the system matrix with the function based method, it is recommended to use main_PET.m and follow its general workflow. This means that you should input the necessary parameters to the initial section (scanner properties, sinogram/raw data properties, etc.). It is also recommended to set options.only_system_matrix = true. Precompute step should be computed if you have set options.precompute_lor = true. The reconstruction step (reconstructions_main) should be ignored. CT data does not have a similar example, but should nevertheless work similarly.

After the reconstruction step, there is a specific section for the system matrix creation and (optionally) to separately compute any selected algorithms. It is important to run the custom_prior_prepass phase before the system matrix is created.

The subset of the system matrix can then be obtained with A = observation_matrix_formation(options, subsetNumber) or the entire system matrix with A = observation_matrix_formation(options). This matrix behaves as any other (sparse) matrix in MATLAB/Octave. Any selected corrections (attenuation, normalization, scatter) are applied to the system matrix automatically.

Since the system matrix creation uses implementation 1, it is HIGHLY recommended to set options.precompute_lor = true. Furthermore, other projectors than improved Siddon are not recommended. PSF blurring needs to be added manually, though you can use PSFKernel function to obtain the Gaussian convolution kernel.