Convex hull, quick algorithm - Printable Version

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