Einsum in geometry

DanPatterson · ‎03-27-2025

Take a square representing a polygon or a closed-loop polyline. Let us call it aoi, perhaps out area of interest.

aoi 
array([[  0.00,   0.00],
       [  0.00,  10.00],
       [ 10.00,  10.00],
       [ 10.00,   0.00],
       [  0.00,   0.00]])
aoi.shape
(5, 2)

Notes:

The first and last points are the same, hence it is a closed geometry
The shape of the poly* feature is 5 rows (0 through 4 inclusive) and 2 columns ( 0 and 1), since python numbering is zero-based

Conceptually we know from out maths, that we just need to take the sequential difference between the points to get the differences in the x's and y's. Whip in a quick euclidean calculation courtesy of the right-angle triangle and other cringeful recollections of things numeric.

So it comes as no surprise that anything that involves rows and columns would be painfully slow to do one calculation at a time. Enter einsum notation and its NumPy representation. (also useful for you Tensor types out there)

diff = aoi[1:] - aoi[:-1]
np.sqrt(np.einsum('ij,ij->i', diff, diff))
array([ 10.00,  10.00,  10.00,  10.00])

Which means multiply the two columns, and summarize, just like

diff
array([[  0.00,  10.00],
       [ 10.00,   0.00],
       [  0.00, -10.00],
       [-10.00,   0.00]])

diff * diff
array([[  0.00,  100.00],
       [ 100.00,   0.00],
       [  0.00,  100.00],
       [ 100.00,   0.00]])

np.sum(diff*diff, axis= 1)  # -- axis 1 means by row, since axis 0 is columns
array([ 100.00,  100.00,  100.00,  100.00])

So you want the hypotenuse of your square? Just use einsum summation on the other axis.

np.sqrt(np.einsum('ij,ij->j', diff, diff))  # -- note summation is by 'j' this time
array([ 14.14,  14.14])

The einsum notation string, just needs to be consistent when referencing the axes to perform actions on

np.sqrt(np.einsum('ab,ab->a', diff, diff))

is the same as

np.sqrt(np.einsum('ij,ij->j', diff, diff))

I use einsum in a lot of my Geo array calculations. For example to determine the area of a polygon "bit" (either an inner or outer ring or parts of multipart polygons).

def _bit_area_(a):
    """Mini e_area, used by `areas` and `centroids`.

    Negative areas are holes.  This is intentionally reversed from
    the `shoelace` formula.
    """
    a = _get_base_(a)
    x0, y1 = (a.T)[:, 1:]  # cross set up as follows
    x1, y0 = (a.T)[:, :-1]
    e0 = np.einsum('...i,...i->...i', x0, y0)  # 2024-03-28 modified
    e1 = np.einsum('...i,...i->...i', x1, y1)
    return np.sum((e0 - e1) * 0.5)

This essentially sets up the 2D crossproduct for a modified 'shoelace' formula.

How about determining the distance from one point to a bunch of other points?

I will use 'aoi' as my other points and define a point 'p'.

def _e_2d_(a, p):
    """Return distance from a point to an array of points"""
    diff = a - p[None, :]
    return np.sqrt(np.einsum('ij,ij->i', diff, diff))

p = np.array([3., 4.])
_e_2d_(aoi, p)
array([  5.00,   6.71,   9.22,   8.06,   5.00])

And to end this off, how about a translation and rotation of a polygon? A little use of einsum to simplify the matrix calculations.

def _r_(a, cent, angle, clockwise):
    """Rotate by `angle` in degrees, about the center."""
    angle = np.radians(angle)
    if clockwise:
        angle = -angle
    c, s = np.cos(angle), np.sin(angle)
    R = np.array(((c, -s), (s, c)))
    return np.einsum('ij,jk->ik', a - cent, R) + cent

Using the function above with these inputs

_r_(aoi, cent=[1., 1.], angle=45., clockwise=True)
 
array([[  1.00,  -0.41],
       [ -6.07,   6.66],
       [  1.00,  13.73],
       [  8.07,   6.66],
       [  1.00,  -0.41]])

yields the above coordinates which looks like this

Enough for now.