Geometry puzzler

[Alex Martelli]

"Paul Winkler" <slinkp23 at yahoo.com> wrote in message
news:3A073A51.AE3691A3 at yahoo.com...
> Frédéric van der Plancke wrote:
> >
> > The comp.graphics.algorithms FAQ answers your geometry questions too.
> > FvdP.
>
> Ooh, that's a good one. Thanks all!
> I'll check out those books too.

Yep, very neat collection of formulas, thanks Frédéric!

Translating those formulas into Python is a simple and
interesting exercise.  For example, representing points
as pairs (2-element tuples) of coordinates, and segments
as pairs of points, one could write:

def intersect( ((Ax,Ay),(Bx,By)), ((Cx,Cy),(Dx,Dy)) ):
    denom = (Bx-Ax)*(Dy-Cy)-(By-Ay)*(Dx-Cx)
    if not denom: return None    # collinear segments
    r = ((Ay-Cy)*(Dx-Cx)-(Ax-Cx)*(Dy-Cy)) / denom
    s = ((Ay-Cy)*(Bx-Ax)-(Ax-Cx)*(By-Ay)) / denom
    if 0<=r<=1 and 0<=s<=1:
        return Ax+r*(Bx-Ax), Ay+r*(By-Ay)

I like the way this exploits explicit-tuple-unpacking
(much like FP language's more extensive use of similar
pattern-matching ideas) to let me call this in terms
of either coordinate-level, point-level, or segment-
level objects, while letting all computations be a
straightforward transcription of the formulas... e.g.:

>>> origin=(0.0,0.0)
>>> diag_ne_sw=((1.0,0.0), (0.0,1.0))
>>> for i in range(-3,4):
 factor = i/10.0
 print i, geom.intersect( (origin,(1.0+factor,1.0-factor)), diag_ne_sw)


-3 (0.34999999999999998, 0.65000000000000002)
-2 (0.40000000000000002, 0.59999999999999998)
-1 (0.45000000000000001, 0.55000000000000004)
0 (0.5, 0.5)
1 (0.55000000000000004, 0.45000000000000001)
2 (0.59999999999999998, 0.40000000000000002)
3 (0.65000000000000002, 0.34999999999999998)
>>>


Of course, intersect as written is doing more work than
it needs to -- take it as a good, simple exercise in
factoring out common expressions...:

def intersect( ((Ax,Ay),(Bx,By)), ((Cx,Cy),(Dx,Dy)) ):
    deXba = Bx-Ax
    deYba = By-Ay
    deYdc = Dy-Cy
    deXdc = Dx-Cx
    denom = deXba*deYdc - deYba*deXdc
    if not denom: return None    # collinear segments
    deYac = Ay-Cy
    deXac = Ax-Cx
    r = (deYac*deXdc - deXac*deYdc) / denom
    if r<0 or r>1: return None
    s = (deYac*deXba - deXac*deYba) / denom
    if s<0 or s>1: return None
    return Ax+r*deXba, Ay+r*deYba

...more verbose, yes, but no less clear -- maybe more,
as it emphasizes the symmetry of the computations --
and it avoids some recomputing (or computing of s
that is not needed if r is outside of range...).  Next
obvious optimization step is avoiding the division by
denom when possible -- but it's trickier, to take
proper account of the signs in play, to do so without
impairing clarity (and risk introducing a bug...).
Neat classroom exercise!


The point-in-polygon formula is also interesting:
representing a polygon by the sequence of its
vertices, we get:

def pnpoly(polygon_vertices, (x,y)):
    inside = 0
    xpj, ypj = polygon_vertices[-1]
    for xpi, ypi in polygon_vertices:
        if (ypi <= y < ypj or ypj <= y < ypi) and (
            x < (xpj-xpi)*(y-ypi)/(ypj-ypi) + xpi):
                inside = not inside
        xpj, ypj = xpi, ypi
    return inside

Here, optimization is not much of an issue, but
clarity and simplicity are always worth attention.
The naming of xpj, etc, is doubtful, for example --
I've chosen it for parallelism with xp[j] &tc in
the original C source.  Maybe taking the slight hit
(in performance terms) of an embedded function...:

    xpj, ypj = polygon_vertices[-1]
    def between(y1, y, y2):
        return y1<=y<y2 or y2<=y<y1
    for xpi, ypi in polygon_vertices:
        if between(ypi, y, ypj) and (

would be worth it for enhanced clarity (coming
just from the fact that we _name_ the test we
are performing).  And again, tuple pack/unpack
and the Python for statement show to good use.

Recoding the second part of the condition in
more symmetrical terms, making it into...:

            (x-xpi)*(ypj-ypi) < (xpj-xpi)*(y-ypi)):

might be nice for both clarity and performance,
BUT it's only correct if ypj-ypi is > 0 -- the
same kind of sign-issue we met in intersect...!


In other words, this seems to be mostly about
style and idiom issues... good things to ponder!-)



Alex