*(profits as yet unrealized)

This post will discuss the how and a little bit of the why for my latest project, a Python package called hexfft, which implements two types of FFT for data sampled on a hexagonal 2D grid.

Geometry of regular hexagonal sampling

There are several interesting advantages to this sampling geometry. For regular hexagonal sampling, each point has 6 equidistant nearest neighbors, compared with 4 in square sampling. This scheme of tiling the 2D plane reduces grid edge effects and also for more flexible path-finding and boundary definitions. An oft-cited figure is that for a circularly bandlimited signal in the plane, 13.4% fewer samples are required to fully specify the spectrum, as a hexagonal bandregion more tightly encloses a circle than a square one. I will leave the details of the advantages of hexagonal sampling to a future discussion, as the primary purpose of this post is introduce the basic concepts used in the hexfft package.

Containing a circular region with hexagonal vs square sampling

Honeycomb-like structures are found throughout nature and science. Besides being highly aesthetically pleasing, this layout is used (or evolves naturally) when a maximally efficient coverage of the 2D plane is required. Two examples are the layout of electrodes on a probe used to detect electrical signals from neurons in the mouse brain, and the arrangement of optical fibers used to image galaxies (see below).

Honeycomb beehive structure. Courtesy of Nautilus magazine

Layout of cone cells in the human retina. Courtesy of Coleman, Sonya & Scotney, Bryan & Gardiner, B.. (2009). Design of Feature Extraction Operators for use on Biologically Motivated Hexagonal Image Structures. Proceedings of the 11th IAPR Conference on Machine Vision Applications, MVA 2009.

Neuropixels 1.0 neuroelectrophysiology probe. Courtesy of the Allen Institute

Layout of optical fibers in an astronomical detector. Courtesy of the Sloan Digital Sky Survey

The seminal paper on processing hexagonally sampled images is Mersereau’s 1979 paper, The Processing of Hexagonally Sampled Two-Dimensional Signals. Mersereau introduces an FFT algorithm for hexagonally sampled data on a hexagon-shaped (and hexagonally periodic in the plane) region of support. This was followed up by Erhardt’s 1993 Hexagonal Fast Fourier Transform with Rectangular output, which extends the method to regions with rectangular periodicity. These two algorithms are implemented in the hexfft library. Rummelt and Birdsong discovered and patented a coordinate system in which the hexagonal Fourier kernel is separable. There is also work covering transforms on hexagonal (and other) lattices by the algebraic signal processing community. A 2005 book, Hexagonal Image Processing, also covers the topic and introduces an array of fascinating and exotic hexagonal coordinate systems.

As I became curious about this topic, I was surprised to find a lack of ready-to-use implementations of hexagonal FFTs. The closest I could find was HexagDLy and the accompanying publication by Steppa and Holch (2019) which implements convolutions on hexagonal grids with an eye towards neural network processing of astronomical data. If hexfft’s niche has already been filled somewhere, I can only blame my own lack of Google savvy.

I would like to add the additional disclaimer that my implementations of the below algorithms in hexfft are likely far from optimal, and my intent with project is primarily (auto)didactic.

The mathematics of the hexagonal FFT

The exposition below closely follows sections 1.2 through 1.4 of Multidimensional Signal Processing, 2nd ed. by Dudgeon and Mersereau.

Coordinate system

Many coordinate systems for hexagonal grids have been proposed, both in the literature above and in the programming community. Mersereau and Ehrhardt both use an oblique coordinate system where the y-axis is slanted 30 degrees from its usual position. In other words, we use the lattice on $\mathbb{R}^2$ generated by the vectors $(1, 0)$ and $(-\frac{1}{2}, \frac{\sqrt{3}}{2})$. In this system, the gridpoints are laid out in a regular hexagonal pattern with a spacing of 1 between all nearest neighbors.

Oblique coordinate system

The plot above shows an 8x8 “square” in oblique coordinates (defined by ${x < 8, y < 8}$), corresponding to a parallelogram on the plane.

Sampling and periodicity

