Welcome, guest | Sign In | My Account | Store | Cart
#!/usr/bin/env python

"""convexhull.py

Calculate the convex hull of a set of n 2D-points in O(n log n) time.  
Taken from Berg et al., Computational Geometry, Springer-Verlag, 1997.
Prints output as EPS file.

When run from the command line it generates a random set of points
inside a square of given length and finds the convex hull for those,
printing the result as an EPS file.

Usage:

    convexhull.py <numPoints> <squareLength> <outFile>

Dinu C. Gherman
"""



import sys, string, random


######################################################################
# Helpers
######################################################################

def _myDet(p, q, r):
   
"""Calc. determinant of a special matrix with three 2D points.

    The sign, "
-" or "+", determines the side, right or left,
    respectivly, on which the point r lies, when measured against
    a directed vector from p to q.
    """


   
# We use Sarrus' Rule to calculate the determinant.
   
# (could also use the Numeric package...)
    sum1
= q[0]*r[1] + p[0]*q[1] + r[0]*p[1]
    sum2
= q[0]*p[1] + r[0]*q[1] + p[0]*r[1]

   
return sum1 - sum2


def _isRightTurn((p, q, r)):
   
"Do the vectors pq:qr form a right turn, or not?"

   
assert p != q and q != r and p != r
           
   
if _myDet(p, q, r) < 0:
       
return 1
   
else:
       
return 0


def _isPointInPolygon(r, P):
   
"Is point r inside a given polygon P?"

   
# We assume the polygon is a list of points, listed clockwise!
   
for i in xrange(len(P[:-1])):
        p
, q = P[i], P[i+1]
       
if not _isRightTurn((p, q, r)):
           
return 0 # Out!        

   
return 1 # It's within!


def _makeRandomData(numPoints=10, sqrLength=100, addCornerPoints=0):
   
"Generate a list of random points within a square."
   
   
# Fill a square with random points.
    min
, max = 0, sqrLength
    P
= []
   
for i in xrange(numPoints):
        rand
= random.randint
        x
= rand(min+1, max-1)
        y
= rand(min+1, max-1)
        P
.append((x, y))

   
# Add some "outmost" corner points.
   
if addCornerPoints != 0:
        P
= P + [(min, min), (max, max), (min, max), (max, min)]

   
return P


######################################################################
# Output
######################################################################

epsHeader
= """%%!PS-Adobe-2.0 EPSF-2.0
%%%%BoundingBox: %d %d %d %d

/r 2 def                %% radius

/circle                 %% circle, x, y, r --> -
{
    0 360 arc           %% draw circle
} def

/cross                  %% cross, x, y --> -
{
    0 360 arc           %% draw cross hair
} def

1 setlinewidth          %% thin line
newpath                 %% open page
0 setgray               %% black color

"""


def saveAsEps(P, H, boxSize, path):
   
"Save some points and their convex hull into an EPS file."
   
   
# Save header.
    f
= open(path, 'w')
    f
.write(epsHeader % (0, 0, boxSize, boxSize))

    format
= "%3d %3d"

   
# Save the convex hull as a connected path.
   
if H:
        f
.write("%s moveto\n" % format % H[0])
       
for p in H:
            f
.write("%s lineto\n" % format % p)
        f
.write("%s lineto\n" % format % H[0])
        f
.write("stroke\n\n")

   
# Save the whole list of points as individual dots.
   
for p in P:
        f
.write("%s r circle\n" % format % p)
        f
.write("stroke\n")
           
   
# Save footer.
    f
.write("\nshowpage\n")


######################################################################
# Public interface
######################################################################

def convexHull(P):
   
"Calculate the convex hull of a set of points."

   
# Get a local list copy of the points and sort them lexically.
    points
= map(None, P)
    points
.sort()

   
# Build upper half of the hull.
    upper
= [points[0], points[1]]
   
for p in points[2:]:
        upper
.append(p)
       
while len(upper) > 2 and not _isRightTurn(upper[-3:]):
           
del upper[-2]

   
# Build lower half of the hull.
    points
.reverse()
    lower
= [points[0], points[1]]
   
for p in points[2:]:
        lower
.append(p)
       
while len(lower) > 2 and not _isRightTurn(lower[-3:]):
           
del lower[-2]

   
# Remove duplicates.
   
del lower[0]
   
del lower[-1]

   
# Concatenate both halfs and return.
   
return tuple(upper + lower)


######################################################################
# Test
######################################################################

def test():
    a
= 200
    p
= _makeRandomData(30, a, 0)
    c
= convexHull(p)
    saveAsEps
(p, c, a, file)


######################################################################

if __name__ == '__main__':
   
try:
        numPoints
= string.atoi(sys.argv[1])
        squareLength
= string.atoi(sys.argv[2])
        path
= sys.argv[3]
   
except IndexError:
        numPoints
= 30
        squareLength
= 200
        path
= "sample.eps"

    p
= _makeRandomData(numPoints, squareLength, addCornerPoints=0)
    c
= convexHull(p)
    saveAsEps
(p, c, squareLength, path)

History