Spherical Harmonic Transforms on Ring-based Pixelizations

In two previous articles, I described the details of performing spherical harmonic synthesis and analysis using an equidistant cylindrical projection (ECP) pixelization scheme as the assumed map format. In this article, I’ll describe the generalizations to the algorithms which permit supporting the transform of any “ring-based” pixelization of the sphere.

Introduction¶

The equidistant cylindrical projects (ECP) pixelization is maybe the most intuitive choice of pixelizing the sphere since it corresponds directly to the familiar longitude-latitude system we are taught in school. This made the ECP grid a very convenient choice of pixelization to use when developing the spherical harmonic synthesis and analysis algorithms in previous articles since we could concentrate on the transforms themselves and not have to discuss various details of how to build the pixelization.

In analytical and numerical terms, though, the ECP grid is not an ideal choice. We saw one example of this fact in the discussion of weighted-quadrature analysis, where we observed that moving the colatitude nodes from the uniformly-spaced ECP lines of latitude to the (slightly different) Gauss-Legendre nodes allows for a much more accurate recovery of a map’s spherical harmonic coefficients.

In using the Gauss-Legendre (GL) grid, we were forced to rewrite the synthesis and analysis functions specifically for the GL pixelization since the prior implementations were all specialized for the ECP grid. The purpose of this article is to show that with a little forethought, we can write a pair of analysis and synthesis routines which are largely agnostic to the specific pixelization grid;¹ the only requirement is that the pixelization is compatible with the symmetries that have been exploited to create these fast transforms, and we’ll generically refer to such compatible pixelizations as “ring-based” pixelizations.

The next section will review the necessary assumptions and symmetries that need to be maintained, and then the rest of the article will demonstrate updating the previously-provided Julia codes to implement the more generic transform capabilities.

Necessary Properties of Ring-based Pixelizations¶

As the name “ring-based pixelizations” suggests, the first obvious requirement for supported pixelization schemes is that the format be describable in terms of isolatitude rings. This was the first requirement we imposed in the first article because it permits greatly reducing the number of times the associated Legendre polynomials must be calculated — from (in principle) every pixel if they all have a unique latitude to only once per isolatitude ring. Therefore, the first per-ring parameter we must be told is the colatitude $\theta$ of each ring.

The second requirement is that each isolatitude ring’s pixels meet the necessary conditions for the Fast Fourier Transform (FFT) to be used — namely, the pixels in azimuth are uniformly spaced and span (exactly) a full $2\pi$ radians. These conditions are met when the azimuth coordinates are at locations $\phi_k = \phi_0 + 2\pi k / N_\phi$ (for integer index $k = [0, 1, \ldots, N_\phi - 1]$), and therefore two additional parameters must be given: the number of pixels within the ring $N_\phi$ and the azimuth offset of the first pixel in the ring $\phi_0$.

These three parameters ($\theta$, $N_\phi$, $\phi_0$) are sufficient for implementing the synthesis transform of a single isolatitude ring. For the analysis transform, though, we also require the integration volume element $d\Omega \equiv \sin\theta , d\theta , d\phi$, or rather its finite approximation $\Delta\Omega \approx \sin\theta , \Delta\theta , \Delta\phi$. It is tempting to assume $\Delta\phi = 2\pi / N_\phi$ and so the only remaining unknown is the ring spacing $\Delta\theta$, but doing so implies that each pixel is approximated well by a rectangle.

The trivial case of a 2×4 ECP grid shows how this is a poor assumption — each pixel covers one octant of the surface of the sphere, so the shape of each pixel is best described as a spherical triangle. We know a priori that the surface area of each pixel is $4\pi/8 \approx 1.571$ steradians by geometric construction, but the naive calculation using $\theta = \pi/4$, $\Delta\theta = \pi/2$, and $\Delta\phi = 2\pi/4$ instead calculates $\Delta\Omega = 1.745$ steradians.

The solution to this problem is to take the volume element itself $\Delta\Omega$ as the necessary parameter. Given this approach, we never explicitly calculate the quantity $\sin\theta$, using only $\cos\theta$ as the argument of the associated Legendre polynomials, so we can actually also revise the first parameter to instead be the cosine of the colatitude angle $\cos\theta$.

Together, this gives us 4 parameters which are required to sufficiently describe the geometry of an isolatitude ring so that we can perform spherical harmonic transforms on it:

