Dev/hexfft docs

docs/ref/hexfft.md

@@ -0,0 +1,105 @@
+---
+layout: page
+title: sm::hexfft
+parent: Reference
+nav_order: 11
+permalink: /ref/hexfft/
+---
+# sm::hexfft
+{: .no_toc}
+## The hexagonal fast Fourier transform
+{: .no_toc}
+```c++
+import sm.hexfft;
+```
+
+Module file: [sm/hexfft.cppm](https://github.com/sebsjames/maths/blob/main/sm/hexfft.cppm). Test code:
+[tests/hexfft1](https://github.com/sebsjames/maths/blob/main/tests/hexfft1.cpp)
+
+**Table of Contents**
+
+- TOC
+{:toc}
+
+## Summary
+
+The Hexagonal FFT algorithm, following Nicholas I. Rummelt's PhD thesis *Array set addressing: Enabling efficient hexagonally sampled image processing*, University of Florida, 2010.
+
+With this code, you can present data arranged over a hexgrid (specifically, a `sm::hexgrid<F, hexalign::point_up>`) with an arbitrary boundary, and compute the two dimensional spatial Fourier transform. A hexgrid in the frequency space (a `sm::hexgrid<F, hexalign::flat_up>`) is created alongside a data container with the FFT result in it.
+
+Defined as:
+```c++
+template<typename F = double, bool construct_hg_asa = false>
+struct fft
+```
+Floating point type `F` is used for hexgrid coordinates and as the element type for `std::complex<F>` values. `construct_hg_asa` is a boolean which may be set true to create an optional hexgrid that is useful for debugging, but not required for the forward or inverse transforms.
+
+The algorithm is made fast by splitting the hexgrid into two rectangular grids of alternating rows. For this reason, the input hexgrid must be enclosed by a perfect rectangular hexgrid (with zero-padding of new elements). Alternating rows are placed in two Array Set Addressing (ASA) grids, then the standard two dimensional, rectangular FFT can be applied.
+
+For a practical, visualized example implementation, you can see the [hex_fft](https://github.com/sebsjames/hex_fft) repository.
+
+## Quick usage guide
+
+Create a hexgrid. The hexgrid constructor args are hex-hex distance, grid width and grid 'z' value (usually set to 0).
+
+```c++
+import sm.hexgrid;
+
+sm::hexgrid<float sm::hexalign::point_up> hg(0.01f, 4.0f, 0.0f);
+hg.set_circular_boundary (1.0f); // or any other boundary setting function in hexgrid
+```
+
+Create some data. The order of the data is defined by the hexgrid indexing. Each hexgrid element has a 'vector iterator', `vi` and provides access to the location of the hex.
+```c++
+import sm.vvec;
+
+sm::vvec<float> data (hg.num(), 0.0f);
+for (auto h : hg.hexen) {
+    data[h.vi] = some_function_of (h.x, h.y);
+}
+```
+
+Create an `sm::hexfft::fft` object and perform a forward transform. The result is stored in `hfft.X_hexgrid`, which is a `sm::vvec` of `std::complex<>` values.
+
+```c++
+import sm.hexfft;
+
+sm::hexfft::fft<float> hfft (&hg); // construct and initialize
+hfft.forward (data); // Perform forward FFT transform
+```
+You can modify the values in `X_hexgrid` to make filters. The values in `X_hexgrid` are associated with a frequency hexgrid, `hexfft::fft::hgf`, which is created when hfft is initialized.
+```c++
+for (auto h : hfft.hgf->hexen) { // hgf is a unique_ptr to a hexgrid
+    std::cout << "FFT Frequency " << h.x << ", " << h.y
+              << " has magnitude " << std::real(hfft.X_hexgrid[h.vi]) << std::endl;
+}
+
+```
+After changing values in `hfft.X_hexgrid` (perhaps by masking) you can then inverse transform from frequency space to image space
+
+```c++
+sm::vvec<std::complex<float>> invimg = hfft.inverse();
+```
+The returned data is defined over your original hexgrid, `hg`.
+
+## fft members
+
+### Attributes populated during initialization
+
+`hexfft::asa_rows` and `hexfft::asa_cols` (both `uint32_t` are populated with the dimensions of the two ASA grids.
+
+`hexfft::ri_min` and `hexfft::gi_min` are configured with the `hex::ri`, `hex::gi` of the padded rectangular hexgrid's (a=0, r=0, c=0) corner.
+
+The hexgrid pointer `hexfft::hg` is the point your provide at construction/init.
+
+The hexgrid `hexfft::hgf` is constructed to match `hg` for the frequency space. It has alignment `hexalign::flat_up`.
+
+`hexfft::Uscale` is a scaling factor (computed from image data hexgrid spacing) for the frequency hexgrid, hgf. Uscale allows the user to visualize the frequency hexgrid on a similar size scale to the input data, regardless of the spacing on the input data.
+
+### Data attributes
+
+`hexfft::d0` and `hexfft::d1` hold copies of the input data (or inverse-transformed data) in ASA format. d0 holds even rows, d1, odd rows. They are both `sm::vmat` containers of `std::complex<F>` values.
+
+`hexfft::X0` and `hexfft::X1` are similar containers that hold the (ASA format) result of the forward transformation of `d0` and `d1`.
+
+The data in `X0` and `X1` are manipulated, and then rearranged into a `vvec` of `std::complex<F>` values, `hexfft::X_hexgrid`. `X_hexgrid` is spatially defined by the frequency hexgrid `hexfft::hgf`.

docs/ref/hexgrid.md

@@ -24,50 +24,56 @@ Module file: [sm/hexgrid.cppm](https://github.com/sebsjames/maths/blob/main/sm/h
 
 ## Summary
 
-`sm::hexgrid` is a hexagonal-tiling counterpart to [`sm::grid`](/maths/ref/grid/): rather than using rectangular elements, it lays out a grid of hexagons (each an `sm::hex`) to manage spatial information for an associated computation.
-It was designed for a study of [reaction-diffusion systems across two dimensional domains](https://elifesciences.org/articles/55588). `std::vector` arrays held the system state variables, and the spatial information for each element was managed in the hexgrid. Here, it was useful to use a hexagonal grid, because this made the [computation of the Laplacian easy](https://elifesciences.org/articles/55588#s4).
+`sm::hexgrid` is a counterpart to the Cartesian [`sm::grid`](/maths/ref/grid/). It lays out a grid of hexagons (each an `sm::hex`) to manage spatial information for an associated computation.
+It was originally designed for a study of [reaction-diffusion systems across two dimensional domains](https://elifesciences.org/articles/55588). `std::vector` arrays held the system state variables, and the spatial information for each element was managed in the hexgrid. Use of a hexagonal grid was motivated by the ease of [computation of the Laplacian easy](https://elifesciences.org/articles/55588#s4).
 
-The design of `hexgrid` differs from that of `sm::grid`, being more similar to [`sm::cartgrid`](/maths/ref/cartgrid/). `sm::hexgrid` defines an initial hexagonal grid of hexagons which you can then clip to an arbitrary boundary in exactly the same spirit as `cartgrid` clips its rectangular lattice.
-Where `cartgrid` and `sm::grid` share three boundary/wrap-related enums, `hexgrid` predates that enum-based design; its boundary shape is chosen by which `set_*_boundary` method you call, rather than by setting a `domain_shape` member. `hexgrid` was written before `cartgrid`, which predated `grid`.
+`sm::hexgrid` defines an initial hexagonal grid of hexagons, built outwards ring-by-ring from a single centre hex until it reaches a requested diameter.
+You can use the grid with this default hexagonal shape, or clip it down to have another boundary shape (such as a circle, ellipse, rectangle or parallelogram). You can also set a boundary you supply yourself as a closed [`sm::bezcurvepath`](https://github.com/sebsjames/maths/blob/main/sm/bezcurvepath.cppm) or list of points.
+Clipping discards every hex outside the boundary and re-links the neighbour relationships of those that remain.
 
+`sm::hexgrid` is a templated class with a coordinate type `F` (a floating point type) and a `sm::hexalign` template parameter `A` which defines whether the hexagonal lattice is arranged in an orientation for which individual hexes have their 'points up' or their 'flats up'.
 
-`sm::hexgrid` always starts out as a  filled hexagon of hexagonal elements, built outwards ring-by-ring from a single centre hex until it reaches a requested diameter. As with `cartgrid`, you can leave the grid as this default hexagonal shape, or clip it down to an arbitrary boundary; a circle, ellipse, rectangle, parallelogram, or a boundary you supply yourself as a closed [`sm::bezcurvepath`](https://github.com/sebsjames/maths/blob/main/sm/bezcurvepath.cppm) or list of points. Clipping discards every hex outside the boundary and re-links the neighbour relationships of what remains.
+![Two hexagonal grids](https://github.com/sebsjames/maths/blob/main/docs/images/hexgrids.png?raw=true)
+*point_up and flat_up hexgrids. The hex elements of the point_up grid have their points up; the overall hexagonal shape is opposite.*
 
-`sm::hexgrid` is a non-templated class with coordinates of type `float` and indices of type `std::uint32_t` (sometimes cast to `std::int32_t`). It does not derive from, or share any types with, `sm::grid`; it does, however, share a very similar design and method-naming convention with `sm::cartgrid`. `hexgrid` was designed first, then `cartgrid` was coded up using the same ideas. (The `set_boundary`/`set_boundary_only`/`set_boundary_on_outer_edge`/`get_region` family, and the flat `d_*` cache-vector convention, are essentially the same functions applied to hexes instead of rects.)
+Indices have type `std::uint32_t` (sometimes cast to `std::int32_t`).
+
+It does not derive from, or share any types with, `sm::grid`; it does, however, share a very similar design and method-naming convention with `sm::cartgrid`.
 
 Defined as:
 ```c++
 export namespace sm
 {
-    class alignas(8) hexgrid
+    template<typename F, sm::hexalign A = sm::hexalign::point_up> requires std::is_floating_point_v<F>
+    struct alignas(8) hexgrid
     {
         // ...
         std::list<hex> hexen;
 ```
 
 ### Hex coordinates
 
-Each `sm::hex` (defined in `sm/hex.cppm`, and re-exported by `sm.hexgrid`, so `import sm.hexgrid;` is enough to use it) stores an axial coordinate `{ri, gi, bi}` alongside its Cartesian `{x, y, z}` position, a 32-bit `flags` word (`HEX_IS_BOUNDARY`, `HEX_INSIDE_BOUNDARY`, `HEX_INSIDE_DOMAIN`, `HEX_IS_REGION_BOUNDARY`, `HEX_INSIDE_REGION`, plus 16 bits reserved for your own use as `HEX_USER_FLAG_0`..`HEX_USER_FLAG_15`), and six neighbour iterators (`ne`, `nne`, `nnw`, `nw`, `nsw`, `nse`; see [Neighbours](#neighbours-in-the-six-hex-directions)). The hexes are 'point-up', spaced `d` apart within a row and `v = d * sqrt(3)/2` apart between rows.
+Each `sm::hex` (defined in `sm/hex.cppm`, and re-exported by `sm.hexgrid`, so `import sm.hexgrid;` is enough to use it) stores an axial coordinate `{ri, gi, bi}` alongside its Cartesian `{x, y, z}` position, a 32-bit `flags` word (`HEX_IS_BOUNDARY`, `HEX_INSIDE_BOUNDARY`, `HEX_INSIDE_DOMAIN`, `HEX_IS_REGION_BOUNDARY`, `HEX_INSIDE_REGION`, plus 16 bits reserved for your own use as `HEX_USER_FLAG_0`..`HEX_USER_FLAG_15`), and six neighbour iterators (`n0`, `n1`, `n2`, `n3`, `n4`, `n5`; see [Neighbours](#neighbours-in-the-six-hex-directions)). The hexes are 'point-up' by default, spaced `d` apart within a row and `v = d * sqrt(3)/2` apart between rows.
 
 ## Create a hexgrid
 
 ```c++
-sm::hexgrid hg (0.01f, 3.0f, 0.0f); // d (hex spacing), x_span (diameter), z (layer)
+sm::hexgrid<float> hg (0.01f, 3.0f, 0.0f); // d (hex spacing), x_span (diameter), z (layer)
 ```
 This builds a full hexagon of hexes with hex-to-hex spacing `d` and a horizontal diameter of approximately `x_span`. `init (d_, x_span_, z_)` re-runs the same construction on an existing `hexgrid`, and the default constructor `hexgrid()` leaves `d = x_span = 1.0f` but, like `cartgrid`'s default constructor, does not build the grid for you.
 
 ## Setting a boundary
 
 The convenience methods compute the boundary points for a given shape and clip the grid to them in one call:
 ```c++
-sm::hexgrid hg (0.01f, 3.0f, 0.0f);
+sm::hexgrid<float> hg (0.01f, 3.0f, 0.0f);
 hg.set_circular_boundary (0.6f);         // radius 0.6, centred at the origin by default
 std::cout << "Number of hexes in grid: " << hg.num() << std::endl;
 ```
 `set_elliptical_boundary`, `set_rectangular_boundary` and `set_parallelogram_boundary` work the same way for their respective shapes. For an arbitrary shape, supply a closed Bezier path:
 ```c++
 sm::bezcurvepath<float, 3> bound = /* ... four curve segments forming a closed loop ... */;
-auto hgrid = std::make_unique<sm::hexgrid> (0.02f, 4.0f, 0.0f);
+auto hgrid = std::make_unique<sm::hexgrid<float>> (0.02f, 4.0f, 0.0f);
 hgrid->set_boundary (bound);
 std::cout << "Number of hexes is: " << hgrid->num() << std::endl;
 ```
@@ -81,6 +87,13 @@ And if you just want to *mark* a boundary for inspection without discarding any
 
 `get_boundary()` returns a copy of the current boundary hexes, and `compute_distance_to_boundary()` fills each hex's `dist_to_boundary` (`0` on the boundary itself, `-100.0f` for any hex outside the boundary, otherwise the distance to the nearest boundary hex).
 
+There is a special `set_rectangular_boundary` function for making Array Set Addressing (ASA) compatible grids:
+```c++
+void set_rectangular_boundary (const std::uint32_t n_x, const std::uint32_t n_y,
+                               const std::uint32_t x0 = 0u, const std::uint32_t y0 = 0u)
+```
+This takes integer arguments for the number of rows and columns, and carefully arranges the rows so that the resulting grids for `point_up` and `flat_up` hexgrids can be used, respectively, as the image lattice and corresponding frequency lattice in a hexagonal FFT computation.
+
 ### Temporary regions
 
 As with `cartgrid`, `get_region` (given a Bezier path or point vector) and `get_hexagonal_region` (given a centre hex index and radius) let you mark and retrieve a sub-set of hexes without discarding anything from the grid; useful for querying, e.g., "which hexes fall within this circle" while leaving the grid itself untouched:
@@ -103,25 +116,29 @@ auto at_axial = hg.find_hex_at ({ 2, -1, 0 }); // {ri, gi, bi}
 
 ## Neighbours in the six hex directions
 
-Each hex has up to six neighbours; East, North-East, North-West, West, South-West and South-East; reachable either through the hex object's own iterators, or via a flat domain index:
+Each hex has up to six neighbours, referred to with iterators `hex<>::n0` to `hex<>::n5`.
+
+In a `point_up` lattice, these are East, North-East, North-West, West, South-West and South-East neighbours.
+
+In a `flat_up` lattice they are North-East, North, North-West, South-West, South and South-East neigbours.
+
+Nieghbour existence can be tested with `has_n0` to `has_n5` methods, and accessed via hex iterators or `hexgrid::n0()` to `hexgrid::n5()` methods.
+
 ```c++
 auto hi = hg.hexen.begin();
-if (hi->has_ne()) {
-    auto east_neighbour = hi->ne; // std::list<hex>::iterator
+if (hi->has_n0()) {
+    auto east_neighbour = hi->n0; // std::list<hex>::iterator
 }
 // or, given a flat domain index `di` (hex::di):
-if (hg.has_ne (di)) {
-    std::int32_t east_di = hg.ne (di); // -1 if there's no such neighbour
+if (hg.has_n0 (di)) {
+    std::int32_t east_di = hg.n0 (di); // -1 if there's no such neighbour
 }
 ```
-**Note:** the domain-index `has_ne`/`has_nw`/`has_nne`/`has_nnw`/`has_nse`/`has_nsw` functions return `std::int32_t`, not `bool`, even though they behave as a boolean presence check (they evaluate to `0` or `1`).
-
-## Wrapping
-
-There's no general wrap enum for `hexgrid`; the only wrapping support is `set_parallelogram_wrap (bool on_r, bool on_g)`, which re-wires the neighbour links at the edges of a parallelogram-shaped domain to point at the opposite edge. **At present it only supports wrapping both axes together**; it throws `std::runtime_error` unless both `on_r` and `on_g` are `true`.
 
 ## Convolution, resampling and shifting data
 
+Methods in the namespace `sm::algo::hexgrid` (module file [algo_hexgrid.cppm](https://github.com/sebsjames/maths/blob/main/sm/algo_hexgrid.cppm)) provide hexgrid-compatible algorithms.
+
 `convolve` performs a 2D convolution of per-hex data against a kernel defined on a second `hexgrid` (which must share the same `d`), walking neighbour links rather than assuming a fixed array stride, so it works correctly on boundary-clipped domains. `resample_image` Gaussian-resamples a rectangular pixel image onto the hex centres, much like the equivalent methods in `sm::grid` and `sm::cartgrid`.
 
 `shiftdata` translates per-hex data by an arbitrary Cartesian vector, splitting the shift into whole hex-hops (following neighbour links, so any wrapping you've set up is respected) plus a sub-hex remainder distributed by exact hex-overlap-area weighting:
@@ -131,6 +148,8 @@ bool ok = hg.shiftdata (image_data, sm::vec<float, 2>{ 0.003f, -0.001f });
 ```
 It returns `false` (leaving `image_data` unmodified) if the overlap geometry couldn't be resolved for the given shift. The `compute_hex_overlap`/`compute_overlap_*`/`setup_hexoverlap_geometry` methods it relies on are public, but are internal machinery for `shiftdata`; you shouldn't normally need to call them directly.
 
+There are also some masking functions in algo_hexgrid.cppm.
+
 ## Geometry
 
 ```c++
@@ -146,7 +165,7 @@ float area = hg.get_hex_area(); // area of one hex
 
 ## Saving and loading
 
-HDF5 persistence lives in a separate module, `sm.hexgrid.hdf` (which re-exports `sm.hexgrid`, so importing it gives you everything above too):
+HDF5 persistence lives in a separate module, `sm.hexgrid.hdf`.
 ```c++
 import sm.hexgrid.hdf;
 
@@ -155,6 +174,6 @@ sm::hexgrid_save (hg, "myhexgrid.h5");
 sm::hexgrid hg2;
 sm::hexgrid_load (hg2, "myhexgrid.h5");
 ```
-Loading reconstructs each hex's six neighbour relationships by matching saved indices against the freshly-loaded hex list; an O(n²) operation for large grids; and throws `std::runtime_error` if any expected neighbour can't be matched. The boundary curve itself (as a `bezcurvepath`) is not saved; only the resulting hex positions, flags and neighbour relationships are.
+Loading a hexgrid reconstructs each hex's six neighbour relationships by matching saved indices against the freshly-loaded hex list; an O(n²) operation for large grids.  It throws `std::runtime_error` if any expected neighbour can't be matched. The boundary curve itself (as a `bezcurvepath`) is not saved; only the resulting hex positions, flags and neighbour relationships are.
 
 *This page was authored with AI, based on human written code in hexgrid.cppm and reviewed by Seb James.*

docs/ref/hexyhisto.md

@@ -31,7 +31,7 @@ Module file: [sm/hexyhisto.cppm](https://github.com/sebsjames/maths/blob/main/sm
 Here's an example where we create a circular `sm::hexgrid`, then
 
 ```c++
-sm::hexgrid hg (0.1f, 2.0f, 0.0f);
+sm::hexgrid<float> hg (0.1f, 2.0f, 0.0f);
 hg.set_circular_boundary (0.5f);
 
 sm::vvec<sm::vec<float>> data; // sm::vec<float> defaults to 3 elements: {x, y, flag}
@@ -52,6 +52,6 @@ T total = hh.datacount;         // how many input points were actually counted
 sm::vvec<T> counts = hh.counts; // raw count per hex, indexed by each hex's vi
 sm::vvec<T> proportions = hh.proportions; // counts, normalized to sum to 1
 ```
-`proportions` is exactly what you'd plot on the `sm::hexgrid` to visualize the histogram as a density map.
+`proportions` is exactly what you'd plot on the `sm::hexgrid<>` to visualize the histogram as a density map.
 
 *This page was authored with AI, based on human written code in hexyhisto.cppm and reviewed by Seb James*