Get oriented bounding box for a set of points

Share your advanced PureBasic knowledge/code with the community.
Post Reply
Seymour Clufley
Addict
Addict
Posts: 1265
Joined: Wed Feb 28, 2007 9:13 am
Location: London

Get oriented bounding box for a set of points

Post by Seymour Clufley »

This is for getting the smallest rectangle, of any orientation, that can possibly contain a set of points. It is also called the "arbitrarily oriented minimum bounding box" or the "minimum area enclosing rectangle". This is something I need fairly often, and maybe other people will find uses for it too.

Note: The point set passed to the procedure must be a convex polygon, or the convex hull of a point set. To get the latter, you can use this code kindly provided by Said.

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



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 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 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



Procedure.b ComputeOrientedBoundingBox(chpnts.i,Array chpnt.PointD(1),Array bb.PointD(1)) ; "ch" = convex hull
  
  ComputePolygonCentroid(chpnts,chpnt(),@cntrd.PointD)
  
  Dim testchpnt.PointD(chpnts)
  smallestmass.d = 99999999
  For p1 = 1 To chpnts
    p2=p1+1
    If p1=chpnts : p2=1 : EndIf
    deg.d = DegreeAngleBetweenTwoPoints(@chpnt(p1),@chpnt(p2))
    
    minx.d = 9999999
    miny.d = 9999999
    maxx.d = -9999999
    maxy.d = -9999999
    For p = 1 To chpnts
      DegreeRotatePointAroundPoint(@chpnt(p),@cntrd,-deg,@testchpnt(p))
      DefeatThis(minx,testchpnt(p)\x)
      DefeatThis(miny,testchpnt(p)\y)
      BeatThis(maxx,testchpnt(p)\x)
      BeatThis(maxy,testchpnt(p)\y)
    Next p
    For p = 1 To chpnts
      testchpnt(p)\x - minx
      testchpnt(p)\y - miny
    Next p
    maxx-minx
    maxy-miny
    mass.d = maxx*maxy
    If mass<smallestmass
      smallestmass = mass
      smallestmassanchor = p1
    EndIf
  Next p1
  
  
  p2=smallestmassanchor+1
  If smallestmassanchor=chpnts : p2=1 : EndIf
  deg.d = DegreeAngleBetweenTwoPoints(@chpnt(smallestmassanchor),@chpnt(p2))
  
  minx.d = 9999999
  miny.d = 9999999
  maxx.d = -9999999
  maxy.d = -9999999
  For p = 1 To chpnts
    DegreeRotatePointAroundPoint(@chpnt(p),@cntrd,-deg,@testchpnt(p))
    DefeatThis(minx,testchpnt(p)\x)
    DefeatThis(miny,testchpnt(p)\y)
    BeatThis(maxx,testchpnt(p)\x)
    BeatThis(maxy,testchpnt(p)\y)
  Next p
  For p = 1 To chpnts
    testchpnt(p)\x - minx
    testchpnt(p)\y - miny
  Next p
  maxx-minx
  maxy-miny
  
  
  Dim tempbb.PointD(4)
  Dim bb(4)
  tempbb(1)\x = minx
  tempbb(1)\y = miny
  tempbb(2)\x = minx+maxx
  tempbb(2)\y = miny
  tempbb(3)\x = minx+maxx
  tempbb(3)\y = miny+maxy
  tempbb(4)\x = minx
  tempbb(4)\y = miny+maxy
  For a = 1 To 4
    DegreeRotatePointAroundPoint(@tempbb(a),@cntrd,deg,@bb(a))
  Next a
  
  ProcedureReturn #True
EndProcedure
Example usage:

Code: Select all

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



