0

As above I am trying to implement a Graham Scan Convex Hull algorithm but I am having trouble with the stack appending too many vertices. The points are read from a .dat file here

For reading points my function is as follows:

def readDataPts(filename, N):
    """Reads the first N lines of data from the input file
          and returns a list of N tuples
          [(x0,y0), (x1, y1), ...]
    """
    count = 0
    points = []
    listPts = open(filename,"r")
    lines = listPts.readlines()
    for line in lines:
        if count < N:
            point_list = line.split()
            count += 1
            for i in range(0,len(point_list)-1):
                points.append((float(point_list[i]),float(point_list[i+1])))

    return points

My graham scan is as follows:

def theta(pointA,pointB):
    dx = pointB[0] - pointA[0]
    dy = pointB[1] - pointA[1]
    if abs(dx) < 0.1**5 and abs(dy) < 0.1**5:
        t = 0
    else:
        t = dy/(abs(dx) + abs(dy))
    if dx < 0:
        t = 2 - t
    elif dy < 0:
        t = 4 + t
    return t*90



def grahamscan(listPts):
    """Returns the convex hull vertices computed using the
         Graham-scan algorithm as a list of 'h' tuples
         [(u0,v0), (u1,v1), ...]  

    """
    listPts.sort(key=lambda x: x[1]) 
    p0 = listPts[0]
    angles = []
    for each in listPts:
        angle = theta(p0,each)
        angles.append((angle,each))
    angles.sort(key=lambda angle: angle[0])
    stack = []
    for i in range(0,3):
        stack.append(angles[i][1])
    for i in range(3, len(angles)):
        while not (isCCW(stack[-2],stack[-1],angles[i][1])):
            stack.pop()
        stack.append(angles[i][1])
    merge_error = stack[-1]
    #stack.remove(merge_error)
    #stack.insert(0,merge_error)
    return stack #stack becomes track of convex hull

def lineFn(ptA, ptB, ptC):
    """Given three points, the function finds the value which could be used to determine which sides the third point lies"""
    val1 = (ptB[0]-ptA[0])*(ptC[1]-ptA[1])
    val2 = (ptB[1]-ptA[1])*(ptC[0]-ptA[0])
    ans = val1 - val2
    return ans 

def isCCW(ptA, ptB, ptC):
    """Return True if the third point is on the left side of the line from ptA to ptB and False otherwise"""    
    ans = lineFn(ptA, ptB, ptC) > 0
    return ans

When I run it using the data set taking the first 50 lines as input it produces the stack:

[(599.4, 400.8), (599.0, 514.4), (594.5, 583.9), (550.1, 598.5), (463.3, 597.2), (409.2, 572.5), (406.0, 425.9), (407.3, 410.2), (416.3, 405.3), (485.2, 400.9)]

but it should produce (in this order):

[(599.4, 400.8), (594.5, 583.9), (550.1, 598.5), (472.6, 596.1), (454.2, 589.4), (410.8, 564.2), (416.3, 405.3), (487.7, 401.5)]

Any ideas?

1 Answers1

0

Angle sorting should be made against some extremal reference point (for example, the most bottom and left), that is undoubtedly included in convex hull. But your implementation uses the first point of list as reference.

Wiki excerpt:

swap points[1] with the point with the lowest y-coordinate

MBo
  • 77,366
  • 5
  • 53
  • 86
  • I normally have this `listPts.sort(key=lambda x: x[1])` in the code as my first line to make the first vertex the one with the minimum y-value. –  May 29 '17 at 05:31
  • OK. And comparing equal angles should account for distance. And are you sure that your pseudo angle function properly compares angles in all cases? – MBo May 29 '17 at 05:42
  • I've been told there are no coincident points within the file so the angles should never be the same. –  May 29 '17 at 06:24
  • A(0,0), B(1,1),C(2,2) are not coincident, but directions AB=AC – MBo May 29 '17 at 06:38
  • I don't think there are any 0 vectors either. –  May 29 '17 at 07:04
  • This is just example. Points are not coincident, they are collinear .Your isCCW function ignores zero (collinear) result (and perhaps initial sorting too) – MBo May 29 '17 at 07:47
  • Ok So if they are collinear it is neither clockwise or anticlockwise as all the points lie on the same straight line? –  May 29 '17 at 08:09
  • so a line like `if ans == 0:` ? –  May 29 '17 at 08:24
  • What should happen when it does equal 0? –  May 29 '17 at 08:33
  • I can not test your implementation, but in general all collinear points might stay in convex hull. Note that you can debug code on small data set. For example, 9 points - square corners, middles of edges and center of square. Center must not be in hull. Middles could - if implementation does not remove them. – MBo May 29 '17 at 08:45