1. The cosine of the ring’s colatitude angle: $\cos\theta$
2. The azimuth angle of the first pixel in the ring: $\phi_0$
3. The number of pixels equally spaced in azimuth within the ring: $N_\phi$
4. The surface area of each pixel within the ring: $\Delta\Omega$

What remains is to describe how multiple rings are assembled into a map that covers the sphere, and here again the ECP grid proves a useful case study. The ECP grid maps well to a matrix (2D array, where colatitude increases across rows and azimuth across columns), but there is an inherent ambiguity in whether the ordering of the matrix in memory is row-major or column-major. If we choose one — say column-major order, consistent with Julia’s matrix representation — then the contents of a given isolatitude ring are non-contiguous in memory. One option would be to require that the transpose of the ECP map matrix be used (so that a ring is a contiguous vector of elements), but another more flexible option is to instead make the stride $s$ of the array an additional ring description parameter. When paired with offset of the first pixel in the ring $o$, we have sufficient information to completely describe an ECP grid in either row- or column-major order.

As one final optimization, though, recall that if the pixel grid is symmetric over the equator, then rings are best considered in pairs, since the associated Legendre polynomials are the same up to a negative sign for colatitude angles $\theta$ and $\pi - \theta$. If we replace the single parameter $o$ with a tuple of parameters $(o_1, o_2)$, we can inform our implementation explicitly whether the ring-pair optimization applies (such as for a north-south pair of rings) or not (like for the singular equatorial row) by encoding either two or one valid indices³ in the tuple, respectively.

Finally, we add these three additional parameters to the above list:

5. Offset to the first pixel in the northern- and southern-hemisphere rings, respectively (if valid): $(o_1, o_2)$
6. The stride across the map array between successive pixels within the ring: $s$

In the following section, we’ll translate this list of parameters into Julia data structures.

Ring-based Map Description Structures¶

Having identified all of the necessary information in the previous section, the structure definition to describe each ring is relatively straight-foward:

struct RingInfo
    offsets::Tuple{Int,Int}
    stride::Int
    nϕ::Int
    cosθ::Float64
    ϕ₀::Float64
    ΔΩ::Float64
end

We then must collect a set of ring descriptions that together describe a particular map pixelization. Here, we choose⁴ a simple container type that holds the Vector{RingInfo}:

struct MapRings
    rings::Vector{RingInfo}
end

With all of the types necessary to define a ring-based map pixelization in place, the next obvious step is to actually instantiate a particular pixelization. Sticking with our standard case, the description of the ECP grid is given by:

"""
    maprings = ecp_maprings(nθ, nϕ)

Constructs a `MapRings` structure describing an equidistant cylindrical projection (ECP) grid of dimensions
`nθ × nϕ` pixels in column-major ordering.
"""
function ecp_maprings(nθ, nϕ)
    nθh = (nθ + 1) ÷ 2  # number of rings in N. hemisphere + equator (if nθ is odd)
    ϕ₀ = 2π / 2nϕ  # constant for all rings
    ΔθΔϕ = (2π / nϕ) * (π / nθ)  # constant portion of ΔΩ = sinθ * ΔθΔϕ
    stride = nθ  # for column-major ordering
    rings = Vector{RingInfo}(undef, nθh)
    for ii in 1:nθh
        o₁ = ii           # northern-ring offset
        o₂ = nθ - ii + 1  # southern-ring offset
        offsets = (o₁, o₁ == o₂ ? 0 : o₂)  # for equator of odd nθ, disable southern ring
        θ_π = (ii - 0.5) / nθ  # θ/π for the ring
        sθ, cθ = sincospi(θ_π)
        ΔΩ = sθ * ΔθΔϕ
        rings[ii] = RingInfo(offsets, stride, nϕ, cθ, ϕ₀, ΔΩ)
    end
    return MapRings(rings)
end

On first glance, this seems to require more effort to describe such a simple grid, but remember that the goal is to describe a variety of pixelization systems with a shared infrastructure.

A good example of the simplification this extra abstraction affords is to also define the ring structure for the Gauss-Legendre (GL) variant grid which was used at the end of Part II. Recall that the GL pixelization is the same as ECP in the azimuthal ($\phi$) direction, but the locations of the rings in $\theta$ are instead the Gauss-Legendre nodes. As an added benefit, note that the associated quadrature weights are effectively the quantity $\sin\theta,\Delta\theta$, so they naturally fit into our ring description without requiring any new parameters.

The implementation that gives a description of the GL grid is:

using FastGaussQuadrature: gausslegendre

