Spherical harmonics

This notebook demonstrates how to use a set of spherical harmonics functions from pyshtools to approximate pointclouds in 3D.

import napari_stress
import vedo
import napari
import matplotlib.pyplot as plt

Let’s create some sample surface data which we would like to approximate with a spherical harmonics expansion (SHE), for instance an ellipsoid:

data = vedo.shapes.Ellipsoid()
viewer = napari.Viewer(ndisplay=3)
viewer.add_points(data.points(), size=0.05, face_color='orange', name='Raw')
napari.utils.nbscreenshot(viewer, canvas_only=True)

WARNING: QWindowsWindow::setGeometry: Unable to set geometry 3258x2028+12+126 (frame: 3290x2116-4+54) on QWidgetWindow/"_QtMainWindowClassWindow" on "\\.\DISPLAY1". Resulting geometry: 7684x2108+20+167 (frame: 7716x2196+4+95) margins: 16, 72, 16, 16 minimum size: 474x554 MINMAXINFO maxSize=0,0 maxpos=0,0 mintrack=1454,1750 maxtrack=0,0)
WARNING: QOpenGLFramebufferObject: Framebuffer incomplete attachment.
WARNING: QOpenGLFramebufferObject: Framebuffer incomplete attachment.
WARNING: QOpenGLFramebufferObject: Framebuffer incomplete attachment.
WARNING: QOpenGLFramebufferObject: Framebuffer incomplete attachment.
WARNING: QOpenGLFramebufferObject: Framebuffer incomplete, missing attachment.

Expansion

Now we use the napari_stress.fit_spherical_harmonics() function to approximate this pointcloud. The important parameter is the max_degree parameter. Higher values for this parameter will tell the function to include spherical harmonic functions of higher orders to approximate the pointcloud and arrive at a better estimation.

You’ll see that the approximation of the ellipsoid becomes better with higher order.

fitted_points = napari_stress.fit_spherical_harmonics(data.points(), max_degree=5)
fitted_points[1]['size'] = 0.05
viewer.add_points(fitted_points[0], **fitted_points[1])
napari.utils.nbscreenshot(viewer, canvas_only=True)

fitted_points = napari_stress.fit_spherical_harmonics(data.points(), max_degree=15)
fitted_points[1]['size'] = 0.05
viewer.add_points(fitted_points[0], **fitted_points[1])
napari.utils.nbscreenshot(viewer, canvas_only=True)

Quantification

You may have noticed that the data that is returned by the fit_spherical_harmonics() function returns more than simply the fitted pointcloud, but a number of additional features, among which are some interesting data, most notably the fit residues and the coefficients of the sperical harmonics function:

fitted_points[1]['features'].keys()

dict_keys(['error'])

fitted_points[1]['metadata'].keys()

dict_keys(['spherical_harmonics_coefficients', 'spherical_harmonics_implementation'])

As for the error, this is encoded in the colorscale of the points in the viewer screenshots above and shows the euclidian distance between the approximated point and the respective input point. We could plot the fit residues as a histogram:

hist = plt.hist(fitted_points[1]['features']['error'], 50)
plt.xlabel('Fit residues')
plt.ylabel('Occurences')

Text(0, 0.5, 'Occurences')

As for the coefficients, this is a bit more complicated.

Mathematical background: A spherical harmonics approximation (a.k.a. expansion) \(f(\theta,\phi)\) can be written as a superposition of multiple single spherical harmonics functions \(Y_l^m(\theta\phi)\) of different degree and order:

\(f(\theta,\phi) = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} f_l^m Y_l^m(\theta\phi)\)

whereas the coefficients \(f_l^m\) determmine how much contribution of which degree and order is neede to achieve the best approximation of the input pointcloud. Currently, napari-stress performs the sperhical harmonics expansion separately for the \(x\), \(y\) and \(z\) direction, hence you’ll receive three distinct shape spectra:

fitted_points[1]['metadata']['spherical_harmonics_coefficients'].shape

(3, 16, 16)

fig, axes = plt.subplots(ncols=3, figsize=(15,5))
axes[0].imshow(fitted_points[1]['metadata']['spherical_harmonics_coefficients'][0])
axes[1].imshow(fitted_points[1]['metadata']['spherical_harmonics_coefficients'][1])
axes[2].imshow(fitted_points[1]['metadata']['spherical_harmonics_coefficients'][2])

axes[0].set_title('X')
axes[1].set_title('Y')
axes[2].set_title('Z')

Text(0.5, 1.0, 'Z')