**Modern tomography**

More modern variations of tomography involve gathering projection data from multiple directions and feeding the data into a tomographic reconstruction software algorithm processed by a computer.Different types of signal acquisition can be used in similar calculation algorithms in order to create a tomographic image. With current (2005) technology, tomograms are derived using several different physical phenomena listed in the following table.

Physical phenomenon | Type of tomogram |
---|---|

X-rays | CT |

gamma rays | SPECT |

radio-frequency waves | MRI |

electron-positron annihilation | PET |

electrons | Electron tomography or 3D TEM |

ions | atom probe |

**1. Acoustic tomography**

Assume the acoustical property of a moving medium is described by the speed of sound, denoted by , and the velocity of the medium is . Here we assume all functions are well-behaved.

Then the problem turns out to be an **eikonal equation**.

,

This equation is called **Jacobi-Hamilton equation**. And it can be reduced to characteristic system of differential equations. The system of DE can be solved by **shooting method.**

**Remark**:

For **eikonal equation:**

, where .

There are some numerical approaches for it.

- Fast Marching Method
- Fast Sweeping Method

For the first one, it uses **Dijkstra’s algorithm**, due to the fact that information only flows outward from the seeding area.

For the second one, it is attributed to ZHAO HONGKAI, from UCI. The algorithm is simple and only uses** FDM** and **upwind** scheme. However, this method is not operational for all kinds of boundary value problems. Here is a example for it.

function fastsweep(N) %It is for solving the sample eikonal equation |Du|=1 in a square [0,1]X[0,1]. h=1/N; for i=1:N+1 for j=1:N+1 u(i,j)=0; end end for i=2:N for j=2:N u(i,j)=1000^i+100^j; end end for k=1:N for i=2:N for j=2:N ux=min(u(i-1,j),u(i+1,j)); uy=min(u(i,j-1),u(i,j+1)); if abs(h^2-(ux-uy)^2)<10^(-2) u_tmp(i,j)=min(ux,uy)+h; else u_tmp(i,j)=0.5*(ux+uy+sqrt(2*h^2-(ux-uy)^2)); end u(i,j)=min(u(i,j),u_tmp(i,j)); end end for i=N:2 for j=2:N ux=min(u(i-1,j),u(i+1,j)); uy=min(u(i,j-1),u(i,j+1)); if abs(h^2-(ux-uy)^2)<10^(-2) u_tmp(i,j)=min(ux,uy)+h; else u_tmp(i,j)=0.5*(ux+uy+sqrt(2*h^2-(ux-uy)^2)); end u(i,j)=min(u(i,j),u_tmp(i,j)); end end for i=N:2 for j=2:N ux=min(u(i-1,j),u(i+1,j)); uy=min(u(i,j-1),u(i,j+1)); if abs(h^2-(ux-uy)^2)<10^(-2) u_tmp(i,j)=min(ux,uy)+h; else u_tmp(i,j)=0.5*(ux+uy+sqrt(2*h^2-(ux-uy)^2)); end u(i,j)=min(u(i,j),u_tmp(i,j)); end end for j=N:2 for i=2:N ux=min(u(i-1,j),u(i+1,j)); uy=min(u(i,j-1),u(i,j+1)); if abs(h^2-(ux-uy)^2)<10^(-2) u_tmp(i,j)=min(ux,uy)+h; else u_tmp(i,j)=0.5*(ux+uy+sqrt(2*h^2-(ux-uy)^2)); end u(i,j)=min(u(i,j),u_tmp(i,j)); end end end surf(u,'FaceColor','yellow','EdgeColor','none') camlight left; lighting phong end Result:

>> fastsweep(100)

Here is a picture for the computing result.

However, there is a singularity in this example, thus there will be some propagation error around the singularity. We decided to sweep the domain many times to reduce the error. In ZHAO’s paper, the cost is at , here I used a scheme to achieve the convergence.

Improved code and algorithm for this problem:

function fastsweep3(N) tic h=1/N; for i=1:N+1 for j=1:N+1 u(i,j)=0; end end for i=2:N for j=2:N u(i,j)=1000^i+100^j; end end for k=1:floor(N/2) j=1+k; for i=1+k:N+1-k ux=min(u(i-1,j),u(i+1,j)); uy=min(u(i,j-1),u(i,j+1)); if abs(h^2-(ux-uy)^2)<10^(-2) u_tmp(i,j)=min(ux,uy)+h; else u_tmp(i,j)=0.5*(ux+uy+sqrt(2*h^2-(ux-uy)^2)); end u(i,j)=min(u(i,j),u_tmp(i,j)); end i=N+1-k; for j=1+k:N+1-k ux=min(u(i-1,j),u(i+1,j)); uy=min(u(i,j-1),u(i,j+1)); if abs(h^2-(ux-uy)^2)<10^(-2) u_tmp(i,j)=min(ux,uy)+h; else u_tmp(i,j)=0.5*(ux+uy+sqrt(2*h^2-(ux-uy)^2)); end u(i,j)=min(u(i,j),u_tmp(i,j)); end j=N+1-k; for i=N+1-k:-1:1+k ux=min(u(i-1,j),u(i+1,j)); uy=min(u(i,j-1),u(i,j+1)); if abs(h^2-(ux-uy)^2)<10^(-2) u_tmp(i,j)=min(ux,uy)+h; else u_tmp(i,j)=0.5*(ux+uy+sqrt(2*h^2-(ux-uy)^2)); end u(i,j)=min(u(i,j),u_tmp(i,j)); end i=1+k; for j=N+1-k:-1:1+k ux=min(u(i-1,j),u(i+1,j)); uy=min(u(i,j-1),u(i,j+1)); if abs(h^2-(ux-uy)^2)<10^(-2) u_tmp(i,j)=min(ux,uy)+h; else u_tmp(i,j)=0.5*(ux+uy+sqrt(2*h^2-(ux-uy)^2)); end u(i,j)=min(u(i,j),u_tmp(i,j)); end end %surf(u,'FaceColor','yellow','EdgeColor','none') %camlight left; lighting phong toc end

This time I used cost to produce he mesh. stands for the total grid points, boundary included. Though on the boundary there are only grid points.

**End of remark.**

2. **Radon Transformation**

For an arbitrary function , its Radon transform is defined as the integral of along a straight line :

,