Sunday, July 3, 2011

pyeuclid, vector math, and polygon offset

UPDATE:  My code here has a bug.  Mr. Rafsanjani corrected it in the comments.  Please refer to them below.  Apologies - this blog is for beginners and for learning, so I've left the original code and the comment as is.  CBT 21NOV2014

I had previously done a series of posts on polygon offset intended as a practical guide for accomplishing the task quickly.  Some kind folks pointed out that I was making things considerably harder than necessary by using trigonometric functions when vector math would be easier and less error prone.

A coworker lent me his Java code that does polygon offset.  I translated it into python (2.5) using the pyeuclid module:

import euclid as eu
import copy

OFFSET = 0.15

# coordinates
             # PT 1
MONASTERY = [(1.1, 0.75),
             # PT 2
             (1.2, 1.95),
             . . .
             # PT 21
             (1.1, 0.75)]

def scaleadd(origin, offset, vectorx):
    """
    From a vector representing the origin,
    a scalar offset, and a vector, returns
    a Vector3 object representing a point 

    offset from the origin.

    (Multiply vectorx by offset and add to origin.)
    """
    multx = vectorx * offset
    return multx + origin

def getinsetpoint(pt1, pt2, pt3):
    """
    Given three points that form a corner (pt1, pt2, pt3),
    returns a point offset distance OFFSET to the right
    of the path formed by pt1-pt2-pt3.

    pt1, pt2, and pt3 are two tuples.

    Returns a Vector3 object.
    """
    origin = eu.Vector3(pt2[0], pt2[1], 0.0)
    v1 = eu.Vector3(pt1[0] - pt2[0], 

                    pt1[1] - pt2[1], 0.0)
    v1.normalize()
    v2 = eu.Vector3(pt3[0] - pt2[0], 

                    pt3[1] - pt2[1], 0.0)
    v2.normalize()
    v3 = copy.copy(v1)
    v1 = v1.cross(v2)
    v3 += v2
    if v1.z < 0.0:
        retval = scaleadd(origin, -OFFSET, v3)
    else:
        retval = scaleadd(origin, OFFSET, v3)
    return retval

polyinset = []
lenpolygon = len(MONASTERY)
i = 0
poly = MONASTERY
while i < lenpolygon - 2:
    polyinset.append(getinsetpoint(poly[i], 

                 poly[i + 1], poly[i + 2]))
    i += 1
polyinset.append(getinsetpoint(poly[-2], 

             poly[0], poly[1]))
polyinset.append(getinsetpoint(poly[0], 

             poly[1], poly[2]))

The result:

7 comments:

  1. Hey thanks heaps for writing this up. You should have seen the crazy quadrant atan2 path I was going down!

    ReplyDelete
  2. Could you link in your post the Java code your coworker passed you on?

    ReplyDelete
  3. Miguel,

    I'm sorry. I no longer have access to the Java code.

    Carl T.

    ReplyDelete
  4. Dear Carl,

    Many thanks for this nice code.
    I think there is a small bug:

    In "getinsetpoint", the vector v3 should be normalized before passing to "scaleadd".

    Furthermore, the final offset is not as the prescribed OFFSET and the angle between vectors should be taken into account.

    A possible solution could be:

    def getinsetpoint(pt1, pt2, pt3,offset):
    """
    Given three points that form a corner (pt1, pt2, pt3),
    returns a point offset distance OFFSET to the right
    of the path formed by pt1-pt2-pt3.

    pt1, pt2, and pt3 are two tuples.

    Returns a Vector3 object.
    """
    origin = eu.Vector3(pt2[0], pt2[1], 0.0)
    v1 = eu.Vector3(pt1[0] - pt2[0],
    pt1[1] - pt2[1], 0.0)
    v1.normalize()

    v2 = eu.Vector3(pt3[0] - pt2[0],
    pt3[1] - pt2[1], 0.0)
    v2.normalize()

    v3 = copy.copy(v1)
    v1 = v1.cross(v2)
    v3 += v2
    v3.normalize()

    cs=v3.dot(v2)

    a1=cs*v2
    a2=v3-a1

    if cs>0:
    alpha=sqrt(a2.magnitude_squared())
    else:
    alpha=-sqrt(a2.magnitude_squared())

    if v1.z < 0.0:
    retval = scaleadd(origin, -offset/alpha, v3)
    else:
    retval = scaleadd(origin, offset/alpha, v3)
    return retval

    ReplyDelete
  5. @Ahmad Thank you for stopping by and for reviewing the code. I don't have gnuplot handy at the moment to see what would happen if I instituted the changes you suggest. I think what I did originally is just inspect the plot and felt that it looked OK. What I need to do is measure the actual offset in the code and report it. I probably won't get to this immediately, but I intend to take a second look at this.

    Again, thank you for the rigor with which you reviewed this code. Carl T.

    ReplyDelete
  6. Mr. Rafsanjani,

    I finally got around to updating this with some data that illustrates the problem you pointed out:

    http://pyright.blogspot.com/2014/11/polygon-offset-with-pyeuclid-revisited.html

    Thanks so much for your help. Your attention to detail and willingness to point me in the right direction are appreciated. CBT

    ReplyDelete