FUNCTION LinearHoughLink,xpts,ypts,THRESHOLD=thresh,NLINES=nlines,$
	RHOBINS=N,THETABINS=M,HOUGH_MATRIX=A1,FINAL_MATRIX=A,H_MATRIX=H,$
	NHD_RAD=nrem,VERBOSE=verbose,MIN_LENGTH=minLen,MAX_LENGTH=maxLen,$
	MAXSTEP=maxStep,RHOMAX=rhomax,RHOMIN=rhomin
;+
;RESULT=LinearHoughLink(xpts,ypts) finds the endpoints of linse in
;sets of points whose coordinates are given by xpts,ypts. The funcion
;uses the Hough transform and is controlled by a number of keyword
;parameters.
;
;RESULT is an array of anonymous structures with the following tags
;and initial values:
;
;Index:-1L			Each line has a unique index [1,2,...,nlines]
;Xends:[0.0,0.0]	X coordinates of beginning and ending points
;Yends:[0.0,0.0]	Y coordinates of beginning and ending points
;THETA:0.0			The THETA of the normal representation
;RHO:0.0			The RHO of the normal representation
;Length:0.0			The length of the line
;Location:0L}		The location in the Hough Parameter matrix
;
;EXAMPLE
;Suppose that IMG is a binary image. An image B of boundary
;points that describe object boundaries can be found by erosion
;with a structuring element S. The line descriptions of the
;boundaries can be found and plotted as an overlay on the
;image.with the following code.
;
;S=REPLICATE(1,3,3)
;B=image AND NOT ERODE(IMG,S)
;imSiz=SIZE(B)
;WINDOW,/FREE,XSIZE=imSiz[1],YSIZE=imSiz[2],TITLE='Edges'
;lines=LinearHoughLink(x,y)
;nlines=MAX(lines.Index) ;Find the number of lines found.
;FOR k=0,nlines-1 DO PLOTS,lines[k].Xends,lines[k].Yends
;
;DESCRIPTION AND KEYWORDS
;The algorithm utilizes the normal form of a straight line
;given by rho=x*Cos(theta)+y*Sin(theta). (theta,rho) space
;is divided into a grid of size NxM over the interval -90<theta<90
;and rhomin<rho<rhomax. The range of theta is fixed and the
;range of rho is found internally by a pass through the data.
;
;KEYWORDS: HOUGH_MATRIX, RHOBINS, THETABINS, RHOMAX
;The (theta,rho) cells correspond to a matrix A, which
;accumulates the count of the number of parametric curves in
;(theta,rho) space that pass through each cell. The number of
;cells can be set through the keywords RHOBINS and THETABINS.
;The default values are RHOBINS=100, THETABINS=100. The
;accumulated A matrix can be accessed via the keyword
;HOUGH_MATRIX. A[i,j] gives the number of lines that have
;theta and rho values that correspond to cell A[i,j].
;The range of rho may be given by the keyword RHOMAX.The
;range is (-RHOMAX,RHOMAX). The default is to find the maximum value
;in the data.
;
;KEYWORD: H_MATRIX
;A second array H is an array of pointers. Each pointer contains
;an array with the set of original (X,Y) points that contributed
;to the parametric curves that are counted in the corresponding cell of A.
;The original points that contributed to the parameter count A[i,j] is
;given by *H[i,j]. H is accessed with the keyword H_MATRIX. The
;list of points for cell [i,j] can be found with the following
;statements, which extract the interlaced coordinate values into
;lists xv and yv.
;	kp=i+j*THETABINS
;	kx=INDGEN(A[kp])*2+1
;	ky=kx+1
;	xv=(*H[kp])[kx]
;	yv=(*H[kp])[ky]
;
;KEYWORDS: THRESHOLD, NLINES
;The search for lines progresses by cycling through the A matrix
;and using a current maximum value to identify a parameter cell.
;The search continues until it finds NLINES or there are no
;remaining cell counts that exceed THRESHOLD. The default value
;of NLINES=1 and the default THRESHOLD=0.9*MAX(A). The actual
;number of lines returned is N_ELEMENTS(Result), which may be less
;than NLINES.
;
;KEYWORDS: NHD_RAD, FINAL_MATRIX, MIN_LENGTH, MAX_LENGTH, MAXSTEP
;The search algorithm provides some parameters that facilitate
;weeding out unwanted lines. Each time a cell of A is used then
;that element is set to a value of -2 and the cells in a
;neighborhood of size NHD_RAD in each direction are set to -1. This
;marks the central cell and the neighborhood and prevents any A-cells
;in the neighborhood from being used. The default value is NHD_RAD=1.
;
;The search also rejects lines that are shorter than MIN_LENGTH or
;longer than MAX_LENGTH. The default values are MIN_LENGTH=0 and
;MAX_LENGTH=1000.
;
;Some lines have breaks. When this happens it is possible to keep
;just the longest unbroken section. That section will be kept if it
;fits the minimum and maximum length conditions, The size of the
;gap between line sections that is acceptable without breaking the
;line is given by the keyword MAXSTEP. The default value is large,
;MAXSTEP=1000.
;
;The FINAL_MATRIX keyword can be used to access the last parameter
;count matrix used by the algorithm just before it exited. This is
;useful in examining the parameter combinations that have
;contributed to the lines that were returned. A visual presentation
;can be constructed as in the code below.
;
;Assume that "FINAL_MATRIX=A1" is part of the function call.
;
;WINDOW,/FREE,TITLE='FINAL MATRIX'
;thetaval=(INDGEN(7)-3)*30
;xtks=FINDGEN(THETABINS)/(THETABINS-1)*180-90
;ytks=INDGEN(RHOBINS)
;CONTOUR,A1,xtks,ytks,LEVELS=INDGEN(15)-2,XRANGE=[-90,90],YRANGE=[0,rhobins],$
;	XTICKS=6,XTICKV=thetaval,XTICKLEN=1,YTICKLEN=1
;
;KEYWORD VERBOSE
;Setting VERBOSE=1 produces a printout that summarizes each line as the
;search progresses. VERBOSE=2 adds the list of coordinate points to the
;printout for each line. The default is VERBOSE=0.
;
;HISTORY
;H. Rhody, September, 1999 for Digital Image Processing class.
;-

