Tomography Note 1

Modern tomography

refer to wikipedia
 

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 c(r), and the velocity of the medium is v(r). Here we assume all functions are well-behaved.

Then the problem turns out to be an eikonal equation.

\displaystyle\frac{c_0}{c(r)}-\frac{v(r)\cdot \nabla \zeta(r)}{c(r)}=|\nabla \zeta(r)|,

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:

|\nabla u(x)|=f(x), where f(x)\ge 0.

There are some numerical approaches for it.

  1. Fast Marching Method
  2. 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 O(n), here I used a O(n^{1.5}) 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 O(N) cost to produce he mesh. N stands for the total grid points, boundary included. Though on the boundary there are only O(\sqrt{N}) grid points.

End of remark.

2. Radon Transformation

For an arbitrary function m(x, y), its Radon transform is defined as the integral of m(x, y) along a straight line \Gamma(r,\phi):

\displaystyle\mathcal{R}m(r, \phi) =\int_{\Gamma(r,\phi)}m(x, y) \mathrm{d}s,

About YiMin

This is just a nerd PhD student of Math@UT Austin.