Find the points which form a convex hull from a set of arbitrary two dimensional points.

For example, given the points (16,3), (12,17), (0,6), (-4,-6), (16,6), (16,-7), (16,-3), (17,-4), (5,19), (19,-8), (3,16), (12,13), (3,-4), (17,5), (-3,15), (-3,-9), (0,11), (-9,-3), (-4,-2) and (12,10) the convex hull would be (-9,-3), (-3,-9), (19,-8), (17,5), (12,17), (5,19) and (-3,15).




See also
 Convex Hull (youtube)
 http://www.geeksforgeeks.org/convex-hull-set-2-graham-scan/


 
#include "crt.bi"
Screen 20
Window (-20,-20)-(30,30)
 
Type Point
    As Single x,y
    As Long done
End Type
 
#macro rotate(pivot,p,a,scale)
Type<Point>(scale*(Cos(a*.0174533)*(p.x-pivot.x)-Sin(a*.0174533)*(p.y-pivot.y))+pivot.x, _
scale*(Sin(a*.0174533)*(p.x-pivot.x)+Cos(a*.0174533)*(p.y-pivot.y))+pivot.y)
#endmacro
 
 
Dim As Point p(1 To ...)={(16,3),(12,17),(0,6),(-4,-6),(16,6),(16,-7),(16,-3),(17,-4),(5,19), _
(19,-8),(3,16),(12,13),(3,-4),(17,5),(-3,15),(-3,-9),(0,11),(-9,-3),(-4,-2),(12,10)}
 
 
Function south(p() As Point,Byref idx As Long) As Point
    Dim As Point s=Type(0,100)
    For n As Long=Lbound(p) To Ubound(p)
        Circle(p(n).x,p(n).y),.2,7,,,,f 
        If s.y>p(n).y Then s=p(n):idx=n 
    Next n
    Return s
End Function
 
Function segment_distance(lx1 As Single, _
                          ly1 As Single, _
                          lx2 As Single, _
                          ly2 As Single, _
                          px As Single,_
                          py As Single, _
                          Byref ox As Single=0,_
                          Byref oy As Single=0) As Single
    Dim As Single M1,M2,C1,C2,B
    B=(Lx2-Lx1):If B=0 Then B=1e-20
    M2=(Ly2-Ly1)/B:If M2=0 Then M2=1e-20
    M1=-1/M2
    C1=py-M1*px
    C2=(Ly1*Lx2-Lx1*Ly2)/B
    Var L1=((px-lx1)*(px-lx1)+(py-ly1)*(py-ly1)),L2=((px-lx2)*(px-lx2)+(py-ly2)*(py-ly2))
    Var a=((lx1-lx2)*(lx1-lx2) + (ly1-ly2)*(ly1-ly2))
    Var a1=a+L1
    Var a2=a+L2
    Var f1=a1>L2,f2=a2>L1
    If f1 Xor f2 Then
        Var d1=((px-Lx1)*(px-Lx1)+(py-Ly1)*(py-Ly1))
        Var d2=((px-Lx2)*(px-Lx2)+(py-Ly2)*(py-Ly2))
        If d1<d2 Then Ox=Lx1:Oy=Ly1 : Return Sqr(d1) Else  Ox=Lx2:Oy=Ly2:Return Sqr(d2)
    End If
    Var M=M1-M2:If M=0 Then M=1e-20
    Ox=(C2-C1)/(M1-M2)
    Oy=(M1*C2-M2*C1)/M
    Return Sqr((px-Ox)*(px-Ox)+(py-Oy)*(py-Oy))
End Function
 
 
Dim As Long idx
Var s= south(p(),idx)
p(idx).done=1
Redim As Point ans(1 To 1)
ans(1)=s
Dim As Point e=s
e.x=1000
Dim As Long count=1
Dim As Single z
Circle(s.x,s.y),.4,5
 
Do
    z+=.05
    Var pt=rotate(s,e,z,1)
    For n As Long=Lbound(p) To Ubound(p)
        If segment_distance(s.x,s.y,pt.x,pt.y,p(n).x,p(n).y)<.05  Then
            s=p(n)
            If p(n).done=0 Then
                count+=1
                Redim Preserve ans(1 To count)
                ans(count)=p(n)
                p(n).done=1
            End If
        End If
        Circle(s.x,s.y),.4,5
    Next n
Loop Until z>360
 
For n As Long=Lbound(ans) To Ubound(ans)
    printf (!"(%2.0f , %2.0f )\n", ans(n).x, ans(n).y)
Next
Sleep
 