Objectives:
1 Review recursive rules in Douglas algorithm (Peucker reading)
2 Show how it is implemented in code (without recursion)
3 Application of the distance-from-line calculation
- Trend line: null hypothesis -
- deviations from trend less than tolerance are insignificant.
- Trend lines constructed between all pairs found significant (iterative process)
Results with different tolerances reflect
different amounts of generalization...
Douglas's version (REDUCE) of the original paper:
Douglas David H. and Peucker, Thomas 1974: Algorithms for the
reduction of the number of points required to represent a digitized
line or its caricature, Canadian Cartographer 11, 112-122.
SUBROUTINE REDUCT (POINTS, N, KATCH, IP, DMIN, SHORT)
C APPLY PEUCKER-DOUGLAS LINE REDUCTION ALGORITHM.
C ADAPTED FOR THE FLECS LANGUAGE FROM PREVIOUS FORTRAN
C VERSIONS BY DOUGLAS AND CHRISMAN.
C
C POINTS -- ARRAY OF COORDINATES (X,Y,...X,Y)
C N ------- THE LENGTH OF POINTS
C KATCH --- WORKSPACE MUST BE AT LEAST N/2 LONG
C IP ------ NUMBER OF POINTS (RETURNED)
C DMIN ---- GENERALIZATION TOLERANCE
C SHORT --- LENGTH OF LINE AT MACHINE PRECISION
C
C R D WHITE, LAB FOR COMPUTER GRAPHICS, SEPTEMBER 1976
C
DIMENSION POINTS (N), KATCH (1)
LOGICAL STACKD
DATA BIG /9E28/
C INITIALIZE THE STACK WITH NODES (N.B.
Commentary)
I = 1
IP = 1
NP = N / 2
JP = NP
KATCH (IP) = 1 accept first point
M = N - 1
KATCH (JP) = M Prime stack with last
STACKD = .TRUE.
C REPEAT UNTIL EVERYTHING IS UNSTACKED
10 IF(.NOT.(STACKD)) GO TO 140 REPEAT
UNTIL done...
C GET TOP OF STACK Each trend line:
I to J
J = KATCH (JP)
C IF THE PROPOSED TREND LINE HAS MORE THAN TWO POINTS
DMAX = 0.0
IF(.NOT.(J .GT. I + 2)) GO TO 90 escape
if nothing between
C COMPUTE PARAMETERS OF LINE
A = (POINTS (J + 1) - POINTS (I + 1))
B = (POINTS (I) - POINTS (J))
C = POINTS (J) * POINTS (I + 1)
C = (C - POINTS (I) * POINTS (J + 1))
D = SQRT (A * A + B * B)
K = I + 2 Set K for LOOPs below
C LOCATE MAXIMUM DEVIANT
IF(.NOT.(D .GT. SHORT)) GO TO 50
C NORMAL CASE distance J, I > SHORT
20 IF(K .GE. J) GO TO 40 LOOP: K from
I to J
DI = A * POINTS(K) + B * POINTS(K + 1)
DI = ABS ((DI + C) / D) distance from
line to K
IF(.NOT.(DI .GT. DMAX)) GO TO 30
DMAX = DI save new max distance
IJ = K save point as IJ
30 K = K + 2 Next K (xy =2)
GO TO 20
40 GO TO 80
C VERY SHORT TREND LINE ELSE (dist < short)
50 IF(K .GE. J) GO TO 70 same loop
as normal case
DX = (POINTS (I) - POINTS (K))
DY = (POINTS (I + 1) - POINTS (K + 1))
DI = SQRT (DX * DX + DY * DY) distance
from I
IF(.NOT.(DI .GT. DMAX)) GO TO 60
DMAX = DI
IJ = K
60 K = K + 2
GO TO 50
70 CONTINUE
80 CONTINUE
C TEST DEVIANT AGAINST MINIMUM SIGNIFICANT DEVIATION
90 IF(.NOT.(DMAX .GT. DMIN)) GO TO 100
C PUSH STACK This accepts point IJ
JP = JP - 1new trends will be
KATCH (JP) = IJ I to IJ; IJ to J
GO TO 130
C POP STACK
100 IP = IP + 1 nothing between I,J
JP = JP + 1 J moves into final list
KATCH (IP) = J
IF(.NOT.(J .GE. M)) GO TO 110
STACKD = .FALSE. done when last point
(M) popped
GO TO 120
110 I = J if anything left on stack,
old J becomes I
120 CONTINUE
130 GO TO 10
C
C NOW CRUNCH OUT UNNEEDED POINTS superfluous?
140 J=1
DO 150 I=1,IP
K=KATCH(I)
POINTS(J)=POINTS(K)
POINTS(J+1)=POINTS(K+1)
J=J+2
150 CONTINUE
RETURN
END
maintains all lines (no topological changes)
can produce intersections with adjacent features (eg. West Virginia
crossing into PA)
can accentuate zig-zags...
Index from here:
Version of 25 January 1999
