Tuesday, November 25, 2014

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.


No comments:

Post a Comment