I've just ported this over from bmax because I needed it for something, so here you go:
This code will find the convex hull of a set of 2d points. The convex hull is the smallest convex polygon which contains all the given points.
The function quickhull returns a list of points which are the vertices of the hull, in clockwise order.
Class Point
Field x#,y#
Method New(_x#,_y#)
x=_x
y=_y
End
End
'returns True if p1 and p2 are on the same side of the line a->b
Function sameside(p1x#,p1y#,p2x#,p2y#,ax#,ay#,bx#,by#)
Local cp1# = (bx-ax)*(p1y-ay)-(p1x-ax)*(by-ay)
Local cp2# = (bx-ax)*(p2y-ay)-(p2x-ax)*(by-ay)
If cp1*cp2 >= 0 Then Return True
End Function
'Clever little trick for telling if a point is inside a given triangle
'If for each pair of points AB in the triangle, P is on the same side of AB as
'the other point in the triangle, then P is in the triangle.
Function pointintriangle(px#,py#,ax#,ay#,bx#,by#,cx#,cy#)
Return (sameside(px,py,ax,ay,bx,by,cx,cy) And sameside(px,py,bx,by,ax,ay,cx,cy) And sameside(px,py,cx,cy,ax,ay,bx,by))
End Function
'Quickhull function - call this one with a set of points.
Function quickhull:List<Point>(s:List<Point>)
If s.Count()<=3 Return s
Local l:Point=Null
Local r:Point=Null
Local p:Point
'get left and right bounds
For p=Eachin s
If l=Null
l=p
Elseif p.x<l.x
l=p
Endif
If r=Null
r=p
Elseif p.x>r.x
r=p
Endif
Next
Local an#=ATan2(r.y-l.y,r.x-l.x)
Local rx#=Cos(an)
Local ry#=Sin(an)
Local sx#=Cos(an+90)
Local sy#=Sin(an+90)
Local s1:=New List<Point>
Local s2:=New List<Point>
For p=Eachin s
If p<>l And p<>r
Local mu#=(l.y-p.y+(ry/rx)*(p.x-l.x))/(sy-sx*ry/rx)
If mu<0
s1.AddLast p
Elseif mu>0
s2.AddLast p
Endif
Endif
Next
Local out1:=findhull(s1,l,r)
Local out2:=findhull(s2,r,l)
Local out:= New List<Point>
out.AddLast l
If out1
For p=Eachin out1
out.AddLast p
Next
Endif
out.AddLast r
If out2
For p=Eachin out2
out.AddLast p
Next
Endif
Return out
End Function
'Findhull helper function - you never need to call this
Function findhull:List<Point>(sk:List<Point>,p:Point,q:Point)
If sk.Count()=0 Return Null
Local c:Point=Null
Local out:=New List<Point>
Local maxdist#=-1
Local an#=ATan2(q.y-p.y,q.x-p.x)
Local rx#=Cos(an)
Local ry#=Sin(an)
Local sx#=-ry
Local sy#=rx
Local mu#
Local tp:Point
For tp=Eachin sk
If tp<>p And tp<>q
mu=Abs((p.y-tp.y+(ry/rx)*(tp.x-p.x))/(sy-sx*ry/rx))
If maxdist=-1 Or mu>maxdist
c=tp
maxdist=mu
Endif
Endif
Next
an=ATan2(c.y-p.y,c.x-p.x)
rx=Cos(an)
ry=Sin(an)
sx=Cos(an+90)
sy=Sin(an+90)
Local s1:=New List<Point>
Local s2:=New List<Point>
For tp=Eachin sk
If tp<>c
If Not pointintriangle(tp.x,tp.y,p.x,p.y,q.x,q.y,c.x,c.y)
mu=(p.y-tp.y+(ry/rx)*(tp.x-p.x))/(sy-sx*ry/rx)
If mu<0
s1.AddLast tp
Elseif mu>0
s2.AddLast tp
Endif
Endif
Endif
Next
Local out1:=findhull(s1,p,c)
Local out2:=findhull(s2,c,q)
If out1
For p=Eachin out1
out.AddLast p
Next
Endif
out.AddLast c
If out2
For p=Eachin out2
out.AddLast p
Next
Endif
Return out
End Function
|
