Polygon Offset Using Vector Math in IronPython

The other day I saw a something retweeted by @leppie (I think) about an experimental hyper-fast vector math driven 3D engine for the dot Net Framework.  This led me to investigate whether there is a default implementation of vector math in the dot Net Framework.  As it turns out, there is.

This is of interest because (I think) this would make IronPython the only Python implementation that has vector math included without having to install a third party library.  Java has a utils.Vector object, but it has nothing to do with vector math (it's a specialized array).  You do need to use the dot Net Framework instead of standard Python modules, but if you're running IronPython, you should have access to that anyway.

The whole, or at least a big part of the idea of running a Python implementation against the dotNet Framework is that you can leverage the power of that big library collection with a language that's fairly dense, easy, and doesn't require compilation.

This was pretty easy on Windows.  The only confusing part is that there are two namespaces in dot Net called System.Windows.  You want the one that references the WindowsBase dll.  This is the one that has our Vector object in it.

The code (including the plotting by Gnuplot - I had to download the Windows version; I did leave out the monastery.py file with the original shape points in it; also, the writetofile.py file is almost exactly like the one from the previous post except that for a Vector object, the x and y names are capitalized):

# vecipy.py

"""
Polygon offset problem using
dot Net Framework.
"""

import clr

WINX = 'WindowsBase'

clr.AddReference(WINX)

from System.Windows import Vector

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 Vector object representing a point
    offset from the origin.

    (Multiply vectorx by offset and add to origin.)
    """
    # Multiply method that takes scalar and Vector.
    multx = Vector.Multiply(vectorx, offset)
    return Vector.Add(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 Vector object.
    """
    origin = Vector(*pt2)
    v1 = Vector(pt1[0] - pt2[0], pt1[1] - pt2[1])
    v1.Normalize()
   
    v2 = Vector(pt3[0] - pt2[0], pt3[1] - pt2[1])
    v2.Normalize()
   
    v3 = copy.copy(v1)

    v1 = Vector.CrossProduct(v1, v2)

    v3 = Vector.Add(v3, v2)
    v3.Normalize()
   
    # In dotNet - Vector.Multiply is overloaded.
    # When it gets two Vector objects as arguments
    # it returns a dot product.
    cs = Vector.Multiply(v3, v2)
   
    # Again multiplication is overloaded.
    # Here it gets a scalar and a Vector
    # as arguments.
    a1 = Vector.Multiply(cs, v2)
    a2 = Vector.Subtraction(v3, a1)
   
    if cs > 0:
        alpha = math.sqrt(a2.LengthSquared)
    else:
        alpha =- math.sqrt(a2.LengthSquared)
   
    if v1 < 0.0:
        return scaleadd(origin, -1.0 * OFFSET/alpha, v3)
    else:
        return scaleadd(origin, OFFSET/alpha, v3)

def generatepoints():
    """
    Create list of offset points
    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


# writetofile.py


"""
Write vector points to file.

Show in gnuplot.
"""

import vecipy 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.
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)
    i += 1
f.close()

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


The result (shown in previous post):

I run OpenBSD on my laptop at home.  So I would be using mono in my cross-platform experiment.

Microsoft just recently (Fall 2014) announced the open sourcing of the dotNet Framework and cross platform capability for it.  The mono project responded very positively to this announcement.  I would imagine this as being good news for IronPython too.

OpenBSD has a package for mono.  From there, I just needed to download the IronPython binaries and run mono against them, or so I thought . . .

As it turns out, my script kept crashing on the overloaded Vector.Multiply method - NotImplementedError.  I tried to research things, wasn't having any luck, and brute forced the problem by wrapping the method in a class in C# class I called vecx:

Note (26NOV2014):  I hacked this C# module up a bit too quickly and didn't have performance or elegance in mind.  If you declare those Multiply methods as static you can save yourself the trouble of instantiating a new instance of the class each time you want to call them.  In fact, you can do the same thing with all the Vector methods you want to use (Add, CrossProduct, etc.).  I was just too hurried and too lazy.  CBT

using System;

public class vecx
{

  public System.Windows.Vector vectorx;

  public vecx()
  {
    System.Windows.Vector vectorx = new System.Windows.Vector(0.0, 0.0);
    this.vectorx = vectorx;
  }

  public vecx(double x, double y)
  {
    System.Windows.Vector vectorx = new System.Windows.Vector(x, y);
    this.vectorx = vectorx;
  }

  public Double Multiply(System.Windows.Vector a, System.Windows.Vector b)
  {
    return System.Windows.Vector.Multiply(a, b);
  }

  public System.Windows.Vector Multiply(Double a, System.Windows.Vector b)
  {
    return System.Windows.Vector.Multiply(a, b);
  }

  public System.Windows.Vector Multiply(System.Windows.Vector a, Double b)
  {
    return System.Windows.Vector.Multiply(a, b);
  }

}



The command line (your paths will probably be different) text for compiling this under mono was:

$ mcs -r:/usr/local/lib/mono/4.5/WindowsBase.dll -target:library vecx.cs

The code using this faux Vector class was a little bit different (and hackish):

"""
Polygon offset problem using
dot Net Framework.

Modified for use with mono.
"""

import clr

# Hacked C# module.
VECX = '/home/carl/vectormath/IronPython/mono/vecx.dll'

clr.AddReference(VECX)

import vecx

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 Vector object representing a point
    offset from the origin.

    (Multiply vectorx by offset and add to origin.)
    """
    # Generic vector for use of Vector type.
    vecgeneric = vecx().vectorx

    # Multiply method that takes scalar and Vector.
    # Using cs module compiled to dll for Multiply
    # methods in mono.
    multx = vecx().Multiply(vectorx, offset)
    return vecgeneric.Add(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 Vector object.
    """
    # Generic vector for use of type.
    vecgeneric = vecx().vectorx

    origin = vecx(*pt2).vectorx
    v1 = vecx(pt1[0] - pt2[0], pt1[1] - pt2[1]).vectorx
    v1.Normalize()
   
    v2 = vecx(pt3[0] - pt2[0], pt3[1] - pt2[1]).vectorx
    v2.Normalize()
   
    v3 = copy.copy(v1)

    v1 = vecgeneric.CrossProduct(v1, v2)

    v3 = vecgeneric.Add(v3, v2)
    v3.Normalize()
   
    # In dotNet - Vector.Multiply is overloaded.
    # When it gets two Vector objects as arguments
    # it returns a dot product.
    # Using cs module compiled to dll for Multiply
    # methods in mono.
    cs = vecx().Multiply(v3, v2)
   
    # Again multiplication is overloaded.
    # Here it gets a scalar and a Vector
    # as arguments.
    # Using cs module compiled to dll for Multiply
    # methods in mono.
    a1 = vecx().Multiply(cs, v2)
    a2 = vecgeneric.Subtract(v3, a1)
   
    if cs > 0:
        alpha = math.sqrt(a2.LengthSquared)
    else:
        alpha =- math.sqrt(a2.LengthSquared)
   
    if v1 < 0.0:
        return scaleadd(origin, -1.0 * OFFSET/alpha, v3)
    else:
        return scaleadd(origin, OFFSET/alpha, v3)

def generatepoints():
    """
    Create list of offset points
    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

Any port in a storm or whatever it takes, as they say.

Thanks again to Mr. Rafsanjani whom I referenced in my previous post.  His methodology and detection of a former bug got me back on track.

And thank you for stopping by.

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.