Polygon Offset With pyeuclid Revisited

A few years back I did two or three posts on polygon offset.  It was a learning experience that I never quite completed to my satisfaction.  A kind visitor to my last post on the subject, Mr. Ahmad Rafsanjani, actually rewrote some of my code in a comment.  I gave him a polite weasel answer thanking him, but dropped the effort and never felt quite right about it.

Well, as the saying goes, better late than never.  He was quite correct in his assessment, but my understanding of vector math was not strong enough to prove this to myself.  I was visually inspecting the results, and, given what I was dealing with at the time, they seemed OK.

Here is the picture we're trying to get (this is with Mr. Rafsanjani's code, but the difference with mine and the original code, although wrong, is not great):

In order to nail down the discrepancy in my original code, I inserted some print statements with a lot of numeric precision (28 digits to the right of the decimal) in the output:

$ more points
1.2231671842700024832595318003      1.7024195134850139687898717966
2.1231671842700023944416898303      1.7024195134850139687898717966
2.2768328157299975167404681997      2.5475804865149860312101282034
1.6635803619063778135966913396      2.5493839809701555054743948858
1.7364196380936223196300716154      3.3506160190298444057077631442
2.5205825797292722434406186949      3.3529128986463621053815131745
2.6794174202707274901058553951      4.1470871013536383387076966756
2.1360193516544989655869812850      4.1228847880778562995374159073

(etc.)

The numbers highlighted in yellow are mismatches in the Y-coordinates of points of the inset offset polygon - each pair of Y coordinates should represent lines parallel to the X axis; in other words, they should be equal.  I have a bug.

Contrast that with the numbers yielded by Mr. Rafsanjani's code:

$ more points
1.2251864530113494300422871675      1.7000000000000001776356839400
2.1251864530113491191798402724      1.7000000000000001776356839400
2.2797319075568038826418160170      2.5499999999999998223643160600
1.6642549229616445671808833140      2.5499999999999998223643160600
1.7369821956889173186766583967      3.3500000000000000888178419700
2.5229705854077835169846366625      3.3500000000000000888178419700
2.6829705854077836590931838145      4.1500000000000003552713678801
2.1880983342360056376207921858      4.1500000000000003552713678801
2.6780983342360054066944030637      4.8499999999999996447286321199
3.1219016657639944156699129962      4.8499999999999996447286321199
(etc.)

Much better.  Lines that are supposed to be perfectly parallel to the X axis are, at least to 28 decimal places precision and the limits of my platform and the C Python interpreter, parallel to the X axis.  For what I am doing, I can more than live with that.

I've included Mr. Rafsanjani's comments in the code.  My modifications to his code were mainly for the purpose of printing some things out and organizing the polygon offset part of this exercise into a module.

I've made a separate main script for gnuplot.  After not looking at everything for three years I realized I had forgotten everything I ever knew about gnuplot and wanted to record it this time.  The file with the 20 points for the shape (monastery.py) is available on request.

Here is the main pyeuclid/polygon offset part of the code (rafsanjanicorrection.py):

"""
Polygon offset problem using
pyeuclid and incorporating corrections
made by Ahmed Rafsanjani.
"""

# Mr. Rafsanjani's comments:

# 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:

import euclid as eu
import math
import copy

import monastery as pic

OFFSET = 0.15

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
    v3.normalize()
   
    cs = v3.dot(v2)
   
    a1 = cs * v2
    a2 = v3 - a1
   
    if cs > 0:
        alpha = math.sqrt(a2.magnitude_squared())
    else:
        alpha =- math.sqrt(a2.magnitude_squared())
   
    if v1.z < 0.0:
        return scaleadd(origin, -1.0 * OFFSET/alpha, v3)
    else:
        return scaleadd(origin, OFFSET/alpha, v3)

def generatepoints():
    """
    Create list of offset points
    (pyeuclid.Vector3 objects) for
    points inset from polygon.

    Return list.
    """
    polyinset = []
    lenpolygon = len(pic.MONASTERY)
    i = 0
    poly = pic.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]))

    return polyinset

The file that prints stuff out and summons gnuplot (writtofile.py):

"""
Write vector points to file.

Show in gnuplot.
"""

# import blogpost as vecx
import rafsanjanicorrection as vecx
import os

# We're using gnuplot.
# It doesn't like commas, so
# we'll use whitespace (6).
FMT = '{0:30.28f}      {1:30.28f}'
FILEX = 'points'
ORIGSHAPE = 'originalshape'

PLOTCMD = 'set xrange[0.0:6.0]\n'
PLOTCMD += 'set yrange[0.0:6.0]\n'
PLOTCMD += 'plot "{0:s}" with lines lt rgb "red" lw 4, '
PLOTCMD += '"{1:s}" with lines lt rgb "blue" lw 4'
GNUPLOTFILE = 'plotfile'
GNUPLOT = 'gnuplot -p {:s}'.format(GNUPLOTFILE)

pts = vecx.generatepoints()
f = open(FILEX, 'w')
i = 1
for ptx in pts:
    print('Printing point {0:d} . . .'.format(i))
    print >> f, FMT.format(ptx.x, ptx.y)
    i += 1
f.close()
# Plot original as well.
# XXX - repetetive - make function.
i = 0
f = open(ORIGSHAPE, 'w')
for ptx in vecx.pic.MONASTERY:
    print('Printing point {0:d} of original shape . . .'.format(i))
    print >> f, FMT.format(ptx[0], ptx[1])
    i += 1
f.close()

f = open(GNUPLOTFILE, 'w')
print >> f, PLOTCMD.format(ORIGSHAPE, FILEX)
f.close()
os.system(GNUPLOT)

pyeuclid, to the best of my knowledge, runs only in Python 2.7 at the moment.  In any case, I got an error on the Python 3.4 install with setup.py so I stuck with 2.7.

Thanks to Mr. Rafsanjani for his help with this and for the rest of you for stopping by.