```python
from __future__ import print_function, division

# Basic python packages
import json
import math
import itertools

# External plotting and analysis tools
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

# Internal plotting and analysis
from omnutils.extends.matplotlib import *
screen_style()
```


```python
from matplotlib.patches import RegularPolygon
from matplotlib.collections import PatchCollection
from matplotlib.colors import Normalize
```

Hex grid
===================
*Seth R Johnson*

*November 12, 2015*



```python
apothem = 1.0 # inner radius
offset = 0.0 # set to 0.0 for flat-topped ("rhex"), 0.5 for pointy-topped (standard KENO "hex")
num_sides = 6
offset = math.fmod(num_sides * 3 + 4 * offset, 4.0) / 4.0
```


```python
hextype = "flattop" if offset == 0.5 else "pointytop"
```


```python
# Face normals and displacement vectors
normals = np.empty((num_sides, 2))
displacements = np.empty((num_sides, 2))
for n in range(num_sides):
    theta = 2 * math.pi * (n + offset) / num_sides
    normals[n] = [math.cos(theta), math.sin(theta)]
    displacements[n] = normals[n] * apothem
```


```python
normals[:3, :]
```




    array([[ 8.66025404e-01,  5.00000000e-01],
           [ 6.12323400e-17,  1.00000000e+00],
           [-8.66025404e-01,  5.00000000e-01]])



With hexes, the unit basis vectors are along the normals through the zeroth and first faces. It's easier to scale the basis vectors by twice the apothem (distance from the center to a face), so that the distance between adjacent hex center points is always 1 in the hex coordinate system:
$$ \vec{u} = 2a \cos \frac{\pi \delta}{3} \vec{i} + 2a \sin \frac{\pi \delta}{3} \vec{j} $$
and 
$$ \vec{v} = 2a \cos \frac{\pi (1 + \delta)}{3} \vec{i} + 2a \sin \frac{\pi (1 + \delta)}{3} \vec{j} $$
Here, $\delta$ is the "net" offset (different from the user-input offset): $\tfrac{1}{2}$ for flat-topped hexes (in KENO, "rhexagon"), and 0 for pointy-topped hexes (in KENO, "hexagon").




```python
bases = 2 * displacements[0:2,:]
bases
```


```python
outer_radius = apothem / math.cos(math.pi / 6)
poly_orient = math.pi * offset
```


```python
1/math.cos(math.pi / 6)
```


```python
fig, ax = plt.subplots()

# Plot hex
patch = RegularPolygon((0., 0.), num_sides, radius=outer_radius, orientation=poly_orient,
                       facecolor='none', edgecolor='k')
ax.set_xlim(-2.2 * apothem, 2.2 * apothem)
ax.set_ylim(-1.5 * apothem, 2.5 * apothem)
ax.set_aspect('equal')
ax.add_patch(patch)
ax.plot(0, 0, 'ko')

# Plot faces
ax.plot(displacements[:,0], displacements[:,1], 'ro')
for n in range(num_sides):
    ax.text(1.25 * displacements[n,0], 1.25 * displacements[n,1], str(n),
           ha='center', va='center')
    
# Plot bases
for (i, lab) in enumerate("uv"):
    basis = bases[i,:]
    ax.arrow(0, 0, basis[0], basis[1], head_width=.1 * apothem, head_length=.15 * apothem, fc='k')
    ax.text(1.1 * basis[0], 1.1 * basis[1], "${}$".format(lab))
# Additional basis function
if 1:
    basis = 2 * displacements[2,:]
    lab = 'w'
    c = (.75,.75,.75)
    ax.arrow(0, 0, basis[0], basis[1], head_width=.1 * apothem, head_length=.15 * apothem, fc=c, color=c)
    ax.text(1.1 * basis[0], 1.1 * basis[1], "${}$".format(lab), color=c)
    
ax.grid()
fig.savefig("../doc/hex-{}-bases.pdf".format(hextype))
```

