PET reconstruction system matrix

DLL to generate a PET system matrix for image reconstruction

Description

This include a MEX-DLL file which is an implementation of R L Siddon's refined method by Jacob for the calculation of ray path integral in PET/CT system image reconstruction. This is based on Line Intersection Length model to calculate probability matrix elements. Inputs are the LOR number and a structure defining the system geometry such as pixel numbers and widths. Help is included on use in a file SYS_MAT.m or the CPP file.

• MATLAB release: MATLAB 7.0.1 (R14SP1)
• Other requirements: Created for Windows and may require Image Processing toolbox commands in *.m file.

[Munir Ahmad]

I would like to detial the use of this function a bit. This function is of a first hand implementation of the method to calculate system amtrix elements in a PET system. It make certain assumptions at this time evalute the matrix elements. It assumes that the detector does not extend out of the reconstruction circle as it cannot deal with the zero matri value (an error of division by zero may occur), though it is not difficult to modify the .CPP file to add this functionality. I will add this functionality in the next version. Also, another method is to calculate intersected area of the projection ray and the pixel, that function will also follow this one later. It is faster if you assume detector width equal to that of the reonstruction pixel width. It works in 2D only at the moment.

[Munir Ahmad]

Please note that this function can be used in iterative reconstruction methods, such as maximum likelihood (ML EM) to evaluate system matrix elements.

Please contact me in case you have any questions at [e-mail]

[Siarhei]

Dear Munir,
Could you please clarify what all parameters do mean?

In case if it is needed to re scale system matrix, what should be changed: i.e. physical volume of image should be constant. Therefore we can change, h(v)pixels -number of pixels in image dimensions, or physical pixel sizes.

[Siarhei]

Hello,
I think I have found a bug in code.

the left down corner of the object should have coordinates:

bx = - (hp - 1)*dx/ 2.0;
by= - (vp - 1)*dy/ 2.0;

instead of

bx = - (hp/2.0 - dx/ 2.0);
b = - (vp/2.0 - dy/ 2.0);

[hamed moradi]