PureBasic Forums - English

http://www.purebasic.com
https://www.purebasic.fr/english/

Get oriented bounding box for a set of points

https://www.purebasic.fr/english/viewtopic.php?t=68639

Page 1 of 1

Get oriented bounding box for a set of points

Posted: Wed Jun 14, 2017 9:06 am
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)

Re: Get oriented bounding box for a set of points

Posted: Wed Jun 14, 2017 2:17 pm
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
All times are UTC+01:00
Page 1 of 1

Powered by phpBB® Forum Software © phpBB Limited