FUNCTION points_near_line,a0,b0,a1,b1,pts,delta
;+
;R=points_near_line(x0,y0,x1,y1,P,delta) or
;R=points_near_line(V,pts,delta) returns the
;points in the array P that are near a line described
;by the endpoints [x0,y0], [x1,y1]. The distance threshold
;is set by delta. The endpoints may be supplied as
;the four-parameters x0,y0,x1,y1 or a vector V=[x0,y0,x1,y1].
;
;P must be a [2,n] array of n points. R is the index list
;of the points in P that are near the line. The points
;may be selected by Q=P[*,R]
;
;H. Rhody
;July, 2004
;-

IF (n_params() EQ 3) THEN BEGIN
	x0=a0 & y0=b0 & x1=a1
	delta=X1
	pts=Y0
	y1=x0[3]
	x1=x0[2]
	y0=x0[1]
	x0=x0[0]
ENDIF ELSE BEGIN
	x0=a0 & y0=b0 & x1=a1 & y1=b1
ENDELSE

;Find a unit vector in the direction of the line

u=[x1-x0,y1-y0]/sqrt((float(x1)-float(x0))^2+(float(y1)-float(y0))^2)

;Orthogonal unit vector
v=Transpose([u[1],-u[0]])

;Find a vector from [x0,y0] to each point in the list.
P=[pts[0,*]-x0,pts[1,*]-y0]

;Project P onto V
Q=P##v
;Find the orthogonal component
S=[Q*P[0,*],Q*P[1,*]]

;Find the length of each vector in the Q list.
L=sqrt(total(S^2,1))

;Find the vectors in L that are shorter than delta
R=where(L LE delta)

ptestx=Pts[0,R]
ptesty=Pts[1,R]
xmin=x0<x1 & xmax=x0>x1 & ymin=y0<y1 & ymax=y0>y1
ip=where(xmin LE ptestx AND xmax GE ptestx $
	AND ymin LE ptesty AND ymax GE ptesty)
RI=R[ip]
return,RI
END