Skip to content

Computing the forward and or backward projections

Jump to bottom
villekf edited this page Mar 23, 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

For CT data, only the class based functionality is supported. Examples are available in walnut2D_CT_main.m, walnut_CT_main.m and Inveon_CT_main.m, where the first one is a pure 2D (fan-beam) example while the latter two are 3D (cone-beam) examples (it is possible to convert the latter two to 2D as well). In all cases, it is assumed that the source is initially at the bottom of the image and rotates counter-clockwise.

There are three different ways to construct the class object. The recommended method is to use one of the example main-files as a basis and construct the class object by using the options struct: A = forwardBackwardProjectCT(options). Both the improved Siddon and volume-based ray tracers are available with this method as well as automatic PSF if selected. With this method, it is also possible to easily affect the initial position of the source by using the options.offangle value.

Alternatively, it is possible to construct the class object without the options struct, either by using the projection angles or by inputting the detector coordinates for every measurement point. In the first case (angles), the constructor is A = forwardBackwardProjectCT(angles, ySize, xSize, Nx, Ny, Nz, FOVa_x, axial_fov, dPitch, sourceToDetector, sourceToCRot, horizontalOffset, bedOffset, subsets, implementation, deviceN). angles are the angles corresponding to each projection, ySize is the number of rows in a single projection image, xSize the number of columns, Nx, Ny, Nz are dimensions of the estimated image, FOVa_x is the size of the transaxial FOV, axial_fov is the size of the axial FOV, dPitch is the crystal pitch (size), sourceToDetector is the distance from source to the detector (mm), sourceToCRot the distance from source to the center of rotation (mm), horizontalOffset is an optional value that can be used to specify possible offset of the source from the center of rotation, bedOffset is an optional value that describes the offsets of each bed position, subsets is the optional number of subsets, implementation is the used implementation (3, OpenCL, and 4, OpenMP, are supported), deviceN is the platform used when using implementation 3.

The alternative way, is to use user-computed detector coordinates. The constructor is then: A = forwardBackwardProjectCT(x,y,z, Nx, Ny, Nz, FOVa_x, axial_fov, subsets, implementation, deviceN), where x, y and z are the detector coordinates in the transaxial and axial directions.

In each of the three cases, the class object can simply be used as A * x to compute forward projection and A' * y to compute backprojection. Subsets can be set with A.subset = subset_number. Alternatively, you can use the class object in function: forwardProject(A, f, subset_number) and backwardProject(A, y, subset_number).

Extracting the system matrix

The class object can also be used to form the system matrix itself, either the entire system matrix or a subset. However, to form the system matrix, the class object has to be created with implementation 1 (e.g. options.implementation = 1). The matrix itself is obtained with B = formMatrix(A, subset_number), where subset_number can be omitted if the entire system matrix is required.

Clone this wiki locally