A lattice on $\mathbb{R}^2$ can be expressed as the integer linear combinations of any basis vectors $\mathbf{v_1}$ and $\mathbf{v_2}$. If a 2D discrete signal is composed of samples taken on a regular lattice, we can relate the discrete function $x: \mathbb{Z}^{2} \rightarrow \mathbb{C}$ and the underlying continuous function $\tilde{x}: \mathbb{R}^2 \rightarrow \mathbb{C}$ by defining

\[\begin{align*} V = \big( \mathbf{v_1} & | \mathbf{v_2} \big) \end{align*}\] \[\begin{align*} x(\mathbf{n}) := \tilde{x}(V\mathbf{n}) \end{align*}\]

For example, if a signal is sampled on a rectilinear grid with a spacing of 0.5 in the $x$ direction and $1$ in the $y$ direction, then

\[\begin{align*} V = \begin{pmatrix} 0.5 & 0 \\ 0 & 1 \end{pmatrix} \end{align*}\]

The value of $x$ at the point $(1, 1)$ on the lattice is then given by

\[\begin{align*} x(1, 1) = \tilde{x}(0.5, 1) \end{align*}\]

We refer to $V$ as the sampling matrix. For a general hexagonal lattice with spacings $d_1$ and $d_2$ in the $x$ and $y$ directions respectively,

\[\begin{align*} V = \begin{pmatrix} d_1 & -\frac{d_1}{2} \\ 0 & d_2 \end{pmatrix} \end{align*}\]

If $d_1 = d_2 = 1$, then this reduces to

\[\begin{align*} \begin{pmatrix} 1 & -0.5 \\ 0 & 1 \end{pmatrix} \end{align*}\]

The hexshow() function in hexfft.plot uses this grid to display hexagonal arrays as there is a convenient correspondence between array indices and physical distance between points. For regular hexagonal spacing, where the distance between each point and its six neighbors is 1,

\[\begin{align*} V = \begin{pmatrix} 1 & -1/2\\ 0 &\sqrt{3}/2 \end{pmatrix} \end{align*}\]

Discrete Fourier Transforms on Lattices

A signal $x: \mathbb{Z}^{2} \rightarrow \mathbb{C}$ is periodic with periodicity matrix $P$ if for any integer vector $\mathbf{r}$,

\[\begin{align*} x(\mathbf{n}) = x(\mathbf{n} + P\mathbf{r}) \end{align*}\]

For a periodic signal regardless of sampling strategy, we use the following definition of the 2D discrete Fourier transform pair:

\[\begin{align} X(k_1, k_2) = \sum_{n_1} \sum_{n_2} x(n_1, n_2) \exp\big(-2\pi i \mathbf{k}^T P^{-1} \mathbf{n}\big) \\ x(n_1, n_2) = \frac{1}{\det P}\sum_{k_1} \sum_{k_2} X(k_1, k_2) \exp\big(2\pi i \mathbf{k}^T P^{-1} \mathbf{n}\big) \end{align}\]

The normalizing element $\frac{1}{\det P}$ corresponds to the area of one periodic unit of the signal. It is easy to see that this form reduces to the usual DFT and IDFT for square sampling on an $N$ by $N$ grid,

\[\begin{align*} P = N\begin{pmatrix}1 & 0\\ 0 & 1\end{pmatrix} \end{align*}\]

This form is agnostic to the specific geometry of the lattice (the sampling matrix $V$ does not come into the equations). It is crucial, however, for interpreting the output of the Fourier transform.

Mersereau’s hexagonal FFT

In The Processing of Hexagonally Sampled Two-Dimensional Signals, Mersereau introduces a hexagonal region of support parameterized by an integer $N$. This region always has area $3N^{2}$. In oblique coordinates $n_1, n_2$ the region is the set

\[\begin{align*} R_H (N) = \\ {n_1, n_2 : 0 \leq n_1 < 2N, \\ 0 \leq n_2 < 2N \\ -N \leq n_1 - n_2 < N\\} \end{align*}\]

