Adjacency ... numpy snippets

[DanPatterson]

Before I forget.

A brief introduction to array broadcasting in numpy.

# ---

# -- Two polygons
a 
array([[  0.00,  10.00],
       [  1.00,  13.00],
       [  3.00,  14.00],
       [  2.50,  13.00],
       [  1.00,  11.50],
       [  0.00,  10.00]])
b
array([[  1.00,  11.50],
       [  2.50,  13.00],
       [  4.00,  12.50],
       [  3.00,  11.00],
       [  2.00,  11.00],
       [  1.00,  11.50]])

Which look like this.

The code is somewhat simple but requires some explanation.

The code is basically one line with the more important documentation.

def _adjacent_(a, b):
    """Check adjacency between 2 polygon shapes on their edges.

    Parameters
    ----------
    a, b : ndarrays
        The arrays of coordinates to check for polygon adjacency.
        The duplicate first/last point is removed so that a count will not
        flag a point as being a meeting point between polygons.

    Note
    ----
    Adjacency is defined as meeting at least 1 point.  Further checks will be
    needed to assess whether two shapes meet at 1 point or 2 non-consequtive
    points.

    """
    s = np.sum((a[:, None] == b).all(-1).any(-1))
    if s > 0:
        return True
    return False

Lets break down line 18.  

a[:, None] == b  # -- `a` is broadcast so that each segment can be compared to `b`
array([[[0, 0],
        [0, 0],
        [0, 0],
        [0, 0],
        [0, 0],
        [0, 0]],

       [[1, 0],  # -- one match
        [0, 1],  # -- one match
        [0, 0],
        [0, 0],
        [0, 0],
        [1, 0]],  # -- one match

       [[0, 0],
        [0, 0],
        [0, 0],
        [1, 0],  # -- one match
        [0, 0],
        [0, 0]],

       [[0, 0],
        [1, 1],  # -- both match --
        [0, 0],
        [0, 0],
        [0, 0],
        [0, 0]],

       [[1, 1],  # -- both match
        [0, 0],
        [0, 0],
        [0, 0],
        [0, 0],
        [1, 1]],  # -- both match

       [[0, 0],
        [0, 0],
        [0, 0],
        [0, 0],
        [0, 0],
        [0, 0]]])

# -- next ... reduce the dimensions by 1

(a[:, None] == b).all(-1)

array([[0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0],
       [1, 0, 0, 0, 0, 1],
       [0, 0, 0, 0, 0, 0]])

# -- next ... and one again.  Note -1, means the last axis

(a[:, None] == b).all(-1).any(-1)
array([0, 0, 0, 1, 1, 0])

# -- finally, the comparison has been reduces from a 3 dimensional array, to 2D to 1D

# -- the result

np.sum((a[:, None] == b).all(-1).any(-1))
2

# two points which `a` meets `b` at points 3 and 4 (good old `nonzero`)

np.nonzero((a[:, None] == b).all(-1).any(-1))[0]
array([3, 4], dtype=int64)

# if we reverse `b` and `a`

np.nonzero((b[:, None] == a).all(-1).any(-1))[0]

array([0, 1, 5], dtype=int64)  # -- 0 and 5 are the start and end point of `b`

Of course, with numpy, this can be scaled up to determine adjacency amongst more than two polygons.  But the simplicity of array `broadcasting` is nice.