"""
    maprings = gl_maprings(nθ, nϕ)

Constructs a `MapRings` structure describing an Gauss-Legendre (GL) grid of dimensions `nθ × nϕ` pixels in
column-major ordering.
"""
function gl_maprings(nθ, nϕ)
    nθh = (nθ + 1) ÷ 2  # number of rings in N. hemisphere + equator (if nθ is odd)
    ϕ₀, Δϕ = 2π / 2nϕ, 2π / nϕ  # constant for all rings
    cθ, sθΔθ = gausslegendre(nθ) # z = cos(θ) and corresponding weights
    stride = nθ  # for column-major ordering
    rings = Vector{RingInfo}(undef, nθh)
    for ii in 1:nθh
        o₁ = ii           # northern-ring offset
        o₂ = nθ - ii + 1  # southern-ring offset
        offsets = (o₁, o₁ == o₂ ? 0 : o₂)  # for equator of odd nθ, disable southern ring
        # N.B. z ∈ [-1, 1] and θ ∈ [0, π] run in opposite directions, and cos(π-θ) == -z
        rings[ii] = RingInfo(offsets, stride, nϕ, -cθ[ii], ϕ₀, sθΔθ[ii] * Δϕ)
    end
    return MapRings(rings)
end

The differences in implementation between ecp_maprings and gl_maprings is rather minimal. Explicit calculation of the coordinates $\theta$ and pixel areas $\Delta\Omega$ are replaced by the outputs from the gausslegendre() function, but otherwise the structure of the rings is highly similar.

Generic Ring-based Synthesis and Analysis¶

We will implement the generic synthesis and analysis functions by taking the ECP-specific functions written in the previous two articles and modifying them to use the MapRings and RingInfo structures. For convenience, the synthesize_ecp and analyze_ecp functions are included within the following hidden section, with minor clarifications and cleanups applied.

using AssociatedLegendrePolynomials: λlm!
using FFTW

"""
    (i′, isconj, isnyq) alias_index(len::Int, i::Int)

Calculates aliased index `i′` given original index `i` and Fourier axis length 'len'.
Also indicates whether conjugation (`isconj`) or twice-real symmetries (`isnyq`) are needed.
"""
@inline function alias_index(len::Int, i::Int)
    isconj, isnyq = false, false
    nyq = max(1, len ÷ 2)
    i < nyq && return (i, isconj, isnyq)
    i = mod(i, len)
    if i > nyq
        i = len - i
        isconj = true
    elseif i == 0 || (iseven(len) && i == nyq)
        isnyq = true
    end
    return (i, isconj, isnyq)
end

"""
    coeffs′ = alias_coeffs(coeffs::NTuple, isconj::Bool, isnyq::Bool)

Applies conjugation (`isconj`) or twice-real (`isnyq`) transformation, if required, in mode aliasing to the
harmonic coefficients `coeffs`.
"""
@inline function alias_coeffs(coeffs::NTuple, isconj::Bool, isnyq::Bool)
    return ifelse(isnyq, complex.(2 .* real.(coeffs)),
                  ifelse(isconj, conj.(coeffs), coeffs))
end

function synthesize_ecp(alms::Matrix{C}, nθ, nϕ) where {C<:Complex}
    R = real(C)
    lmax, mmax = size(alms) .- 1
    nϕr = nϕ ÷ 2 + 1 # real-symmetric FFT's Nyquist length (index)
    nθh = (nθ + 1) ÷ 2 # number of rings in northern hemisphere

    # pixel grid definition for ECP
    θr = range(R(π)/2nθ, R(π), step = R(π)/nθ)
    ϕ₀ = R(π) / nϕ

    ecp = Matrix{R}(undef, nθ, nϕ)
    Λ = Matrix{R}(undef, lmax + 1, mmax + 1)
    λ₁ = zeros(C, nϕr) # phase factors for northern ring
    λ₂ = zeros(C, nϕr) # phase factors for southern ring

    for j in 1:nθh
        j′ = nθ - j + 1
        θ = θr[j]
        λlm!(Λ, lmax, mmax, cos(θ))

        for m in 0:mmax
            acc₁, acc₂ = zero(C), zero(C)
            # Σ_{ℓ = m}^{lmax}
            for ℓ in m:lmax
                term = alms[ℓ+1,m+1] * Λ[ℓ+1,m+1]
                acc₁ += term
                acc₂ += isodd(ℓ + m) ? -term : term
            end
            acc₁, acc₂ = (acc₁, acc₂) .* cis(m * ϕ₀) # ring offset rotation
            i, isconj, isnyq = alias_index(nϕ, m)
            acc₁, acc₂ = alias_coeffs((acc₁, acc₂), isconj, isnyq)
            λ₁[i+1] += acc₁
            λ₂[i+1] += acc₂
        end

        ecp[j,:] = brfft(λ₁, nϕ)
        λ₁ .= zero(C)
        ecp[j′,:] = brfft(λ₂, nϕ)
        λ₂ .= zero(C)
    end
    return ecp