;Set up defaults
IF KEYWORD_SET(nlines) EQ 0 THEN nlines=1
IF N_ELEMENTS(N) EQ 0 THEN N=100
IF N_ELEMENTS(M) EQ 0 THEN M=100
IF KEYWORD_SET(minLen) EQ 0 THEN minLen=0
IF KEYWORD_SET(maxLen) EQ 0 THEN maxLen=1000
IF KEYWORD_SET(maxStep) EQ 0 THEN maxStep=1000
IF KEYWORD_SET(verbose) EQ 0 THEN verbose=0

;Set up the neighborhood that is eliminated around each
;accepted point in the parameter count matrix during the
;search process.
IF KEYWORD_SET(nrem) EQ 0 THEN nrem=1
meshgrid,2*nrem+1,2*nrem+1,xNbrVec,yNbrVec
xNbrVec=xNbrVec-nrem
yNbrVec=yNbrVec-nrem

;Set up the parameter count and point array matrices.
A=INTARR(N,M)
H=PTRARR(N,M,/ALLOCATE_HEAP)
FOR I=0,n-1 DO FOR J=0,M-1 DO *H[I,J]=[0]

;Define an array of line structures to be filled in by the search.
;The index fields are initialized to -1 to facilitate finding
;the number of lines actually located should the search
;fall short.
lines=REPLICATE({Index:-1L,Xends:[0.0,0.0],Yends:[0.0,0.0],$
	THETA:0.0,RHO:0.0,Length:0.0,Location:0L},nlines)

