Spherical harmonic transforms

Spherical harmonic transforms for MW and symmetric MW samplings are created via the factory function harmonic_transform. The factory function takes at least one parameter, lmax.

from tensossht import harmonic_transform
transform = harmonic_transform(lmax=8, spin=0, dtype=tf.float64)

Optionally, we can also specify lmin, mmin, and mmax, though little testing has been done apart from mmin=0. Finally, the spin can be specified either with spin=0 for single-spin transforms, or with smin=-1 and smax=1 for multi-spin transforms. The latter case is another way to create a transform for real image-space signals. for this specific case:

rtransform = harmonic_transform(lmax=8, mmin=0, spin=0, dtype=tf.float64)

The underlying floating point can be specified via the dtype keyword argument. It defaults to tf.float32 (just like any other tensorflow operation). Furthermore, the image-space sampling scheme can be specified via the sampling keyword. It defaults to “MW”.

MW Sampling in harmonic-space and image-space

The transform works on a fixed basis in harmonic space, fully specifid by lmax, lmin, mmax, mmin. The order of the basis functions can be assessed via the two properties:

>>> transform.llabels
<tf.Tensor: shape=(64,), dtype=int32, numpy=
array([0, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4,
       4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6,
       6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7],
      dtype=int32)>

>>> transform.mlabels
<tf.Tensor: shape=(64,), dtype=int32, numpy=
array([ 0, -1,  0,  1, -2, -1,  0,  1,  2, -3, -2, -1,  0,  1,  2,  3, -4,
       -3, -2, -1,  0,  1,  2,  3,  4, -5, -4, -3, -2, -1,  0,  1,  2,  3,
        4,  5, -6, -5, -4, -3, -2, -1,  0,  1,  2,  3,  4,  5,  6, -7, -6,
       -5, -4, -3, -2, -1,  0,  1,  2,  3,  4,  5,  6,  7], dtype=int32)>

The transform for real signals only contains \(m \geq 0\) components:

>>> rtransform.mlabels
<tf.Tensor: shape=(36,), dtype=int32, numpy=
array([0, 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 5, 0,
       1, 2, 3, 4, 5, 6, 0, 1, 2, 3, 4, 5, 6, 7], dtype=int32)>

The sampling can be queried for \(\theta\) and \(\phi\):

>>> transform.thetas.numpy().round(3)
array([0.209, 0.628, 1.047, 1.466, 1.885, 2.304, 2.723, 3.142])

>>> transform.phis.numpy().round(3)
array([0.   , 0.419, 0.838, 1.257, 1.676, 2.094, 2.513, 2.932, 3.351,
       3.77 , 4.189, 4.608, 5.027, 5.445, 5.864])

A \(\theta\) by \(\phi\) grid is also provided for convenience:

>>> transform.grid.shape
TensorShape([2, 8, 15])

As well as list of \((\theta, \phi)\) points:

>>> transform.points.shape
TensorShape([120, 2])

Forward and Inverse Transforms with MW Sampling

The transforms themselves are called via forward() and inverse(). For instance, we can create a set of two random images and transform them back and forth.

from pytest import approx
nthetas = transform.lmax
nphis = 2 * transform.lmax - 1
assert (int(nthetas), int(nphis)) == transform.sampling.grid.shape[1:]
images = tf.complex(
    tf.random.uniform((2, nthetas, nphis), dtype=transform.real_dtype),
    tf.random.uniform((2, nthetas, nphis), dtype=transform.real_dtype),
)
representable_images = transform.inverse(transform.forward(images))
assert tf.reduce_any(representable_images != 0)

A first back-and-forth removes from the images components that cannot be represented with the transform’s harmonic-space basis. The second back-and-forth, however, should yield exactly the same result as the first, except for numerical noise.

toandfro = transform.inverse(transform.forward(representable_images))
assert toandfro.numpy() == approx(representable_images.numpy())

Similarly, for a real image-space signal:

images = tf.random.uniform((2, nthetas, nphis), dtype=transform.real_dtype)
representable_images = rtransform.inverse(rtransform.forward(images))
toandfro = rtransform.inverse(rtransform.forward(representable_images))
assert toandfro.numpy() == approx(representable_images.numpy())

Note that both forward() and inverse() accept keywords specifying which dimension corresponds to \(\theta\), \(\phi\) and to the harmonic space coefficients. By default, and to avoid extra transpose operations, \(\theta\) and the harmonic-space coefficients should be stored along the last dimension, and \(\phi\) on the penultimate.

Forward and Inverse Transforms with MWSS Sampling

MW symmetric sampling is a more widespread sampling for 360 images. It requires slightly more memory and the transforms require slightly more compute. However, both samplings are quite similar. A transform with MWSS sampling can be created as follows:

from tensossht import harmonic_transform
transform = harmonic_transform(lmax=5, spin=0, sampling="MWSS", dtype=tf.float64)

The sampling scheme is available through the transform:

>>> transform.sampling.thetas.numpy().round(3)
array([0.   , 0.628, 1.257, 1.885, 2.513, 3.142])

>>> transform.sampling.phis.numpy().round(3)
array([0.   , 0.628, 1.257, 1.885, 2.513, 3.142, 3.77 , 4.398, 5.027,
       5.655])

Let’s create a full set of spherical harmonics with MWSS sampling. For illustration purposes, the \(\phi\) and \(\theta\) dimensions are ordered in a less orthodox (and less efficient) manner:

