Smallest Enclosing Circle / Minimum Bounding Circle
Posted: Fri Aug 14, 2020 9:02 pm
This is for finding the smallest circle that will enclose a set of arbitrary points. Derived from this tutorial, it utilises functions that I have posted here and here. For simplicity, I'm just including all of the necessary code below.
Demo code (press SPACE to cycle through and ESCAPE to quit):
Code: Select all
;-
;- GENERIC FUNCTIONS
Macro CopyPoint(cpp1,cpp2)
cpp2\x=cpp1\x
cpp2\y=cpp1\y
EndMacro
Macro SwapPoints(spp1,spp2)
slave.PointD
slave\x = spp1\x
spp1\x = spp2\x
spp2\x = slave\x
slave\y = spp1\y
spp1\y = spp2\y
spp2\y = slave\y
EndMacro
Structure PointD
x.d
y.d
EndStructure
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.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)
ProcedureReturn Sqr(Pow(*a\x - *b\x, 2) + Pow(*a\y - *b\y, 2))
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
;-
;- CONVEX HULL CODE
Structure CHPointD
x.d
y.d
d.d
EndStructure
Procedure.d ccw(*p1.CHPointD, *p2.CHPointD, *p3.CHPointD)
ProcedureReturn ((*p2\x - *p1\x)*(*p3\y - *p1\y)) - ((*p2\y - *p1\y)*(*p3\x - *p1\x))
EndProcedure
Procedure.i ComputeConvexHull(pnts.i,Array opnt.PointD(1),Array chpnt.PointD(1))
Dim tempchpnt.CHPointD(pnts)
For a = 1 To pnts
CopyPoint(opnt(a),tempchpnt(a))
Next
SortStructuredArray(tempchpnt(), #PB_Sort_Descending, OffsetOf(CHPointD\y), #PB_Double, 1, pnts)
For a = 2 To pnts
tempchpnt(a)\d = ATan2(tempchpnt(1)\y - tempchpnt(a)\y, tempchpnt(1)\x - tempchpnt(a)\x)
Next
SortStructuredArray(tempchpnt(), #PB_Sort_Descending, OffsetOf(CHPointD\d), #PB_Double, 2, pnts)
tempchpnt(0)\x=tempchpnt(pnts)\x : tempchpnt(0)\y=tempchpnt(pnts)\y
M = 1
For i = 2 To pnts
While ccw(tempchpnt(M-1), tempchpnt(M), tempchpnt(i)) <= 0
If M > 1
M - 1
Continue
Else
If i = pnts
Break
Else
i + 1
EndIf
EndIf
Wend
M + 1
SwapPoints( tempchpnt(M) , tempchpnt(i) )
Next i
ReDim chpnt(M)
For p = 1 To M
CopyPoint(tempchpnt(p),chpnt(p))
Next p
ProcedureReturn M
EndProcedure
;-
;- "THREE POINTS" CODE
Structure CircleD
c.PointD
r.d
EndStructure
Procedure.b FindCircleThroughThreePoints_Centre(*b.PointD, *c.PointD,*I.PointD)
B.d = *b\x * *b\x + *b\y * *b\y
C.d = *c\x * *c\x + *c\y * *c\y
D.d = *b\x * *c\y - *b\y * *c\x
*I\x = (*c\y * B - *b\y * C) / (2 * D)
*I\y = (*b\x * C - *c\x * B) / (2 * D)
EndProcedure
Procedure.b FindCircleThroughThreePoints(*A.PointD, *B.PointD, *C.PointD, *cir.CircleD)
test1.PointD
test1\x = *B\x - *A\x
test1\y = *B\y - *A\y
test2.PointD
test2\x = *C\x - *A\x
test2\y = *C\y - *A\y
FindCircleThroughThreePoints_Centre(@test1, @test2, *cir\c)
*cir\c\x + *A\x
*cir\c\y + *A\y
*cir\r = DistanceBetweenTwoPoints(*cir\c,*A)
EndProcedure
;-
;- "SMALLEST ENCLOSING CIRCLE" CODE
Procedure.b CircleContainsPoint(*c.CircleD,*p.PointD)
If DistanceBetweenTwoPoints(*c\c,*p) > *c\r
ProcedureReturn #False
EndIf
ProcedureReturn #True
EndProcedure
Procedure.b CircleContainsPoints(*c.CircleD,pnts.i,Array pnt.PointD(1))
For p = 1 To pnts
If DistanceBetweenTwoPoints(*c\c,@pnt(p)) > *c\r
ProcedureReturn #False
EndIf
Next p
ProcedureReturn #True
EndProcedure
Macro CopyCircle(cc1,cc2)
CopyPoint(cc1\c,cc2\c)
cc2\r = cc1\r
EndMacro
Procedure.b ComputeMinimumBoundingCircle(pnts.i,Array pnt.PointD(1),*c.CircleD)
; get the convex hull
Dim chpnt.PointD(1)
chpnts.i = ComputeConvexHull(pnts,pnt(),chpnt())
*c\r = 9999999
; look at pairs of hull points...
For i = 1 To chpnts-1
For j = i+1 To chpnts
; find the circle through these two points
test.CircleD
MidPoint(@chpnt(i),@chpnt(j),@test\c)
test\r = DistanceBetweenTwoPoints(@test\c,@chpnt(i))
; see if this circle would be an improvement
If test\r < *c\r
; see if this circle encloses all of the points
If CircleContainsPoints(@test,chpnts,chpnt())
; save this solution
CopyCircle(test,*c)
EndIf
EndIf
Next j
Next i
; look at triples of hull points...
For i = 0 To chpnts-2
For j = i+1 To chpnts-1
For k = j+1 To chpnts
; find the circle through these three points
FindCircleThroughThreePoints(@chpnt(i),@chpnt(j),@chpnt(k),@test.CircleD)
; see if this circle would be an improvement
If test\r < *c\r
; see if this circle encloses all the points
If CircleContainsPoints(@test,chpnts,chpnt())
; save this solution
CopyCircle(test,*c)
EndIf
EndIf
Next k
Next j
Next i
EndProcedure
Code: Select all
iw = 1600
ih = 900
img = CreateImage(#PB_Any,iw,ih)
win = OpenWindow(#PB_Any,0,0,iw,ih,"Minimum bounding circle",#PB_Window_ScreenCentered)
imgad = ImageGadget(#PB_Any,0,0,iw,ih,ImageID(img))
space = 5
AddKeyboardShortcut(win,#PB_Shortcut_Space,space)
esc = 6
AddKeyboardShortcut(win,#PB_Shortcut_Escape,esc)
Repeat
pnts = Random(40,4)
Dim pnt.PointD(pnts)
For p = 1 To pnts
pnt(p)\x = Random(iw*0.7,iw*0.3)
pnt(p)\y = Random(ih*0.7,ih*0.3)
Next p
ClearDebugOutput()
ComputeMinimumBoundingCircle(pnts,pnt(),@cir.CircleD)
img = CreateImage(#PB_Any,iw,ih,24,#Black)
StartDrawing(ImageOutput(img))
Circle(cir\c\x,cir\c\y,cir\r,#Blue)
For p = 1 To pnts
Circle(pnt(p)\x,pnt(p)\y,5,#Yellow)
Next p
StopDrawing()
SetGadgetState(imgad,ImageID(img))
Repeat : Until WaitWindowEvent(5)=#PB_Event_Menu
If EventMenu()=esc : Break : EndIf
ForEver