; Construct the theta vector of length N and the associated
; index vector.
theta=FINDGEN(N)*!PI/(N-1)-!pi/2 ; -pi/2<=theta<=pi/2 in N steps
kt=INDGEN(N) ;Vector of theta indices.

; Construct vectors for the trig functions so we only have to
; do the trig stuff once.
costheta=COS(theta)
sintheta=SIN(theta)

; If the lengths of the coordinate vectors are not equal,
; then use the smaller number of points.
NPTS=MIN([N_ELEMENTS(xpts),N_ELEMENTS(ypts)])
xpts=xpts[0:NPTS-1] & ypts=ypts[0:NPTS-1]

; Find the limits of rho. This is needed to quantize the
; parameter space
IF KEYWORD_SET(rhomax) EQ 0 THEN $
rhomax= SQRT(MAX(xpts^2+ypts^2)) ELSE rhomax=FLOAT(rhomax)
rhomin=-rhomax
deltaRho=(rhomax-rhomin)/(M-1)

IF KEYWORD_SET(VERBOSE) THEN $
PRINT,FORMAT='("Rho limits: [",f6.1,",",f6.1,"] deltaRho=",g8.4)',$
	[rhomin,rhomax,deltaRho]

; Construct the count in parameter space, for each point
; p=(x,y) provided in the function call.
FOR p=0L,NPTS-1 DO BEGIN
	; Calculate the rho(theta) vector for each point
	rho=xpts[p]*costheta+ypts[p]*sintheta
	;Quantize the theta rho vector so that the indexes
	;have acceptable values in the range (0,M-1).
	kr=LONG((rho-rhomin)/deltaRho)
	;Calculate the array index vector
	k=kt+kr*N	; Get the array index of each (theta,rho) pair
	A[k]=A[k]+1	; INCREMENT the counting matrix at each point.
	; Put the data points in the correct coordinate arrays in
	; the H matrix
	FOR I=0L,N_ELEMENTS(K)-1 DO BEGIN
		*H[K[I]]=[*H[K[I]],[xpts[p],ypts[p]]]
	ENDFOR
ENDFOR

;Preserve the original counting matrix to return if it is requested
;with the HOUGH_MATRIX keyword.
A1=A


;==================================================
;Conduct the search procedure. The HISTOGRAM function
;with the REVERSE_INDICES keyword is used with each
;iteration to find the next peak. This seems to be
;more efficient than other methods.
;==================================================

