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:

2D geometry + infinite and zero slopes

My last post dealt with a similar subject, but it didn't quite tell the whole story.  No sooner had I posted it and attempted to implement it in a project, I ran into a problem.  Granted, it was on a Java GUI in Java code (not a Python project) - but the concept is universal.  When is a floating point number small enough that another mathematical operation on it will take it to zero or to the realm of infinity?

Python (I'm using CPython 2.5) is pretty forgiving in this sense.

>>> 1.0/1.0e-308
1e+308
>>> 1.0/1.0e-309
inf


There comes a point when even a representable number equals zero.

>>> 0.0 == 1.0e-10000
True
 
With some trigonometric functions from the math module:


>>> math.atan2(1.0, 0.0000000000000001)
1.5707963267948966
>>> 2 * _
3.1415926535897931
>>> _ == math.pi
True


>>> math.cos(1.0000000000000000001e-250)
1.0


The equation for the slope of a line in the Cartesian plane is (y2 - y1)/(x2 - x1) where y2, y1, x2, and x1 are the x and y coordinates of the enpoints of the line.


One way I've handled the infinite slope problem in the past is by setting a value as infinity:

# one million

>>> INFINITY = 1.0e+6
>>> slope = 1000000.0
>>> slope == INFINITY
True
This works sometimes, but it's not totally safe from error.


What I ended up doing this time was checking for zero or near zero values in the numerator and denominator of the slope equation.


NEARZERO = 1.0e-6
if math.abs(y2 - y1) < NEARZERO:
    # set y2 == y1 or vice versa
    # and proceed with the calculations
    # for a zero slope
if math.abs(x2 - x1) < NEARZERO:
    # set x2 == x1 or vice versa
    # and proceed with the calculations
    # for an infinite slope

This solved my problem in the sense that it gave reasonable results that served the problem domain well.

Lastly, you could use a similar method to check for duplicate points.  If both numerator and denominator of the slope equation are near zero, the line segment is too small and the point should be done away with:


if (math.abs(y2 - y1) < NEARZERO and
    math.abs(x2 - x1) < NEARZERO):
    # duplicate point
    # delete it


Disclaimer:  I am not a floating point math expert.  There may be better solutions out there.  These are a couple practical ones that have worked for me.




 

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.

jython + joda-time

joda-time is a date-time library written in Java.  It has advantages over Java's built in time utilities in terms of power and ease of use.  It has some functionality that Python's own datetime library does not.  Here is a sampling from within jython:


[OpenJDK Client VM (Sun Microsystems Inc.)] on java1.7.0-internal
Type "help", "copyright", "credits" or "license" for more information.
>>> import sys
>>> # make sure the jar is in the CLASSPATH
>>> sys.path.append('/home/carl/Downloads/joda-time-1.6.2/joda-time-1.6.2.jar')
>>> from org.joda import time as joda
>>> rightnow = joda.DateTime()
>>> rightnow
2011-03-24T19:54:45.544-06:00
>>> rightnow.getZone()
America/Edmonton
>>> # I'm actually in New Mexico - close enough
>>> rightnow.dayOfWeek().getAsText()
u'Thursday'
>>> rightnow.monthOfYear().getAsText()
u'March'
>>> rightnow.getDayOfMonth()
24
>>> # days are one indexed along with months

>>> # a date from the past
>>> declofindependence = joda.DateTime(1776, 7, 4, 0, 0, 0, 0)



The Period object allows you to count down to dates in the future of count from dates in the past in terms of years, months, and days:

>>> americasage = joda.Period(declofindependence, rightnow)
>>> americasage.getYears()
234
>>> americasage.getMonths()
8
>>> americasage.getWeeks()
2
>>> americasage.getDays()
6


joda-time can handle dates far into the future as well as those in the far distant past:

>>> farfuture = joda.DateTime(10000, 1, 1, 1, 0, 0, 0)
>>> farfuture.getYear()
10000
>>> farfuture
10000-01-01T01:00:00.000-07:00
>>>

Like most things Java, joda-time can be a bit more verbose than its Python equivalent.  Still, the ability to get dates beyond 9999 and some of the functionality may make it well worth the trouble.

WKT (Well Known Text) + JTS Topology + jython

The last few posts have dealt with the JTS Topology library.  This one will wrap up that topic with a method for converting JTS's Geometry objects to and from text, WKT.

WKT is a standard of the Open Geospatial Consortium (OGC).  Objects expressed in WKT are not much different than their string representations in JTS.  Below is a code example of working back and forth between JTS and WKT.


Jython 2.5.2 (Release_2_5_2:7206, Mar 2 2011, 23:12:06)
[OpenJDK Client VM (Sun Microsystems Inc.)] on java1.7.0-internal
Type "help", "copyright", "credits" or "license" for more information.
>>> import sys
>>> sys.path.append('/home/carl/Downloads/jts/lib/jts-1.11.jar')
>>> from com.vividsolutions.jts.geom import GeometryFactory
>>> from com.vividsolutions.jts.geom import Coordinate
>>> from com.vividsolutions.jts.geom import LinearRing
>>> from com.vividsolutions.jts.io import WKTWriter
>>> from com.vividsolutions.jts.io import WKTReader
>>> coord1 = Coordinate(0, 0)
>>> coord2 = Coordinate(0, 1)
>>> coord3 = Coordinate(1, 1)
>>> coord4 = Coordinate(1, 0)
>>> coord5 = Coordinate(0, 0)
>>> coords = [coord1, coord2, coord3, coord4, coord5]
>>> from com.vividsolutions.jts.geom.impl import CoordinateArraySequence
>>> cas = CoordinateArraySequence(coords)
>>> gf = GeometryFactory()
>>> lr = LinearRing(cas, gf)
>>> asquare = gf.createPolygon(lr, None)
>>> f = open('wkttest.xml', 'w')
>>> wktwriter = WKTWriter()
>>> # this is just my own personal preference
>>> # it is inefficient for storage, but
>>> # sometimes easier on the eyes
>>> wktwriter.setMaxCoordinatesPerLine(1)
>>> text = wktwriter.writeFormatted(asquare)
>>> f.write(text)
>>> f.close()
>>> print text
POLYGON ((0 0,
  0 1,
  1 1,
  1 0,
  0 0))
>>> wktreader = WKTReader()
>>> f = open('wkttest.xml', 'r')
>>> text = f.read()
>>> asquare2 = wktreader.read(text)
>>> asquare2
POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))


