Reconstruction of 2D PET data with Monte Carlo generated natural pixels

In discrete detector PET, natural pixels are image basis functions calculated from responses of detector pairs. By using reconstruction with natural pixels the discretization of the object into a predefined grid can be avoided. Instead of trying to find the object x from system matrix A and measured data p by solving Ax = p, one solves Mq = AA^Tq = p. The backprojection of the final q results in the estimated image of the object x. In previous work natural pixels were calculated by intersections of strip functions or combinations of strip functions. However, in PET systems the detector response is not correctly described by a strip function. In PET the detector response is better described by a Gaussian tube. It is quite difficult to calculate the intersections of all possible lines of response (LOR) in the scanner. This paper proposes an easy and efficient way to generate the matrix M directly by Monte Carlo simulation. Using the natural pixel matrix M has several advantages over using the conventional system matrix A. The discretisation of the object into pixels, which requires an angular dependent projector and backprojector, is avoided. Due to rotational symmetry in the PET scanner the matrix M is block circulant and only the first block row needs to be stored. The block circulant property of the matrix allows use of fast direct methods to calculate the solution. Iterative methods to solve the system are preferred because of numerical stability.