Minimum area bounding box

Seymour Clufley · Post by **Seymour Clufley** » Sun Aug 30, 2020 7:12 pm

This code finds the smallest bounding box at any angle, returning both the bounding box points and the degree angle at which they were found.

This box is also known as the minimum area enclosing rectangle.

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

Code: Select all

Macro DefeatThis(a,b)
  If a>b
      a=b
  EndIf
EndMacro
Macro BeatThis(a,b)
  If b>a
      a=b
  EndIf
EndMacro




Structure PointD
  x.d
  y.d
EndStructure

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




Procedure.b ComputePolygonCentroid(points.i,Array vertex.PointD(1),*centroid.PointD)
  signedArea.d = 0
  x0 = 0; Current vertex X
  y0 = 0; Current vertex Y
  x1 = 0; Next vertex X
  y1 = 0; Next vertex Y
  a = 0 ;  Partial signed area
  
  ; for all points except last
  For i = 1 To points-1
    x0 = vertex(i)\x
    y0 = vertex(i)\y
    x1 = vertex(i+1)\x
    y1 = vertex(i+1)\y
    a = x0*y1 - x1*y0
    signedArea + a
    *centroid\x + (x0 + x1)*a
    *centroid\y + (y0 + y1)*a
  Next
  
  ; do last vertex
  x0 = vertex(points)\x
  y0 = vertex(points)\y
  x1 = vertex(1)\x
  y1 = vertex(1)\y
  a = x0*y1 - x1*y0
  signedArea + a
  *centroid\x + (x0 + x1)*a
  *centroid\y + (y0 + y1)*a
  
  signedArea * 0.5
  *centroid\x / (6*signedArea)
  *centroid\y / (6*signedArea)
EndProcedure




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), TypeOf(CHPointD\y), 1, pnts)
  
  ; polar angle between all points with pnt(1) ... stored in \d
  For a = 2 To pnts
    tempchpnt(a)\d = ATan2(tempchpnt(1)\y - tempchpnt(a)\y, tempchpnt(1)\x - tempchpnt(a)\x)
  Next
  
  ;sort points by polar angle With pnt(1)
  SortStructuredArray(tempchpnt(), #PB_Sort_Descending, OffsetOf(CHPointD\d), TypeOf(CHPointD\d), 2, pnts)
  
  ; We want pnt[0] to be a sentinel point that will stop the loop
  tempchpnt(0)\x=tempchpnt(pnts)\x : tempchpnt(0)\y=tempchpnt(pnts)\y
  
  ; M will denote the number of points on the convex hull
  M = 1
  For i = 2 To pnts
    ; Find next valid point on convex hull
    While ccw(tempchpnt(M-1), tempchpnt(M), tempchpnt(i)) <= 0
      If M > 1
        M - 1
        Continue
        ; All points are collinear
      Else
        If i = pnts
          Break
        Else
          i + 1
        EndIf
      EndIf
    Wend
    
    ; Update M and Swap pnt[i] to the correct place
    M + 1
    SwapPoints( tempchpnt(M) , tempchpnt(i) )
  Next i
  
  
  Dim chpnt(M)
  For p = 1 To M
    CopyPoint(tempchpnt(p),chpnt(p))
  Next p
  
  ProcedureReturn M
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.b DegreeRotatePointAroundPoint(*p.PointD,*rc.PointD,degrees.d,*np.PointD)
  radia.d = Radian(degrees)
  *np\x = *rc\x + ( Cos(radia) * (*p\x - *rc\x) - Sin(radia) * (*p\y - *rc\y) )
  *np\y = *rc\y + ( Sin(radia) * (*p\x - *rc\x) + Cos(radia) * (*p\y - *rc\y) )
EndProcedure




Procedure.d ComputeSmallestBoundingBox(pnts.i,Array pnt.PointD(1),Array perimrect.PointD(1))
  
  Dim hullpnt.PointD(pnts)
  chpnts.i = ComputeConvexHull(pnts,pnt(),hullpnt())
  
  ComputePolygonCentroid(chpnts,hullpnt(),@cd.PointD)
  
  smallest_area.d = 99999999
  smallest_deg.d = 0
  ReDim perimrect(4)
  
  For s = 1 To chpnts
    s2=s+1 : If s=chpnts : s2=1 : EndIf
    deg.d = DegreeAngleBetweenTwoPoints(@hullpnt(s),@hullpnt(s2))
    
    xmin.d = 99999999
    ymin.d = 99999999
    xmax.d = -99999999
    ymax.d = -99999999
    
    ; rotate all points...
    Dim rotpnt.PointD(chpnts)
    For rs = 1 To chpnts
      DegreeRotatePointAroundPoint(@hullpnt(rs),@cd,-deg,@rotpnt(rs))
      DefeatThis(xmin,rotpnt(rs)\x)
      DefeatThis(ymin,rotpnt(rs)\y)
      BeatThis(xmax,rotpnt(rs)\x)
      BeatThis(ymax,rotpnt(rs)\y)
    Next rs
    
    this_area.d = (xmax-xmin) * (ymax-ymin)
    If this_area<smallest_area
      smallest_area = this_area
      smallest_deg = deg
      
      perimrect(1)\x=xmin : perimrect(1)\y=ymin
      perimrect(2)\x=xmax : perimrect(2)\y=ymin
      perimrect(3)\x=xmax : perimrect(3)\y=ymax
      perimrect(4)\x=xmin : perimrect(4)\y=ymax
    EndIf
  Next s
  
  For r = 1 To 4
    DegreeRotatePointAroundPoint(@perimrect(r),@cd,smallest_deg,@np.PointD)
    CopyPoint(np,perimrect(r))
  Next r
  ProcedureReturn smallest_deg
  
EndProcedure




iw = 1600
ih = 800
win = OpenWindow(#PB_Any,0,0,iw,ih,"Minimum Perimeter Rectangle",#PB_Window_ScreenCentered)
cnv = CanvasGadget(#PB_Any,0,0,iw,ih)
space=5
AddKeyboardShortcut(win,#PB_Shortcut_Space,space)
esc=6
AddKeyboardShortcut(win,#PB_Shortcut_Escape,esc)


Repeat
  pnts = Random(20,5)
  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.8,ih*0.2)
  Next p
  
  Dim rotpnt.PointD(4)
  deg.d = ComputeSmallestBoundingBox(pnts,pnt(),rotpnt())
  rpnts=4
  
  
  StartDrawing(CanvasOutput(cnv))
  Box(0,0,OutputWidth(),OutputHeight(),#Black)
  
  For p = 1 To pnts
    Circle(pnt(p)\x,pnt(p)\y,5,#Red)
  Next p
  
  For p = 1 To rpnts
    p2=p+1 : If p=rpnts : p2=1 : EndIf
    LineXY(rotpnt(p)\x,rotpnt(p)\y,rotpnt(p2)\x,rotpnt(p2)\y,#Green)
  Next p
  
  StopDrawing()
  
  Repeat : Until WaitWindowEvent(5)=#PB_Event_Menu
  If EventMenu()=esc : Break : EndIf
ForEver