It's fairly easy to store JTS Geometry objects in human readable (WKT) format and reuse them.  

There is a good bit more to WKT than what's shown here.  This was just a simple example for getting started.


WKB (Well Known Binary) is also available in JTS - this is used for efficient storage of Geometry objects, mainly in databases.


JTS also has a WKTFileReader class which can read in multiple geometries from a file.

Polygon Buffering with JTS Topology and jython

Previously, I had done a post regarding polygon offset.  It turns out that JTS Topology has polygon buffering capability.  I wanted to try this out and compare it to my results from the polygon offset post.

As with most things with a library like JTS, the software does most of the work for you:

import sys
# path to jar into CLASSPATH
sys.path.append('/home/carl/Downloads/jts/lib/jts-1.11.jar')

from com.vividsolutions.jts.geom import Coordinate
from com.vividsolutions.jts.geom import GeometryFactory
# LinearRing is for the creation of the Polygon
from com.vividsolutions.jts.geom import LinearRing
# CoordinateSequenceArray is for the 
# creation of the LinearRing
from com.vividsolutions.jts.geom.impl import CoordinateArraySequence
# CAP_ROUND is for the Polygon.buffer operation
from com.vividsolutions.jts.operation.buffer.BufferOp import CAP_ROUND

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

gf = GeometryFactory()
coords = [Coordinate(*coord) for coord in MONASTERY]
cas = CoordinateArraySequence(coords)
lr = LinearRing(cas, gf)
polyx = gf.createPolygon(lr, None)

# offset of 0.15 inward
# rounded edges
# 10 points per quarter circle of rounded edges
polyinset = polyx.buffer(-0.15, 10, CAP_ROUND)

The result is shown below:

The rounded edges give an entirely different look to the shape. My version is shown below:

I'm glad I did the egg the way I did it (it looks more realistic).  Nonetheless, the polygon buffering capability of JTS is a useful tool, particularly for calculating offsets and distances for scientific or geographic purposes.


Notes:  Shapely is a Python package (CPython) with much of the same functionality as JTS.

Clipping Voronoi Polygons with JTS Topology and jython

Following up on the last post on Voronoi diagrams and the JTS topology library, I will show clipping of the Voronoi polygons with another polygon (a polygonal bounding region) in this post.

