It is currently Wed Oct 21, 2020 3:27 am

All times are UTC + 1 hour




Post new topic Reply to topic  [ 1 post ] 
Author Message
 Post subject: Minimum area bounding box
PostPosted: Sun Aug 30, 2020 7:12 pm 
Offline
Addict
Addict

Joined: Wed Feb 28, 2007 9:13 am
Posts: 1063
Location: London
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:
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

_________________
JACK WEBB: "Coding in C is like sculpting a statue using only sandpaper. You can do it, but the result wouldn't be any better. So why bother? Just use the right tools and get the job done."


Top
 Profile  
Reply with quote  
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 1 post ] 

All times are UTC + 1 hour


Who is online

Users browsing this forum: No registered users and 16 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum

Search for:
Jump to:  

 


Powered by phpBB © 2008 phpBB Group
subSilver+ theme by Canver Software, sponsor Sanal Modifiye