(Note correspondence between faces of adjacent hexes: `(n + 3) % 6`. We will need to check surface senses, however, to determine whether they're orthogonal or not (because external senses will change). Erf.)

The basis functions allow us to convert from the hexagonal coordinate system $u \vec{u} + v \vec{v}$ to the Cartesian coordinate system $x \vec{i} + y \vec{j}$. We simply evaluate the basis vectors given the rotation offset $\delta$.

For flat-topped hexes, the "hex coordinates" $(u,v)$ relate to the cartesian coordinates $(x, y)$ by the transform
$$
  \begin{bmatrix}
     x \\
     y
  \end{bmatrix} 
   =
   2a
  \begin{bmatrix}
    \tfrac{\sqrt{3}}{2} & \tfrac{1}{2} \\
    0 & 1
   \end{bmatrix} 
   \begin{bmatrix}
     u \\
     v
   \end{bmatrix} 
$$

For pointy-topped hexes ($\delta = \tfrac{1}{2}$), the transformation is 
$$
  \begin{bmatrix}
     x \\
     y
  \end{bmatrix} 
   =
   2a
  \begin{bmatrix}
    1 & 0 \\
    \tfrac{1}{2} & \tfrac{\sqrt{3}}{2}
   \end{bmatrix} 
   \begin{bmatrix}
     u \\
     v
   \end{bmatrix} 
$$


```python
def retick(ax, dx=2):
    (lo, hi) = ax.get_xlim()
    ax.set_xticks(np.arange(math.ceil(lo), math.floor(hi + dx), dx))
    (lo, hi) = ax.get_ylim()
    ax.set_yticks(np.arange(math.ceil(lo), math.floor(hi + dx), dx))
```


```python
fig, ax = plt.subplots()

patches = []
(umax, vmax) = (3, 5)
for (u, v) in itertools.product(range(umax), range(vmax)):
    uv = np.array([u, v])
    xy = uv.dot(bases) # multiply by the column vector (u,v)
    patches.append(RegularPolygon(xy, num_sides, outer_radius, poly_orient))
    ax.text(xy[0], xy[1], "$({},{})$".format(u, v), ha='center', va='center')

ax.set_xlim(-1.2 * outer_radius, xy[0] + 1.2 * outer_radius)
ax.set_ylim(-1.2 * outer_radius, xy[1] + 1.2 * outer_radius)
ax.set_aspect('equal')
retick(ax)
ax.add_collection(PatchCollection(patches, facecolor='none', edgecolor='k'))
ax.grid()
fig.savefig("../doc/hex-{}-coords-{:d}x{:d}.pdf".format(hextype, umax, vmax))
```

## Convert from Cartesian to $(u,v)$

To convert from Cartesian to grid coordinates, we invert the basis functions: $$
  \begin{bmatrix}
     u \\
     v
  \end{bmatrix} 
   =
   B^{-1}
   \begin{bmatrix}
     x \\
     y
   \end{bmatrix} 
   $$


```python
inv_bases = np.linalg.inv(bases)
inv_bases
```


```python
def xy_to_uv(x, y):
    xy = np.array([x,y])
    return xy.dot(inv_bases)
```

This gives us decimal coordinates in the $(u,v)$ basis system


```python
print(xy_to_uv(2.0, 0))
print(xy_to_uv(1.0, 2.0))
print(xy_to_uv(8.0, 10.0))
```


```python
x = np.linspace(-2 * apothem, 2 * apothem, 65)
y = x
(X, Y) = np.meshgrid(x, y)
XY  = np.stack([X,Y])
```


```python
UV = np.empty((2, y.size, x.size))
for (j, i) in itertools.product(range(y.size), range(x.size)):
    UV[:,j,i] = XY[:,j,i].dot(inv_bases)
U = UV[0,:,:]
V = UV[1,:,:]
```


```python
norm = Normalize(vmin=-1, vmax=1.)
cmap = 'rb_linear'

def plot_hex_field(ax, val):
    pc = ax.pcolormesh(X, Y, val, norm=norm, cmap=cmap)
    patches = []
    points = []
    for (u, v) in itertools.product(range(-1,2), range(-1,2)):
        uv = np.array([u, v])
        xy = uv.dot(bases) # multiply by the column vector (u,v)
        patches.append(RegularPolygon(xy, num_sides, outer_radius, poly_orient))
        patch = RegularPolygon((0., 0.), num_sides, radius=outer_radius, orientation=poly_orient,
                           facecolor='none', edgecolor='k')
        points.append(xy)
    ax.add_collection(PatchCollection(patches, facecolor='none', edgecolor='k'))
    points = np.array(points)
    ax.plot(points[:,0], points[:,1], 'ko')
    ax.set_aspect('equal')
    ax.set_xlim(x[0], x[-1])
    ax.set_ylim(y[0], y[-1])
    retick(ax, 1.0)
    return pc
```


```python
(fig, (axu, axv)) = plt.subplots(ncols=2, figsize=(8, 3))
pcu = plot_hex_field(axu, U)
axu.set_title('$u$ value')
pcv = plot_hex_field(axv, V)
axv.set_title('$v$ value')
fig.colorbar(pcu)
```

These pseudocolor plots show that our $(x, y) \to (u, v)$ conversion correctly project the hex axes; isovalues are along the $\vec u$ and $\vec v$ normals.

## Third basis function

A solution to rounding these correctly to the integer hex coordinates is to rephrase the coordinate system to add a third "unit vector", drawn in the earlier graph as the grayed-out $w$ along face index 2. It's a simple linear combination of the other two axes, defined so that $u + v + w = 0$:


```python
W = -U - V
```


```python
fig, ax = plt.subplots()
pc = plot_hex_field(ax, W)
axv.set_title('$w$ value')
fig.colorbar(pc)
```


```python
max_uvw = np.maximum.reduce([U,V,W,-U,-V,-W,])

fig, ax = plt.subplots()
pc = plot_hex_field(ax, max_uvw)
ax.contour(X, Y, max_uvw, [0.5, 1.0], colors='g')
fig.colorbar(pc);
```

Next, we round the $(u,v,w)$ coordinates to the nearest whole number:


```python
UVW = np.stack([U, V, W])

rounded = np.round(UVW)
```


```python
fig, axes = plt.subplots(ncols=3, figsize=(12, 2.5))
for (i, ax) in enumerate(axes):
    pc = plot_hex_field(ax, rounded[i,:,:])
fig.colorbar(pc)
```

To be consistent, the rounded $(u,v,w)$ values should add to zero just as the unrounded did. We can see that on the corners of the hex, this is not the case:


```python
summed = rounded[0] + rounded[1] + rounded[2]

fig, ax = plt.subplots()
pc = plot_hex_field(ax, summed)
fig.colorbar(pc)
```


```python
diffs   = np.abs(rounded - UVW)
(du, dv, dw) = (diffs[i] for i in range(3))
u_worst = (du > dw) & (du > dv)
v_worst = ~u_worst & (dv > dw)
w_worst = ~u_worst & ~v_worst
```


```python
(ru, rv, rw) = (rounded[i] for i in range(3))
ru[u_worst] = -rv[u_worst] - rw[u_worst]
rv[v_worst] = -ru[v_worst] - rw[v_worst]
rw[w_worst] = -ru[w_worst] - rv[w_worst]
```


```python
summed = ru + rv + rw

fig, ax = plt.subplots()
pc = plot_hex_field(ax, summed)
fig.colorbar(pc)
```

This shows that the rounded values satisfy the $u + v + w = 0$ relationship. Finally, we look at the rounded $(u,v,w)$ values themselves to see if they give us the correct coordinates:


```python
fig, axes = plt.subplots(ncols=3, figsize=(12, 2.5))
for (label, data, ax) in zip("uvw", (ru, rv, rw), axes):
    pc = plot_hex_field(ax, data)
    ax.set_title(label)
fig.colorbar(pc)
```

## Simpler 'rounding' function

This is based on the "Branchless Method" described in http://www.redblobgames.com/grids/hexagons/more-pixel-to-hex.html#neighbors

> All the pixel to hex conversions I’ve seen use branches or a lookup table. I was mystified when Charles Chambers sent me pixel to hex conversion code that uses floor() five times, but no branches. First, divide x and y by size * sqrt(3); then find q, r with:
>
>     temp = floor(x + sqrt(3) * y + 1)
    q = floor((floor(2*x+1) + temp) / 3);
    r = floor((temp + floor(-x + sqrt(3) * y + 1))/3);

Extracting the embedded expressions, it's clear that the floored expressions are coordinates along the $u$, $v$, and $w$ basis vectors:


```python
uvw_bases = normals[:3] / apothem
uvw_bases
```


```python
A = np.floor(uvw_bases[0,0] * X + uvw_bases[0,1] * Y)
fig, ax = plt.subplots()
plot_hex_field(ax, A)
```


```python
B = np.floor(uvw_bases[1,0] * X + uvw_bases[1,1] * Y)
fig, ax = plt.subplots()
plot_hex_field(ax, B)
```


```python
C = np.floor(uvw_bases[2,0] * X + uvw_bases[2,1] * Y)
fig, ax = plt.subplots()
plot_hex_field(ax, C)
```

Adding these is effectively like 'OR'ing the enclosed regions; dividing by 2 or 3 allows the non-integer regions to expand to fill parts of the hexes. Different combinations select different regions; we want regions that have isovalues inside the hexagon.


```python
fig, axes = plt.subplots(ncols=3, nrows=2, figsize=(10,6))
plot_hex_field(axes[0,0], (A + B) / 3)
plot_hex_field(axes[0,1], (B + C) / 3)
plot_hex_field(axes[0,2], (A + C) / 3)
plot_hex_field(axes[1,0], (A - B) / 3)
plot_hex_field(axes[1,1], (B - C) / 3)
plot_hex_field(axes[1,2], (C - A) / 3)
```

We can translate these rhombi as needed; in this case we want to move the $(C - A)$ choice up to select the third corner of the hex in analog to $(A+B)$ and $(B + C)$...


```python
fig, axes = plt.subplots(ncols=3, nrows=1, figsize=(10,6))
plot_hex_field(axes[0], (A - B + 1) / 3)
plot_hex_field(axes[1], (B - C + 1) / 3)
plot_hex_field(axes[2], (A - C + 1) / 3)
```


```python
fig, axes = plt.subplots(ncols=3, figsize=(10,3))
plot_hex_field(axes[0], (A + B))
plot_hex_field(axes[1], (A + B) / 3)
plot_hex_field(axes[2], np.ceil((A + B) / 3))
```


```python
fig, axes = plt.subplots(ncols=4, figsize=(10,3))
plot_hex_field(axes[0], (B + C))
plot_hex_field(axes[1], (B + C) / 3)
plot_hex_field(axes[2], np.ceil((B + C) / 3))
```


```python
fig, axes = plt.subplots(ncols=3, figsize=(10,3))
plot_hex_field(axes[0], (1 + A - C))
plot_hex_field(axes[1], (1 + A - C) / 3)
plot_hex_field(axes[2], np.floor((1 + A - C) / 3))
```

This is much more straightforward than that cube-and-rounding system of doing things.


```python
uvw_bases = normals[:3]
def xy_to_uv(x, y):
    # Dot basis functions with xy coordinates
    (a, b, c) = np.array(uvw_bases.dot([x, y]) / apothem, dtype=int)
    u = -((-(a + b)) // 3)
    v = (1 + a - c) // 3
    return (u, v)
```


```python
print(xy_to_uv(2.0, 0))
print(xy_to_uv(1.0, 2.0))
print(xy_to_uv(8.0, 10.0))
```


```python
uvw_bases
```


```python
uvw_bases.dot([0, 1])
```

### Shexagonal

Keno's 'hexagonal' and 'triangular' are the same; they generate a rhombus-looking array with pointy tops. The 'shexagonal' option staggers the grid to look like a square layout.


```python
fig, ax = plt.subplots()

patches = []
umax, vmax = (3, 5)
(urangemax, vrangemax) = (umax, vmax)
if offset == 0.0: # pointy top
    urangemax += (vmax - 1) // 2
else:
    vrangemax += (umax - 1) // 2

for (u, v) in itertools.product(range(urangemax), range(vrangemax)):
    uv = np.array([u, v])
    xy = uv.dot(bases)
    fc = 'none'
    if offset == 0.0:
        uprime = u + v // 2 - (vmax - 1) // 2
        vprime = v
    else:
        vprime = v + u // 2 - (umax - 1) // 2
        uprime = u
    if (0 <= uprime < umax) and (0 <= vprime < vmax):
        fc = (0,0,0,.25)
    patches.append(RegularPolygon(xy, num_sides, outer_radius, poly_orient, fc=fc, edgecolor='k'))
    ax.text(xy[0], xy[1], "{},{}".format(uprime, vprime), ha='center', va='center')

ax.set_xlim(-1.2 * outer_radius, xy[0] + 1.2 * outer_radius)
ax.set_ylim(-1.2 * outer_radius, xy[1] + 1.2 * outer_radius)
ax.set_aspect('equal')
retick(ax)
ax.add_collection(PatchCollection(patches, match_original=True))
ax.grid()
fig.savefig("../doc/hex-{}-rect-{:d}x{:d}.pdf".format(hextype, umax, vmax))
```

So we need to map:
\begin{align}
u_0 &= \lfloor \tfrac{v_\text{max} - 1}{2} \rfloor \\
u' &= u + \lfloor \tfrac{v}{2} \rfloor - u_0 \\
v' &= v
\end{align}
for $u = \{0,\ldots,u_\text{max} + u_0\}$, $v = \{0,\ldots, v_\text{max}\})$
in the pointy-topped case and
\begin{align}
v_0 &= \lfloor \tfrac{u_\text{max} - 1}{2} \rfloor \\
u' &= u \\
v' &= v - \lfloor \tfrac{u}{2} \rfloor - v_0
\end{align}
for $u = \{0,\ldots,u_\text{max}\}$, $v = \{0,\ldots, v_\text{max} + v_0\})$
in the flat-topped case.

d700e6b9b0dfcf7ab8f8c4a19312e71c16fc9c7d
ipynb
Jupyter Notebook
doc/.ipynbcheckpoints/hex-grid-checkpoint.ipynb
celeritas-project/orange-port
9aa2d36984a24a02ed6d14688a889d4266f7b1af
Apache-2.0 MIT
doc/.ipynbcheckpoints/hex-grid-checkpoint.ipynb
celeritas-project/orange-port
9aa2d36984a24a02ed6d14688a889d4266f7b1af
Apache-2.0 MIT
doc/.ipynbcheckpoints/hex-grid-checkpoint.ipynb
celeritas-project/orange-port
9aa2d36984a24a02ed6d14688a889d4266f7b1af
Apache-2.0 MIT
Qwen/Qwen-72B
1. YES
2. YES