Message 1 of 2
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report
I am in the process of creating a lisp to help with decimating polylines with many points while minimizing the variation from the original polyline. I am using the Ramer–Douglas–Peucker algorithm. Reference can be found here. I think there is an issue with the recursion part and outputting the new polyline points. I would also appreciate any advice in cleaning up the code in general and simplifying it. Below is the code I have currently written as well as a version I have previously made to work in excel for reference.
; # source: https://karthaus.nl/rdp/
; function DouglasPeucker(PointList[], epsilon)
; # Find the point with the maximum distance
; dmax = 0
; index = 0
; end = length(PointList)
; for i = 2 to (end - 1) {
; d = ShortestDistance(PointList[i], Line(PointList[1], PointList[end]))
; if (d > dmax) {
; index = i
; dmax = d
; }
; }
;
; ResultList[] = empty
;
; # If max distance is greater than epsilon, recursively simplify
; if (dmax > epsilon) {
; # Recursive call
; recResults1[] = DouglasPeucker(PointList[1...index], epsilon)
; recResults2[] = DouglasPeucker(PointList[index...end], epsilon)
;
; # Build the result list
; ResultList[] = {recResults1[1...length(recResults1) - 1], recResults2[1...length(recResults2)]}
; } else {
; ResultList[] = {PointList[1], PointList[end]}
; }
; # Return the result
; return ResultList[]
;
;--------------------------------------------------------------------------------------------------------------------
(defun C:RDP_Decimate (/ ent PL0 PL RowLast epsilon output)
(setq ent (entsel "\nPick POLYLINE to list vertices:"))
(if
(and ent
(wcmatch (cdr (assoc 0 (entget (setq ent (car ent))))) "*POLYLINE") ;all types
) ;and
(setq PL (PolyVert ent)) ;make list
)
;(setq PLz (mapcar 'caddr PL))
;(setq PL (mapcar 'reverse (mapcar 'cdr (mapcar 'reverse PL0)))) ;remove z coordinates
(setq RowLast (- (length PL) 1) ;end = length(PointList) // subtracting 1 because indices start at 0 while length command starts at 1
epsilon (getreal "\nEpsilon:") ;set value of epsilon
)
(setq output (RDP PL epsilon RowLast)) ;save the result
;(princ output) ;print the result (this is an error checking step during development. Comment out at end)
(PLDraw output) ;draw the result
)
;--------------------------------------------------------------------------------------------------------------------
(defun RDP (PLxyz eps end / PLx PLy dmax d index iter ArrRes1 ArrRes2 RecRes1 RecRes2
Res
)
(setq PLx (mapcar 'car PLxyz)
PLy (mapcar 'cadr PLxyz)
dmax 0.0
index 0
iter 1 ;this is assuming that the indices start at 0
)
(while (<= iter (- end 1)) ;for i = 2 to (end - 1)
(setq d (MinDist PLx PLy)) ;d = ShortestDistance(PointList[i], Line(PointList[1], PointList[end]))
(if (> d dmax) ;if (d > dmax) {
(setq index iter ;index = i
dmax d ;dmax = d
)
)
(setq iter (1+ iter))
)
(if (> end index) ;assuming indices start at 0
(progn
;ResultList[] = empty
(setq ArrRes1 (car (ArrCut PLxyz index)) ;ArrRes1=PLxys(1,...,index)
ArrRes2 (cadr (ArrCut PLxyz index)) ;ArrRes2=PLxyz(index,...,end)
;ArrRes2 (member (nth index PLxyz) PLxyz) ;ArrRes2=PLxyz(index,...,end)
;ArrRes1 (reverse (member (car ArrRes2) (reverse PLxyz))) ;ArrRes1=PLxys(1,...,index)
)
(if (> dmax eps) ;if (dmax > epsilon) {
;# Recursive call
(setq RecRes1 (RDP ArrRes1 eps (length ArrRes1)) ;recResults1[] = DouglasPeucker(PointList[1...index], epsilon)
RecRes2 (RDP ArrRes2 eps (length ArrRes2)) ;recResults2[] = DouglasPeucker(PointList[index...end], epsilon)
;# Build the result list. Note that the last term of RecRes1 should be the same as the first term of RecRes2
Res (append RecRes1 (cdr RecRes2)) ;ResultList[] = {recResults1[1...length(recResults1) - 1], recResults2[1...length(recResults2)]}
) ;}
(setq Res (list (car PLxyz) (last PLxyz))) ;else {ResultList[] = {PointList[1], PointList[end]}
) ;}
)
(setq Res PLxyz)
)
;sort the result list from least to greatest x values
;(vl-sort Res
; (function
; (lambda (e1 e2)
; (< (car e1) (car e2))
; )
; )
;)
) ;return Res
;--------------------------------------------------------------------------------------------------------------------
; PolyVert
; returns vertices list for any type of polyline (in WCS coordinates)
; arg POLY: polyline to list (ename or vla-object)
(defun PolyVert (POLY / par pt1 lst)
(vl-load-com)
(setq par (if (vlax-curve-isClosed POLY)
(vlax-curve-getEndParam POLY) ; else
(1+ (vlax-curve-getEndParam POLY))
)
)
(while (setq pt1 (vlax-curve-getPointAtParam POLY (setq par (- par 1))))
(setq lst (cons pt1 lst))
)
) ;; return lst
;--------------------------------------------------------------------------------------------------------------------
; If the line passes through two points P1=(x1,y1) and P2=(x2,y2) then the distance of (x0,y0) from the line is:
; MinDist(P1,P2,(x0,y0))=Abs((y2-y1)*x0-(x2-x1)*y0+x2*y1-y2*x1)/sqrt((y2-y1)^2+(x2-x1)^2)
(defun MinDist (ArgX ArgY / y0 x0 y1 x1 y2 x2)
(setq y0 (nth iter ArgY)
x0 (nth iter ArgX)
y1 (car ArgY)
x1 (car ArgX)
y2 (last ArgY)
x2 (last ArgX)
)
; MinDist(P1,P2,(x0,y0))=Abs((y2-y1)*x0-(x2-x1)*y0+x2*y1-y2*x1)/sqrt((y2-y1)^2+(x2-x1)^2)
; Abs((y2-y1)*x0-(x2-x1)*y0+x2*y1-y2*x1)->(abs (+ (- (* (- y2 y1) x0) (* (- x2 x1) y0) (* y2 x1)) (* x2 y1)))
; sqrt((y2-y1)^2+(x2-x1)^2)-> (+ (expt (- y2 y1) 2) (expt (- x2 x1) 2))
(setq d (/
(abs (+ (- (* (- y2 y1) x0) (* (- x2 x1) y0) (* y2 x1)) (* x2 y1)))
(+ (expt (- y2 y1) 2) (expt (- x2 x1) 2))
)
)
)
;--------------------------------------------------------------------------------------------------------------------
;draw the polyline
(defun PLDraw (listpts / x)
(command "_pline") ;start pline command
(while (= (getvar 'cmdactive) 1) ;while a command is running loop (it is because we haven't entered any points)
(repeat (setq x (length listpts)) ;repeat loop (sets x to length of your list) repeat x times
(command (nth (setq x (- x 1)) listpts)) ;because this is still inside the pline command it evaluates the expression (nth item (x-1) of list listpts
) ;end repeat ;length gives you number of items in list but list is 0 based so first item is 0 and last is (x-1)
(command "") ;end pline command with a return ("")
) ;end while
)
;--------------------------------------------------------------------------------------------------------------------
;Cuts an array at a specified index and returns a list of 2 arrays that start and end at the given index respectively
;Note that indices start at 0 when choosing the CutIndex
(defun ArrCut (Arr CutIndex / Arr1 Arr2 ArrList)
(setq Arr2 (member (nth CutIndex Arr) Arr) ;Arr2=PLxyz(index,...,end)
Arr1 (reverse (member (car Arr2) (reverse Arr))) ;ArrRes1=PLxys(1,...,index)
)
(setq ArrList (list Arr1 Arr2))
)
Sub CallRamerDouglasPeucker()
Dim pointList As Variant
Dim rowCount As Integer
Dim result As Variant
Dim epsilon As Double
Call ClearOutput
pointList = Sheets("Full").Range("Full")
rowCount = UBound(pointList)
epsilon = Sheets("Full").Range("epsilon")
result = RamerDouglasPeucker(pointList, epsilon, rowCount)
WriteResult (result)
End Sub
Function RamerDouglasPeucker(pointList As Variant, epsilon As Double, rowCount As Integer) As Variant
Dim dMax As Double
Dim Index As Integer
Dim d As Double
Dim arrResults1 As Variant
Dim arrResults2 As Variant
Dim resultList As Variant
Dim recResults1 As Variant
Dim recResults2 As Variant
' Find the point with the maximum distance
dMax = 0
Index = 0
For i = 2 To (rowCount - 1)
d = Abs((pointList(rowCount, 1) - pointList(1, 1)) * (pointList(1, 2) - pointList(i, 2)) - (pointList(1, 1) - pointList(i, 1)) * (pointList(rowCount, 2) - pointList(1, 2))) / Sqr((pointList(rowCount, 1) - pointList(1, 1)) ^ 2 + (pointList(rowCount, 2) - pointList(1, 2)) ^ 2)
If d > dMax Then
Index = i
dMax = d
End If
Next i
'Testing if can stop cut going to index 0
If Index > 1 Then
arrResults1 = Cut_Array(pointList, 1, Index)
arrResults2 = Cut_Array(pointList, Index, rowCount)
' If max distance is greater than epsilon, recursively simplify
If (dMax > epsilon) Then
' Recursive call
recResults1 = RamerDouglasPeucker(arrResults1, epsilon, UBound(arrResults1))
recResults2 = RamerDouglasPeucker(arrResults2, epsilon, UBound(arrResults2))
' Build the result list
resultList = Join_Array(recResults1, recResults2)
Else
ReDim resultList(1 To 2, 1 To 2)
resultList(1, 1) = pointList(1, 1)
resultList(1, 2) = pointList(1, 2)
resultList(2, 1) = pointList(rowCount, 1)
resultList(2, 2) = pointList(rowCount, 2)
End If
Else
resultList = pointList
End If
' Return the result
RamerDouglasPeucker = resultList
End Function
Function Cut_Array(arr As Variant, arrStart As Integer, arrEnd As Integer) As Variant
Dim resultList As Variant
ReDim resultList(1 To (arrEnd - arrStart) + 1, 1 To 2)
For i = arrStart To arrEnd
For j = 1 To 2
resultList((i - arrStart) + 1, j) = arr(i, j)
Next j
Next i
Cut_Array = resultList
End Function
Function Join_Array(arr1 As Variant, arr2 As Variant) As Variant
Dim resultList As Variant
Dim arr1Length As Integer
Dim arr2Length As Integer
arr1Length = UBound(arr1) - 1
arr2Length = UBound(arr2)
newArrLength = arr1Length + arr2Length
ReDim resultList(1 To newArrLength, 1 To 2)
For i = 1 To arr1Length
For j = 1 To 2
resultList(i, j) = arr1(i, j)
Next j
Next i
For i = 1 To arr2Length
For j = 1 To 2
resultList(i + arr1Length, j) = arr2(i, j)
Next j
Next i
Join_Array = resultList
End Function
Sub WriteResult(result As Variant)
Dim Filtered As ListObject
Set Filtered = Sheets("Filtered").ListObjects("Filtered")
'Copy information loop
For i = 1 To UBound(result)
For j = 1 To UBound(result, 2)
Filtered.Range.Cells(i + 1, j).Value = result(i, j)
Next j
Next i
End Sub
Sub ClearOutput()
Dim OutputTable As ListObject
Dim StartRow As Integer
Set OutputTable = Sheets("Filtered").ListObjects("Filtered")
StartRow = OutputTable.Range.Cells(1, 1).Row + 1
If Not OutputTable.InsertRowRange Is Nothing Then
'Pass
Else
Sheets("Filtered").Rows(StartRow & ":" & (OutputTable.DataBodyRange.Rows.Count + StartRow)).Delete
End If
End Sub
Solved! Go to Solution.
