Computing the convex hull for a set of points

Just starting out? Need help? Post your questions and find answers here.
Seymour Clufley
Addict
Addict
Posts: 1265
Joined: Wed Feb 28, 2007 9:13 am
Location: London

Computing the convex hull for a set of points

Post 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)
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."
said
Enthusiast
Enthusiast
Posts: 342
Joined: Thu Apr 14, 2011 6:07 pm

Re: Computing the convex hull for a set of points

Post 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)
Seymour Clufley
Addict
Addict
Posts: 1265
Joined: Wed Feb 28, 2007 9:13 am
Location: London

Re: Computing the convex hull for a set of points

Post by Seymour Clufley »

Thanks very much! :)
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."
User avatar
Comtois
Addict
Addict
Posts: 1431
Joined: Tue Aug 19, 2003 11:36 am
Location: Doubs - France

Re: Computing the convex hull for a set of points

Post by Comtois »

Nice :)

My version using Jarvis algo (animation added for fun) :
http://www.purebasic.fr/english/viewtop ... 12&t=14049
Please correct my english
http://purebasic.developpez.com/
User avatar
Demivec
Addict
Addict
Posts: 4270
Joined: Mon Jul 25, 2005 3:51 pm
Location: Utah, USA

Re: Computing the convex hull for a set of points

Post by Demivec »

Also for reference, JHPJHP has included solutions for this and other image problems in his demos for his PureBasic interface to OpenCV.
Post Reply