end

function analyze_ecp(map::Matrix{R}, lmax::Integer, mmax::Integer = lmax) where {R<:Real}
    C = complex(R)
    nθ, nϕ = size(map)
    nϕr = nϕ ÷ 2 + 1 # real-symmetric FFT's Nyquist length (index)
    nθh = (nθ + 1) ÷ 2 # number of rings in northern hemisphere

    # pixel grid definition for ECP
    θr = range(R(π)/2nθ, R(π), step = R(π)/nθ)
    ϕ₀ = R(π) / nϕ
    ΔΩ = 2R(π)^2 / (nθ * nϕ)

    alms = fill(zero(C), lmax + 1, mmax + 1)
    Λ = zeros(R, lmax + 1, mmax + 1)

    for j in 1:nθh
        j′ = nθ - j + 1 # southern ring index
        sθ, cθ = sincos(θr[j])
        λlm!(Λ, lmax, mmax, cθ) # λ_ℓ^m(cos θ) factors

        h₁ = rfft(map[j,:])  # phase factors for northern ring
        h₂ = rfft(map[j′,:]) # phase factors for southern ring

        sθΔΩ = sθ * ΔΩ
        for m in 0:mmax
            i, isconj, isnyq = alias_index(nϕ, m)
            a₁, a₂ = alias_coeffs((h₁[i+1], h₂[i+1]), isconj, isnyq)
            a₁, a₂ = (a₁, a₂) .* (sθΔΩ * cis(m * -ϕ₀))
            # Σ_{ℓ = m}^{lmax}
            for ℓ in m:lmax
                c = isodd(ℓ+m) ? a₁ - a₂ : a₁ + a₂
                alms[ℓ+1,m+1] += c * Λ[ℓ+1,m+1]
            end
        end
    end
    return alms
end

For brevity, we’ll walk through implementing a new synthesize function and simply present the new analyze function at the end since the pair largely mirror each other (and the previous two articles have already described their differences in detail).

In the original ECP implementation, we knew a priori the size of the map and the required size of the phase accumulation buffers trivially from the given nθ and nϕ arguments (nθ*nϕ and nϕ÷2+1, respectively). For our new generalized implementation, we must instead infer these from the description of the map rings. Since this is a common operation shared between synthesize and analyze, we’ll extract this pre-processing⁵ step into the helper function _mapring_counts:

"""
    (; npix, nϕmax) = _mapring_counts(maprings::MapRings)

Counts total number of pixels (`npix`) described by `maprings` and the longest observed ring (`nϕmax`).
"""
function _mapring_counts(maprings::MapRings)
    nϕmax, npix = 0, 0
    for ring in maprings.rings
        mult = mapreduce(!=(0), +, ring.offsets)  # single or paired rings?
        npix += mult * ring.nϕ                    # add 1 or 2 times Nϕ to total count
        nϕmax = max(nϕmax, ring.nϕ)               # keep track of longest seen ring
    end
    return (; npix, nϕmax)
end

With a new helper function in hand, we are now ready to implement our new synthesize function. The function signature takes only the MapRings description rather than explicit ECP dimensions, and we then allocate the return map and temporary phase-accumulation buffers based on what _mapring_counts returns:

function synthesize(maprings::MapRings, alms::Matrix{C}) where {C<:Complex}
    R = real(C)
    lmax, mmax = size(alms) .- 1
    (; npix, nϕmax) = _mapring_counts(maprings)
    nϕr = nϕmax ÷ 2 + 1  # (maximum) real-symmetric FFT's Nyquist length (index)

    mapbuf = Vector{R}(undef, npix)
    Λ = Matrix{R}(undef, lmax + 1, mmax + 1)
    λ₁ = zeros(C, nϕr) # phase factors for northern ring
    λ₂ = zeros(C, nϕr) # phase factors for southern ring

Most of the core loop is largely (structurally) unchanged, with only minor renaming of a couple of variables as long as we unpack the necessary constants at the start of the loop.

    for ring in maprings.rings
        nϕ, ϕ₀, cθ = ring.nϕ, ring.ϕ₀, ring.cosθ