The plot below is what we have currently, a Voronoi diagram (with points) clipped along a rectangular boundary:

Shown below is what I'd like to accomplish - clipping the red polygons with the blue surrounding one:

Here is the code, continued from the last post:

# because JTS can handle holes and "donuts", it
# has a LinearRing structure for constructing
# polygons
from com.vividsolutions.jts.geom import LinearRing
# the LinearRing requires a CoordinateArraySequence
# (it will not accept a list)
from com.vividsolutions.jts.geom.impl import CoordinateArraySequence


boundingcoords = [(10, 0), (2, 0),
                  (-2, 4), (2, 12),
                  (10, 6), (10, 0)]
boundingcoords = [Coordinate(*coord) for coord in boundingcoords]
casboundingcoords = CoordinateArraySequence(boundingcoords)
lr = LinearRing(casboundingcoords, geomfactx)
# None here fills in the spot where any inner rings
# would be in the polygon
boundingpoly = geomfactx.createPolygon(lr, None)
 
newpolygons = []
numpolygons = diagramx.getNumGeometries()
for numx in range(numpolygons):
    newpolygons.append(diagramx.getGeometryN(
        numx).intersection(boundingpoly))

And the result:

Note:  it is possible to get more than one polygon as the result of each intersection.  Through visual inspection I determined there would only be one polygon per intersection with the bounding area and structured my list (newpolygons) accordingly.

jython, JTS Topology, and Voronoi Diagrams

Last time I barely scratched the surface of the JTS library.  As it turns out, through ignorance I made a few things harder than they needed to be.  The GeometryFactory class is quite versatile and has a number of methods for creating objects from simple Coordinate objects.  In the example below, I've managed to create a clipped Voronoi diagram without too much code or difficulty.

import sys
# this line will be different depending on where
# the jts jar is stored and may be unnecessary if
# the jar is in your CLASSPATH
sys.path.append('/home/carl/Downloads/jts/lib/jts-1.11.jar')

from com.vividsolutions.jts.geom import Coordinate
from com.vividsolutions.jts.geom import GeometryFactory
from com.vividsolutions.jts.triangulate import VoronoiDiagramBuilder
from com.vividsolutions.jts.geom import Envelope



# first some coordinates
coord1 = Coordinate(5, 7)
coord2 = Coordinate(3, 6)
coord3 = Coordinate(1, 4)
coord4 = Coordinate(2, 5)
coord5 = Coordinate(4, 2)
coord6 = Coordinate(7, 2)
coord7 = Coordinate(5, 5)
coord8 = Coordinate(2, 6)
coord5 = Coordinate(4, 2)
coords = [coord1, coord2,
          coord3, coord4,
          coord5, coord6,
          coord7, coord8]



# the handy GeometryFactory
geomfactx = GeometryFactory()

# a MultiPoint Geometry object
mp = geomfactx.createMultiPoint(coords)

# the diagram builder class
vdb = VoronoiDiagramBuilder()

# it would be nice to clip the
# diagram with a rectangle;
# this is what the Envelope
# object is for
lowerleft = Coordinate(0, 0)
upperright = Coordinate(8, 8)
env = Envelope(lowerleft, upperright)
vdb.setClipEnvelope(env)


# load our sites (points) from the
# MultiPoint object
vdb.setSites(mp)


# this is where the builder object clips the
# the polygons with the Envelope rectangle
diagramx = vdb.getDiagram(geomfactx)

print 'Number of polygons = %d' % diagramx.getNumGeometries()
polygonx = diagramx.getGeometryN(0)
print 'One polygon:'
for coordsetx in polygonx.getCoordinates():
    print coordsetx



Output:


Number of polygons = 8
One polygon:
(-5.0, -4.0, NaN)
(-5.0, 8.25, NaN)
(0.5, 5.5, NaN)
(2.7, 3.3, NaN)
(-2.166666666666666, -4.0, NaN)
(-5.0, -4.0, NaN)


There's more that can be done here (for instance, clipping the polygons with a non-rectangular shape).  As was true of the last post, this is just a small fraction of what the JTS Topology library can do.

jython + the JTS geometry library

My favorite library for dealing with two dimensional points, lines, and polygons is Polygon by Jörg Rädler.  In a Java/jython environment, though, that's not available.  This led me to try my hand at JTS.  

