| Re: backprojection [message #52180 is a reply to message #52177] |
Wed, 17 January 2007 08:32   |
Wox
Messages: 184
Registered: August 2006
|
Senior Member |
|
|
On 17 Jan 2007 07:00:22 -0800, "TimLS" <tlsmith@qinetiq.com> wrote:
> Is there anyone out there in the medical imaging world who has
> some code or who can give me a clue as to what I need to do?
Try the code below. The function 'createsinogram' makes a sinogram.
testtomo uses the radon function for reconstruction as Mr. Fanning
suggested.
Some things you have to supply to radon:
1. projcen: the center of projection (i.e. the projection of the
center of rotation)
2. angleinc: rotation increment
The first might be tricky. For info on how to extract the center of
projection from a sinogram:
http://www-fp.mcs.anl.gov/xray-cmt/rivers/tutorial.html
They have also some IDL code there, including filters on the sinogram
to improve the reconstructed image.
Remark: Notice that in the example, the center of projection is not in
the middle. The reconstruction is not very good because the object
falls out of the detection area during rotation. Using filter (as
metioned above) may improve images like that.
Anyway, hope this helps.
function createsinogram
; 1. Create an image with rings plus random noise:
N=250
M=250
R0=(N<M)*0.1
R1=(N<M)*0.12
xc=N/2.+20
yc=M/2.
x = (LINDGEN(N,M) MOD N)
y = (LINDGEN(N,M)/N)
radius = SQRT((x-50)^2 + (y-80)^2)
realslice = (radius GT R0) AND (radius LT R1)
radius = SQRT((x-40)^2 + (y-100)^2)
realslice OR= (radius GT R0) AND (radius LT R1)
radius = SQRT((x-60)^2 + (y-100)^2)
R0=0
realslice OR= (radius GT R0) AND (radius LT R1)
realslice += RANDOMU(seed,N,M)
realslice[xc-1:xc+1,yc-1:yc+1] = 2*max(realslice)
tvscl,realslice,0, 0, /NORMAL
; 2. Create sinogram
sinogram = RADON(realslice,xmin=-xc,ymin=-yc,drho=1,NRHO=N,RMIN=-xc)
tvscl,sinogram,0, 0.5, /NORMAL
return,sinogram
end;function createsinogram
pro testtomo
sinogram=createsinogram()
s=size(sinogram)
N=s[2]
NAngles=s[1]
projcen=N/2.+20
angleinc=180./NAngles
anglestart=0.
angleinc*=!pi/180
anglestart*=!pi/180
angles=anglestart+angleinc*findgen(NAngles)
rho=findgen(N)-projcen
tomogram=radon(sinogram,/BACKPROJECT,rho=rho,theta=angles,/L INEAR,nx=N,ny=N,xmin=-projcen)
tvscl,tomogram,0.5, 0, /NORMAL
end
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
