Identifying the "bendy points" on a path

Just starting out? Need help? Post your questions and find answers here.
Seymour Clufley
Addict
Addict
Posts: 1233
Joined: Wed Feb 28, 2007 9:13 am
Location: London

Identifying the "bendy points" on a path

Post by Seymour Clufley »

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:
Image

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

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)
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.
JACK WEBB: "Coding in C is like sculpting a statue using only sandpaper. You can do it, but the result wouldn't be any better. So why bother? Just use the right tools and get the job done."
User avatar
skywalk
Addict
Addict
Posts: 3999
Joined: Wed Dec 23, 2009 10:14 pm
Location: Boston, MA

Re: Identifying the "bendy points" on a path

Post by skywalk »

Fun problem. Remembering way back in calculus, you find inflection points("bendy") with the 2nd derivative.
You need a derivative(dY/dX) procedure and pass in your path coordinates. Then pass the resulting 1st derivative coordinates to get the 2nd derivative. Then scan the 2nd derivative for peaks.

The problems with discrete derivatives are the need for continuous values. The incremental differences of each coordinate point(dY/dX) will blow up if there are 2 identical X's. So, I guess you would start a new path when that occurs.
The nice thing about standards is there are so many to choose from. ~ Andrew Tanenbaum
User avatar
Michael Vogel
Addict
Addict
Posts: 2677
Joined: Thu Feb 09, 2006 11:27 pm
Contact:

Re: Identifying the "bendy points" on a path

Post by Michael Vogel »

Maybe something like this could be an approach...

Step 1 - get line data (black and white image only)...

Code: Select all

#Dot=#Black
#Paper=#White

Global sx,sy
Global x,y
Global w,h

Global Dim nx(8)
Global Dim ny(8)

Global NewList Poly.Point()

For i=1 To 8
	Read.i nx(i)
	Read.i ny(i)
Next i

DataSection
	Data.i 0,1, 0,-1, 1,0, -1,0, -1,-1, 1,-1, 1,1, -1,1
EndDataSection

Procedure GetStart()

	StartDrawing(ImageOutput(0))

	z=1
	x=w
	While z And x
		x-1
		y=h
		While z And y
			y-1
			If Point(x,y)=#Dot
				z=0
				Break
			EndIf
		Wend
	Wend

	StopDrawing()

EndProcedure
Procedure GetNeighbour(n)

	Protected i

	For i=1 To 8
		If Point(x+nx(i),y+ny(i))=#Dot
			ProcedureReturn i
		EndIf
	Next i

	ProcedureReturn #False

