| Re: backprojection [message #52177 is a reply to message #52176] |
Thu, 18 January 2007 00:40   |
TimLS
Messages: 5
Registered: August 2006
|
Junior Member |
|
|
I got the radon transform to work for me. The link to Mark River's
web-site was very interesting and useful.
Thanks, Tim
Wox wrote:
> 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
