Projection reconstruction and the Radon Transform 2

The Radon Transform

The forward Radon transform is to convert a 2-D object into a set of projects within the plane.

\(\(Radon[object_, n_, fov_] :=
         Table[Integrate[
             object\ \ DiracDelta[
                 m\  - \ x\ Cos[θ]\  -                    y\ Sin[θ]], {x, \(-fov\), fov}, {y, \(-fov\),
               fov}], 
							{m, \(-2\)\ fov/\((n -
                     1)\), \(+2\)\ fov/\((n - 1)\), 2\ fov/n}, {θ, 0,
             Pi, Pi/\((n - 1)\)}];\)\)

The above double integral is very slow to evaluate, and so we reduce it to a line integral along the line defined by the delta function.

Radon2[object_, x_, y_, n_, fov_] :=
       Table[NIntegrate[object[x, y], {yp, -fov, fov},
              {PrecisionGoal -> 4}], {xp, -fov, fov,
             2*fov/(n - 1)}, {θ, 0, Pi, Pi/(n - 1)}]  // N

Define a simple test object

object1[x_, y_] :=
       If[x^2  +  y^2  < 256, 1, 10^(-6)]  // N;

Plot3D[
       object1[x, y], {x, -64, 64}, {y, -64, 64}, {
         PlotRange -> All, PlotPoints -> {64, 64}}]

Robject1 =
       Radon2[object1, xp Cos[θ] - yp Sin[θ],
         yp Cos[θ]  +  xp Sin[θ], 64, 64];

ListPlot3D[Robject1]

Filtered Back Projection

Fdata = Fourier[Robject1];

ListPlot3D[Re[Fdata], {PlotRange -> All}]

\(\(\(Filt =
       Table[If[x <= 32\  && y <= 32, \ Sqrt[x^2 + y^2],
           
			
           If[x <= 32\  && y > 32, \ Sqrt[x^2\  + \ \((65 - y)\)^2],
             
				
             If[x > 32\  && \ y <= 32, \ Sqrt[\((65 - x)\)^2 + y^2],
               
					
               If[x > 32\  && \ y > 32,                  Sqrt[\((65 - x)\)^2 + \((65 - y)\)^2]]]]], {x, 0, 63}, {y, 0,
           63}]; \)
\)\)

ListPlot3D[Filt]

FiltFdata = Fdata*Filt;

ListPlot3D[Re[FiltFdata], {PlotRange -> All}]

Filtdata = Fourier[FiltFdata];

ListPlot3D[Re[Filtdata], {PlotRange -> All}]

Back Projection of Filtered

BFlimited[x_, y_, n_] :=
         If[x^2 + y^2 > 32^2, 0,
           1/(2 Pi) Sum[Transpose[
                     Re[Filtdata]][[m*64 + 1]][[Floor[
                       x Cos[m*Pi] + y Sin[m Pi]] + 33]], {m, 0, 1 - 1/n,
                 1/n}]];

Plot3D[
       BFlimited[x, y, 4], {x, -32, 32}, {y, -32,
         32}, {PlotRange -> All, PlotPoints -> {64, 
64}}]

DensityPlot[
       BFlimited[x, y, 4], {x, -32, 32}, {y, -32,
         32}, {PlotRange -> All, PlotPoints -> {64, 64},
         Mesh -> False}]

Image =
         Table[BFlimited[x, y, 64], {x, -32, 32}, {y, -32,
             32}];

ListPlot3D[Image, {PlotRange -> All}]

ListDensityPlot[
       Image, {PlotRange -> {0, 500}, Mesh -> False}]