Pozitron-Emission Tomography Reconstruction ( A computer graphics view )

1. Pozitron-Emission TomographyReconstruction(A computer graphics view) Szirmay-Kalos László

2. Scalar fields 3D tér pontjaiban egy skalár érték v(x,y,z) v(x,y,z) z y x • Skalár: hőmérséklet, sűrűség • nyomás, potenciál, … • Származás: Euler-i szimuláció, • Rekonstrukció (tomográfia) tárolás: 3D textúra vagy „voxel tömb”

3. Visualization of scalar fields • Megjelenítés fényszóró anyag (participating media) analógiáját felhasználva (belsejébe belelátunk) • Adott szintfelület kiemelése (külsőt lehámozzuk) Kép szintézis Transzfer függvény grad v v Kép Optikai paraméterek Sűrűség + deriváltak

4. Fényszóró közeg ds Hatáskeresztmetszet, alias kioltási tényező ·ds = P(ütközés) • Albedo a: az elnyelődés valószínűsége feltéve, hogy az ütközés bekövetkezett • Fekete test: albedo= 0 A=1

5. ds Sugársűrűség változása s L(s+ds) L(s) // Kiszóródás+abszorbció // Emisszió // Beszóródás L(s+ds)= L(s) – L(s)·(s)·ds + Le(s)·ds + (s)·a(s)·ds·f(‘,)Li(‘)d‘ dL(s)/ds= –L(s)·(s)+Le(s)+Linscatter(s) Megoldás fényelnyelő közegre (emisszió és beszóródás nincs): L(s)= L(0)·exp(–s (s)ds)

6. Tomography Mediso NanoPETTM/CT Absorption Emission P N N P e+ e- N P N L(s)= L(0)·exp(–(s)ds) (s)ds = – log L(s)/L(0) L(s)=Le(s)ds

7. Reconstruction

8. FBP = Filtered backprojection Impulse response: w(x,y)=w(r)  1/r i(x,y)=(x,y) Measurement + back projection Circle of radius R: rd • w(x,y)dxdy=2Rw(r)rdrd= • =2Rw(r)rdr R dr r Correction in frequency domain o(x,y)=i(x,y)w(x,y) FxFy o[x,y] = FxFy i[x,y]FxFyw[x,y] Ramp filtering 1/|| FxFy i[x,y] = FxFy o[x,y] ||

9. Algebric back projection Lin egyenlet (V<D): d= Av, d = [d1,..., dD] v = [v1,..., vV] Moore féle pseudo-inverz: AT d= ATA v v = (ATA)-1AT d =A+d (D×V) (D) (V) v1 v1 d4 v3 v4 d3 d2 (V×D) (D) (V×V) d1 = A11v1 +A13v3

10. Iteratív séma Skalármező Projekció: vonalintegrálok Összevetés a mért értékekkel Skalármező korrekciója

11. Light photons • Relativistic mass is small: E = mc2 = hf • Photon’s energy (wavelength) does not change upon elastic collision • Absorption probability is energy dependent • Wavelengts can be handled independently e- e-

12. Gamma photons e- e- • Relativistic mass is significant • Photon’s energy (wavelength) changes upon collision • Absorption probability is energy dependent • Wavelengts cannot be handled independently

13. Compton scattering and the Klein-Nishina phase function x E’ Compton formula: scattered photon incident photon  1 φ  z collision E Klein-Nishina differential cross section (extinction coefficient): y Cross section E/mec2 C(1+cos2) Electron density P(E, ) = … E/mec2 u Rejection sampling Table driven sampling

14. Gamma photon transport • Compton scattering + • Klein-Nishina formula • = New • Direction • Energy (freq) Detector grid

15. Why is it important? Density distribution from CT Source distribution ??? Measured detector response

16. Iterative reconstruction Source estimation Density distribution from CT Compute detector response Estimated detector response Source distribution ??? Measured detector response Source correction

17. PET Intensity: x e- e+ LOR: y

18. Finite-element approximation x(v) NxV bV (v) • Higher order filter • more accurate • more sensitive to noise • Non-CC grid? 1 box 1 tent b1 b1 b2 b2 b3 b3 constant Tri-linear

19. BCC and FCC grids CC BCC FCC

20. System matrix voxel • Problems: • Huge matrix: 108×107 • A matrix element is an infinite-dimensional integral • Matrix element depends on the measured object (CT) • Compute matrix elements on-the-fly (like in radiosity) 107 voxel intensities: 108 LOR responses: Probabilities that a photon pair emitted in voxel V is detected by LOR L

