Sunday, March 27, 2011

Polygon Offset Revisited

I had previously written a post concerning polygon offset that was incomplete in its scope.  Specifically, it dodged the problem of horizontal and vertical lines in polygons.

Since then I've had to confront this problem again.  Shown below is the code I'm using for the more complete solution.  It's a bit verbose, but it works.

import math

def calcoffsetpoint(pt1, pt2, offset):
    """
    Get a point offset from the line
    segment pt1-pt2 distance "offset".

    Return two tuple of coordinates.
    """
    theta = math.atan2(pt2[1] - pt1[1],
                       pt2[0] - pt1[0])
    theta += math.pi/2.0
    return (pt1[0] - math.cos(theta) * offset,
            pt1[1] - math.sin(theta) * offset)

def getoffsetintercept(pt1, pt2, m, offset):
    """
    From points pt1 and pt2 defining a line
    in the Cartesian plane, the slope of the
    line m, and an offset distance,
    calculates the y intercept of
    the new line offset from the original.
    """
    x, y = calcoffsetpoint(pt1, pt2, offset)
    return y - m * x

def getpt(pt1, pt2, pt3, offset):
    """
    Gets intersection point of the two
    lines defined by pt1, pt2, and pt3;
    offset is the distance to offset
    the point from the polygon.

    Valid for lines with slopes other
    than zero or infinity.

    Returns a two tuple of coordinates.
    """
    # get first offset intercept
    m = (pt2[1] - pt1[1])/(pt2[0] - pt1[0])
    boffset = getoffsetintercept(pt1, pt2, m, offset)
    # get second offset intercept
    mprime = (pt3[1] - pt2[1])/(pt3[0] - pt2[0])
    boffsetprime = getoffsetintercept(pt2, pt3, mprime, offset)
    # get intersection of two offset lines
    newx = (boffsetprime - boffset)/(m - mprime)
    newy = m * newx + boffset
    return newx, newy

def getslopeandintercept(pt1, pt2, offset):
    """
    Gets the slope and the intercept of the
    offset line.
    Result returned as a two tuple.
    """
    m = (pt2[1] - pt1[1])/(pt2[0] - pt1[0])
    b = getoffsetintercept(pt1, pt2, m, offset)
    return m, b

def getoffsetcornerpoint(pt1, pt2, pt3, offset):
    """
    Gets intersection point of the two
    lines defined by pt1, pt2, and pt3;
    offset is the distance to offset
    the point from the polygon.

    Returns a two tuple of coordinates.
    """
    # starting out with horizontal line
    if (pt2[1] - pt1[1]) == 0.0:
        ycoord = pt1[1] - math.cos(math.atan2(0.0, pt2[0] - pt1[0])) * offset
        # a vertical line follows
        if (pt3[0] - pt2[0]) == 0.0:
            xcoord = pt2[0] + math.sin(math.atan2(pt3[1] - pt2[1], 0.0)) * offset
        # a sloped line follows
        else:
            m, offsetintercept = getslopeandintercept(pt2, pt3, offset)
            # calculate for x with ycoord
            xcoord = (ycoord - offsetintercept)/m
    # starting out with a vertical line
    if (pt2[0] - pt1[0]) == 0.0:
        xcoord = pt1[0] + math.sin(math.atan2(pt2[1] - pt1[1], 0.0)) * offset
        # a horizontal line follows
        if (pt3[1] - pt2[1]) == 0.0:
            ycoord = pt2[1] - math.cos(math.atan2(0.0, pt3[0] - pt2[0])) * offset
        # a sloped line follows
        else:
            m, offsetintercept = getslopeandintercept(pt2, pt3, offset)
            # calculate for y with xcoord
            ycoord = m * xcoord + offsetintercept
    # starting out with sloped line
    if (pt2[1] - pt1[1]) != 0.0 and (pt2[0] - pt1[0]) != 0.0:
        # if second line is horizontal
        if (pt3[1] - pt2[1]) == 0.0:
            ycoord = pt2[1] - math.cos(math.atan2(0.0, pt3[0] - pt2[0])) * offset
            m, offsetintercept = getslopeandintercept(pt1, pt2, offset)
            # calculate for x with y coord
            xcoord = (ycoord - offsetintercept)/m
        # if second line is vertical
        elif (pt3[0] - pt2[0]) == 0.0:
            xcoord = pt2[0] + math.sin(math.atan2(pt3[1] - pt2[1], 0.0)) * offset
            m, offsetintercept = getslopeandintercept(pt1, pt2, offset)
            # solve for y with x coordinate
            ycoord = m * xcoord + offsetintercept
        # if both lines are sloped
        else:
            xcoord, ycoord = getpt(pt1, pt2, pt3, offset)
    return xcoord, ycoord

def offsetpolygon(polyx, offset):
    """
    Offsets a clockwise list of coordinates
    polyx distance offset to the inside of
    the polygon.
    Returns list of offset points.
    """
    polyy = []
    # need three points at a time
    for counter in range(0, len(polyx) - 3):
        # get first offset intercept
        pt = getoffsetcornerpoint(polyx[counter],
                   polyx[counter + 1],
                   polyx[counter + 2],
                   offset)
        # append new point to polyy
        polyy.append(pt)
    # last three points
    pt = getoffsetcornerpoint(polyx[-3], polyx[-2], polyx[-1], offset)
    polyy.append(pt)
    pt = getoffsetcornerpoint(polyx[-2], polyx[-1], polyx[0], offset)
    polyy.append(pt)
    pt = getoffsetcornerpoint(polyx[-1], polyx[0], polyx[1], offset)
    polyy.append(pt)
    return polyy
 

I used this polygon to test it.


This is a fairly simple geometric problem, but it can be confusing.  Hopefully this will save someone the trouble of working through a solution.

Notes:

1) If you're using Python 2.x, feed your coordinates in as floats.  This will ensure the absence of integer division, which often yields zeros, which in turn yields ZeroDivision errors.

2) The polygon fed to the offsetpolygon function should not have a closing point.  The closing point will cause ZeroDivision errors.  Instead, append the point before drawing the polygon in gnuplot or some other plotting tool.

3 comments:

  1. Thanks for your ideas.
    Still have a devide by zero when putting a polygon point on a straight line (ok, bit exceptional case).

    I continued on your code and i am looking for code to have full polygon offset algorithm.
    Vertex events are solved, split events i am still searching for and working on...

    ReplyDelete
  2. This is a very interesting algorithm but as you mentioned it has some problems.
    A more robust one can be found here:

    http://pyright.blogspot.ch/2011/07/pyeuclid-vector-math-and-polygon-offset.html

    ReplyDelete
  3. many thanks i have converted it into Javascript and its working like anything.. thank you once again..

    ReplyDelete