EndProcedure
Procedure GetLine()

	Protected n,i

	StartDrawing(ImageOutput(0))

	For i=0 To 1
		x=sx
		y=sy


		Repeat
			n=GetNeighbour(n)
			If n
				Plot(x,y,#Green)
				x+nx(n)
				y+ny(n)
				AddElement(Poly())
				Poly()\x=x
				Poly()\y=y
			EndIf
		Until n=0

	Next i

	StopDrawing()


EndProcedure

Procedure main()

	LoadImage(0,"test.bmp")

	w=ImageWidth(0)
	h=ImageHeight(0)

	GetStart()
	sx=x
	sy=y

	GetLine()

	OpenWindow(0,0,0,w,h,":)")
	CanvasGadget(0,0,0,w,h)

	StartDrawing(CanvasOutput(0))
	DrawImage(ImageID(0),0,0)
	Circle(sx,sy,3,#Red)
	StopDrawing()

	Repeat : Until WaitWindowEvent()=#PB_Event_CloseWindow


EndProcedure

main()
Step 2 Reducing lines (lines, not curves)...
See here...
User avatar
Derren
Enthusiast
Enthusiast
Posts: 313
Joined: Sat Jul 23, 2011 1:13 am
Location: Germany

Re: Identifying the "bendy points" on a path

Post by Derren »

I don't get how you identify bends.
It's not a mathematical "twist", where a curve changes from a left-turn into a right-turn.

the first red point (top left) would not identify such a twist.
Yet, why are there 3 points in the "lagoon", but there's not point between the 1st and 2nd red dot? There's quite a big "bend" there, I would say.

Do you want want to use Bezier-curves to aproximate the arbitrary curve?
User avatar
Michael Vogel
Addict
Addict
Posts: 2677
Joined: Thu Feb 09, 2006 11:27 pm
Contact:

Re: Identifying the "bendy points" on a path

Post by Michael Vogel »

Here's a quick demo for 'step 2' (use keys 'b' and 'n' to modify the number of points)...

Code: Select all

DataSection
	ShapeData:
	Data.w 435, 588,443, 583,444, 578,445, 573,446, 568,449, 563,452, 558,453, 553,456, 548,459, 543,462, 539,467, 535,472, 531,477, 527,482, 523,487, 518,489, 514,494, 509,497, 504,500, 499,503, 495,508, 490,510, 485,512, 480,515, 476,520, 471,521, 466,524, 461,527, 457,532, 452,534, 447,538, 442,540, 437,543, 433,548, 428,551, 424,556, 419,558, 415,563, 410,563, 405,565, 400,566, 395,567, 390,569, 385,570, 380,572, 375,574, 370,576, 365,578, 360,580, 355,582, 350,584, 345,586, 340,587, 335,588, 330,588, 325,589, 320,589, 315,589, 310,589, 305,590, 300,590, 295,591, 290,593, 285,594, 280,595, 275,597, 270,599, 265,600, 260,602, 255,604, 250,604, 245,605, 240,606, 235,606, 230,606, 225,606, 220,606, 215,605, 210,604, 205,600, 200,598, 195,597, 190,594, 185,590, 181,585, 180,580, 182,575, 185,570, 190,566, 195,562, 200,558, 205,554, 207,549, 208,544, 211,539, 216,535, 221,531, 226,527, 231,523, 236,519, 241,515, 246,512, 251,510, 256,508, 261,505, 266,501, 271,499, 276,497, 281,495, 286,493, 291,492, 296,491, 301,489, 306,488, 311,486, 316,484, 321,482, 326,480, 331,478, 336,476, 341,474, 346,471, 351,469, 356,468, 361,466, 366,464, 371,461, 376,459, 381,456, 386,452, 391,448, 396,446, 401,443, 406,440, 411,438, 416,434, 421,431, 426,427, 431,423, 436,419, 441,415, 446,411, 449,406, 453,401, 458,397, 463,393, 468,389, 473,385, 478,381, 483,378, 488,375, 493,371, 498,368, 503,366, 508,364, 513,360, 518,356, 522,351, 524,346, 526,341, 528,336, 528,331, 529,326, 530,321, 531,316, 531,311, 531,306, 531,301, 531,296, 531,291, 530,286, 529,281, 528,276, 528,271, 528,266, 527,261, 526,256, 526,251, 525,246, 524,241, 523,236, 521,231, 518,226, 515,221, 511,216, 506,214, 501,212, 496,209, 491,207, 486,206, 481,203, 476,200, 471,197, 466,195, 461,193, 456,190, 451,189, 446,186, 441,183, 436,179, 431,175, 426,171, 421,167, 417,162, 412,158, 407,154, 402,151, 397,148, 392,145, 387,141, 382,138, 377,134, 372,131, 367,128, 362,126, 357,124, 352,123, 347,121, 342,119, 337,118, 332,117, 327,116, 322,116, 317,115, 312,114, 307,112, 302,110, 297,109, 292,106, 287,105, 282,102, 277,101, 272,99, 267,98, 262,96, 257,95, 252,95, 247,95, 242,95, 237,95, 232,95, 227,95, 222,96, 217,99, 213,104, 211,109, 211,114, 211,119, 211,124, 211,129, 212,134, 213,139, 213,144, 214,149, 214,154, 216,159, 218,164, 222,169, 224,174, 227,179, 231,184, 232,189, 237,193, 239,198, 243,203, 246,208, 251,208, 256,211, 259,216, 264,219, 268,224, 272,229, 276,234, 279,239, 282,244, 283,249, 283,254, 283,259, 283,264, 282,269, 281,274, 279,279, 277,284, 276,289, 275,294, 273,299, 271,304, 269,309, 267,314, 263,319, 260,324, 256,329, 252,334, 249,339, 247,344, 243,349, 239,354, 234,356, 230,361, 225,362, 220,364, 215,366, 210,367, 205,370, 200,371, 195,373, 190,374, 185,375, 180,375, 175,376, 170,376, 165,375, 160,374, 155,372, 150,369, 145,366, 140,363, 135,360, 130,359, 125,356, 120,352, 115,348, 110,344, 106,339, 101,335, 97,330, 95,325, 94,320, 94,315, 95,310, 97,305, 99,300, 101,295, 103,290, 106,285, 109,280, 114,277, 119,273, 121,268, 126,264, 127,259, 130,254, 132,249, 134,244, 136,239, 138,234, 140,229, 141,224, 143,219, 144,214, 146,209, 148,204, 150,199, 153,194, 154,189, 155,184, 156,179, 157,174, 157,169, 157,164, 155,159, 153,154, 151,149, 149,144, 145,139, 141,134, 137,129, 133,124, 129,119, 124,116, 119,114, 114,114, 109,114, 104,114, 99,116, 94,119, 89,122, 85,127, 80,129, 76,134, 73,139, 69,144, 65,149, 62,154, 59,159, 56,164, 53,169, 50,174, 47,179, 44,184, 41,189, 39,194, 37,199, 34,204, 32,209, 30,214, 29,219, 29,224, 29,229, 29,234, 28,239, 28,244, 26,249, 25,254, 24,259, 23,264, 22,269, 21,274, 21,279, 21,284, 21,289, 21,294, 22,299, 22,304, 22,309, 23,314, 23,319, 24,324, 25,329, 25,334, 25,339, 25,344, 25,349, 24,354, 24,359, 23,364, 23,369, 21,374, 18,379, 17,384, 16,389, 16,394, 16,399, 16,404, 17,406
EndDataSection

Global x_=GetSystemMetrics_(#SM_CXFULLSCREEN)
Global y_=GetSystemMetrics_(#SM_CYFULLSCREEN)

Global xkm.f,ykm.f,xdy.f
Global xxx.f,rul,rux.s

Structure PointType
	x.f
	y.f
EndStructure

Global Dim P.PointType(10000)
Global Dim Q.PointType(10000)

Global np,nq
Global epsilon.f

Global dotinfo
Global geopos.s

Procedure LinePlot(x1,y1,x2,y2,color1,color2,size=1)

	LineXY(x1,y1,x2,y2,color1)
	If x1>=0 And y1>=0 And x1<x_ And y1<y_
		Box(x1-size>>1,y1-size>>1,size,size,color2)
	EndIf
	If x2>=0 And y2>=0 And x2<x_ And y2<y_
		Box(x2-size>>1,y2-size>>1,size,size,color2)
	EndIf

EndProcedure
Procedure.d Distance(px.d,py.d, ax.d,ay.d, bx.d,by.d)

	Protected x0.d, y0.d, l.d, t.d

	x0 = bx-ax
	y0 = ay-by
	l = x0*x0 + y0*y0

	If l=0
		x0 = bx-px
		y0 = by-py
		ProcedureReturn Sqr(x0*x0+y0*y0)
	EndIf

	t = ((bx-px)*y0+(by-py)*x0)/l
	If t<0.0
		t = -t
	EndIf

	ProcedureReturn t*Sqr(l)

EndProcedure
Procedure ShapeRead()

	Protected xmax,xmin
	Protected ymax,ymin
	Protected xmid,ymid

	xmin=99999
	ymin=99999
	xmax=-99999
	ymax=-99999

	np=0
	Restore ShapeData

	Read.w np
	For i=1 To np
		With P(i)
			Read.w \x
			Read.w \y

			If \x>xmax
				xmax=\x
			EndIf
			If \x<xmin
				xmin=\x
			EndIf
			If \y>ymax
				ymax=\y
			EndIf
			If \y<ymin
				ymin=\y
			EndIf

		EndWith

	Next i

	xmid=(xmin+xmax)/2
	ymid=(ymin+ymax)/2

	#GlobePrecision=			2<<28
	#GlobeAdaption=#GlobePrecision>>15

	ykm=111142.56676698782
	xkm=Cos(ymid*180.0*#GlobeAdaption/#GlobePrecision*#PI/180)*ykm
	xdy=xkm/ykm
	xxx=(xmax-xmin)*xkm*180.0*#GlobeAdaption/#GlobePrecision/x_

EndProcedure
Procedure ShapeMaxima()

	Global minx.f=8e8
	Global miny.f=8e8
	Global maxx.f=-8e8
	Global maxy.f=-8e8

	For i=1 To np
		If p(i)\x<minx
			minx=p(i)\x
		EndIf
		If p(i)\x>maxx
			maxx=p(i)\x
		EndIf
		If p(i)\y<miny
			miny=p(i)\y
		EndIf
		If p(i)\y>maxy
			maxy=p(i)\y
		EndIf
	Next i

	#Border=32

	maxx-minx
	maxy-miny
	maxx/(x_-#Border)*xdy
	maxy/(y_-#Border)
	maxx=1/maxx
	maxy=1/maxy

	If maxx>maxy
		maxx=maxy
	Else
		maxy=maxx
	EndIf

	minx=minx*maxx-#Border>>1
	miny=miny*maxy-#Border>>1

EndProcedure
Procedure ShapeDraw(Array o.PointType(1),n,color,dots=0,ox=0,oy=0)

	Protected x,y,v
	Protected xt,yt,tt,t.s

	o(n+1)\x=o(1)\x
	o(n+1)\y=o(1)\y

	v=dots>>1


	For i=1 To n-1
		x=o(i)\x*maxx-minx+ox
		y=o(i)\y*maxy-miny+oy
		x*xdy
		LinePlot(x,y,(o(i+1)\x*maxx-minx+ox)*xdy,o(i+1)\y*maxy-miny+oy,color,#Blue)
		If dots
			Box(x-v,y-v,dots,dots,#Black)
		EndIf

		If i=1
			DrawingMode(#PB_2DDrawing_Outlined)
			Box(x-6,y-6,12,12,#Blue)
			DrawingMode(#PB_2DDrawing_Transparent|#PB_2DDrawing_Outlined)
		EndIf
	Next i

EndProcedure
Procedure ShapeAddPoint(n)

	nq+1
	q(nq)\x = p(n)\x
	q(nq)\y = p(n)\y

EndProcedure
Procedure ShapeOptimize(first,last)

	Protected i.i,b.i
	Protected di.f,db.f,x.f,y.f

	If last > first + 1
		x = p(last)\x - p(first)\x
		y = p(last)\y - p(first)\y

		b = first+1
		db = Distance(p(b)\x,p(b)\y,p(first)\x,p(first)\y,p(last)\x,p(last)\y)
		i = b + 1
		While i < last
			di = Distance(p(i)\x,p(i)\y,p(first)\x,p(first)\y,p(last)\x,p(last)\y)
			If di >db
				b = i
				db = di
			EndIf
			i = i + 1
		Wend

		If db >= epsilon
			ShapeOptimize(first,b)
			ShapeAddPoint(b)
			ShapeOptimize(b,last)
		EndIf
	EndIf

EndProcedure
Procedure ShapeSimplify(tolerance.f)

	epsilon=tolerance
	nq=0

	ShapeAddPoint(1)
	ShapeOptimize(1,np)
	ShapeAddPoint(np)

EndProcedure
Procedure ShapeDoAll(tolerance.f)

	Protected s.s
	Protected n

	LoadFont(0,"Arial",7)
	StartDrawing(CanvasOutput(0))
	DrawingFont(FontID(0))
	DrawingMode(#PB_2DDrawing_Transparent)

	Box(0,0,x_,y_,#White)
	DrawingMode(#PB_2DDrawing_Transparent|#PB_2DDrawing_Outlined)
	ShapeDraw(p(),np,#Gray,0)
	DrawText(0,0,Str(np)+" Dots",#Gray)

	ShapeSimplify(tolerance)
	ShapeDraw(q(),nq,#Red,2,0)

	n=nq
	If nq>1 And q(1)\x=q(n)\x And q(1)\y=q(n)\y
		n-1
		s="*"
	EndIf

	DrawText(65,0,Str(nq)+s+" Dots, Tolerance = "+Str(tolerance),#Red)
	DrawText(200,0,"Tolerance: [+]  N   B  [-]",$404080)
	StopDrawing()

EndProcedure

Procedure Main(tolerance=5)

	ShapeRead()
	ShapeMaxima()

	OpenWindow(0,0,0,x_,y_,"World",#PB_Window_BorderLess)
	CanvasGadget(0,0,0,x_,y_)

	ShapeDoAll(tolerance)

	Repeat
		If WaitWindowEvent()=#WM_CHAR
			key=EventwParam()
			Select key

			Case #ESC
				End

			Case 'n','N'
				tolerance+1
				ShapeDoAll(tolerance)

			Case 'b','B'
				tolerance-1
				ShapeDoAll(tolerance)


			EndSelect
		EndIf

	ForEver

EndProcedure

Main()
Seymour Clufley
Addict
Addict
Posts: 1233
Joined: Wed Feb 28, 2007 9:13 am
Location: London

Re: Identifying the "bendy points" on a path

Post by Seymour Clufley »

Michael, thank you so much! That's very helpful indeed.

I have adapted your code and integrated it into my own, for anybody who is interested:

Code: Select all

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
EndMacro
Macro BeatThis(a,b)
  If b>a
      a=b
  EndIf
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


Macro CopyPoint(p1,p2)
  p2\x = p1\x
  p2\y = p1\y
EndMacro

Structure PointD
  x.d
  y.d
EndStructure

Structure PointI
  x.i
  y.i
EndStructure



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.d VertexDivergence(*p.PointI, *a.PointI, *b.PointI)
  
  Protected x0.d, y0.d, l.d, t.d
  
  x0 = *b\x-*a\x
  y0 = *a\y-*b\y
  l = x0*x0 + y0*y0
  
  If l=0
    x0 = *b\x-*p\x
    y0 = *b\y-*p\y
    ProcedureReturn Sqr(x0*x0+y0*y0)
  EndIf
  
  t = ((*b\x-*p\x)*y0+(*b\y-*p\y)*x0)/l
  If t<0.0
    t = -t
  EndIf
  
  ProcedureReturn t*Sqr(l)
  
EndProcedure


Structure SimplifyStructure
  simple_level.d
  spnts.i
  Array spnt.PointI(1)
EndStructure


Macro ShapeAddPoint(simplify,orig)
  simplify\spnts+1
  CopyPoint(orig,simplify\spnt(simplify\spnts))
EndMacro


Procedure ShapeOptimize(pnts.i,Array pnt.PointI(1),*simplify.SimplifyStructure,first.i,last.i)
  
  Protected x.d,y.d
  
  If last > first + 1
    x = pnt(last)\x - pnt(first)\x
    y = pnt(last)\y - pnt(first)\y
    
    b = first+1
    db.d = VertexDivergence(@pnt(b),@pnt(first),@pnt(last))
    i = b + 1
    While i < last
      di.d = VertexDivergence(@pnt(i),@pnt(first),@pnt(last))
      If di >db
        b = i
        db = di
      EndIf
      i = i + 1
    Wend
    
    If db >= *simplify\simple_level
      ShapeOptimize(pnts,pnt(),*simplify,first,b)
      ShapeAddPoint(*simplify,pnt(b))
      ShapeOptimize(pnts,pnt(),*simplify,b,last)
    EndIf
  EndIf
  
EndProcedure


Procedure.b GetPathBends(pnts.i,Array pnt.PointI(1),*bnd.SimplifyStructure)
  With *bnd
    \spnts = 0
    Dim \spnt(pnts)
    ShapeAddPoint(*bnd,pnt(1))
    ShapeOptimize(pnts,pnt(),*bnd,1,pnts)
    ShapeAddPoint(*bnd,pnt(pnts))
  EndWith
EndProcedure




Procedure DrawPolygon(Array o.PointI(1),pnts.i,clr.i,draw_lines.b=#True,draw_points.b=#False)
  
  For i = 1 To pnts-1
    If draw_lines
      LineXY(o(i)\x,o(i)\y,o(i+1)\x,o(i+1)\y,clr)
    EndIf
    If draw_points
      Circle(o(i)\x,o(i)\y,3,clr)
    EndIf
  Next i
  If draw_points
    Circle(o(pnts)\x,o(pnts)\y,3,clr)
  EndIf
  
EndProcedure


Procedure ShapeDoAll(cnv.i,pnts.i,Array pnt.PointI(1),level.d)
  
  LoadFont(0,"Arial",7)
  StartDrawing(CanvasOutput(cnv))
  DrawingFont(FontID(0))
  Box(0,0,OutputWidth(),OutputHeight(),RGB(255,255,192))
  DrawingMode(#PB_2DDrawing_Transparent|#PB_2DDrawing_Outlined)
  DrawText(0,0,Str(pnts)+" points",#Gray)
  DrawingMode(#PB_2DDrawing_Transparent)
  
  simplify.SimplifyStructure
  
  With simplify
    \simple_level = level
    GetPathBends(pnts,pnt(),@simplify)
    
    DrawPolygon(pnt(),pnts,#Black)
    Circle(pnt(1)\x,pnt(1)\y,5,#Green)
    Circle(pnt(pnts)\x,pnt(pnts)\y,5,#Red)
    DrawPolygon(\spnt(),\spnts,#Gray,#False,#True)
    
    DrawText(60,0,"Simplification: [-] B "+Str(\simple_level)+"  N [+]",$404080)
    DrawText(200,0,"Simple points: "+Str(\spnts),#Red)
    StopDrawing()
  EndWith
  
  FreeFont(0)
  
EndProcedure





pathpnts.i = 0
Dim pathpnt.PointI(10000)



simplify_level.d = 10

iw.i = 1400;GetSystemMetrics_(#SM_CXFULLSCREEN)
ih.i = 700;GetSystemMetrics_(#SM_CYFULLSCREEN)
win = OpenWindow(#PB_Any,0,0,iw,ih,"World",#PB_Window_ScreenCentered|#PB_Window_BorderLess)
cnv = CanvasGadget(#PB_Any,0,0,iw,ih)
StartDrawing(CanvasOutput(cnv))
Box(0,0,OutputWidth(),OutputHeight(),RGB(255,255,192))
StopDrawing()

Repeat
  we = WaitWindowEvent()
  Select we
    Case #WM_CHAR
      Select EventwParam()
        Case #ESC
          End
        Case 'n','N'
          simplify_level+10
          ShapeDoAll(cnv,pathpnts,pathpnt(),simplify_level)
        Case 'b','B'
          simplify_level-10
          ShapeDoAll(cnv,pathpnts,pathpnt(),simplify_level)
      EndSelect
      
    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()
              pathpnts = 0
              
            Case #PB_EventType_LeftButtonUp
              tracking = #False
              pathpnts = JoinUpDots(pathpnts,pathpnt())
              ShapeDoAll(cnv,pathpnts,pathpnt(),simplify_level)
              
            Case #PB_EventType_MouseMove
              If tracking
                x = GetGadgetAttribute(cnv,#PB_Canvas_MouseX)
                y = GetGadgetAttribute(cnv,#PB_Canvas_MouseY)
                If pathpnts=0 Or x<>pathpnt(pathpnts)\x Or y<>pathpnt(pathpnts)\y
                  ;Debug Str(x)+","+Str(y)
                  pathpnts+1
                  pathpnt(pathpnts)\x = x
                  pathpnt(pathpnts)\y = y
                  
                  StartDrawing(CanvasOutput(cnv))
                  Circle(pathpnt(pathpnts)\x,pathpnt(pathpnts)\y,3,#Black)
                  StopDrawing()
                EndIf
              EndIf
          EndSelect
      EndSelect
      
  EndSelect
  Delay(10)
ForEver
JACK WEBB: "Coding in C is like sculpting a statue using only sandpaper. You can do it, but the result wouldn't be any better. So why bother? Just use the right tools and get the job done."
User avatar
Andre
PureBasic Team
PureBasic Team
Posts: 2058
Joined: Fri Apr 25, 2003 6:14 pm
Location: Germany (Saxony, Deutscheinsiedel)
Contact:

Re: Identifying the "bendy points" on a path

Post by Andre »

Very nice, thank you! :D
Bye,
...André
(PureBasicTeam::Docs & Support - PureArea.net | Order:: PureBasic | PureVisionXP)
Post Reply