Macro StandardReportingWindowStuff(win)
   escapekey = 1
   returnkey = 2
   AddKeyboardShortcut(win,#PB_Shortcut_Escape,escapekey)
   AddKeyboardShortcut(win,#PB_Shortcut_Return,returnkey)
   
   MoveWindowToTop(win,#True)
   
   
   Repeat
      we = WindowEvent()
      If we
         If we=#PB_Event_Menu
            Break
         EndIf
      Else
         Delay(10)
      EndIf
   ForEver
   
   CloseWindow(win)
   
EndMacro

Procedure.b MoveWindowToTop(win.i,stayontop.b)
   If Not WindowID(win) : MessageRequester("Error","CAN'T FIND WINDOW ID") : ProcedureReturn #False : EndIf
   If stayontop
      SetWindowPos_(WindowID(win),#HWND_TOPMOST,0,0,0,0,#SWP_NOMOVE|#SWP_NOSIZE)    ; gets it on top to stay 
   Else
      SetWindowPos_(WindowID(win),#HWND_NOTOPMOST,0,0,0,0,#SWP_NOMOVE|#SWP_NOSIZE)  ; now it is on top, allow it to go behind
   EndIf
   ProcedureReturn #True
EndProcedure

Procedure.b ResizeToFitInsideFrame(iw.d,ih.d,fw.d,fh.d,*ns.PointD)
   
   *ns\x = iw
   *ns\y = ih
   
   While *ns\x>fw Or *ns\y>fh
     *ns\x * 0.99
     *ns\y * 0.99
   Wend
   ProcedureReturn #True
   
EndProcedure

Procedure RI(img.i,title.s="",bgclr.i=#Green,fit_on_screen.b=#True,timelimit.i=0,free.b=#False)
  
  If Not IsImage(img)
    MessageRequester("Error","PROC RI (ReportImage):"+c13+"img is not an image!")
    ProcedureReturn
  EndIf
  
  iw.d = ImageWidth(img)
  ih.d = ImageHeight(img)
  
  simg = CreateImage(#PB_Any,iw,ih,32,bgclr)
  StartDrawing(ImageOutput(simg))
  DrawAlphaImage(ImageID(img),0,0)
  StopDrawing()
  
  If fit_on_screen
    d = ExamineDesktops()
    ResizeToFitInsideFrame(iw,ih,DesktopWidth(0),DesktopHeight(0),@ns.PointD)
    iw = ns\x
    ih = ns\y
    ResizeImage(simg,iw,ih)
  EndIf
  
  If title=""
    title = "Report Image"
  EndIf
  
  win = OpenWindow(#PB_Any,0,0,iw,ih,title,#PB_Window_BorderLess|#PB_Window_ScreenCentered)
  
  imgad = ImageGadget(#PB_Any,0,0,iw,ih,ImageID(simg))
  
  StandardReportingWindowStuff(win)
  
  FreeImage(simg)
  
EndProcedure



Macro XYDrawable(x,y)
  Bool(x=>0 And y=>0 And x<(OutputWidth()-1) And y<(OutputHeight()-1))
EndMacro

Procedure.b LineAA(*p1.PointD,*p2.PointD,clr.i,Thickness.d=1)
	
	width.d = *p2\x-*p1\x
	hight.d = *p2\y-*p1\y
	Protected SensX, SensY, n, nn, Epaisseur.d, x2.d, y2.d, Couleur_Fond.l, Application.d, Distance.d
	; On mets la droite toujours dans le même sens pour l'analyse
	; La sauvegarde du sens permettra de dessiner la droite ensuite dans le bon sens
	If Width >= 0
		SensX = 1
	Else
		SensX = -1
		Width = - Width
	EndIf
	If Hight >= 0
		SensY = 1
	Else
		SensY = -1
		Hight = - Hight
	EndIf
	
	
	; Demi épaisseur de la ligne
	Epaisseur.d = Thickness / 2
	
	; calcul pour le changement de repère qui permet de connaitre l'épaisseur du trait et de gérer l'AA
	Distance.d = Sqr(Width * Width + Hight * Hight)
	CosAngle.d = Width / Distance
	SinAngle.d = -Sin(ACos(CosAngle))
	
	; Dessin de la ligne
	For n = -Thickness To Width + Thickness
		For nn = -Thickness To Hight + Thickness
			dx = *p1\x + n * SensX
			dy = *p1\y + nn * SensY
			If XYDrawable(dx,dy)
				
				; changement de base
				; les y représentent l'épaisseur de la ligne
				x2 = n * CosAngle - nn * SinAngle
				y2 = Abs(n * SinAngle + nn * CosAngle)
				
				;If XYDrawable(x2,y2)
				If y2 <= Epaisseur + 0.5
					Application =  0.5 + Epaisseur - y2
					If Application > 1
						Application = 1
					EndIf
					If x2 > -1 And x2 < Distance + 1
						If x2 < 0
							Application * (1 + x2)
						ElseIf x2 > Distance
							Application * (1 - x2 + Distance)
						EndIf
					Else
						Application = 0
					EndIf
					If Application > 0
						Plot(dx,dy,clr)
					EndIf
				EndIf
			EndIf
		Next
	Next
	
EndProcedure



iw.d = 1600
ih.d = 900
img = CreateImage(#PB_Any,iw,ih,32)
cnt.PointD
cnt\x = iw/2
cnt\y = ih/2
Dim bb.PointD(4)
Repeat
  pnts.i = 3+Random(12)
  Dim pnt.PointD(pnts)
  For a = 1 To pnts
    DegreeCoordsFromPoint(@cnt,360/pnts*a,ih/(2.5+Random(4)),@pnt(a))
  Next a
  ComputeOrientedBoundingBox(pnts,pnt(),bb())
  StartDrawing(ImageOutput(img))
  Box(0,0,iw,ih,#Black)
  For a1 = 1 To pnts
    a2=a1+1 : If a1=pnts : a2=1 : EndIf
    LineAA(@pnt(a1),@pnt(a2),#White,4)
  Next a1
  For a1 = 1 To 4
    a2=a1+1 : If a1=4 : a2=1 : EndIf
    LineAA(@bb(a1),@bb(a2),#Red,1)
  Next a1
  StopDrawing()
  RI(img)
Until GetAsyncKeyState_(#VK_ESCAPE)
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
said
Enthusiast
Enthusiast
Posts: 342
Joined: Thu Apr 14, 2011 6:07 pm

Re: Get oriented bounding box for a set of points

Post by said »

.... you can use this code kindly provided by Said
Thanks but ... this is not true :o It is your code and i just helped you fixing few issues

Thanks for sharing, not sure i have a use for a bounding box but it is always good to have something like this (many problems do share common tricks/tasks ...) :D
Top
Post Reply

Return to “Tricks 'n' Tips”