Split a concave polygon into convex polygons
Posted: Thu Aug 27, 2020 8:41 am
A fairly niche job, I admit, but you never know when somebody may need to do this exact thing.
The procedure can be set to use the first available split, but by default, it looks for the most "equal" split, so that the two resulting convex polygons will be as similar in area to each other as possible. This takes more time, obviously.
Because my code sometimes "over-splits", I have written a procedure for merging contiguous polygons (as long as the resulting polygon will be convex). This procedure is useful for other purposes.
Also included are procedures for calculating the area of a polygon, and identifying its type (concave, convex, triangle, etc.).
Demo code (press SPACE to cycle through and ESCAPE to quit):
The procedure can be set to use the first available split, but by default, it looks for the most "equal" split, so that the two resulting convex polygons will be as similar in area to each other as possible. This takes more time, obviously.
Because my code sometimes "over-splits", I have written a procedure for merging contiguous polygons (as long as the resulting polygon will be convex). This procedure is useful for other purposes.
Also included are procedures for calculating the area of a polygon, and identifying its type (concave, convex, triangle, etc.).
Code: Select all
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
Macro CopyPoint(cpp1,cpp2)
cpp2\x=cpp1\x
cpp2\y=cpp1\y
EndMacro
Structure LineD
pnt1.PointD
pnt2.PointD
deg.d
length.d
EndStructure
Macro CopyLine(cpl1,cpl2)
CopyPoint(cpl1\pnt1,cpl2\pnt1)
CopyPoint(cpl1\pnt2,cpl2\pnt2)
EndMacro
Structure PolygonStructure
points.i
Array pnt.PointD(3)
EndStructure
Macro PointsIdentical(p1,p2)
Bool(p1\x=p2\x And p1\y=p2\y)
EndMacro
Macro PointsVirtuallyIdentical(p1,p2)
Bool( Difference(p1\x,p2\x)<1 And Difference(p1\y,p2\y)<1 )
EndMacro
Macro LinesVirtuallyIdentical(l1,l2)
((PointsVirtuallyIdentical(l1\pnt1,l2\pnt1) And PointsVirtuallyIdentical(l1\pnt2,l2\pnt2) ) Or ( PointsVirtuallyIdentical(l1\pnt1,l2\pnt2) And PointsVirtuallyIdentical(l1\pnt2,l2\pnt1) ) )
EndMacro
Macro LinesIdentical(l1,l2)
((PointsIdentical(l1\pnt1,l2\pnt1) And PointsIdentical(l1\pnt2,l2\pnt2) ) Or ( PointsIdentical(l1\pnt1,l2\pnt2) And PointsIdentical(l1\pnt2,l2\pnt1) ) )
EndMacro
;-
;- GEOMETRY STUFF
Procedure.d DoubleCrossProductD(*a.PointD,*b.PointD,*c.PointD)
dx1.d = *b\x-*a\x
dy1.d = *b\y-*a\y
dx2.d = *c\x-*b\x
dy2.d = *c\y-*b\y
ProcedureReturn (dx1*dy2) - (dy1*dx2)
EndProcedure
Procedure.d Angle2D(*p1.PointD,*p2.PointD)
theta1.d = ATan2(*p1\y,*p1\x)
theta2.d = ATan2(*p2\y,*p2\x)
dtheta.d = theta2 - theta1
TWOPI.d = #PI*2
While dtheta > #PI
dtheta - TWOPI
Wend
While dtheta < -#PI
dtheta + TWOPI
Wend
ProcedureReturn dtheta
EndProcedure
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
Procedure.d DistanceBetweenTwoPoints(*a.PointD,*b.PointD)
xdif.d = Difference(*a\x,*b\x)
ydif.d = Difference(*a\y,*b\y)
ProcedureReturn Sqr((xdif*xdif)+(ydif*ydif))
EndProcedure
Procedure.b RadianCoordsFromPoint(*base.PointD,radia.d,distance.d,*c.PointD)
*c\x = (Cos(radia)*distance)+*base\x
*c\y = (Sin(radia)*distance)+*base\y
EndProcedure
Procedure.b DegreeCoordsFromPoint(*base.PointD,degrees.d,distance.d,*c.PointD)
RadianCoordsFromPoint(*base,Radian(degrees),distance,*c)
EndProcedure
Procedure.b MidPoint(*pnt1.PointD,*pnt2.PointD,*mp.PointD)
deg.d = DegreeAngleBetweenTwoPoints(*pnt1,*pnt2)
dist.d = DistanceBetweenTwoPoints(*pnt1,*pnt2)
DegreeCoordsFromPoint(*pnt1,deg,dist/2,*mp)
EndProcedure
;-
;- POLYGON STUF
Enumeration 0 ; polygon types
#Polygon_NonTriangulable
#Polygon_Convex
#Polygon_Concave
#Polygon_Triangle
EndEnumeration
Procedure.s PolygonTypeName(tp.b)
Select tp
Case #Polygon_NonTriangulable
ProcedureReturn "non-triangulable"
Case #Polygon_Convex
ProcedureReturn "convex"
Case #Polygon_Concave
ProcedureReturn "concave"
Case #Polygon_Triangle
ProcedureReturn "triangle"
EndSelect
ProcedureReturn "unknown!"
EndProcedure
Procedure.b IdentifyPolygonType(points.i,Array pnt.PointD(1))
;R("IPT points: "+Str(points))
Select points
Case 0, 1, 2
ProcedureReturn #Polygon_NonTriangulable
Case 3
ProcedureReturn #Polygon_Triangle
EndSelect
hasneg.b
haspos.b
zcrossproduct.d = DoubleCrossProductD(@pnt(1),@pnt(2),@pnt(3))
If zcrossproduct>0
haspos=#True
Else
If zcrossproduct<0
hasneg=#True
EndIf
EndIf
For p = 2 To points
a = p
b = p+1
c = p+2
If p=points-1
c=1
Else
If p=points
b=1
c=2
EndIf
EndIf
zcrossproduct.d = DoubleCrossProductD(@pnt(a),@pnt(b),@pnt(c))
If zcrossproduct>0
haspos=#True
Else
If zcrossproduct<0
hasneg=#True
EndIf
EndIf
Next p
;R("HasPos: "+ByteTruth(haspos)+c13+"HasNeg: "+ByteTruth(hasneg))
If haspos And hasneg
ProcedureReturn #Polygon_Concave
Else
ProcedureReturn #Polygon_Convex
EndIf
EndProcedure
Procedure.b EnsurePolygonIsPolyline(*pg.PolygonStructure)
If PointsIdentical(*pg\pnt(1),*pg\pnt(*pg\points))
;R("IS NOT POLYLINE")
*pg\points-1
ReDim *pg\pnt(*pg\points)
EndIf
EndProcedure
Procedure.i PointDIsInside(pnts.i,Array pnt.PointD(1),*p.PointD)
angle.d = 0
p1.PointD
p2.PointD
For i = 1 To pnts
tp1 = i
If i<pnts
tp2 = i+1
Else
tp2 = 1
EndIf
p1\x = pnt(tp1)\x - *p\x
p1\y = pnt(tp1)\y - *p\y
p2\x = pnt(tp2)\x - *p\x
p2\y = pnt(tp2)\y - *p\y
angle + Angle2D(@p1,@p2)
Next
If Abs(angle) < #PI
ProcedureReturn #False
Else
ProcedureReturn #True
EndIf
EndProcedure
Procedure.i PointDIsInsidePolygon(*g.PolygonStructure,*p.PointD)
ProcedureReturn PointDIsInside(*g\points,*g\pnt(),*p)
EndProcedure
Procedure.b PolygonFromWithinPolygon(*opg.PolygonStructure,pa.i,pb.i,*npg.PolygonStructure,ltr.s)
*npg\points = 0
Dim *npg\pnt(*opg\points*2)
narr.s = ""
If pb>pa
For a = pa To pb
*npg\points+1
narr + Str(a)+c32
CopyPoint(*opg\pnt(a),*npg\pnt(*npg\points))
Next a
Else
For a = pa To *opg\points
*npg\points+1
narr + Str(a)+c32
CopyPoint(*opg\pnt(a),*npg\pnt(*npg\points))
Next a
For a = 1 To pb
If a=1
If PointsIdentical(*opg\pnt(1),*opg\pnt(*opg\points))
Continue
EndIf
EndIf
*npg\points+1
narr + Str(a)+c32
CopyPoint(*opg\pnt(a),*npg\pnt(*npg\points))
Next a
EndIf
ReDim *npg\pnt(*npg\points)
EndProcedure
Procedure.d ComputePolygonArea(pnts.i,Array pnt.PointD(1))
area.d = 0
For p = 1 To pnts
p2=p+1 : If p=pnts : p2=1 : EndIf
area + (pnt(p)\x*pnt(p2)\y - pnt(p)\y*pnt(p2)\x)
Next
ProcedureReturn Abs(area/2)
EndProcedure
;-
Procedure.i SplitConcavePolygonIntoConvexPolygons(*pg.PolygonStructure,Array cvx.PolygonStructure(1),use_first_match.b)
If IdentifyPolygonType(*pg\points,*pg\pnt())<>#Polygon_Concave
ProcedureReturn #False
EndIf
CopyStructure(*pg,unencpg.PolygonStructure,PolygonStructure)
EnsurePolygonIsPolyline(@unencpg)
smallest_dif.d = 9999999
Dim cvx(2)
With unencpg
; sides, for intersection testing
;Dim side.LineD(*pg\points)
For p1 = 1 To \points
; test each other point
For p2 = 1 To \points-1
If p1=p2 Or Difference(p1,p2)<2 : Continue : EndIf
Midpoint(@\pnt(p1),@\pnt(p2),@mp.PointD)
If Not PointDIsInsidePolygon(@unencpg,@mp)
Continue
EndIf
; test the two polygons created by this hypothetical split...
PolygonFromWithinPolygon(@unencpg,p1,p2,@test1.PolygonStructure,"A")
type1.b = IdentifyPolygonType(test1\points,test1\pnt())
PolygonFromWithinPolygon(@unencpg,p2,p1,@test2.PolygonStructure,"B")
type2.b = IdentifyPolygonType(test2\points,test2\pnt())
;Debug "A: "+PolygonTypeName(type1)
;Debug "B: "+PolygonTypeName(type2)
If type1=#Polygon_Convex Or type1=#Polygon_Triangle
If type2=#Polygon_Convex Or type2=#Polygon_Triangle
;R("SUCCESS!")
If use_first_match
CopyStructure(@test1,@cvx(1),PolygonStructure)
CopyStructure(@test2,@cvx(2),PolygonStructure)
ProcedureReturn 2
EndIf
area1.d = ComputePolygonArea(test1\points,test1\pnt())
area2.d = ComputePolygonArea(test2\points,test2\pnt())
dif.d = Difference(area1,area2)
If dif<smallest_dif
smallest_dif = dif
CopyStructure(@test1,@cvx(1),PolygonStructure)
CopyStructure(@test2,@cvx(2),PolygonStructure)
EndIf
EndIf
EndIf
Next p2
Next p1
If smallest_dif<>9999999
ProcedureReturn 2
EndIf
; it hasn't been able to split it into two convex polygons. recursion is necessary
For p1 = 1 To \points
; test each other point
For p2 = 1 To \points-1
If p1=p2 Or Difference(p1,p2)<2 : Continue : EndIf
Midpoint(@\pnt(p1),@\pnt(p2),@mp.PointD)
If Not PointDIsInsidePolygon(@unencpg,@mp)
Continue
EndIf
; test the two polygons created by this hypothetical split...
PolygonFromWithinPolygon(@unencpg,p1,p2,@test1.PolygonStructure,"A")
type1.b = IdentifyPolygonType(test1\points,test1\pnt())
If type1<>#Polygon_Convex And type1<>#Polygon_Triangle
Continue
EndIf
CopyStructure(@test1,@cvx(1),PolygonStructure)
PolygonFromWithinPolygon(@unencpg,p2,p1,@ccave.PolygonStructure,"B")
; this second polygon will be concave
Dim ncvx.PolygonStructure(ccave\points)
ncvxs = SplitConcavePolygonIntoConvexPolygons(@ccave,ncvx(),use_first_match)
;R("NCVXS: "+Str(ncvxs))
ReDim cvx(1+ncvxs)
c=1
For nc = 1 To ncvxs
c+1
CopyStructure(@ncvx(nc),@cvx(c),PolygonStructure)
Next nc
ProcedureReturn 1+ncvxs
Next p2
Next p1
ProcedureReturn 0
EndWith
EndProcedure
#d1="|"
#d2="~"
Procedure.b IdentifyPolygonNeighbours(Map quad.PolygonStructure(),Map neighbour.s())
ClearMap(neighbour())
NewMap tempquad.PolygonStructure()
CopyMap(quad(),tempquad())
ForEach quad()
qkey.s = MapKey(quad())
Dim qline.LineD(quad(qkey)\points)
For qp = 1 To quad(qkey)\points
qp2=qp+1 : If qp=quad(qkey)\points : qp2=1 : EndIf
CopyPoint(quad(qkey)\pnt(qp),qline(qp)\pnt1) : CopyPoint(quad(qkey)\pnt(qp2),qline(qp)\pnt2)
Next qp
ForEach tempquad()
tkey.s = MapKey(tempquad())
If tkey=qkey : Continue : EndIf
For tp = 1 To tempquad(tkey)\points
tp2=tp+1 : If tp=tempquad(tkey)\points : tp2=1 : EndIf
tline.LineD : CopyPoint(tempquad(tkey)\pnt(tp),tline\pnt1) : CopyPoint(tempquad(tkey)\pnt(tp2),tline\pnt2)
For qp = 1 To quad(qkey)\points
If LinesVirtuallyIdentical(tline,qline(qp))
;Debug "MATCH"
neighbour(qkey) + Str(qp)+#d2+tkey+#d2+Str(tp)+#d2+#d1
EndIf
Next qp
Next tp
Next
Next
FreeMap(tempquad())
ProcedureReturn #True
EndProcedure
Procedure.b MergeTheseContiguousPolygons(Map pg.PolygonStructure(),*mgp.PolygonStructure)
Dim templinebucket.LineD(100)
lns = 0
ForEach pg()
key.s = MapKey(pg())
For s = 1 To pg()\points
lns+1
If lns>ArraySize(templinebucket()) : ReDim templinebucket(lns+100) : EndIf
s2=s+1 : If s=pg()\points : s2=1 : EndIf
CopyPoint(pg()\pnt(s),templinebucket(lns)\pnt1)
CopyPoint(pg()\pnt(s2),templinebucket(lns)\pnt2)
Next s
Next
ReDim templinebucket(lns)
;Debug "LNS: "+Str(lns)
NewList linebucket.LineD()
For a = 1 To lns
duplicate_found.b = #False
For b = 1 To lns
If b=a : Continue : EndIf
If LinesIdentical(templinebucket(a),templinebucket(b))
;Debug "dupe: "+Str(a)+" and "+Str(b)
duplicate_found = #True
Break
EndIf
Next b
If Not duplicate_found
AddElement(linebucket())
CopyLine(templinebucket(a),linebucket())
EndIf
Next a
FreeArray(templinebucket())
*mgp\points = ListSize(linebucket())*2
Dim *mgp\pnt(*mgp\points)
FirstElement(linebucket())
p=1
CopyPoint(linebucket()\pnt1,*mgp\pnt(p))
p=2
CopyPoint(linebucket()\pnt2,*mgp\pnt(p))
DeleteElement(linebucket())
Repeat
ForEach linebucket()
If PointsIdentical(*mgp\pnt(p),linebucket()\pnt1)
p+1
CopyPoint(linebucket()\pnt2,*mgp\pnt(p))
DeleteElement(linebucket())
Break
Else
If PointsIdentical(*mgp\pnt(p),linebucket()\pnt2)
p+1
CopyPoint(linebucket()\pnt1,*mgp\pnt(p))
DeleteElement(linebucket())
Break
EndIf
EndIf
Next
Until p=*mgp\points Or ListSize(linebucket())=0
*mgp\points = p
If PointsIdentical(*mgp\pnt(1),*mgp\pnt(*mgp\points))
*mgp\points-1
EndIf
ProcedureReturn #True
EndProcedure
Procedure.i MergeContiguousPolygonsIntoConvexPolygons(pgs.i,Array pg.PolygonStructure(1),Array cvx.PolygonStructure(1))
NewMap pgmap.PolygonStructure()
NewMap neighbour.s()
NewMap mergemap.PolygonStructure()
NewMap tested_combination.b()
For p = 1 To pgs
key.s = Str(p)
CopyStructure(@pg(p),@pgmap(key),PolygonStructure)
Next p
count = pgs
Repeat
mergings.b = #False
IdentifyPolygonNeighbours(pgmap(),neighbour())
keyarr.s = ""
ForEach pgmap()
keyarr + MapKey(pgmap())+#d1
Next
ks = CountString(keyarr,#d1)
For k = 1 To ks
key.s = StringField(keyarr,k,#d1)
For n = 1 To CountString(neighbour(key),#d1)
this_one.s = StringField(neighbour(key),n,#d1)
nkey.s = StringField(this_one,2,#d2)
combkey.s = key+#d1+nkey
If FindMapElement(tested_combination(),combkey)
Continue
EndIf
tested_combination(combkey) = #True
ClearMap(mergemap())
CopyStructure(@pgmap(key),@mergemap(key),PolygonStructure)
CopyStructure(@pgmap(nkey),@mergemap(nkey),PolygonStructure)
MergeTheseContiguousPolygons(mergemap(),@mg.PolygonStructure)
mgtype.b = IdentifyPolygonType(mg\points,mg\pnt())
If mgtype=#Polygon_Convex Or mgtype=#Polygon_Triangle
mergings = #True
DeleteMapElement(pgmap(),nkey)
DeleteMapElement(pgmap(),key)
count+1
CopyStructure(@mg,@pgmap(Str(count)),PolygonStructure)
Break 2
EndIf
Next n
Next
Until Not mergings
ReDim cvx(MapSize(pgmap()))
p=0
ForEach pgmap()
p+1
CopyStructure(@pgmap(),@cvx(p),PolygonStructure)
Next
ProcedureReturn p
EndProcedure
Code: Select all
;-
;- DEMO CODE
Procedure.b FindViableFillpointInsidePolygon(pnts.i,Array pnt.PointD(1),*fp.PointD,xlimit.i,ylimit.i)
*fp\x=-1 : *fp\y=-1
Select pnts
Case 1, 2
ProcedureReturn #False
Case 3
MidPoint(@pnt(1),@pnt(2),@test.PointD)
MidPoint(@test,@pnt(3),*fp)
Default
For p = 3 To pnts
MidPoint(@pnt(1),@pnt(p),@test.PointD)
If test\x>0 And test\x<xlimit And test\y>0 And test\y<ylimit
If PointDIsInside(pnts,pnt(),@test)=#True
CopyPoint(test,*fp)
ProcedureReturn #True
EndIf
EndIf
Next p
EndSelect
ProcedureReturn #False
EndProcedure
Procedure.b DrawPolygon(*pg.PolygonStructure,clr.i,do_fill.b=#False,xoff.i=0,yoff.i=0)
With *pg
For p = 1 To \points
p2=p+1 : If p=\points : p2=1 : EndIf
LineXY(\pnt(p)\x+xoff,\pnt(p)\y+yoff,\pnt(p2)\x+xoff,\pnt(p2)\y+yoff,clr)
Next p
If do_fill
FindViableFillpointInsidePolygon(*pg\points,*pg\pnt(),@fp.PointD,OutputWidth(),OutputHeight())
;MidPoint(\pnt(1),\pnt(3),@fp.PointD)
FillArea(fp\x+xoff,fp\y+yoff,-1,clr)
;Circle(fp\x+xoff,fp\y+yoff,3,RGBA(0,0,0,255))
EndIf
EndWith
EndProcedure
;-
;- DEMO
iw = 1920
ih = 1000
cnt.PointD : cnt\x=iw/2 : cnt\y=ih/2
win = OpenWindow(#PB_Any,0,0,iw,ih,"",#PB_Window_ScreenCentered|#PB_Window_BorderLess)
space = 5
AddKeyboardShortcut(win,#PB_Shortcut_Space,space)
esc = 6
AddKeyboardShortcut(win,#PB_Shortcut_Escape,esc)
cnv = CanvasGadget(#PB_Any,0,0,iw,ih)
Repeat
; random polygon
Repeat
pg.PolygonStructure
pg\points = Random(15,5)
Dim pg\pnt(pg\points)
offset.d = Random(90)
For p = 1 To pg\points
dist = Random(ih/4,20)
switch = 1-switch
If switch
dist * 2
EndIf
DegreeCoordsFromPoint(@cnt,offset+(360/pg\points*p),dist,@pg\pnt(p))
Next p
Until IdentifyPolygonType(pg\points,pg\pnt())=#Polygon_Concave
;DrawPolygon(@pg,RGBA(255,0,0,255))
Dim cvx.PolygonStructure(2)
cvxs.i = SplitConcavePolygonIntoConvexPolygons(@pg,cvx(),#False)
;R("FINAL CVXS: "+Str(cvxs))
If cvxs>0
Dim ncvx.PolygonStructure(cvxs)
cvxs.i = MergeContiguousPolygonsIntoConvexPolygons(cvxs,cvx(),ncvx())
StartDrawing(CanvasOutput(cnv))
Box(0,0,OutputWidth(),OutputHeight(),RGBA(50,50,50,255))
For c = 1 To cvxs
DrawPolygon(@ncvx(c),RGBA(Random(255),Random(255),Random(255),255),#True)
Next c
StopDrawing()
Else
MessageRequester("ERROR","COULDN'T DO IT!")
EndIf
Repeat : Until WaitWindowEvent(5)=#PB_Event_Menu
If EventMenu()=esc : Break : EndIf
ForEver