The accumulation of the ring phases (which mixes the provided harmonic coefficients with the associated Legendre polynomials) is then almost completely unchanged:

        λlm!(Λ, lmax, mmax, cθ)

        for m in 0:mmax
            acc₁, acc₂ = zero(C), zero(C)
            # Σ_{ℓ = m}^{lmax}
            for ℓ in m:lmax
                term = alms[ℓ+1,m+1] * Λ[ℓ+1,m+1]
                acc₁ += term
                acc₂ += isodd(ℓ + m) ? -term : term
            end
            acc₁, acc₂ = (acc₁, acc₂) .* cis(m * ϕ₀) # ring offset rotation
            i, isconj, isnyq = alias_index(nϕ, m)
            acc₁, acc₂ = alias_coeffs((acc₁, acc₂), isconj, isnyq)
            λ₁[i+1] += acc₁
            λ₂[i+1] += acc₂
        end

The end of the loop changes a little more substantially, though:

1. Either of the offsets (nominally representing the northern- and southern-hemisphere rings) may be disabled, so the final fast Fourier transform (FFT) and assignment to the map must only occur when the offset is non-zero.
2. The phases buffers λ₁ and λ₂ are sized for the longest ring in the map, but each ring may be shorter (i.e. nϕ < nϕmax). The (inverse) FFT requires the buffer be sized appropriately for the actual nϕ number of destinations pixels, so we must be careful to pass along only the correctly-sized portion of the buffer.
3. The destination pixels in the map buffer may have non-unit stride and start at arbitrary offsets, so we must scatter to the correct set of destination pixel elements.

        stride = ring.stride
        o₁, o₂ = ring.offsets
        nϕh = nϕ ÷ 2 + 1
        if o₁ != 0
            idx = range(o₁, step = stride, length = nϕ)
            @views mapbuf[idx] = brfft(λ₁[1:nϕh], nϕ)
        end
        if o₂ != 0
            idx = range(o₂, step = stride, length = nϕ)
            @views mapbuf[idx] = brfft(λ₂[1:nϕh], nϕ)
        end
        λ₁ .= zero(C)
        λ₂ .= zero(C)
    end

    return mapbuf
end  # function synthesize

That completes the generalization of the synthesize function, and updating the corresponding analyze function proceeds very similarly. The hidden code section below presents both functions in full.

function synthesize(maprings::MapRings, alms::Matrix{C}) where {C<:Complex}
    R = real(C)
    lmax, mmax = size(alms) .- 1
    (; npix, nϕmax) = _mapring_counts(maprings)
    nϕr = nϕmax ÷ 2 + 1  # (maximum) real-symmetric FFT's Nyquist length (index)

    mapbuf = Vector{R}(undef, npix)
    Λ = Matrix{R}(undef, lmax + 1, mmax + 1)
    λ₁ = zeros(C, nϕr) # phase factors for northern ring
    λ₂ = zeros(C, nϕr) # phase factors for southern ring

    for ring in maprings.rings
        nϕ, ϕ₀, cθ = ring.nϕ, ring.ϕ₀, ring.cosθ

        λlm!(Λ, lmax, mmax, cθ)

        for m in 0:mmax
            acc₁, acc₂ = zero(C), zero(C)
            # Σ_{ℓ = m}^{lmax}
            for ℓ in m:lmax
                term = alms[ℓ+1,m+1] * Λ[ℓ+1,m+1]
                acc₁ += term
                acc₂ += isodd(ℓ + m) ? -term : term
            end
            acc₁, acc₂ = (acc₁, acc₂) .* cis(m * ϕ₀) # ring offset rotation
            i, isconj, isnyq = alias_index(nϕ, m)
            acc₁, acc₂ = alias_coeffs((acc₁, acc₂), isconj, isnyq)
            λ₁[i+1] += acc₁
            λ₂[i+1] += acc₂
        end

        stride = ring.stride
        o₁, o₂ = ring.offsets
        nϕh = nϕ ÷ 2 + 1
        if o₁ != 0
            idx = range(o₁, step = stride, length = nϕ)
            @views mapbuf[idx] = brfft(λ₁[1:nϕh], nϕ)
        end
        if o₂ != 0
            idx = range(o₂, step = stride, length = nϕ)
            @views mapbuf[idx] = brfft(λ₂[1:nϕh], nϕ)
        end
        λ₁ .= zero(C)
        λ₂ .= zero(C)
    end

    return mapbuf