linesFound=0
IF KEYWORD_SET(thresh) EQ 0 THEN thresh=0.9*MAX(A)
REPEAT BEGIN
	HA=MAX(A,kp) & kp=kp[0]; Get the index kp of the peak
	xp=kp MOD N	;The horizontal coordinate of the peak
	yp=kp/N		;The vertical coordinate of the peak
	tnbrs=xp+xNbrVec ;The coordinates of the neighbors
	rnbrs=yp+yNbrVec
	;Check that the neighbors are within the array boundaries.
	inbrs=WHERE(tnbrs GE 0 AND tnbrs LT N AND rnbrs GE 0 AND rnbrs LT M)
	nbrs=tnbrs[inbrs]+N*rnbrs[inbrs]; Indices of valid neighbors
	;The neighbors will be removed later if the cell is accepted.

	;Get the points that contributed to A. These are extracted from H.
	kx=INDGEN(A[kp])*2+1
	ky=kx+1
	xv=(*H[kp])[kx]
	yv=(*H[kp])[ky]
	;Calculate some initial spacing properties of the list of points.
	xmin=MIN(xv) & xmax=MAX(xv) & dx=xmax-xmin
	ymin=MIN(yv) & ymax=MAX(yv) & dy=ymax-ymin

	;Sort the points and find the spacing between
	;successive points in the sequence. Eliminate
	;those points that are too far from the others
	IF dx GE dy THEN BEGIN
		iv=SORT(xv)
		xv=xv[iv] & yv=yv[iv]
		xs=[xv[0],xv] & ys=[yv[0],yv]
		dv=sqrt((xv-xs)^2+(yv-ys)^2)
	ENDIF ELSE BEGIN
		iv=SORT(yv)
		xv=xv[iv] & yv=yv[iv]
		xs=[xv[0],xv] & ys=[yv[0],yv]
		dv=sqrt((xv-xs)^2+(yv-ys)^2)
	ENDELSE
	;This is where we look for gaps and select the longest
	;line section.
	jv=WHERE(dv GT maxStep)
	IF jv[0] GE 0 THEN BEGIN ;Point spacing too big somewhere
		jv1=[0,jv,A[kp]]
		jv2=[0,0,jv]
		jv3=jv1-jv2
		JMAX=MAX(jv3,jv4)
		jv4=jv4[0]
		jv5=INDGEN(jv1[jv4]-jv1[jv4-1])+jv1[jv4-1]
		;jv5 is the list of indices of points in the longest
		;section. Now select just those points from xv,yv.
		xv=xv[jv5] & yv=yv[jv5]
	ENDIF	;End check on point spacing

	;Calculate final spacing properties of the list of points.
	xmin=MIN(xv) & xmax=MAX(xv) & dx=xmax-xmin
	ymin=MIN(yv) & ymax=MAX(yv) & dy=ymax-ymin
	;Calculate the length of the line.
	linelen=SQRT(dx^2+dy^2)
	;Check to see if the line length is in the acceptable range.
	;Otherwise skip where the region of A is eliminated and then
	;repeat the search.

	;Remember the value at the peak then set that array value to\
	;zero.
	peak=A[kp]
	A[kp]=0

	;Check that the line length satisfies the conditions.
	IF linelen GE minLen AND linelen LE maxLen THEN BEGIN

	; Find the rho and theta for the Hough index kp.
	kt1=kp MOD N
	kr=kp/N
	th=theta[kt1]
	r=rhomin+kr*(rhomax-rhomin)/(M-1)
	;Look for the endpoints of the line. These are the values
	;returned to describe the line.
	IF dx GE dy THEN BEGIN
		x0=xmin & x1=xmax
		y0=(r-x0*cos(th))/sin(th)
		y1=(r-x1*cos(th))/sin(th)
	ENDIF ELSE BEGIN
		y0=ymin & y1=ymax
		x0=(r-y0*sin(th))/cos(th)
		x1=(r-y1*sin(th))/cos(th)
	ENDELSE
	linesFound=linesFound+1
	IF verbose GE 1 THEN BEGIN
	PRINT,' '
	PRINT,'===================================================='
	PRINT,'Line Number ',LinesFound
	PRINT,FORMAT='("Hough Index=",i7,"  Number of points=",i7)',[kp,peak]
	PRINT,'Theta of Line=',180*th/!pi
	PRINT,'Rho of Line=',r
	PRINT,'Endpoints: ',x0,y0,x1,y1
	PRINT,'Length: ',linelen
	ENDIF
	IF verbose EQ 2 THEN BEGIN
		PRINT,'X Coordinates=',xv
		PRINT,'Y Coordinates=',yv
	ENDIF

	;We like the line. Keep the data in the structure.
	lines[linesFound-1].Index=linesFound
	lines[linesFound-1].Xends=[x0,x1]
	lines[linesFound-1].Yends=[y0,y1]
	lines[linesFound-1].THETA=180*th/!pi
	lines[linesFound-1].RHO=r
	lines[linesFound-1].Length=linelen
	lines[linesFound-1].Location=kp
	;Set the neighborhood values to -1 and the point value
	;to -2 to mark them. These values will not be used
	;again in the search.
	A[nbrs]=-1
	A[kp]=-2
ENDIF

ENDREP UNTIL linesFound GE nlines OR peak LT thresh
IF linesFound GT 0 THEN $
RETURN,lines[0:linesFound-1] ELSE RETURN,-1
END