1

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.

Amro
  • 123,847
  • 25
  • 243
  • 454
Nick X Tsui
  • 2,737
  • 6
  • 39
  • 73
  • 1
    Well, Matlab does matrix multiplication and transpose off-the-shelf...: `P = A*M`, `M = A.'*P` – Luis Mendo Jan 02 '14 at 21:54
  • No, I mean how to implement the projection operator, such that by multiplying the operator with an image, I can get the projection. – Nick X Tsui Jan 02 '14 at 22:05
  • 2
    The entry p(i,j) is the sensor data of sensor i if there was an image with pixel j set to 1 and all other pixels set to zero. You need to flatten the image and projection to vectors, so 8x8 gets 64x1. The projection matrix depends on your setup (parallel, fanbeam, helical scan, sensor grid, etc.), so you need to be more specific. – mars Jan 02 '14 at 23:00
  • @mars Hi Mars, you are absolutely right. I have no idea how the image should be arranged, or vectorized. I think for most of the times, the image is vectorized, but I am not sure in this case. For the projection pattern, I would say fan beam with equiangular detector setup. Do you have any idea how I should implement this? Thanks a lot in advance. – Nick X Tsui Jan 03 '14 at 17:45
  • Are you asking about [something like this](http://www.mathworks.co.uk/matlabcentral/fileexchange/34608-ct-reconstruction-package) for example, or did I misunderstand? – Roger Rowland Jan 10 '14 at 15:19
  • @RogerRowland Yes something like this. However, the projection method in this package (myRadon.m) did not employ the distance-driven methods, so it might have high-frequency artifacts involved. Also I don't think it actually implemented a projection operator whose conjugate transpose does the backprojection job. – Nick X Tsui Jan 13 '14 at 15:08
  • This question could use a small scale example. That way it is most clear what you have and what exactl you need. – Dennis Jaheruddin Jan 14 '14 at 15:00
  • 4
    This question is only for people who know what it is. If you don't, then you don't have to write any comments here. – Nick X Tsui Jan 14 '14 at 16:31
  • 1
    Note your notation that "$P = AM$" and "$M = A'P$" is incorrect, and in particular you won't get the original image back just by applying the transpose operator. In general matrix terms, if $P = AM$ and $Q = A'P$, then $Q = A' A M$. But $Q=M$ if and only if $A$ is an orthogonal matrix. In general projection is not such a matrix, so $Q \neq M$. People have spent entire careers working on reconstruction methods, essentially just to invert $A$. – Kevin Holt Jan 20 '16 at 21:45

4 Answers4

7

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.

Glorfindel
  • 21,988
  • 13
  • 81
  • 109
phyrox
  • 2,423
  • 15
  • 23
  • Thanks a lot. I am not expecting to find example code here, so this although does not give the direct answer, but I think it is very informational and closer, so I think I will take this one as the answer. – Nick X Tsui Jan 15 '14 at 17:19
  • 2
    Thanks man! (if you find other useful example code, please post it here so everybody will be able to use it) – phyrox Jan 16 '14 at 09:24
5

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.

lp251
  • 166
  • 3
  • Hi Luke, thanks a lot for your comments. I can certainly understand the projection of the boundaries (basically linear interpolation) of pixels or detectors onto some axis, but I am not sure about the overlap kernel. If there is a matrix (I know it is probably too large to implement in reality), what on earth the entries in the matrix? I suppose something like the length of the segment within each pixel? You mentioned the kernels are well described by some pseudo codes in some papers, can you provide some links? Thanks a lot! – Nick X Tsui Jan 20 '14 at 16:09
  • The matrix you need is the same one that arises for discretized integration. For the forward projection operation, the overlap kernel computes j-th detector value d_j as the sum of the pixel values, p_i, multiplied by some weights, w_i. So you can write d_j as the inner product between a vector containing the weights and a vector containing the pixel values. Your matrix would be formed by stacking the corresponding weighting vectors (as row vectors). But, obviously, this isn't the way to go. – lp251 Jan 20 '14 at 18:36
  • @NickXTsui As for references, look at this paper: [Branchless distance-driven projection](http://proceedings.spiedigitallibrary.org/proceeding.aspx?articleid=728852), in particular, at Figure 2b. The focus of that paper is on a modified version of the overlap kernel, but the first part describes the original version in good detail. – lp251 Jan 20 '14 at 18:37
2

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:

and for 3D DBT:

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.

rvimieiro
  • 1,043
  • 1
  • 10
  • 17
1

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.

Matt J
  • 1,127
  • 7
  • 15
  • Thanks a lot for your reply. So can I ask how distance driven is implemented? Can it be prototyped in matlab? I am wondering how I can fill in each element in the system matrix for distance driven. Really stuck here. Thanks. – Nick X Tsui Jan 14 '14 at 16:29
  • 2
    Here is [a paper by the developers of distance-driven projection](http://iopscience.iop.org/0031-9155/49/11/024/). I haven't read it, but it, and the papers it references, should be what you need to code your own. Once you have both forward and back proj in function form, I don't think you actually need to build the matrix, although you could do so with FUNC2MAT, referenced in my answer. – Matt J Jan 14 '14 at 16:39
  • Thanks a lot. One more question, when you say "Distance-driven projection is not implemented in stock MATLAB", you mean the case in industry? Or you mean the distance -driven cannot be coded into matlab due to its computational complexity, even for a simple case (ex. small image)? – Nick X Tsui Jan 14 '14 at 16:42
  • I mean The Mathworks has not implemented and packaged it for you in standard MATLAB or any of its toolboxes. You would have to write the MATLAB code for it yourself. – Matt J Jan 14 '14 at 16:44
  • I would also add that if, as you say, you are working with small images and simple cases, there really is not much reason to be pursuing distance-driven projection. The point of the distance-driven algorithm is to achieve speed when processing real-life data sizes. That's not something you should be caring about for small examples. Your main problem is really just finding a representation of the projection function, e.g., radon() or fanbeam(), that is transposable, and that you can readily do with FUNC2MAT. – Matt J Jan 15 '14 at 16:12
  • Hi Matt: Thanks for the followup, but I'd really like to stick with distance driven since it is the state of art for now. To start, small image is easy to analyze. Once the algorithm is clear to me, I can move forward with large scale. So I really don't see any conflict in this logic. – Nick X Tsui Jan 15 '14 at 16:48