Convex hull, quick algorithm

[Marcus]

I need to calculate and display the convex hull of 2D polygons in the editor for that 3d thing, so I implemented a quick (as in "divide and conquer") algorithm for it. I remember writing a quick algorithm for the convex hull of a set of random points in 3D at the university, which is a bit trickier.

If you ever need to calculate the convex hull for a polygon or a bunch of points, you can use this

Code:

' Convex hull test
' ----------------

set window "test", 640, 480
set redraw off

do
    set color 0, 0, 0
    cls
    ' Create 100 random points.
    points = []
    for i = 1 to 100
        a = rnd()*2*PI
        r = rnd(200)
        points[sizeof(points)] = 320 + cos(a)*r
        points[sizeof(points)] = 240 + sin(a)*r
    next
    ' Draw points.
    set color 255, 255, 255
    pcount = sizeof(points)/2
    for i = 0 to pcount - 1
        draw pixel points[i*2], points[i*2 + 1]
    next
    ' Display and wait one second.
    redraw
    wait 1000
    ' Calculate convex hull and display it.
    hull = ConvexHull(points)
    draw poly hull
   
    set caret width(primary)/2, 4
    center "Press space bar to continue ..."
   
    redraw
    while not keydown(KEY_SPACE, true)  fwait 60
loop

' ConvexHull
' ----------
' Calculate the convex hull for a set of points in an array [x0, y0, x1, y1 ..], return as a polygon
' in the same format.
function ConvexHull(points)
    pcount = sizeof(points)/2
    xminIndex = 0; xmaxIndex = 0
    for i = 1 to pcount - 1
        x = points[i*2]
        if x < points[xminIndex*2] xminIndex = i
        if x > points[xmaxIndex*2] xmaxIndex = i
    next
    lpoints = []; rpoints = []
    for i = 0 to pcount - 1
        if i <> xminIndex and i <> xmaxIndex
            cp = CrossProd(points[xminIndex*2], points[xminIndex*2 + 1],
                    points[xmaxIndex*2], points[xmaxIndex*2 + 1],
                    points[i*2], points[i*2 + 1])
            if cp < 0
                lpoints[sizeof(lpoints)] = points[i*2]
                lpoints[sizeof(lpoints)] = points[i*2 + 1]
            elseif cp > 0
                rpoints[sizeof(rpoints)] = points[i*2]
                rpoints[sizeof(rpoints)] = points[i*2 + 1]
            endif
        endif
    next
    hull = []
    hull[sizeof(hull)] = points[xminIndex*2]
    hull[sizeof(hull)] = points[xminIndex*2 + 1]
    HullRec(hull, lpoints, points[xminIndex*2], points[xminIndex*2 + 1], points[xmaxIndex*2], points[xmaxIndex*2 + 1])
    hull[sizeof(hull)] = points[xmaxIndex*2]
    hull[sizeof(hull)] = points[xmaxIndex*2 + 1]
    HullRec(hull, rpoints, points[xmaxIndex*2], points[xmaxIndex*2 + 1], points[xminIndex*2], points[xminIndex*2 + 1])
   
    return hull
   
    function HullRec(hull, points, x0, y0, x1, y1)
        if not sizeof(points)  return
        pcount = sizeof(points)/2
        maxd = 0; maxdIndex = -1
        for i = 0 to pcount - 1
            x = points[i*2]; y = points[i*2 + 1]
            d = LinePointDist(x0, y0, x1, y1, x, y)
            if d > maxd
                maxd = d
                maxdIndex = i
            endif
        next
           
        npoints = []
        for i = 0 to pcount - 1
            if i <> maxdIndex
                x = points[i*2]; y = points[i*2 + 1]
                cp = CrossProd(x0, y0, points[maxdIndex*2], points[maxdIndex*2 + 1], x, y)
                if cp < 0
                    npoints[sizeof(npoints)] = x
                    npoints[sizeof(npoints)] = y
                endif
            endif
        next
        if sizeof(npoints)
            HullRec(hull, npoints, x0, y0, points[maxdIndex*2], points[maxdIndex*2 + 1])
        endif
        hull[sizeof(hull)] = points[maxdIndex*2]
        hull[sizeof(hull)] = points[maxdIndex*2 + 1]
       
        npoints = []
        for i = 0 to pcount - 1
            if i <> maxdIndex
                x = points[i*2]; y = points[i*2 + 1]
                cp = CrossProd(points[maxdIndex*2], points[maxdIndex*2 + 1], x1, y1, x, y)
                if cp < 0
                    npoints[sizeof(npoints)] = x
                    npoints[sizeof(npoints)] = y
                endif
            endif
        next
        if sizeof(npoints)
            HullRec(hull, npoints, points[maxdIndex*2], points[maxdIndex*2 + 1], x1, y1)
        endif
    endfunc
   
    function CrossProd(ax, ay, bx, by, cx, cy)
        return (bx - ax)*(cy - ay) - (by - ay)*(cx - ax)
    endfunc
   
    function LinePointDist(x0, y0, x1, y1, px, py)
        dx = x1 - x0; dy = y1 - y0; d = dx*dx + dy*dy
        if d <= 0  return 0
        return |(x1 - x0)*(py - y0) - (px - x0)*(y1 - y0)|/sqr(d)
    endfunc   
endfunc