Tutorial: your own filtration in Rust

This tutorial builds persistent homology for a filtration phlite has never seen: the cubical complex of a 3D greyscale volume (a CT scan, a simulation, a 3D image), with the lower-star filtration: a cell appears with its brightest voxel. It starts with the least you can implement, and each step adds a piece that unlocks an optimisation.

When is a lazy oracle the right tool?

phlite’s oracles never store the matrix, and that pays off when most of the matrix never needs touching, not merely when it wouldn’t fit in memory:

• Most pairs are zero-length apparent pairs. In a lower-star filtration nearly every cell is paired with a neighbour that appears at the same moment; in Vietoris–Rips nearly every simplex is. Such cells needn’t be reduced, recorded, or even listed (apparent pairs).
• The top dimension needn’t be enumerated. Cohomology only meets the top cells as rows, produced on demand.
• Columns are cheap to compute from a compact label, so recomputing them costs less than storing them.

Vietoris–Rips complexes often fit in memory, and storing them is still the wrong choice: they are huge, sparse and almost entirely apparent pairs. The opposite case is a complex where many columns need genuine reduction: dense boundaries, or cubical complexes of pure noise. There the work is in adding columns, and PHAT’s tricks (dense indices, pivot columns, storing R) win. For those, build the complex explicitly with the phlite Python package (which applies all of them), or use PHAT.

A smooth 3D volume is in the lazy camp: in the example below, 94.5% of the cells are in zero-length apparent pairs.

The code is two examples in the phlite crate, sharing the volume and its geometry:

cargo run --release -p phlite --example volume_basic                       # steps 1–4
cargo run --release -p phlite --example volume_lazy --features parallel    # steps 5–7

The complex

Voxels are vertices. Between two adjacent voxels is an edge, between four a square, between eight a cube. It is convenient to put every cell on a grid twice as fine as the volume: a coordinate that is odd spans two voxels, so a cell’s dimension is its number of odd coordinates. Here is the same idea in 2D:

In 2D: a vertex where both coordinates are even, an edge where one is odd, a square where both are. The volume adds a third coordinate.

Numbering the grid point by point gives every cell a u64. A cell’s faces are one step either way along each odd coordinate, and its cofaces one step either way along each even coordinate, so both are computed from the number alone. Nothing about the complex is stored but the voxels:

/// Voxels are vertices. The cells sit at the points of a grid of side `2 size - 1`, at coordinates `(x, y, z)`:
/// a coordinate that is odd spans two voxels, so a cell's dimension is its number of odd coordinates (a vertex, an
/// edge, a square or a cube). Cells are numbered `(x · side + y) · side + z`.
impl Volume {
    pub fn side(&self) -> u64 {
        2 * self.size as u64 - 1
    }

    pub fn n_cells(&self) -> u64 {
        self.side().pow(3)
    }

    fn coordinates(&self, cell: u64) -> [u64; 3] {
        let side = self.side();
        [cell / (side * side), (cell / side) % side, cell % side]
    }

    fn number(&self, [x, y, z]: [u64; 3]) -> u64 {
        let side = self.side();
        (x * side + y) * side + z
    }

    pub fn dimension(&self, cell: u64) -> usize {
        self.coordinates(cell).iter().filter(|&&c| c % 2 == 1).count()
    }

    /// The faces of a cell: one step either way along each odd coordinate (at most six).
    pub fn faces(&self, cell: u64) -> impl Iterator<Item = u64> {
        let coords = self.coordinates(cell);
        let mut faces = [None; 6];
        for axis in 0..3 {
            if coords[axis] % 2 == 1 {
                for (k, step) in [-1i64, 1].into_iter().enumerate() {
                    let mut face = coords;
                    face[axis] = (face[axis] as i64 + step) as u64;
                    faces[2 * axis + k] = Some(self.number(face));
                }
            }
        }
        faces.into_iter().flatten()
    }

    /// The cofaces of a cell: one step either way along each even coordinate, inside the grid (at most six), in
    /// increasing order of their numbers: the steps down (along x, then y, then z), then the steps up (z, y, x).
    pub fn cofaces(&self, cell: u64) -> impl Iterator<Item = u64> {
        let coords = self.coordinates(cell);
        let side = self.side();
        let mut cofaces = [None; 6];
        for axis in 0..3 {
            if coords[axis] % 2 == 0 {
                if coords[axis] > 0 {
                    let mut down = coords;
                    down[axis] -= 1;
                    cofaces[axis] = Some(self.number(down));
                }
                if coords[axis] + 1 < side {
                    let mut up = coords;
                    up[axis] += 1;
                    cofaces[5 - axis] = Some(self.number(up));
                }
            }
        }
        cofaces.into_iter().flatten()
    }

    /// The lower-star filtration: a cell appears with its brightest voxel (of the up to eight it spans).
    pub fn value(&self, cell: u64) -> f32 {
        let [x, y, z] = self.coordinates(cell);
        let n = self.size as u64;
        let mut value = f32::NEG_INFINITY;
        for i in [x / 2, x.div_ceil(2)] {
            for j in [y / 2, y.div_ceil(2)] {
                for k in [z / 2, z.div_ceil(2)] {
                    value = value.max(self.voxels[((i * n + j) * n + k) as usize]);
                }
            }
        }
        value
    }
}

The cofaces are produced in increasing order of their numbers, which step 5 puts to use.