Joel Gales - 2006-09-12

Logged In: YES
user_id=958613

I'm not sure if this helps but there is a RADON function
implemented in GDL.

Joel

------------------------------------------------------------

PRO test_radon

DEVICE, DECOMPOSED=0

;Create an image with a ring plus random noise:
x = (LINDGEN(128,128) MOD 128) - 63.5
y = (LINDGEN(128,128)/128) - 63.5
radius = SQRT(x^2 + y^2)
array = (radius GT 40) AND (radius LT 50)
array = array + RANDOMU(seed,128,128)

;Create display window, set graphics properties:
WINDOW, XSIZE=440,YSIZE=700, TITLE='Radon Example'
!P.BACKGROUND = 255 ; white
!P.COLOR = 0 ; black
!P.FONT=2
ERASE

XYOUTS, .05, .94, 'Ring and Random Pixels', /NORMAL
;Display the image. 255b changes black values to white:
TVSCL, 255b - array, .05, .75, /NORMAL

;Calculate and display the Radon transform:
XYOUTS, .05, .70, 'Radon Transform', /NORMAL
result = RADON(array, RHO=rho, THETA=theta)
TVSCL, 255b - result, .08, .32, /NORMAL
PLOT, theta, rho, /NODATA, /NOERASE, $
POSITION=[0.08,0.32, 1, 0.68], $
XSTYLE=9,YSTYLE=9,XTITLE='Theta', YTITLE='R'

;For simplicity in this example, remove everything except
;the two stripes. A better (and more complicated) method
would
;be to choose a threshold for the result at each value of
theta,
;perhaps based on the average of the result over the theta
;dimension.
result[*,0:55] = 0
result[*,73:181] = 0
result[*,199:*] = 0

;Find the Radon backprojection and display the output:
XYOUTS, .05, .26, 'Radon Backprojection', /NORMAL
backproject = RADON(result, /BACKPROJECT, RHO=rho,
THETA=theta)
TVSCL, 255b - backproject, .05, .07, /NORMAL