Greg Corradini has a post that highlights JTS use with jython and geoscript.  His overview was very helpful, but I wanted to get to know JTS more thoroughly.

The first thing I came up against when using JTS was the documentation.  It's all javadocs.  This isn't necessarily a bad thing, but, along with the library's java-centric design, it takes a little getting used to.  If concepts like abstract class, iterface, extending, and implementing don't immediately ring a bell to the user, use of the javadocs will make them clear over time.

That shouldn't discourage us, though.  Let's fire up the interpreter and try something simple.

Jython 2.5.2 (Release_2_5_2:7206, Mar 2 2011, 23:12:06)
[OpenJDK Client VM (Sun Microsystems Inc.)] on java1.7.0-internal
Type "help", "copyright", "credits" or "license" for more information.
>>> import sys
>>> sys.path.append('/home/carl/Downloads/jts/lib/jts-1.11.jar')

I usually store things I'm testing out in my Downloads folder (I'm using KDE on OpenBSD).  This may not be the best practice.  In any case, to import from jts I need access to the jar in my sys.path.

Making a Point object in JTS in a bit complicated.  There are Coordinate objects and there are Point objects.  The Coordinate is a more basic class with less functionality than the Point class.  Point inherits from the abstract class Geometry.  This leads it to have a number of methods, some meaningless, like convexHull.  Shown below is the process required to make a Point object.

>>> from com.vividsolutions.jts.geom import Point
>>> from com.vividsolutions.jts.geom import GeometryFactory
>>> from com.vividsolutions.jts.geom.impl import CoordinateArraySequence
>>> # get CoordinateArraySequence of size one
>>> cas = CoordinateArraySequence(1)
>>> factoryx = GeometryFactory()
>>> ptx = Point(cas, factoryx)

We have a Point object.  A couple explanatory notes:

1) the documentation for the default GeometryFactory reads thus:  "Constructs a GeometryFactory that generates Geometries having a floating PrecisionModel and a spatial-reference ID of 0."  For what I'm doing this is fine.

2) The preferred constructor for the Point object takes a CoordinateArraySequence of length one instead of a single Coordinate object.  This is analogous to a character in Python being a string of length one.

We don't know yet what our values are for our Point's coordinates.  Let's have a look.

>>> ptx
POINT (0 0)

OK, zero-zero.  Let's assign something to those coordinates away from the origin.

>>> ptx.getCoordinates()[0].x = 44
>>> ptx.getCoordinates()[0].y = 22
>>> ptx
POINT (44 22)

Even though the Point has only one coordinate, we still need to specify that with the index 0.

Well, that was a lot of work for just one point.  Let's try something more efficient and fun.

>>> cas2 = CoordinateArraySequence(6)
>>> somecoordinates = [(1, 5), (7, 14), (22, 44), (36, 12), (19, 1), (12, 4)]
>>> for numx in range(len(somecoordinates)):
...     cas2.getCoordinate(numx).x = somecoordinates[numx][0]
...     cas2.getCoordinate(numx).y = somecoordinates[numx][1]
...
>>> cas2
((1.0, 5.0, NaN), (7.0, 14.0, NaN), (22.0, 44.0, NaN), (36.0, 12.0, NaN), (19.0, 1.0, NaN), (12.0, 4.0, NaN))
>>> from com.vividsolutions.jts.geom import LineString
>>> ls = LineString(cas2, factoryx)
>>> ch = ls.convexHull()
>>> ch
POLYGON ((19 1, 1 5, 22 44, 36 12, 19 1))
>>> 

LineString is one of the classes that can be instantiated from a CoordinateArraySequence.  I went with that mainly for ease of use.

There is a great deal more to JTS.  This post was mainly for getting started with the library.

jython + Google Guava's base.CaseFormat

I was trying to research Google Collections for use in jython and stumbled upon Solid Craft's blog and his description of Guava base's CaseFormat utility.  Where I work we're a Java shop, and use of camelCase infiltrates all code, even jython code.  CaseFormat provides functionality for changing these variable names to a number of different formats.

