Faster way to calculate hexagon grid coordinates

Question

I'm using the following procedure to calculate hexagonal polygon coordinates of a given radius for a square grid of a given extent (lower left --> upper right):

def calc_polygons(startx, starty, endx, endy, radius):
    sl = (2 * radius) * math.tan(math.pi / 6)

    # calculate coordinates of the hexagon points
    p = sl * 0.5
    b = sl * math.cos(math.radians(30))
    w = b * 2
    h = 2 * sl


    origx = startx
    origy = starty

    # offsets for moving along and up rows
    xoffset = b
    yoffset = 3 * p

    polygons = []
    row = 1
    counter = 0

    while starty < endy:
        if row % 2 == 0:
            startx = origx + xoffset
        else:
            startx = origx
        while startx < endx:
            p1x = startx
            p1y = starty + p
            p2x = startx
            p2y = starty + (3 * p)
            p3x = startx + b
            p3y = starty + h
            p4x = startx + w
            p4y = starty + (3 * p)
            p5x = startx + w
            p5y = starty + p
            p6x = startx + b
            p6y = starty
            poly = [
                (p1x, p1y),
                (p2x, p2y),
                (p3x, p3y),
                (p4x, p4y),
                (p5x, p5y),
                (p6x, p6y),
                (p1x, p1y)]
            polygons.append(poly)
            counter += 1
            startx += w
        starty += yoffset
        row += 1
    return polygons

This works well for polygons into the millions, but quickly slows down (and takes up very large amounts of memory) for large grids. I'm wondering whether there's a way to optimise this, possibly by zipping together numpy arrays of vertices that have been calculated based on the extents, and removing the loops altogether – my geometry isn't good enough for this, however, so any suggestions for improvements are welcome.

Su

This ought to be migrated to codereview. – Micronesia

Actually it's a pretty good question about hex grids. Is the ordering of polygons within the list relevant to you? – Fra

@MariaZverina Not necessarily, no. – Su

Urschrei, are you interested in finding one enclosing hexagon, belonging to the grid, that contains a specific P(x,y), without computing the whole grid in advance? Further, asking everyone, is there an established method to do so? thank you in advance – Colonic

@Colonic No, I'm calculating the grid in order to perform some spatial operations – this isn't a game-dev question; this might be useful to you: redblobgames.com/grids/hexagons/#pixel-to-hex – Su

You totally don't need this, in order to convert Cartesian to hexagonal coordinates (offset, axial or cube). See the redblobgames.com/grids/hexagons/#coordinates tutorial that urschrei. Coordinate conversion is like ~15 lines of code. No memory required. – Angwantibo

Answer 1

Decompose the problem into the regular square grids (not-contiguous). One list will contain all shifted hexes (i.e. the even rows) and the other will contain unshifted (straight) rows.

def calc_polygons_new(startx, starty, endx, endy, radius):
    sl = (2 * radius) * math.tan(math.pi / 6)

    # calculate coordinates of the hexagon points
    p = sl * 0.5
    b = sl * math.cos(math.radians(30))
    w = b * 2
    h = 2 * sl


    # offsets for moving along and up rows
    xoffset = b
    yoffset = 3 * p

    row = 1

    shifted_xs = []
    straight_xs = []
    shifted_ys = []
    straight_ys = []

    while startx < endx:
        xs = [startx, startx, startx + b, startx + w, startx + w, startx + b, startx]
        straight_xs.append(xs)
        shifted_xs.append([xoffset + x for x in xs])
        startx += w

    while starty < endy:
        ys = [starty + p, starty + (3 * p), starty + h, starty + (3 * p), starty + p, starty, starty + p]
        (straight_ys if row % 2 else shifted_ys).append(ys)
        starty += yoffset
        row += 1

    polygons = [zip(xs, ys) for xs in shifted_xs for ys in shifted_ys] + [zip(xs, ys) for xs in straight_xs for ys in straight_ys]
    return polygons

As you predicted, zipping results in much faster performance, especially for larger grids. On my laptop I saw 3x speedup when calculating 30 hexagon grid - 10x speed for 2900 hexagon grid.

>>> from timeit import Timer
>>> t_old = Timer('calc_polygons_orig(1, 1, 100, 100, 10)', 'from hexagons import calc_polygons_orig')
>>> t_new = Timer('calc_polygons_new(1, 1, 100, 100, 10)', 'from hexagons import calc_polygons_new')
>>> t_old.timeit(20000)
9.23395299911499
>>> t_new.timeit(20000)
3.12791109085083
>>> t_old_large = Timer('calc_polygons_orig(1, 1, 1000, 1000, 10)', 'from hexagons import calc_polygons_orig')
>>> t_new_large = Timer('calc_polygons_new(1, 1, 1000, 1000, 10)', 'from hexagons import calc_polygons_new')
>>> t_old_large.timeit(200)
9.09613299369812
>>> t_new_large.timeit(200)
0.7804560661315918

There might be an opportunity to create an iterator rather than the list to save memory. Depends on how your code is using the list of the polygons.

Fra

Passer by, but it may be good to state what each coordinate variable represents. It is not clear from sl, p, b, w, h. – Urbai

Answer 2

Here is a solution that doesn't require any looping. It creates grid of 50x50 hexagons:

coord_x, coord_y = np.meshgrid(50, 50, sparse=False, indexing='xy')

ratio = np.sqrt(3) / 2
coord_y = coord_y * ratio          # Condense the coordinates along Y-axes
coord_x = coord_x.astype(np.float)
coord_x[1::2, :] += 0.5            # Shift every other row of the grid
coord_x = coord_x.reshape(-1, 1)   # Flatten the grid matrices into [2500, 1] arrays
coord_y = coord_y.reshape(-1, 1)

radius = 5                         # Inflate each hexagon to the required radius
coord_x *= radius 
coord_y *= radius

This is a code snippet from my hexalattice package in python, that does that for you (+ advanced options of grid rotation, size and plotting)

Here you can find additional examples and link to the repo: LINK

Goldsberry