Smallest Enclosing Circle / Minimum Bounding Circle

Seymour Clufley · Post by **Seymour Clufley** » 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.

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

Demo code (press SPACE to cycle through and ESCAPE to quit):

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

Seymour Clufley · Post by **Seymour Clufley** » Sat Aug 15, 2020 1:58 am

I've just realised that you can also do it without using the convex hull. This eliminates quite a lot of the code above but will be a bit slower in execution since all of the points are examined. Here is the alternative version of the procedure:

Code: Select all

Procedure.b ComputeMinimumBoundingCircle(pnts.i,Array pnt.PointD(1),*c.CircleD)
 
  *c\r = 9999999
 
  ; look at pairs of points...
  For i = 1 To pnts-1
    For j = i+1 To pnts
      ; find the circle through these two points
      test.CircleD
      MidPoint(@pnt(i),@pnt(j),@test\c)
      test\r = DistanceBetweenTwoPoints(@test\c,@pnt(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,pnts,pnt())
          ; save this solution
          CopyCircle(test,*c)
        EndIf
      EndIf
    Next j
  Next i
 
  ; look at triples of points...
  For i = 0 To pnts-2
    For j = i+1 To pnts-1
      For k = j+1 To pnts
        ; find the circle through these three points
        FindCircleThroughThreePoints(@pnt(i),@pnt(j),@pnt(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,pnts,pnt())
            ; save this solution
            CopyCircle(test,*c)
          EndIf
        EndIf
      Next k
    Next j
  Next i
 
EndProcedure

PureBasic Forums - English

Smallest Enclosing Circle / Minimum Bounding Circle

Smallest Enclosing Circle / Minimum Bounding Circle

Re: Smallest Enclosing Circle / Minimum Bounding Circle