$ /usr/local/jdk-1.7.0/bin/java -jar jython.jar
Jython 2.5.1 (Release_2_5_1:6813, Sep 26 2009, 13:47:54)
[OpenJDK Client VM (Sun Microsystems Inc.)] on java1.7.0-internal
Type "help", "copyright", "credits" or "license" for more information.
>>> import sys
>>> # get downloaded jar into sys.path
>>> sys.path.append('/home/carl/Downloads/guava-r08/guava-r08.jar')
>>> from com.google.common.base import CaseFormat as cf
>>> dir(cf)
['LOWER_CAMEL', 'LOWER_HYPHEN', 'LOWER_UNDERSCORE', 'UPPER_CAMEL', . . .]
>>> varname = 'someVariable'
>>> cf.LOWER_CAMEL.to(cf.LOWER_UNDERSCORE, varname)
u'some_variable'
>>> cf.LOWER_CAMEL.to(cf.UPPER_UNDERSCORE, varname)
u'SOME_VARIABLE'
>>> cf.LOWER_CAMEL.to(cf.UPPER_CAMEL, varname)
u'SomeVariable'
>>>


For cycling through a file full of camel case or underscored variable names, even by hand in the interpreter, this could save some time and typing.
 

jython + java.lang.Character.UnicodeBlock

In my post on Unicode Blocks in regular expressions, I mentioned there wasn't support for Unicode Blocks in CPython regular expressions.  Unicode Blocks are contiguous (by number) sections of the Unicode tables with some commonality among the characters.  There is another Java language feature that can be helpful for those interested in Unicode Blocks:  the java.lang.Character.UnicodeBlock object.  In jython:

Jython 2.5.1 (Release_2_5_1:6813, Sep 26 2009, 13:47:54)
[Java HotSpot(TM) 64-Bit Server VM (Sun Microsystems Inc.)] on java1.6.0_22
Type "help", "copyright", "credits" or "license" for more information.
>>> from java.lang.Character import UnicodeBlock
>>> UnicodeBlock.of(236)
LATIN_1_SUPPLEMENT
>>> UnicodeBlock.of('a')
BASIC_LATIN
>>> UnicodeBlock.of(738)
SPACING_MODIFIER_LETTERS
>>> UnicodeBlock.of(922)
GREEK
>>> UnicodeBlock.of(0xffee)
HALFWIDTH_AND_FULLWIDTH_FORMS

The "of" method can accept either a character or a Unicode numeric identifier as it's argument.  It provides a shorthand method of finding out roughly where a character is in Unicode and what it might represent.

More SVG - trace from bitmap with Inkscape and color editing

I'm attempting to convert some of the pysanky eggs I drew with POV-ray to SVG so that I can scale them freely.  Also, once they're converted, I would like to be able to edit their colors in the same way I did with the Python logo in my last post.

To the best of my knowledge, Inkscape (potrace) is the most readily available tool for converting the bitmap images from POV-ray to SVG.  There are a few things to look out for, but this is still really easy:

1) Remove all reflection and shading from the textures (colors) in the POV-ray scene.  Work on getting even lighting for the whole scene.  Generate the scene without antialiasing (jagged edges are fine, and even desired).

2) Follow the directions here for tracing the SVG drawing.

 
Now to use the standard library's xml.etree.ElementTree package to get at the colors.  Before jumping in, I first removed any references to the old bmp image from POV-ray and checked on the three color values in the SVG file.


The three colors in the image are #fefefe (slightly off white), #000000 (black), and #00fd00 (light green).


# svgchangecolor.py


"""
Manipulation of colors in an SVG image.
"""

from xml.etree import ElementTree as ET

def drilldown(nodex, oldcolor, newcolor):
    """
    Walk node of element tree in search
    of oldcolor.

    Replace oldcolor with newcolor where
    oldcolor occurs.
    """
    childrenx = nodex.getchildren()
    if childrenx:
        for childx in childrenx:
            drilldown(childx, oldcolor, newcolor)
    else:
        itemsx = nodex.items()
        for itemx in itemsx:
            if itemx[1].find(oldcolor) > -1:
                oldcolorstr = itemx[1]
                newcolorstr = oldcolorstr.replace(oldcolor, newcolor)
                nodex.set(itemx[0], newcolorstr)

def changecolor(oldcolor, newcolor, filename):
    """
    Replaces oldcolor with newcolor in an
    SVG file named filename.

    All three arguments are strings.

    oldcolor and newcolor are hex strings
    in the format #xxxxxx.
    """
    # the colors are at the lowest level
    # in the XML (SVG) file
    # drill down until they are reached
    svgobj = ET.parse(filename)
    rootx = svgobj.getroot()
    drilldown(rootx, oldcolor, newcolor)
    svgobj.write(filename)

