Re: Recursive Function Program in IDL [message #46797 is a reply to message #46710] |
Wed, 14 December 2005 11:43  |
b_gom
Messages: 105 Registered: April 2003
|
Senior Member |
|
|
David Fanning wrote:
> Does anyone have a handy recursive function that does something neat?
> Someone is asking, and I don't have time to work on this.
Here is a recursive algorithm for simplifying polyline vertices. See
'poly_simplify.pro' in the IDL user contributed library.
http://www.rsinc.com/codebank/search.asp?FID=307
Brad
;**************************************
pro simplifyDP,tol,vertices,j,k,mk
; This is the Douglas-Peucker recursive simplification routine
; It just marks vertices that are part of the simplified polyline
; for approximating the polyline subchain vertices[j] to vertices[k].
; Input: tol = approximation tolerance
; vertices[] = polyline array of vertex points
; j,k = indices for the subchain vertices[j] to vertices[k]
; Output: mk[] = array of markers matching vertex array vertices[]
if (k le j+1) then return ; there is nothing to simplify
; check for adequate approximation by segment S from vertices[j] to
vertices[k]
maxi = j ; index of vertex farthest from S
maxd2 = 0. ; distance squared of farthest vertex
S = [[vertices[*,j]], [vertices[*,k]]] ; segment from vertices[j] to
vertices[k]
u = S[*,1]-S[*,0] ; segment direction vector
cu = ps_dot(u,u); segment length squared
;test each vertex vertices[i] for max distance from S
;compute using the Feb 2001 Algorithm's dist_Point_to_Segment()
;Note: this works in any dimension (2D, 3D, ...)
;Pb = base of perpendicular from vertices[i] to S
;dv2 = distance vertices[i] to S squared
for i=j+1,k-1 do begin
;compute distance squared
w = vertices[*,i] - S[*,0]
cw = ps_dot(w,u)
if cw le 0 then begin
dv2 = ps_d2(vertices[*,i], S[*,0]);
endif else begin
if cu le cw then begin
dv2 = ps_d2(vertices[*,i], S[*,1])
endif else begin
b = cw / cu;
Pb = S[*,0] + b * u;
dv2 = ps_d2(vertices[*,i], Pb);
endelse
endelse
;test with current max distance squared
if dv2 le maxd2 then continue
;vertices[i] is a new max vertex
maxi = i
maxd2 = dv2
endfor
if (maxd2 gt tol^2) then begin ;// error is worse than the tolerance
; split the polyline at the farthest vertex from S
mk[maxi] = 1 ; mark vertices[maxi] for the simplified polyline
; recursively simplify the two subpolylines at vertices[*,maxi]
simplifyDP, tol, vertices, j, maxi, mk ; // polyline vertices[j] to
vertices[maxi]
simplifyDP, tol, vertices, maxi, k, mk ; // polyline vertices[maxi]
to vertices[k]
endif
; else the approximation is OK, so ignore intermediate vertices
return
end
|
|
|