05-02-2024, 04:24 PM
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
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
