Skip to content

Computing the forward and or backward projections

Jump to bottom
villekf edited this page Mar 4, 2021 · 9 revisions

This page details on how to separately compute the forward (Ax) and/or backward projections (A^T y) when using OMEGA.

Computing the forward and/or backward projections

PET data

When you are manually computing the forward and/or backward projections, it is recommended to use the forward_backward_projections_example.m and modify it according to your specifications.

The initial part works very much like main_PET.m, however, there are no adjustments for the reconstruction algorithms. You still need to input the necessary information (scanner properties, sinogram/raw data properties, corrections, etc.) before proceeding. For options.precompute_lor = true you also need to perform the precomputation phase if it hasn’t been performed before. Note that if you want the measurement data to be precorrected, you need to do that manually. Corrections during reconstruction are applied automatically, though in some cases (such as randoms correction) you need to select the correction data when it is prompted.

For the reconstruction, implementations 1, 3 and 4 are supported with 3 recommended.

After the initial setup phase, there are two different ways to actually compute Ax and A^Ty. The first, recommended, method is by using a MATLAB/Octave class that is demonstrated in the initial examples in forward_backward_projections_example.m, the second method is the old one using functions.

Class based examples

First is a class-based exampled computing OSEM/MLEM estimates (line 821) with the use of the special forward and backward projection operators. The major steps needed for class based use is the creation of the class object with forwardBackwardProject(options) and the load of the measurement data. Sensitivity image can be obtained from the backprojection. PSF is applied automatically if selected.

You can also create the same class object as above without the need for any of the main-files, however, several features are then unavailable (such as corrections). This is achieved by inputting the required variables manually: A = forwardBackwardProject(Nx, Ny, Nz, FOVa_x, axial_fov, diameter, rings, cr_p, cr_pz, cryst_per_block, blocks_per_ring, det_per_ring, Ndist, Nang,NSinos, span, ring_difference, subsets, implementation, use_device, projector_type, tube_width_z, tube_radius, voxel_radius, use_psf, FWHM);

To use your own measurement data, replace load('Cylindrical_PET_example_cylpet_example_sinograms_combined_static_200x168x703_span3.mat','raw_SinM') with your own custom load or simply run the code up to that and manually load your data (you can use loadMeasurementData to load custom data, such as interfiles). You should either name your measurement data as raw_SinM or rename all instances of raw_SinM to the same name as your input measurement data (options.SinM if you use loadMeasurementData). Your input data needs to be oriented the same way as the measurement data used in OMEGA. For sinograms that is thus options.Ndist x options.Nang x options.NSinos. Dynamic data needs to be manually correctly formatted.

The second class-based example computing Conjugate-gradient Least-squares (CGLS) utilizes the operator overloading. This means that after you have created the class (e.g. with A = forwardBackwardProject(options);) you can then simply multiply this with the desired vector (e.g. A * x or A * y). Whether forward or backward projection is computed, is determined automatically from the dimensions of the input vector. Alternative, backprojection can be forced with A' * y.

Other examples include LSQR, SA-WLS (includes TOF example), ATP-WLS, EM-PKMA with TV regularization and ADMM with TV. Of these, EM-PKMA and ADMM provide examples on how to use the built-in priors. Prepass functions must be run for the following priors:

MRP (options = computeOffsets(options);)

Quadratic prior (options = computeWeights(options); and options = quadWeights(options, options.empty_weight);)

Huber prior (options = computeWeights(options); and options = huberWeights(options);)

L-filter (options = computeWeights(options);, options = computeOffsets(options); and options.a_L = lfilter_weights(options.Ndx, options.Ndy, options.Ndz, options.FOVa_x/options.Nx, options.FOVa_y/options.Ny, options.axial_fov/options.Nz, options.oneD_weights);)

FMH (options = computeWeights(options);, options = computeOffsets(options); and options = fmhWeights(options);)

Weighted mean (options = computeWeights(options); and options = weightedWeights(options);)

TV (options = TVPrepass(options);, the following ONLY when options.TVtype = 3: options = computeWeights(options);, options = computeOffsets(options); and options = quadWeights(options, options.empty_weight);)

APLS (options = APLSPrepass(options);)

NLM (options = NLMPrepass(options);).

Furthermore, you can also obtain the (subset of the) system matrix from the same class object. This is demonstrated in the System matrix (OSEM) example between SA-WLS and ATP-WLS.

Function based examples

This is the very last example script in forward_backward_projections_example.m. It is no longer recommended to compute the forward and/or backward projections this way.

You’ll need to run the index_maker.m as specified to get the necessary subset indices. The line [gaussK, options] = PSFKernel(options); is only applicable if PSF reconstruction is used as it creates the PSF kernel gaussK.

To use your own measurement data, replace load('Cylindrical_PET_example_cylpet_example_sinograms_combined_static_200x168x703_span3.mat','raw_SinM') with your own custom load or simply run the code up to that and manually load your data (you can use loadMeasurementData to load custom data, such as interfiles). You should either name your measurement data as raw_SinM or rename all instances of raw_SinM to the same name as your input measurement data (options.SinM if you use loadMeasurementData). Your input data needs to be oriented the same way as the measurement data used in OMEGA. For sinograms that is thus options.Ndist x options.Nang x options.NSinos. Dynamic data needs to be manually correctly formatted.

After the above steps, simply use the functions forward_project and backproject as you like. The example provided uses them the same way as in OSEM, but as they simply implement the operations y = Ax and x = A^Tb, respectively, you can implement any algorithm that uses these matrix-vector operations. An example of PSF reconstruction is included as well. Sensitivity image is obtained with backprojection, but can also be omitted by simply having only one output parameter.

All selected corrections are applied to the forward and/or backward projections. If you wish to use different projectors for the forward or backward projections, you need to manually set options.projector_type between them. E.g. compute forward projection when options.projector_type = 2 and then set options.projector_type = 1 and compute backprojection. This, however, is untested and may especially not work when using implementation 1. If you wish to use different projectors, it is recommended to use implementation 3 and set options.precompute_lor = false.

The optional PSF deblurring phase can also be performed and an example of that is also included.

List-mode reconstruction/custom detectors

There is also a "simplified" version of the forward/backward projection class. This simplified version accepts the detector coordinates for each coincidence as well as the FOV size and image size. The input structure is thus:

A = forwardBackwardProject(x, y, z, Nx, Ny, Nz, FOVtr, axial_fov, subsets, implementation, platform)

where x, y and z are the detector coordinates corresponding to every element of the input measurement data (vector or matrix), Nx/y/z are the image size in x/y/z-directions, FOVtr the FOV size in transaxial direction (length of one side of the square), axial_fov the size of the FOV in axial direction, subsets is an optional number of subsets (can be omitted), implementation the optional implementation that should be used (4 is default, can be omitted) and platform the device platform for implementation 3 (default is 0, can be omitted, ignored when implementation 4 is used). Implementations 3 and 4 are currently supported.

When used like this, however, no corrections or PSF are supported. All these operations need to be done manually by the user. An example of this type of reconstruction is presented in Inveon_PET_main_listmode_example.m.

CT data

Clone this wiki locally