$ python2.5
Python 2.5.4 (r254:67916, Aug  9 2010, 08:57:51)
[GCC 4.2.1 20070719 ] on openbsd4
Type "help", "copyright", "credits" or "license" for more information.
>>> import svgchangecolor as svgcc
>>> svgcc.changecolor('#00fd00', '#ffff00', 'crossdroppullegg.svg')
>>> svgcc.changecolor('#000000', '#00aa00', 'crossdroppullegg.svg')
>>> svgcc.changecolor('#fefefe', '#000000', 'crossdroppullegg.svg')

>>>

 

There are only three elements with color in the file (the dark green one is doing all the "work" of filling in the pattern).  The script is really overkill in this case.  For a drawing with a number of separate paths with the same color, though, this could be very useful.

SVG, XML, and the Python logo

Legal notice:  the terms of use for the Python logo restrict the alteration or bastardization of the logo, including its colors.  I have contacted the Python Software Foundation and gotten permission to publish this blog entry.


The logo restrictions I violate in this post, for purposes of demonstrating working with svg, are:


   1) changing the color of the logo
   2) removal of the trademark symbol


While people are encouraged to use the logo, they are also encouraged to contact the PSF trademarks committee regarding questions of alteration.  Thank you for your consideration in this matter.

As a means of becoming more familiar with SVG (scalable vector graphics) and XML, I wanted to see if I could change the colors of the Python logo in code.  Although this is neither encouraged nor permitted (see notice above), it has been done on occasion in the past; I wanted to see how.

Since SVG is a subset of XML, I'm going to use an XML tool from the standard library (xml.etree.ElementTree) to work with the SVG image.

First, I downloaded the logo from python.org. I'm using the Java based tool Batik/Squiggle to view the SVG images:

To simplify the SVG file I removed the gray parts (text, trademark character, and shadow).  This I did by hand by looking for any gray colors of the hex form #xyxyxy in the file.

Also, I changed the name of the file I was referencing inside the SVG file to pythonlogochanged.svg.  This is important, otherwise the SVG file will continue to read whatever file is named inside it.

 

There are now six main elements to the file:

    1) the two shapes in the logo

    2) the two gradients for color variation across each shape

    3) the two sets of colors for start and stop in the gradation (two yellows and two blues)

It's number 3) that we're interested in.

I started out trying to use minidom from the xml standard library, but found ElementTree to be better suited to the task.  This tutorial from Effbot was helpful.

$ python2.5

Python 2.5.4 (r254:67916, Aug  9 2010, 08:57:51)
[GCC 4.2.1 20070719 ] on openbsd4
Type "help", "copyright", "credits" or "license" for more information.
>>> # open file
...
>>> svgobj = ET.parse('pythonlogochanged.svg')

>>> # get toplevel element

...
>>> rootx = svgobj.getroot()

Now it is time to drill down into the element tree.  I walked through this myself previously, so you'll have to trust that I know where the colors are in the tree.

>>> elementsx = rootx.getchildren()

>>> defsx = elementsx[2]
>>> elementsx = defsx.getchildren()

>>> yellowelement = elementsx[0]

>>> blueelement = elementsx[1]

First we'll look at the yellow color definitions.

>>> yellowchildren = yellowelement.getchildren()

>>> yellowchildren[0].items()

[('style', 'stop-color:#ffd43b;stop-opacity:1'), ('id', 'stop4673'), ('offset', '0')]

>>> yellowchildren[1].items()

[('style', 'stop-color:#ffe873;stop-opacity:1'), ('id', 'stop4675'), ('offset', '1')]

The yellow colors across the gradient are #ffd43b and #ffe873.

Now we'll look at the blue colors.

>>> bluechildren = blueelement.getchildren()

>>> bluechildren[0].items()

[('style', 'stop-color:#5a9fd4;stop-opacity:1'), ('id', 'stop4691'), ('offset', '0')]
>>> bluechildren[1].items()
[('style', 'stop-color:#306998;stop-opacity:1'), ('id', 'stop4693'), ('offset', '1')]

#5a9fd4 and #306998 are the blue colors across the gradient.

Let's try taking the color gradients out first without changing the blue and yellow colors.

>>> bluechildren[1].set('style', 'stop-color:#5a9fd4:stop-opacity:1')

