FPGA Accelerated 3-D Tomography

1. FPGA Accelerated3-D Tomography Richard Dorrance Progress Update: 09/07/12

2. Outline • Introduction to Tomography • Reconstruction Methods • Analytical • Backprojection • Filtered Backprojection • Algebraic • Algebraic Reconstruction Technique (ART) • Simultaneous Iterative Reconstruction Technique (SIRT) • Simultaneous Algebraic Reconstruction Technique (SART) • Modeling Performance of Reconstruction Methods • Future Work

3. Tomography • Cross-sectional imaging technique using transmissionor reflection data from multiple angles • Basis for CAT scan, MRI,PET, SPECT, ET, etc. • Computed Tomography (CT):A form of tomographic reconstruction on computers

4. Cross-Sections by X-Ray Projections • Project X-ray through biological tissue;measure total absorption of ray by tissue • Projection Pθ(t) is the Radontransform of object functionf(x,y): • Total set of projections calledsinogram

5. Shepp-Logan Phantom • Standard test image for tomographic reconstructions

6. Example Image with Projections

7. CT Reconstruction • Restore image from projection data • Inverse Radon transform • Most common algorithm is filtered backprojection • “Smear” each projection over image plane • Accuracy of reconstruction depends on the number of detectors and projection angles Original 4 Angles 16 Angles 64 Angles 256 Angles

8. Analytical Reconstruction Methods • (Filtered) Backprojection Pseudo Code: • Input: sinogram sino(θ, N) • Output: image img(x,y) for each θ filter sino(θ,:) ; only for FBP for each x for each y n = x*cos(θ)+ y*sin(θ) img(x,y)= sino(θ,n)+ img(x,y)

9. Backprojection (Step 1)

10. Backprojection (Step 2)

11. Backprojection (Step 3)

12. Backprojection (Step 4)

13. Backprojection (Step 5)

14. Backprojection vs. Original • Final Step: normalize image power • Divide each pixel by θ·N

15. Note On Filtering No Filtering With Filtering

16. Filtered Backprojection (Step 1)

17. Filtered Backprojection (Step 2)

18. Filtered Backprojection (Step 3)

19. Filtered Backprojection (Step 4)

20. Filtered Backprojection (Step 5)

21. Filtered Backprojection vs. Original

22. Conventional Algebraic Reconstruction Methods

23. Problem Formulation • We want to formulate it as a Linear Inverse Problem: • Where x is a column vector of length N2 representing the pixels of the original image, A is an M by N2 matrix representing the data acquisition process, and b is a column vector of length M representing the measured projection data. • We want to find a solution such that:

24. Notes on the Discretized Image x • The discretized image is denoted by:and by:where x is obtained by stacking the columns of X.

25. Notes on the projection data b • There are a total of d detectors and θ projection angles, so that a total of M = d · θare used. • Then the measured projection data is denoted by:and by:where b is obtained by stacking the columns of B.

26. Notes on the Acquisition Matrix A • The acquisition of projection data b from x is modeled by:where: • ai,j is the contribution of pixel j to projection i. • Also, let:be a column matrix that represents the ith ray which computes the value of the ith projection.

27. Iterative Reconstruction Algorithm • Let x(k) denote the kth estimation of the reconstruction. • Then:where the relaxation factor λ is a scalar.

28. Proof of Convergence [1] • Let • Then

29. Proof of Convergence [2] • If ATA is positive definite and λ is chosen so that the spectral radius of Δ is less than 1, then:and

30. Proof of Convergence [3] • Therefore:

31. # of Projections needed for ART • Reconstruction on a square grid (N×N) with N detectors • Assuming a circular reconstruction region, we can ignore pixels outside this region

32. # of Projections needed for FBP [1] • Reconstructing region with diameter L • Sampling interval is at least:with a maximum frequency of: • Due to polar sampling,the density of samplesdecreases as we gooutward on the polar grid

33. # of Projections needed for FBP [2] • To ensure a sampling rate of at least Δω everywhere:therefore:

34. Matrix Formulation with Normalization • Introduce diagonal matrices V and W: • V: diagonal matrix of theinverse of the row sums • W: diagonal matrix of theinverse of the column sums

35. Reconstruction Methods • Algebraic Reconstruction Technique • Update image after each ray is processed • Simultaneous Iterative Reconstruction Technique • Update image after all rays are processed • Simultaneous Algebraic Reconstruction Technique • Update image after all rays in a single projection angle are processed

36. ART • Image update method: • After each ray is processed • Pseudocode: for k = 1:K for i = 1:M end end

37. ART (Iterations 1-6)

38. SIRT • Image update method: • After all rays are processed • Pseudocode: for k = 1:K end

39. SIRT (Iterations 1-6, λ = 0.5)

40. SART • Image update method: • After all rays in a single projection angle are processed • Pseudocode: for k = 1:K forθ = 1:Θ end end

41. SART (Step 1, Iteration 1, Theta 1)

42. SART (Step 2, Iteration 1, Theta 1)

43. SART (Step 1, Iteration 1, Theta 2)

44. SART (Step 2, Iteration 1, Theta 2)

45. SART (Step 1, Iteration 1, Theta 3)

46. SART (Step 2, Iteration 1, Theta 3)

47. SART (Step 1, Iteration 1, Theta 4)

48. ART (Step 2, Iteration 1, Theta 4)

49. SART (Iterations 1-6)

50. Modeling Performance (CPU, GPU, FPGA) • Write C pseudo code for Matrix-Vector multiplication and Vector-Vector addition • Convert C pseudo code to application specific pseudo code (CPU = x86, GPU = OpenCL/CUDA) • Model latency and throughput of pseudo code given: • CPU architecture: • Cache structure, freq., total # of threads, etc… • Image reconstruction problem: • N, d, θ, A matrix sparsity (α), # of iterations, etc…