Note: Mersereau defines this region more generally with parametrizable height and width. For the purposes of hexfft I restricted this to the hexagonal region defined above

The hexagonal region of support $R_H(N)$ is shown for $N=3$ (shaded in light blue). The periodic extensions of $R_H(3)$ are shown in blue dashed outlines. Red arrows show linear combinations of the periodicity vectors $(3, 6)$ and $(6, 3)$ which point to the locations of the periodic repetitions. These are the columns of periodicity matrix $P$.

Given $N$, it’s easy to see from the figure above that this region fits snugly into an oblique “square” of size $2N$. The function hsupport() from hexfft.utils can automatically compute a mask of this region for a square array of a given size:

A hexagonal region à la Mersereau with $N=8$ embedded in a 16x16 oblique array, generated via hexfft.utils.hsupport()

In oblique coordinates, the periodicity matrix for $R_H(N)$ is given by

\[\begin{align*} P = N\begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} \end{align*}\]

Then

\[\begin{align*} P^{-1} = \frac{1}{3N}\begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix} \end{align*}\]

Note the column vectors are in oblique coordinates.

Using the formula from the previous section, we can compute the hexagonal Fourier kernel:

\[\begin{align*} X(k_1, k_2) = \sum_{n_1} \sum_{n_2} x(n_1, n_2) \exp\bigg(-\frac{2\pi i}{3N}\begin{pmatrix} k_1\\ k_2 \end{pmatrix}^{T} \begin{pmatrix} 2 & -1\\ -1 & 2 \end{pmatrix} \begin{pmatrix} n_1\\ n_2\end{pmatrix} \bigg) \end{align*}\] \[\begin{align*} = \sum_{n_1} \sum_{n_2} x(n_1, n_2) \exp\bigg(-\frac{2\pi i}{3N} \big(2 k_1 n_1 + 2 k_2 n_2 - k_1 n_2 - k_2 n_2\big) \bigg) \end{align*}\] \[\begin{align} = \sum_{n_1} \sum_{n_2} x(n_1, n_2) \exp\bigg(-\pi i \bigg(\frac{(2 k_1 - k_2)(2 n_1 - n_2)}{3N} + \frac{k_2 n_2}{N}\bigg) \bigg) \end{align}\]

Where $n_1, n_2 \in R_H(N)$ and $k_1, k_2$ also $\in R_H(N)$, i.e. the frequency values are sampled on the same gridpoints. This last expression is the form given by Mersereau in the paper. The problem, which Mersereau calls an “insurmountable” difficulty, is that unlike the 2D FFT in rectilinear coordinates, this kernel is not separable. Therefore it is not possible, in this coordinate system, to reduce the 2D transform to two sequences of 1D FFTs. Through a clever resampling of $R_H(N)$ and splitting the problem into smaller hexagonal “slow” DFTs, a speed-up from $(3N^{2})^2 = 9 N^{4}$ to $9N^{2}/4$ complex multiplications is achieved by the FFT described by Mersereau. The process is quite complex and I refer the reader to the original paper for a full description. This transform is implemented via hexfft.fft() and hexfft.ifft() with the periodicity="hex" option.

As mentioned above, a separable kernel has since been found using an alternate coordinate system.

“Rectangular” Hexagonal FFT

Mersereau’s FFT algorithm described above can accurately be called a “hexagonal FFT with hexagonal periodicity” due to the fact that the sampling and periodicity matrix are both hexagonal. While the case of a hexagonal region of support could be useful in representing signals with support on, say, the unit disc, it is not the most general format in which one could find hexagonally sampled data. What we want, perhaps, is a transform which takes in data on a rectangular region of support, but still with hexagonal sampling. In Hexagonal Fast Fourier Transform with Rectangular output, Ehrhardt provides just this. It is worth noting that both the input and output of this method are rectangular, and that the frequency points are arranged hexagonally (“Rectangular output” perhaps implies some sort of re-sampling to rectlinear grid points). It might be strictly correct to call this transform a “hexagonal FFT with rectangular periodicity”, but for brevity I’ll refer to it as a “rectangular” FFT with the understanding that all sampling in this post is assumed to be hexagonal.

