Identifying the "bendy points" on a path
Posted: Wed Mar 14, 2018 9:18 pm
I have an array of XY points, making a continuous path. It is a completely arbitrary path. (For example, let's say the path of a river, with sharp corners, smooth bends, etc. all quite random.)
What I need to do is identify the points along the path where a "bend" occurs.
So, if this is the path:

then I would want the code to identify the following points:

I have written a procedure - GetPathBends - but it doesn't work very well. It works by going along the path and, for every pixel, looking backwards along the path and forwards along the path by #JudgeDistancePixels, getting the angle
Here is a demo. You need to draw a path with the mouse. A procedure will join up the dots to make it a continuous path, and then the GetPathBends procedure will be run and the bends will be illustrated with red circles.
The code does correctly identify the "bendy points", but it also comes up with false positives, and there is still the problem of clusters of bend points, instead of the central one in each cluster getting identified.
Any help people could give me with this would be appreciated.
What I need to do is identify the points along the path where a "bend" occurs.
So, if this is the path:

then I would want the code to identify the following points:

I have written a procedure - GetPathBends - but it doesn't work very well. It works by going along the path and, for every pixel, looking backwards along the path and forwards along the path by #JudgeDistancePixels, getting the angle
Here is a demo. You need to draw a path with the mouse. A procedure will join up the dots to make it a continuous path, and then the GetPathBends procedure will be run and the bends will be illustrated with red circles.
Code: Select all
Macro R(t)
MessageRequester("Report",t,0)
EndMacro
Procedure.d Beat(a.d,b.d)
If a>b
ProcedureReturn a
Else
ProcedureReturn b
EndIf
EndProcedure
Procedure.d Defeat(a.d,b.d)
If a<b
ProcedureReturn a
Else
ProcedureReturn b
EndIf
EndProcedure
Macro DefeatThis(a,b)
If a>b
a=b
EndIf
;a = Defeat(a,b)
EndMacro
Macro BeatThis(a,b)
If b>a
a=b
EndIf
;a = Beat(a,b)
EndMacro
Procedure.d Difference(a.d,b.d)
If a=b
ProcedureReturn 0
EndIf
If a>b
ProcedureReturn a-b
EndIf
ProcedureReturn b-a
EndProcedure
Structure PointD
x.d
y.d
EndStructure
Structure PointI
x.i
y.i
EndStructure
Procedure.d DegreeAngleBetweenTwoPoints(*o.PointD,*b.PointD)
calcAngle.d = Degree(ATan2(*b\x-*o\x,*b\y-*o\y))
If calcAngle<0
calcAngle = 360-Abs(calcAngle)
EndIf
ProcedureReturn calcAngle
EndProcedure
Macro CopyPoint(p1,p2)
p2\x = p1\x
p2\y = p1\y
EndMacro
Macro MoreThanOnePixelAway(x1,y1,x2,y2)
( Difference(x1,x2)>1 Or Difference(y1,y2)>1 )
EndMacro
Procedure.b IPointsOnePixelApart(*p1.PointI,*p2.PointI)
xdif.i = Difference(*p1\x,*p2\x)
ydif.i = Difference(*p1\y,*p2\y)
If xdif=1 And ydif=1
ProcedureReturn #True
EndIf
If xdif=1 And ydif=0
ProcedureReturn #True
EndIf
If xdif=0 And ydif=1
ProcedureReturn #True
EndIf
ProcedureReturn #False
EndProcedure
Procedure.i BresenhamLine(Array pnt.PointI(1),*ip1.PointI,*ip2.PointI)
; Line drawing algorithm in Pascal -> http://www.gamedev.net/
; A general-purpose implementation of Bresenham's line algorithm.
; By Mark Feldman.
bplot_x1.i = *ip1\x
bplot_y1.i = *ip1\y
bplot_x2.i = *ip2\x
bplot_y2.i = *ip2\y
;deltax and deltay
bdx = Abs(bplot_x2 - bplot_x1)
bdy = Abs(bplot_y2 - bplot_y1)
;always 1
bx2 = 1
by2 = 1
;initialize vars based on which is the independent variable
If bdx >= bdy ;x is independent variable
bnp = bdx + 1
bd1 = bdy << 1
bdi = bd1 - bdx
bd2 = (bdy - bdx) << 1
bx1 = 1
by1 = 0
Else ;y is independent variable
bnp = bdy + 1
bd1 = bdx << 1
bdi = bd1 - bdy
bd2 = (bdx - bdy) << 1
bx1 = 0
by1 = 1
EndIf
;make sure x and y move in the right directions
If bplot_x1 > bplot_x2
bx1 = -bx1
bx2 = -bx2
EndIf
If bplot_y1 > bplot_y2
by1 = -by1
by2 = -by2
EndIf
;start drawing at
bxi = bplot_x1
byi = bplot_y1
;draw the pixels
pnts.i = 0
ReDim pnt(bnp)
For bli = 1 To bnp
pnts+1
pnt(pnts)\x = bxi
pnt(pnts)\y = byi
;PlotPixel(bxi, byi, bclr)
If bdi < 0
bdi = bdi + bd1
bxi = bxi + bx1
byi = byi + by1
Else
bdi = bdi + bd2
bxi = bxi + bx2
byi = byi + by2
EndIf
Next bli
ProcedureReturn pnts
EndProcedure
Macro PointsIdentical(p1,p2)
Bool(p1\x=p2\x And p1\y=p2\y)
EndMacro
Procedure.i JoinUpDots(pathpoints,Array pnt.PointI(1))
;R("pathpoints: "+Str(pathpoints))
Dim pxl.PointI(pathpoints*4)
pxls.i = 1
CopyPoint(pnt(1),pxl(pxls))
For p = 2 To pathpoints
If IPointsOnePixelApart(pnt(p-1),pnt(p))
pxls+1
CopyPoint(pnt(p),pxl(pxls))
EndIf
If MoreThanOnePixelAway(pnt(p-1)\x,pnt(p)\x,pnt(p-1)\y,pnt(p)\y)
Dim np.PointI(10)
newpoints = BresenhamLine(np(),pnt(p-1),pnt(p))
ReDim pxl(ArraySize(pxl())+newpoints+10)
For n = 1 To newpoints
pxls+1
CopyPoint(np(n),pxl(pxls))
Next n
Else
pxls+1
CopyPoint(pnt(p),pxl(pxls))
EndIf
Next p
ReDim pnt(pxls)
lp.i = 1
npxls.i = 1
CopyPoint(pxl(lp),pnt(npxls))
For p = 2 To pxls
If PointsIdentical(pxl(lp),pxl(p))
Continue
EndIf
npxls+1
CopyPoint(pxl(p),pnt(npxls))
lp = p
Next p
ReDim pnt(npxls)
ProcedureReturn npxls
;Debug "CORRECTION COMPLETE."
EndProcedure
Procedure.b GetPathBends(pathpoints.i,Array ipathpnt.PointI(1),List bend.i())
; convert from integers to doubles, so that we can do Math operations
Dim pathpnt.PointD(pathpoints)
For p = 1 To pathpoints
CopyPoint(ipathpnt(p),pathpnt(p))
Next p
; for every point along the path, jump backward #JudgeDistance points and forward #JudgeDistance points and compare the two angles (current point being the central one of three points)
; if the difference between the two angles is greater than 50 degrees, mark it as a bend point
Dim rawbend.b(pathpoints)
#JudgeDistance = 50
For p = #JudgeDistance+1 To pathpoints-#JudgeDistance
ang1.d = DegreeAngleBetweenTwoPoints(@pathpnt(p-#JudgeDistance),@pathpnt(p))
ang2.d = DegreeAngleBetweenTwoPoints(@pathpnt(p),@pathpnt(p+#JudgeDistance))
dif.d = Difference(ang1,ang2)
If dif>50
rawbend(p) = #True
EndIf
Next p
; since the code above tends to find clusters of "bendpoints", we need to look at each cluster and calculate the middle point - where the bend actually happens
Dim majorbend.b(pathpoints)
close_together_bends = 0
consecutive_straights = 0
bendchainstart = 0
For p = 1 To pathpoints
If rawbend(p)
consecutive_straights = 0
If Not bendchainstart
bendchainstart = p
EndIf
bendchainend = p
close_together_bends + 1
Else
consecutive_straights + 1
If consecutive_straights>25
;Debug "CTB: "+Str(close_together_bends)
If close_together_bends>2
p = bendchainstart + ((bendchainend-bendchainstart)/2)
;Debug "MAJOR BEND: "+Str(p)
majorbend(p) = #True
close_together_bends = 0
consecutive_straights = 0
bendchainstart = 0
bendchainend = 0
EndIf
EndIf
EndIf
Next p
ClearList(bend())
For p = 1 To pathpoints
If majorbend(p)
AddElement(bend())
bend() = p
EndIf
Next p
EndProcedure
ww = 1000
wh = 700
win = OpenWindow(#PB_Any,0,0,ww,wh,"",#PB_Window_ScreenCentered|#PB_Window_BorderLess)
cnv = CanvasGadget(#PB_Any,0,0,ww,wh)
StartDrawing(CanvasOutput(cnv))
Box(0,0,ww,wh,RGB(255,255,192))
StopDrawing()
pnts=0
Dim pathpnt.PointI(5000)
NewList bend.i()
Repeat
we = WindowEvent()
Select we
Case #PB_Event_Gadget
Select EventGadget()
Case cnv
Select EventType()
Case #PB_EventType_LeftButtonDown
tracking = #True
StartDrawing(CanvasOutput(cnv))
Box(0,0,ww,wh,RGB(255,255,192))
StopDrawing()
Case #PB_EventType_LeftButtonUp
tracking = #False
pnts = JoinUpDots(pnts,pathpnt())
GetPathBends(pnts,pathpnt(),bend())
StartDrawing(CanvasOutput(cnv))
Box(0,0,ww,wh,RGB(255,255,192))
; draw the path
For p = 1 To pnts
Plot(pathpnt(p)\x,pathpnt(p)\y,#Black)
Next p
; illustrate the bends
ForEach bend()
p = bend()
Circle(pathpnt(p)\x,pathpnt(p)\y,5,#Red)
Next
StopDrawing()
pnts=0
Case #PB_EventType_MouseMove
If tracking
x = GetGadgetAttribute(cnv,#PB_Canvas_MouseX)
y = GetGadgetAttribute(cnv,#PB_Canvas_MouseY)
If pnts=0 Or x<>pathpnt(pnts)\x Or y<>pathpnt(pnts)\y
;Debug Str(x)+","+Str(y)
pnts+1
pathpnt(pnts)\x = x
pathpnt(pnts)\y = y
StartDrawing(CanvasOutput(cnv))
Circle(pathpnt(pnts)\x,pathpnt(pnts)\y,3,#Black)
StopDrawing()
EndIf
EndIf
EndSelect
EndSelect
EndSelect
Delay(5)
Until GetAsyncKeyState_(#VK_ESCAPE)
Any help people could give me with this would be appreciated.