I have been working on polygon triangulation. The ideal would be a solution that was guaranteed to work for any polygon - Jonathan Shewchuk's
, for example - but in the absence of that I have attempted the challenge myself. I'm pretty amazed, but it seems to work. I haven't tried it with "real world" polygons yet, just randomly drawn ones, so maybe it will fail with more unpredictable shapes. But as I say it seems to work, so I'll put it up here in case anyone else can use it or improve it. It might save somebody some time in the future. The demo below will keep creating and triangulating polygons - press RETURN to get to the next one, press ESCAPE to abort. If you find any bugs, please let me know.
Code: Select all
;- "REPORTING" FUNCTIONS
Global c13.s = Chr(13)
Macro R(t)
MessageRequester("Report",t,0)
EndMacro
Macro StandardReportingWindowStuff(win)
escapekey = 1
spacekey = 2
returnkey = 3
AddKeyboardShortcut(win,#PB_Shortcut_Escape,escapekey)
AddKeyboardShortcut(win,#PB_Shortcut_Space,spacekey)
AddKeyboardShortcut(win,#PB_Shortcut_Return,returnkey)
Repeat
If WindowEvent()=#PB_Event_Menu : Break : EndIf
Delay(10)
ForEver
CloseWindow(win)
EndMacro
;- GENERAL FUNCTIONS
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.f Difference(a.f,b.f)
If a=b
ProcedureReturn 0
EndIf
If a>b
ProcedureReturn a-b
Else
ProcedureReturn b-a
EndIf
EndProcedure
;- GEOMETRIC FUNCTIONS
Structure PointF
x.f
y.f
EndStructure
Procedure.f DoubleCrossProduct(*a.PointF,*b.PointF,*c.PointF)
dx1.f = *b\x-*a\x
dy1.f = *b\y-*a\y
dx2.f = *c\x-*b\x
dy2.f = *c\y-*b\y
ProcedureReturn dx1*dy2 - dy1*dx2
EndProcedure
Procedure.f DistanceBetweenTwoPoints(*a.PointF,*b.PointF)
xdif.f = Difference(*a\x,*b\x)
ydif.f = Difference(*a\y,*b\y)
ProcedureReturn Sqr((xdif*xdif)+(ydif*ydif))
EndProcedure
Procedure.f DegreeAngleBetweenTwoPoints(*o.PointF,*b.PointF)
x.f = *b\x - *o\x
y.f = *b\y - *o\y
ProcedureReturn Degree(ATan(y/x))
EndProcedure
Procedure.b CoordsFromPoint(*base.PointF,radia.d,distance.d,*c.PointF)
*c\x = (Cos(radia)*distance)+*base\x
*c\y = (Sin(radia)*distance)+*base\y
EndProcedure
Macro DegreeCoordsFromPoint(base,degrees,distance,c)
CoordsFromPoint(base,Radian(degrees),distance,c)
EndMacro
Macro PointsIdentical(p1,p2)
(p1\x=p2\x And p1\y=p2\y)
EndMacro
Enumeration ; polygon types
#Polygon_Convex
#Polygon_Concave
#Polygon_Triangle
EndEnumeration
Procedure.b IdentifyPolygonType(points.i,Array pnt.PointF(1))
If points=3
ProcedureReturn #Polygon_Triangle
EndIf
hasneg.b
haspos.b
zcrossproduct.f = DoubleCrossProduct(@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.f = DoubleCrossProduct(@pnt(a),@pnt(b),@pnt(c))
If zcrossproduct>0
haspos=#True
Else
If zcrossproduct<0
hasneg=#True
EndIf
EndIf
Next p
If haspos And hasneg
ProcedureReturn #Polygon_Concave
Else
ProcedureReturn #Polygon_Convex
EndIf
EndProcedure
Procedure.b GetUsableCentrePoint(points.i,Array pnt.PointF(1),*c.PointF)
br.PointF
For p = 1 To points
BeatThis(br\x,pnt(p)\x)
BeatThis(br\y,pnt(p)\y)
Next p
CopyStructure(@br,@tl.PointF,PointF)
For p = 1 To points
DefeatThis(tl\x,pnt(p)\x)
DefeatThis(tl\y,pnt(p)\y)
Next p
dist.f = DistanceBetweenTwoPoints(@tl,@br)
deg.f = DegreeAngleBetweenTwoPoints(@tl,@br)
DegreeCoordsFromPoint(@tl,deg,dist/2,*c)
EndProcedure
;- TRIANGULATION
Structure TriangleStructure
a.PointF
b.PointF
c.PointF
EndStructure
Structure Triangulation
triangles.i
Array triangle.TriangleStructure(0)
EndStructure
Procedure.b TriangulateConvexPolygon(points.i,Array pnt.PointF(1),*t.Triangulation)
GetUsableCentrePoint(points,pnt(),@cnt.PointF) ; maybe getting the centroid would be better, but this seems to work
*t\triangles = points
Dim *t\triangle(points)
For p = 1 To points
If p=points
p2=1
Else
p2=p+1
EndIf
CopyStructure(@pnt(p),*t\triangle(p)\a,PointF)
CopyStructure(@cnt,*t\triangle(p)\b,PointF)
CopyStructure(@pnt(p2),*t\triangle(p)\c,PointF)
Next p
EndProcedure
Procedure.f TriangleArea(*t.TriangleStructure)
ProcedureReturn 0.5*(-*t\b\y**t\c\x + *t\a\y*(-*t\b\x + *t\c\x) + *t\a\x*(*t\b\y - *t\c\y) + *t\b\x**t\c\y)
EndProcedure
Procedure.b PointIsInsideTriangle(*p.PointF,*t.TriangleStructure)
Area.f = TriangleArea(*t)
If Area=0
R("AREA: "+StrF(Area)+c13+c13+StrF(*t\a\x)+","+StrF(*t\a\y)+c13+StrF(*t\b\x)+","+StrF(*t\b\y)+c13+StrF(*t\c\x)+","+StrF(*t\c\y))
EndIf
s.f = 1/(2*Area)*(*t\a\y**t\c\x - *t\a\x**t\c\y + (*t\c\y - *t\a\y)* *p\x + (*t\a\x - *t\c\x) * *p\y)
t.f = 1/(2*Area)*(*t\a\x**t\b\y - *t\a\y**t\b\x + (*t\a\y - *t\b\y)* *p\x + (*t\b\x - *t\a\x)* *p\y)
If s>0 And t>0 And (1-s-t)>0
ProcedureReturn #True
EndIf
EndProcedure
Procedure.b TriangleCentre(*t.TriangleStructure,*c.PointF)
*c\x = (*t\a\x + *t\b\x + *t\c\x) / 3
*c\y = (*t\a\y + *t\b\y + *t\c\y) / 3
EndProcedure
Macro DrawTriangle(t,clr)
LineXY(t\a\x,t\a\y,t\b\x,t\b\y,clr)
LineXY(t\b\x,t\b\y,t\c\x,t\c\y,clr)
LineXY(t\c\x,t\c\y,t\a\x,t\a\y,clr)
TriangleCentre(@t,@dtcnt.PointF)
FillArea(dtcnt\x,dtcnt\y,-1,clr)
EndMacro
Macro TriangleHasPoint(t,p)
( PointsIdentical(t\a,p) Or PointsIdentical(t\b,p) Or PointsIdentical(t\c,p) )
EndMacro
Structure UsePoint
p.i
valid.b
EndStructure
Procedure.i RemoveDeadPoints(points.i,Array nt.UsePoint(1))
Dim nt2.UsePoint(points)
newusepoints = 0
For p = 1 To points
If nt(p)\valid
newusepoints+1
CopyStructure(@nt(p),@nt2(newusepoints),UsePoint)
EndIf
Next p
CopyArray(nt2(),nt())
ProcedureReturn newusepoints
EndProcedure
Procedure.i UsePointSum(Array u.UsePoint(1))
usepoints = 0
For a = 1 To ArraySize(u())
If u(a)\valid = #True
usepoints+1
EndIf
Next
ProcedureReturn usepoints
EndProcedure
Procedure.b PointsVirtuallyIdentical(*p1.PointF,*p2.PointF)
If Difference( *p1\x,*p2\x)<1 And Difference(*p1\y,*p2\y)<1
ProcedureReturn #True
EndIf
ProcedureReturn #False
EndProcedure
Procedure.f AngleBetweenThreePoints(*c.PointF,*p1.Point,*p2.PointF)
tpnt.PointF
cdist1.f = DistanceBetweenTwoPoints(*c,*p1)
cdeg1.f = DegreeAngleBetweenTwoPoints(*c,*p1)
DegreeCoordsFromPoint(*c,cdeg1,cdist1,@tpnt)
While Not PointsVirtuallyIdentical(*p1,@tpnt)
cdeg1+180
DegreeCoordsFromPoint(*c,cdeg1,cdist1,@tpnt)
Wend
If cdeg1<0
cdeg1+360
Else
While cdeg1>360 : cdeg1-360 : Wend
EndIf
cdist2.f = DistanceBetweenTwoPoints(*c,*p2)
cdeg2.f = DegreeAngleBetweenTwoPoints(*c,*p2)
DegreeCoordsFromPoint(*c,cdeg1,cdist1,@tpnt)
While Not PointsVirtuallyIdentical(*p2,@tpnt)
cdeg2+180
DegreeCoordsFromPoint(*c,cdeg2,cdist2,@tpnt)
Wend
If cdeg2<0
cdeg2+360
Else
While cdeg2>360 : cdeg2-360 : Wend
EndIf
mdeg.f = cdeg2-cdeg1
If mdeg<0 : mdeg+360 : ElseIf mdeg>360 : mdeg-360 : EndIf
ProcedureReturn mdeg
EndProcedure
Procedure.b TriangulateConcavePolygon(points.i,Array pnt.PointF(1),*t.Triangulation)
Dim u.UsePoint(points)
For p = 1 To points
u(p)\p = p
u(p)\valid = #True
Next p
NewList temptri.TriangleStructure()
usum = UsePointSum(u())
While usum>3
For a = 1 To usum
u1 = a
u2 = a+1
u3 = a+2
If a=usum
u2 = 1
u3 = 2
Else
If a=usum-1
u3 = 1
EndIf
EndIf
; now test whether this triangle is an "ear"
ear.b = #True
temp.TriangleStructure
CopyStructure(@pnt(u(u1)\p),@temp\a,PointF)
CopyStructure(@pnt(u(u2)\p),@temp\b,PointF)
CopyStructure(@pnt(u(u3)\p),@temp\c,PointF)
; it must have some area - x coords must differ and y coords must differ
If (temp\a\x=temp\b\x And temp\a\x=temp\c\x) Or (temp\a\y=temp\b\y And temp\a\y=temp\c\y)
ear = #False
EndIf
; it must have some area - no two consecutive points can be identical
If PointsIdentical(temp\a,temp\b) Or PointsIdentical(temp\b,temp\c) Or PointsIdentical(temp\c,temp\a)
R("CRASH. TWO POINTS THAT ARE CONSECUTIVE (AT LEAST AFTER CLIPPING) ARE ALSO IDENTICAL.")
ProcedureReturn #False
ear = #False
EndIf
; it must not be convex
If ear And AngleBetweenThreePoints(@pnt(u(u2)\p),@pnt(u(u3)\p),@pnt(u(u1)\p)) > 180
ear = #False
EndIf
; no other points must be inside it
If ear
For p = 1 To points
If p<>u(u1)\p And p<>u(u2)\p And p<>u(u3)\p
If PointIsInsideTriangle(@pnt(p),@temp)
;Debug "Point #"+Str(p)+" is inside triangle"
ear = #False
Break
EndIf
EndIf
Next p
EndIf
If ear
AddElement(temptri())
CopyStructure(@temp,@temptri(),TriangleStructure)
u(u2)\valid = #False
usum = RemoveDeadPoints(usum,u())
;Break 2
;Break
EndIf
Next a
Wend
;R("END USUM: "+Str(usum))
If usum=3
AddElement(temptri())
founds = 0
For a = 1 To usum
If u(a)\valid
p = u(a)\p
founds+1
Select founds
Case 1
CopyStructure(@pnt(p),@temptri()\a,PointF)
Case 2
CopyStructure(@pnt(p),@temptri()\b,PointF)
Case 3
CopyStructure(@pnt(p),@temptri()\c,PointF)
EndSelect
EndIf
Next a
EndIf
*t\triangles = ListSize(temptri())
Dim *t\triangle(*t\triangles)
a=0
ForEach temptri()
a+1
CopyStructure(@temptri(),*t\triangle(a),TriangleStructure)
Next
FreeList(temptri())
EndProcedure
Procedure.b TriangulatePolygon(points,Array pnt.PointF(1),*t.Triangulation)
Select IdentifyPolygonType(points,pnt())
Case #Polygon_Triangle ; the polygon IS a triangle, so no need to split it up
*t\triangles = 1
Dim *t\triangle(1)
CopyStructure(@pnt(1),*t\triangle(1)\a,PointF)
CopyStructure(@pnt(2),*t\triangle(1)\b,PointF)
CopyStructure(@pnt(3),*t\triangle(1)\c,PointF)
ProcedureReturn #True
Case #Polygon_Concave
ProcedureReturn TriangulateConcavePolygon(points,pnt(),*t)
Case #Polygon_Convex
ProcedureReturn TriangulateConvexPolygon(points,pnt(),*t)
EndSelect
EndProcedure
Procedure.b ThickLine(*p1.PointF,*p2.PointF,clr.i)
LineXY(*p1\x,*p1\y,*p2\x,*p2\y,clr)
LineXY(*p1\x+1,*p1\y,*p2\x+1,*p2\y,clr)
LineXY(*p1\x,*p1\y+1,*p2\x,*p2\y+1,clr)
LineXY(*p1\x+1,*p1\y+1,*p2\x+1,*p2\y+1,clr)
LineXY(*p1\x-1,*p1\y,*p2\x-1,*p2\y,clr)
LineXY(*p1\x,*p1\y-1,*p2\x,*p2\y-1,clr)
LineXY(*p1\x-1,*p1\y-1,*p2\x-1,*p2\y-1,clr)
EndProcedure
ExamineDesktops()
iw.f = DesktopWidth(0)
ih.f = DesktopHeight(0)
win = OpenWindow(#PB_Any,0,0,iw,ih,"Triangulation",#PB_Window_BorderLess)
SetWindowColor(win,#Black)
returnkey = 1
AddKeyboardShortcut(win,#PB_Shortcut_Return,returnkey)
escapekey = 2
AddKeyboardShortcut(win,#PB_Shortcut_Escape,escapekey)
imgad = ImageGadget(#PB_Any,0,0,iw,ih,0)
quit.b = #False
cnt.PointF
cnt\x = iw/2
cnt\y = ih/2
r.f = (ih/2)-5
Repeat
points = 5+Random(15)
;points=30
Dim pnt.PointF(points)
deg.f = 360/(points-1)
For p = 0 To points-1
dist.f = r/100*(20+Random(75))
DegreeCoordsFromPoint(@cnt,deg*p,dist,@pnt(p+1))
Next p
;If IdentifyPolygonType(points,pnt())<>#Polygon_Concave : Continue : EndIf ; use this line to see only concave, convex, etc.
pgimg = CreateImage(#PB_Any,iw,ih,32,#PB_Image_Transparent)
StartDrawing(ImageOutput(pgimg))
DrawingMode(#PB_2DDrawing_AlphaBlend)
outlineclr = RGBA(255,255,255,255)
For p = 1 To points-1
ThickLine(@pnt(p),@pnt(p+1),outlineclr)
Next p
ThickLine(@pnt(points),@pnt(1),outlineclr) ; closer line (back to first point)
StopDrawing()
img = CreateImage(#PB_Any,iw,ih,32)
StartDrawing(ImageOutput(img))
start = ElapsedMilliseconds()
TriangulatePolygon(points,pnt(),@tri.Triangulation)
finish = ElapsedMilliseconds()
If tri\triangles
For t = 1 To tri\triangles
clr = RGB(Random(255),Random(255),Random(255))
DrawTriangle(tri\triangle(t),clr)
Next t
EndIf
DrawAlphaImage(ImageID(pgimg),0,0,255)
For p = 1 To points
Circle(pnt(p)\x,pnt(p)\y,2,#White)
Next p
FreeImage(pgimg)
DrawText(20,20,Str(points)+" vertices",#White)
DrawText(20,36,Str(tri\triangles)+" triangles",#White)
DrawText(20,52,"Processing time: "+Str(finish-start)+"ms",#White)
StopDrawing()
SetGadgetState(imgad,ImageID(img))
Repeat
Delay(10)
If WindowEvent()=#PB_Event_Menu
Select EventMenu()
Case returnkey
Break
Case escapekey
quit = #True
Break
EndSelect
EndIf
ForEver
FreeImage(img)
Until quit
CloseWindow(win)
