Page 1 of 1

Computing the convex hull for a set of points

Posted: Tue Jun 13, 2017 7:33 pm
by Seymour Clufley
Based on the pseudocode here, I am trying to write a convex hull procedure. Unfortunately it isn't working. I think the problem might be on line 124 - I'm not sure what the comment means by "polar angle",or if that sorting is what the rest of the code actually does. Can anyone who is better at maths help out?

Code: Select all

Structure PointD
  x.d
  y.d
EndStructure


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




Procedure.d ccw(*p1.PointD, *p2.PointD, *p3.PointD)
  ProcedureReturn ((*p2\x - *p1\x)*(*p3\y - *p1\y)) - ((*p2\y - *p1\y)*(*p3\x - *p1\x))
EndProcedure

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.i ConvexHull(pnts.i,Array opnt.PointD(1),Array pnt.PointD(1))
  ;let pnt(pnts+1) = the Array of points
  
  Dim pnt.PointD(pnts)
  CopyArray(opnt(),pnt())
  
  lowest_y.d = 99999999
  For a = 1 To pnts
    If pnt(a)\y<lowest_y
      lowest_y = pnt(a)\y
      ly=a
    EndIf
  Next a
  SwapPoints( pnt(1), pnt(ly) ) ; the point with the lowest y-coordinate
  
  ;sort points by polar angle With pnt(1)
  
  ; We want pnt[0] to be a sentinel point that will stop the loop
  ;let pnt(0) = pnt(pnts)
  pnt(0)\x=pnt(pnts)\x : pnt(0)\y=pnt(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(pnt(M-1), pnt(M), pnt(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( pnt(M) , pnt(i) )
  Next i
  ProcedureReturn M
EndProcedure




iw = 1000
ih = 700
img = CreateImage(#PB_Any,iw,ih,32)

pnts.i = 10+Random(20)
Dim pnt.PointD(pnts)
#Margin = 200
For p = 1 To pnts
  pnt(p)\x = #Margin+Random(iw-#Margin-#Margin)
  pnt(p)\y = #Margin+Random(ih-#Margin-#Margin)
Next p

StartDrawing(ImageOutput(img))
For p = 1 To pnts
  Circle(pnt(p)\x,pnt(p)\y,8,RGB(48,48,48))
Next p

Dim chpnt.PointD(1)
chpnts.i = ConvexHull(pnts,pnt(),chpnt())
For p = 0 To chpnts-1
  Circle(chpnt(p)\x,chpnt(p)\y,1,#White)
Next p

StopDrawing()
RI(img)

Re: Computing the convex hull for a set of points

Posted: Tue Jun 13, 2017 10:52 pm
by said
Hi,

You can try this code:

Code: Select all

Structure PointD
  x.d
  y.d
  d.d
EndStructure


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




Procedure.d ccw(*p1.PointD, *p2.PointD, *p3.PointD)
  ProcedureReturn ((*p2\x - *p1\x)*(*p3\y - *p1\y)) - ((*p2\y - *p1\y)*(*p3\x - *p1\x))
EndProcedure

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.i ConvexHull(pnts.i,Array opnt.PointD(1),Array pnt.PointD(1))
  ;let pnt(pnts+1) = the Array of points
  
  Dim pnt.PointD(pnts)
  CopyArray(opnt(),pnt())
  
  SortStructuredArray(pnt(), #PB_Sort_Descending, OffsetOf(PointD\y), TypeOf(PointD\y), 1, pnts)
  
  ; polar angle between all points with pnt(1) ... stored in \d
  For a = 2 To pnts
    pnt(a)\d = ATan2(pnt(1)\y - pnt(a)\y, pnt(1)\x - pnt(a)\x)
  Next
  
  ;sort points by polar angle With pnt(1)
  SortStructuredArray(pnt(), #PB_Sort_Descending, OffsetOf(PointD\d), TypeOf(PointD\d), 2, pnts)
  
  
;   lowest_y.d = 99999999
;   For a = 1 To pnts
;     If pnt(a)\y<lowest_y
;       lowest_y = pnt(a)\y
;       ly=a
;     EndIf
;   Next a
;   SwapPoints( pnt(1), pnt(ly) ) ; the point with the lowest y-coordinate
;   
;   ;sort points by polar angle With pnt(1)
;   
   ; We want pnt[0] to be a sentinel point that will stop the loop
   ;let pnt(0) = pnt(pnts)
   pnt(0)\x=pnt(pnts)\x : pnt(0)\y=pnt(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(pnt(M-1), pnt(M), pnt(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( pnt(M) , pnt(i) )
  Next i
  ProcedureReturn M
EndProcedure




iw = 1000
ih = 700
img = CreateImage(#PB_Any,iw,ih,32)

pnts.i = 10+Random(20)
Dim pnt.PointD(pnts)
#Margin = 200
For p = 1 To pnts
  pnt(p)\x = #Margin+Random(iw-#Margin-#Margin)
  pnt(p)\y = #Margin+Random(ih-#Margin-#Margin)
Next p

StartDrawing(ImageOutput(img))
For p = 1 To pnts
  Circle(pnt(p)\x,pnt(p)\y,2,#White);RGB(48,48,48))
Next p

Dim chpnt.PointD(1)
chpnts.i = ConvexHull(pnts,pnt(),chpnt())

For p = 0 To chpnts-1
  LineXY(chpnt(p)\x,chpnt(p)\y, chpnt(p+1)\x,chpnt(p+1)\y,#Red)
Next p
LineXY(chpnt(chpnts)\x,chpnt(chpnts)\y, chpnt(0)\x,chpnt(0)\y,#Red)

StopDrawing()
RI(img)

Re: Computing the convex hull for a set of points

Posted: Tue Jun 13, 2017 11:50 pm
by Seymour Clufley
Thanks very much! :)

Re: Computing the convex hull for a set of points

Posted: Wed Jun 14, 2017 4:37 pm
by Comtois
Nice :)

My version using Jarvis algo (animation added for fun) :
http://www.purebasic.fr/english/viewtop ... 12&t=14049

Re: Computing the convex hull for a set of points

Posted: Wed Jun 14, 2017 7:27 pm
by Demivec
Also for reference, JHPJHP has included solutions for this and other image problems in his demos for his PureBasic interface to OpenCV.