21. Iterative solver Measures LORs Computed LORs Current voxel intensities ~ ~ ~ ~ Correction Back-projection Forward-projection

22. Back-projectionMaximum likelihood estimation Measurement State: x result: y What is x based on y? Find x that maximizes P(y|x) Bayes estimation: P(x|y)=P(y|x)P(x)/P(y)

23. ~ ~ Back-projection • Measured LOR hits have Poisson-distribution with means of the current estimate: • Joint probability: • Likelihood: logarithm of the joint probability: • Estimate x maximizes the likelihood: +const

24. ~ ~ Parallel implementation issues Ratios: Same relative error in all detector pairs For each LOR L: For each voxel V: • Gathering not shooting! • Forward: voxels to LORs • Backward: LORs to voxels

25. Solution alternatives for gamma photon transport • Physicists’ approach: • Direct mathematical model of a phenomenon • Solution algorithm is the simulation of the nature • Porting to the computer • Performance tuning • Computer graphics approach: • Computer architecture • Algorithm that is efficient on this architecture • Tuning to the problem

26. Physicists’ approach Detector grid

27. Physicists’ approach Detector grid

28. Physicists’ approach Detector grid

29. Physicists’ approach Detector grid

30. Woodcock tracking: Free path length sampling in heterogeneous medium CG approach: x  e-maxx s S max s max

31. Physicists’ approach • Pros • Direct simulation of nature • Rejection sampling can mimic the Klein-Nishina phase function and free path length. • Paths are computed independently, so it scales well on MIMD (multi-CPU = slow and expensive) • Cons • Bad on SIMD (rejection sampling is a loop) • Similar absolute error in detectors, the relative error is 1/n where n is the number of hits. • Paths are computed independently, so it cannot exploit coherence (cannot reuse information) • Random writes

32. Computer graphics approach • SIMD-like algorithm (GPU) • Independent threads (gathering) • “No” conditional statements or variable length loops • Reuses paths • Detector view (gathering) • Same relative error in all detector pairs • Same number of samples in all detector pairs • Detector view (gathering) • Solve the adjoint problem  • Like path tracing instead of photon tracing • Integral transformation (Jacobi determinant)

33. CUDA: implementation platform GPU CPU + Kernel programs: Threads Thread blocks Warps, … Host program Thread block Fast shared memory SIMD Quasi-SIMD

34. Adding two N-element vectors Runs on GPU, called from CPU __global__ void AddVectorGPU( float *C, float *A, float *B, int N ) { int i = blockIdx.x*blockDim.x + threadIdx.x; // thread id if (i < N) C[i] = A[i] + B[i]; } float C[100000], A[100000], B[100000]; int main ( ) { … int N = 100000; … int blockDim = 256; // number of threads per block int gridDim = (N + blockDim – 1) / blockDim; // number of thread blocks AddVectorGPU<<<gridDim, blockDim>>>(C, A, B, N); … } 0 ,…, gridDim.x-1 0 ,…, blockDim.x-1

35. Path reuse! Detector 1 Detector 2

36. Integral transformation D2 dz2 dv: differential volume dw2: solid angle in which dz2 is visible from emission point dA: perpendicular area dl: chord length D1 dw: solid angle in which dz1and dz2are visible from dz1

37. Na voxel Nt Detector 1 Detector 2 Forward projection = =

38. Na ~ ~ voxel Nt Detector 1 Detector 2 Back projection ~ Solid angle of the farther detector: • It can change in time! • OSEM • sample perturbation

39. Detector sampling

40. Rnd Poisson-disk

41. Detector modeling

42. Measurements and phantoms • 18x324 LOR • 1283… 2563voxel

43. Time: NV GeForce GTX 285 GPU utilization: 25%

44. GPU optimization • Registers assigned to a thread (32->64) • Forward and back-projection re-structuring • Forward: • Back:

45. Error: 1283voxels

46. Reconstruction quality N_detline=16, Poisson diszk N_detline=1 ref

47. Single scatter compensation Scattering point = + 3 dimensions Reusing rays: Evaluate many integrals in parallel What to sample? What mimics all integrands? Sample the location of the scattering point

48. accumulated emission accumulated attenuation accumulated attenuation Single scatter compensation Importance sampling:

49. 1. Scattering points General scatter compensation 3. Combination of paths 2. Ray marching from detector to scattering points

50. Scatter compensation