end  # function synthesize

function analyze(maprings::MapRings, mapbuf::Vector{R}, lmax::Integer, mmax::Integer = lmax) where {R<:Real}
    C = complex(R)
    (; npix, nϕmax) = _mapring_counts(maprings)
    nϕr = nϕmax ÷ 2 + 1 # (maximum) real-symmetric FFT's Nyquist length (index)

    alms = fill(zero(C), lmax + 1, mmax + 1)
    Λ = zeros(R, lmax + 1, mmax + 1)
    λ₁ = zeros(C, nϕr) # phase factors for northern ring
    λ₂ = zeros(C, nϕr) # phase factors for southern ring

    for ring in maprings.rings
        nϕ, ϕ₀, cθ, ΔΩ = ring.nϕ, ring.ϕ₀, ring.cosθ, ring.ΔΩ

        λlm!(Λ, lmax, mmax, cθ) # λ_ℓ^m(cos θ) factors

        stride = ring.stride
        o₁, o₂ = ring.offsets
        nϕh = nϕ ÷ 2 + 1
        if o₁ != 0
            idx = range(o₁, step = stride, length = nϕ)
            @views λ₁[1:nϕh] = rfft(mapbuf[idx])
        end
        if o₂ != 0
            idx = range(o₂, step = stride, length = nϕ)
            @views λ₂[1:nϕh] = rfft(mapbuf[idx])
        end

        for m in 0:mmax
            i, isconj, isnyq = alias_index(nϕ, m)
            a₁, a₂ = alias_coeffs((λ₁[i+1], λ₂[i+1]), isconj, isnyq)
            a₁, a₂ = (a₁, a₂) .* (ΔΩ * cis(m * -ϕ₀))
            # Σ_{ℓ = m}^{lmax}
            for ℓ in m:lmax
                c = isodd(ℓ+m) ? a₁ - a₂ : a₁ + a₂
                alms[ℓ+1,m+1] += c * Λ[ℓ+1,m+1]
            end
        end
        λ₁ .= zero(C)
        λ₂ .= zero(C)
    end

    return alms
end  # function analyze

As a repeat of the conclusions from Part II, we can show that the GL grid better recovers the input spherical harmonic coefficients than the naive ECP grid, but this time we don’t require separate implementations of the synthesis and analysis functions.

function draw_realization(lmax)
    alms0 = zeros(ComplexF64, lmax+1, lmax+1)  # C_ℓ ∝ ℓ¯² ~ ⟨|alm|^2⟩
    for ℓ in 1:lmax
        alms0[ℓ+1,1] = randn(rng, Float64) / ℓ
        alms0[ℓ+1,2:ℓ+1] .= randn.(rng, ComplexF64) ./ ℓ
    end
    return alms0
end

function reanalyze(pixelization, alms)
    lmax, mmax = size(alms) .- 1
    return analyze(pixelization, synthesize(pixelization, alms), lmax, mmax)
end

alms0 = draw_realization(100)
alms_ecp = reanalyze(ecp_maprings(150, 300), alms0)
alms_gl  = reanalyze( gl_maprings(150, 300), alms0)

Conclusions¶

This article demonstrates that with only minor tweaks to the equidistant cylindrical projection (ECP) specific spherical harmonic transforms, we can generalize the algorithms to support a wide class of ring-based pixelizations on the sphere. Doing so only requires that each pixelization be able to describe itself with a collection of seven parameters per isolatitude ring.

A generalized algorithm was already demonstrated to be a useful tool — in Part II where we considered both the ECP and Gauss-Legendre (GL) grids, we were forced to copy-paste the ECP synthesize/analyze functions to make only minor edits to account for the difference in pixelizations. With this generalized algorithm, we can simply describe the two different grids and share the same spherical harmonic transforms.

As a concrete example where the difference between grids is less trivial (than ECP vs GL grids where only the $\theta$ coordinates and pixel areas are changed slightly but otherwise well approximated as a square grid on the sphere), the HEALPix pixelization⁶ is a very popular choice in cosmic microwave background (CMB) experiments. The HEALPix grid is constructed such that every pixel has uniform surface area, and as a concequence the number of pixels per ring is non-uniform across latitudes. Despite the HEALPix grid having a very different layout than the more familiar ECP grid, though, the generalized ring-based algorithm presented here is sufficient for performing spherical harmonic transforms on a HEALPix map.