The rectangular Fourier transform is defined on a “rectangle” in oblique coordinates, corresponding to a parallelogram in “physical” space.

\[R_R(N_1, N_2) = \\{ n_1, n_2 : n_1 < N_1, n_2 < N_2 \\}\]

For a reason that will be explained below, we pick the following periodicity matrix:

\[\begin{align*} P = \begin{pmatrix} N_1 & \frac{N_2}{2} \\ 0 & N_2 \end{pmatrix} \end{align*}\]

so that

\[\begin{align*} P^{-1} = \frac{1}{2}\begin{pmatrix} \frac{2}{N_1} & -\frac{1}{N_1} \\ 0 & \frac{1}{N_2} \end{pmatrix} \end{align*}\]

(Note $N_2$, or whichever direction is chosen to be the “row” or y-axis must be even.)

Our rectangular Fourier transform is therefore:

\[\begin{align*} X(k_1, k_2) = \sum_{n_1} \sum_{n_2} x(n_1, n_2) \exp\bigg(-2\pi i \bigg(\frac{(k_1(n_1 - n_2/2)}{N_1} + \frac{k_2 n_2}{N_2}\bigg) \bigg) \end{align*}\]

$n_1, n_2 \in R_R(N_1, N_2)$.

Rectangular region of support for a hexagonal FFT. See text below.

The figure above demonstrates why this periodicity matrix was chosen. A parallelogram in physical space is not the most convenient region of support for real-life signals. What we really want is a rectangular region. The region shaded in light blue above, $R_R(6, 8)$ can be re-sampled horizontally by moving the lower right triangular region onto the left side (shown by the overlap of the blue and pink), creating a rectangular region (in pink) whose periodicity is rectangular. One periodic extension of this re-sampled region is shown outlined in red.

So to take the hexagonal FFT of data on a physically rectangular (rather than coordinate-rectangular, which would physically be a parallelogram) region, the region can be “shifted” onto an equivalent parallelogram-shaped region of support, transformed in this “natural” coordinate system for the hexagonal grid, and then the result shifted back to obtain the FT on a rectangular grid. This is done automatically when the hexfft.fft() is called on an array in offset coordinates. (see this notebook for details).

The rect_shift and rect_unshift functions of hexfft.array shift a signal to/from rectangular vs parallelogram regions of support.

The algorithm for the rectangular FFT is quite straightforward. The kernel above is broken into two sequences of 1D FFTs with a phase shift. The following derivation is directly from Ehrhardt’s paper. Define:

\[\begin{align*} F_1(k_1, n_2) = \sum_{n_1} x(n_1, n_2) \exp(-2\pi i k_1 n_1 / N_1) \\ F_2(k_1, n_2) = F_1(k_1, n_2) \exp(\pi i k_1 n_2 / N_1) \end{align*}\]

Then

\[\begin{align*} X(k_1, k_2) = \sum_{n_1} F_2(k_1, n_2) \exp(-2\pi i k_2 n_2 / N_2) \end{align*}\]

The time complexity is thus $N_1 N_2 \log N_1 N_2$. This transform is implemented via hexfft.fft() and hexfft.ifft() with the periodicity="rect" option.

Conclusion

The above is only a (far from exhaustive) introduction to the basic concepts underlying hexagonal Fourier transforms. I have only barely scratched the surface of experimenting with data processing on hexagonal grids. The hodgepodge of tools in the hexfft package and a few examples of basic operations are demonstrated in the example notebooks. Future directions that come to mind are including support for coordinate systems that leverage hexagonal geometry more than the oblique system, designing discrete filters for hexagonal grids, and writing algorithms that can handle more general sampling strategies and periodicity conditions. More practically, I expect to make improvements to the basic API and documentation of the package. Another important step will be verification on hexagonal grids of classical properties of Fourier transforms such as Parseval’s identity and the shift and convolution theorems (and how those might need to be modified based on this geometry).