CT projection (distance-driven) operator implementation?

Question

I am trying to use MATLAB to implement a CT (computed tomography) projection operator, A, which I think is also referred as "system matrix" often times.

Basically, for a N x N image M, the projection data, P, can be obtained by multiplication of the project operator to the image:

P = AM

and the backprojection procedure can be performed by multiplying the (conjugate) transpose of the projection operator to the projection data:

M = A'P

Anyone has any idea/example/sample code on how to implement matrix A (for example: Radon transform)? I would really like to start with a small size of matrix, say 8 x 8, or 16 x 16, if possible.

My question really is: how to implement the projection operator, such that by multiplying the operator with an image, I can get the projections, and by multiplying the (conjugate) transpose of the operator with the projections, I can get the original image back.

EDIT:

Particularly, I would like to implement distance-driven projector, in which case beam trajectory (parallel, fan, or etc) would not matter. Very simple example (MATLAB preferred) will be the best for me to start.

Answer 1

As far as I'm aware, there are no freely available implementations of the distance-driven projector/backprojector (it is patented). You can, however, code it yourself without too much difficulty.

Start by reading the papers and understanding what the projector is doing. There are only a few key parts that you need:

• Projecting pixel boundaries onto an axis.
• Projecting detector boundaries onto an axis.
• The overlap kernel.

The first two are simple geometry. The overlap kernel is described in good detail (and mostly usable pseudocode) in the papers.

Note that you won't wind up with an actual matrix that does the projection. The system would be too large for all but the tiniest examples. Instead, you should write a function that implements the linear operator corresponding to distance-driven projection.

Answer 2

Distance-driven projection is not implemented in stock MATLAB. For forward projection, there is the fanbeam() and radon() command, depending on what geometry you're looking for. I don't consider fanbeam to be very good. It exhibits nonlinear behavior, as of R2013a, see here for details

As for a matching transpose, there is no function for that either for fanbeam or parallel geometry. Note, iradon and ifanbeam are not operator implementations of the matching transpose. However, you might consider using FUNC2MAT. It will let you convert any linear operator from function form to matrix form and then you can transpose freely.

Answer 3

You have different examples:

• Here there is a Matlab example related to 3d Cone Beam. It can be a good starting point.
• Here you also have another operator
• Here you have a brief explanation of the Distance-Driven Method. So using the first example and the explanation in this book, you can obtain some ideas.

If not, you can always go to the Distance-Driven operator paper and implement it using the first example.

Answer 4

Although there are already a lot of satisfactory answers, I would like to mention that I have implemented the Distance Driven method for 2D Computed Tomography (CT) and 3D Digital Breast Tomosynthesis (DBT) on MATLAB.

Until now, for 2D CT, these codes are available:

• Simple Distance-Driven, base on the original papers [1] and [2],

• Branchless Distance-Driven, for acceleration on GPU, based on the papers [3] and [4],

and for 3D DBT:

• Simple Distance-Driven, based on the book [5].

Note that:

1 - The code for DBT is strictly for limited angle tomography; however it is straightforward to extend to a full rotation angle.

2 - All the codes are implemented for CPU.

Please, report any issue on the codes so we can keep improving it.