Prev: GUIKit: How to set up a "tableHeader" in table
Next: ReplaceAll reloaded
From: cire g on 10 Feb 2010 03:37
Hello Mathgroup,
I am reproducing the sinogram presented in wikipedia:
http://en.wikipedia.org/wiki/Radon_transform

Notice how easy is to do this in mathematica, the code resembles the text.

But my code is terrible slow in computing the radon transform and the
image reconstruction... I will appreciate your insight in how to code it
faster. Where is the bottleneck?

Best regards,
eric

this is my code:

kostmoPhantomF[x_, y_] := UnitBox[(x - 16)/32, (y - 16)/32] + UnitBox[(x
+ 16)/32, (y + 16)/32];
kostmoPhantom = Table[kostmoPhantomF[i, j], {i, -64, 64, 1}, {j, -64,
64, 1}];
Image[kostmoPhantom] (* this is the phantom in the wikipidea page named
after his creator*)

(* radon transform*)
radonTransform[a_, s_, f_] := NIntegrate[f[t Sin[a] + s Cos[a], -t
Cos[a] + s Sin[a]], {t, -Infinity, Infinity}]

Timing[sinogram32p64 = Table[radonTransform[a, r, kostmoPhantomF], {a,
0, Pi, Pi/64}, {r, -50, 50, 100/32}];]
(* results in my system {48.6976, Null} *)

Timing[sinogram32p64 = Table[radonTransform[a, r, kostmoPhantomF], {a,
0., 1.0 Pi, Pi/64.0}, {r, -50., 50., 100/32.0}];]
(*results in my system {23.6514, Null}*)

(* sinogram image*)
Image[sinogram32p64*.01]

(*-------simple backprojection reconstruction----------*)
n = Dimensions[sinogram32p32][[1]];
a = Table[aa, {aa, 0.0, 1.0 Pi, Pi/32.0}]; (* 32 projections *)

ImageReconF[x_, y_] := (1/n)*(Sum[radonTransform[a[[i]], x Cos[a[[i]]] +
y Sin[a[[i]]],kostmoPhantomF], {i, 1, n}])
Timing[ImageRecon = Table[ImageReconF[i, j], {i, -4, 4}, {j, -4, 4}];]
(* results in my system {34.8057, Null}, for only 9x9 pixels! *}
Image[ImageRecon*.005]




 | 
Pages: 1
Prev: GUIKit: How to set up a "tableHeader" in table
Next: ReplaceAll reloaded