>>> yellowchildren[1].set('style', 'stop-color:#ffd43b;stop-opacity:1')

Now the result needs to be written to disk.

>>> svgobj.write('pythonlogochanged.svg')

Still blue and yellow, but significantly different without the gradient.

Now let's make a black silhouette.

>>> bluechildren[0].set('style', 'stop-color:#000000:stop-opacity:1')

>>> bluechildren[1].set('style', 'stop-color:#000000:stop-opacity:1')
>>> yellowchildren[0].set('style', 'stop-color:#000000;stop-opacity:1')
>>> yellowchildren[1].set('style', 'stop-color:#000000;stop-opacity:1')
>>> svgobj.write('pythonlogochanged.svg')

The hardest part of this excercise was finding the appropriate parts of the SVG file to edit.  ElementTree makes it easy from there.

Simple Polygon Offset

I was working on another POV-ray pysanky egg design and had to work with offsetting polygons.  Stack Overflow had a brief description of how this could be done.  I chose to go with the simplest solution: drawing offset lines and connecting their intersection points.

Update:  ΤΖΩΤΖΙΟΥ was kind enough to introduce me to the math.atan2 function, which correctly identifies the quadrant the angle theta falls in, and rewrote my verbose function(s) concisely:

def calcoffsetpoint(pt1, pt2, offset):
    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)

getoffsetintercept gets the b in y = mx + b needed to calculate the new point:

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

The function that gets a single point along the polygon:

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.
    """
    # 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

Lastly, the function that works the polygon offset:

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 = getpt(polyx[counter],
                   polyx[counter + 1],
                   polyx[counter + 2],
                   offset)
        # append new point to polyy
        polyy.append(pt)
    # last three points
    pt = getpt(polyx[-3], polyx[-2], polyx[-1], offset)
    polyy.append(pt)
    pt = getpt(polyx[-2], polyx[-1], polyx[0], offset)
    polyy.append(pt)
    pt = getpt(polyx[-1], polyx[0], polyx[1], offset)
    polyy.append(pt)
    return polyy

Shown below is the shape I was trying to offset.

As the "Simple" in the entry's title suggests, this is a best case scenario:

    1) no degenerate polygons inside or outside
    2) no zero or infinite slopes

Actually, I should have some zero slopes, but I cheated and changed the coordinates.  The differences won't be seen by the naked eye, but are big enough for the computer to handle them.

Here is the egg I was working on; it's a Hutsul design from the Sixty Score of Easter Eggs book:

Thanks for having a look.

jython copy versus clone

We're a Java shop at work.  We do a lot of passing of Java objects back and forth through the jython API to our application.

The other day I had the need to copy a Java object.  It was a flat structure and my instinct was to use the object's clone method, which I did.  This got me thinking about the use of clone versus the copy module in jython and I gave it a look:

Jython 2.5.1 (Release_2_5_1:6813, Sep 26 2009, 13:47:54)
[Java HotSpot(TM) 64-Bit Server VM (Sun Microsystems Inc.)] on java1.6.0_22
Type "help", "copyright", "credits" or "license" for more information.
>>> from java.util import ArrayList
>>> a = ArrayList()
>>> a
[]
>>> a.add(1)
True
>>> a.add([2, 3])
True
>>> import copy
>>> b = a.clone()
>>> b
[1, [2, 3]]
>>> b[1]
[2, 3]
>>> b[1][1]
3
>>> b[1][1] = 22
>>> b
[1, [2, 22]]
>>> a
[1, [2, 22]]
>>> b = copy.deepcopy(a)
>>> b
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
SystemError: Automatic proxy initialization should only occur on proxy classes

OK, clone works similar to copy.copy (shallow copy).  Not sure what's going on with deepcopy, but it doesn't appear to be available in the case of an ArrayList object.

Update:  as noted in the comments, the copy behavior on Java objects is a known bug and is scheduled to be fixed in jython version 2.5.2.  In the meantime, clone() would appear to be a decent option.

Let's test the shallow copy capability of clone():

>>> a = ArrayList()
>>> a.add(1)
True
>>> a.add(2)
True
>>> a.add(3)
True
>>> a
[1, 2, 3]
>>> b = a.clone()
>>> b
[1, 2, 3]
>>> b[2] = 22
>>> b
[1, 2, 22]
>>> a
[1, 2, 3]

Perfect - it works just like copy.copy() would.