from tensossht import spherical_harmonics
grid = transform.sampling.grid
sph = spherical_harmonics(grid[0], grid[1], lmax=transform.lmax)

sph is a three-dimensional array with the first dimension refering to \(theta\), the second to \(phi\), and the second dimension to the spherical harmonics with \(l \leq 4\).

>>> sph.shape
TensorShape([6, 10, 25])

In spherical-harmonic space, we expect the set of images to be transformed to the identity matrix, since each image is a single spherical harmonic:

from pytest import approx
sph_space = transform.forward(sph, theta_dim=0, phi_dim=1).numpy()
assert sph_space == approx(tf.eye(sph.shape[-1]).numpy(), abs=1e-7)

Similarly, for a real image-space signal:

rsph = tf.math.real(sph)
rtransform = transform.real
rsph_space = rtransform.forward(rsph, theta_dim=0, phi_dim=1)
as_complex = transform.forward(
    tf.complex(rsph, tf.constant(0, rtransform.real_dtype)),
    theta_dim=0,
    phi_dim=1,
)
as_complex = tf.transpose(
    tf.boolean_mask(tf.transpose(as_complex), transform.mlabels >= 0)
)
assert rsph_space.numpy() == approx(as_complex.numpy())

We can try the same operations for inverse transforms. If we feed the inverse transform an identity matrix, it should spew out the spherical harmonics under the MWSS sampling scheme:

coefficients = tf.eye(
    transform.lmax * transform.lmax, dtype=transform.complex_dtype
)
real_space = transform.inverse(coefficients, phi_dim=1, theta_dim=0)
assert real_space.numpy() == approx(sph.numpy(), abs=1e-6)

Similarly for real-space signals:

coefficients = (
    tf.eye(rtransform.labels.shape[1], dtype=rtransform.complex_dtype)
    * tf.cast(
        tf.where(rtransform.mlabels == 0, 1.0, 0.5),
        dtype=rtransform.complex_dtype
    )
)
real_space = rtransform.inverse(coefficients, phi_dim=1, theta_dim=0)
rsph = spherical_harmonics(grid[0], grid[1], lmax=rtransform.lmax, mmin=0)
assert real_space.numpy() == approx(tf.math.real(rsph).numpy(), abs=1e-5, rel=1e-5)

API

class tensossht.transforms.transforms.HarmonicTransform(forward, inverse, image_sampling, harmonic_sampling)[source]

Fast transform between image space and spherical harmonic space.

Parameters:

  • forward (Callable) –

  • inverse (Callable) –

  • image_sampling (ImageSamplingBase) –

  • harmonic_sampling (HarmonicSampling) –

forward(images: tf.Tensor, theta_dim: int = -2, phi_dim: int = -1, coeff_dim: int = -1) → tf.Tensor

Transform from image to harmonic space. The image axes theta_dim and phi_dim, as well as the harmonic-coefficient axis coeff_dim can be specified on input. They default to the innermost axes (.i.e. \(\phi\) is contiguous in memory).

inverse(images: tf.Tensor, theta_dim: int = -2, phi_dim: int = -1, coeff_dim: int = -1) → tf.Tensor

Transform from harmonic to image space. The image axes theta_dim and phi_dim, as well as the harmonic-coefficient axis coeff_dim can be specified on input. They default to the innermost axes (.i.e. \(\phi\) is contiguous in memory).

property complex

Transform for complex signals.

If the transform is already for complex signals (\(m_\mathrm{min} \leq 0\)), then this property returns self. Otherwise it returns an equivalent transform with \(m_\mathrm{min} = -l_\mathrm{max}\).

property complex_dtype
property grid
harmonic_sampling

Harmonic-space sampling

property labels
property llabels
property lmax
property lmin
property mlabels
property mmax
property mmin
property ncoeffs
property phis
property points

Image-space points.

property real

Transform for real signals.

If the transform is already for real signals (\(m_\mathrm{min} \geq 0\)), then this property returns self. Otherwise it returns an equivalent transform with \(m_\mathrm{min} = 0\).

property real_dtype
sampling

Image-space sampling

property slabels
property smax
property smin
property thetas
class tensossht.transforms.forward.ForwardTransform(images_to_fmm, fmm_to_gmm, gmm_to_coeffs)[source]
Parameters:

  • images_to_fmm (Callable) –

  • fmm_to_gmm (Callable) –

  • gmm_to_coeffs (GmmToCoeffs) –

__call__(images, theta_dim=-2, phi_dim=-1, coeff_dim=-1)[source]
Parameters:

  • images (Tensor) –

  • theta_dim (int) –

  • phi_dim (int) –

  • coeff_dim (int) –

Return type:

Tensor

class tensossht.transforms.inverse.InverseTransform(coeffs_to_fmm, fft, slicer)[source]
Parameters:

  • coeffs_to_fmm (Callable[[ndarray, HarmonicAxes], ndarray]) –

  • fft (Callable[[ndarray], ndarray]) –

  • slicer (Callable[[ndarray], ndarray]) –

__call__(coefficients, theta_dim=-2, phi_dim=-1, coeff_dim=-1)[source]
Parameters:

  • coefficients (ndarray) –

  • theta_dim (int) –

  • phi_dim (int) –

  • coeff_dim (int) –

Return type:

ndarray