Re: Looking for Hough and/or Radon transform code

• Subject: Re: Looking for Hough and/or Radon transform code
• From: "Bobstrosity" <dejastockwell(at)my-dejanews.com>
• Date: Mon, 17 May 1999 10:12:20 -0600
• Newsgroups: comp.lang.idl-pvwave
• Organization: The University of Western Ontario, London, Ont. Canada
• References: <7heo6b\$mkb\$1@coyote.dres.dnd.ca>
• Xref: news.doit.wisc.edu comp.lang.idl-pvwave:14773

```
John McFee wrote in message <7heo6b\$mkb\$1@coyote.dres.dnd.ca>...
>
>Would anyone have or know the whereabouts of IDL code to do either the
Hough
>transform or any variants of the Radon transform?
>
>Thanks
>
>John McFee

Included below is a function that calculates the radon
transform using linearly interpolation.
Copy it to a file called radon.pro

Below the function is a small piece of test code.
Run radon.pro for an example of how it can be used.
(it will produce three plots, and 3 lines of text output.
This is a preliminary piece of code so, as the
wise one (mike brady) says ``caveat emptor''.

Cheers,
bob stockwell

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
; Copyright (c) 1999-5000, R.G. Stockwell
;       Unauthorized modification prohibited.
;
; DISCLAIMER:  this is from RGS's collection of hackware. It
;       is not documented or optimized for speed. Use at your own risk.
;
;
; NAME:
;
; PURPOSE:
;   this incorporates the linear interpolation scheme
;   takes an image h(x,y) and returns  the radon transform
;   r(p, tau), where  p is the slope from -1 to 1
;   and tau is the y offset along a y axis at x= slength/2
;   ie that is halfway horizontally across the image.
;   if you want to find slopes > 1 or < -1 then rotate the
;   image 90 degrees.
;
; CATEGORY:
;       image pro.
;
; CALLING SEQUENCE:
;
; INPUTS:
;       image:    An 2D image
;
; OUTPUTS:
;       p:     1 D vector of the sampled slopes
;       tau:     1 D vector of the sampled y offsets
;
;
; KEYWORD PARAMETERS:
;      RETURNS EXTENDED INFORMATION, does not affect program execution.
;     all_rays:  a 3 dimensional array (p,tau,ray) holding the ray for each
slope and y offset
;   raylen:    the length of the ray for each ray in the above 3D array.
;      NOTE these parameters are not needed in normal use of the radon
transform.
;     They return the individual rays used in the radon calcuation in
;     you (or I) want to perform a different operation than averaging.
;
; EXAMPLE:  follows the function
;
; REFERENCE:
;       The Analysis of Time Series (Fourth Edition)
;       C. Chatfield
;       ISBN 0-412-31820-2
;
; MODIFICATION HISTORY:
;       Written by:  R.G. Stockwell
;     Research Scientist
;    3380 Mitchell Lane
;    stockwell@co-ra.com
;    early 1999
;-

if  n_params() lt 1 then return,-1
sz = size(image)
if sz(0) ne 2 then return,-1
xl = sz(1)
yl = sz(2)
x = findgen(xl)-xl/2
y = findgen(yl)
; max size of inray is max of indexes

kl=xl
hl=yl
ml=xl
nl=yl
delta_x =(delta_y=1)
p = (findgen(kl)/(kl-1)*kl-(kl)/2)*delta_y/max(abs(x))
tau = findgen(hl)*delta_y
x_min = min(x)
y_min = min(y)
alpha = p * delta_x/delta_y

raylen   =  fltarr(kl,hl)
all_rays =  fltarr(kl,hl,xl > yl)

for k = 0,kl-1 do begin
for h = 0,hl-1 do begin
beta = (p(k)*x_min+tau(h) - y_min)/delta_y
if alpha(k) gt 0 then begin
m_min = 0 > ceil((-beta)/alpha(k))
m_max = ml - 1 < floor((nl-1-beta)/alpha(k))
endif else begin
m_min = 0 > ceil((nl-1-beta)/alpha(k))
m_max = ml - 1 < floor((-beta)/alpha(k))
endelse
sum=0
mv = findgen(m_max-m_min+1)+m_min
nfloatv = (alpha(k)*mv+beta)
wneg = where(nfloatv lt 0.0,wnegcount)
wre = where(abs(nfloatv) lt 0.01,wrecount)
if wrecount gt 0 then nfloatv(wre) = 0
nv = floor(nfloatv)
wv = nfloatv-nv
sum = 0
if max(nfloatv) gt nl-1 then stop ;this is an error in indexing
wno = where(wv eq 0,wnocount)
inray = fltarr(n_elements(mv))
if wnocount gt 0 then begin
inray(wno) = image(mv(wno),nv(wno))
wint = where(wv ne 0,wintcount)
if wintcount gt 0 then begin

inray(wint)=image(mv(wint),nv(wint))*(1-wv)+image(mv(wint),nv(wint)+1)*wv
endif
endif else inray = image(mv,nv)*(1-wv)+image(mv,nv)*wv
raylen(k,h) =n_elements(inray)
all_rays(k,h,0:raylen(k,h)-1) = inray
endfor
endfor

end

;__________________________________________________
;__________________________________________________
; test code here

xl=(yl=100)
image = fltarr(xl,yl)

;slope=-0.8
;offset = 0
slope= (randomu(seed,1))(0)*2-1
offset = (randomu(seed,1))(0)*xl

line = 10
for i = 0,yl-1 do begin
yind = slope*(i-xl/2)+offset
if yind ge 0 and yind lt yl then \$
image(i,yind) = image(i,yind) + line
endfor

image=smooth(image,5)
image = image+randomn(seed,xl,yl)

print
print,'Actual Values for slope and y offset:......... P=',slope,  '    Tau =
: ',offset
print,'Measured slope and offset from radon transform: P=',p(wm(0)),'    Tau
= : ',tau(wm(1))

!P.multi=[0,1,3]
!P.charsize=2
fill=1
contour,image,fill=fill,tit='Image (line + noise)',xtit='X',ytit='Y'
axis,xl/2,/data,yaxis=1