球面谐波实数形式的可视化

关于 SH：球面调和 Spherical Harmonics (SH)

MRtrix3 中使用的 SH 基函数 为：

python可视化程序：

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
# The following import configures Matplotlib for 3D plotting.
from mpl_toolkits.mplot3d import Axes3D
from scipy.special import sph_harm

# Grids of polar and azimuthal angles
theta = np.linspace(0, np.pi, 100)
phi = np.linspace(0, 2*np.pi, 100)
# Create a 2-D meshgrid of (theta, phi) angles.
theta, phi = np.meshgrid(theta, phi)
# Calculate the Cartesian coordinates of each point in the mesh.
xyz = np.array([np.sin(theta) * np.sin(phi),
                np.sin(theta) * np.cos(phi),
                np.cos(theta)])

def plot_Y(ax, el, m):
    """Plot the spherical harmonic of degree el and order m on Axes ax."""

    # NB In SciPy's sph_harm function the azimuthal coordinate, theta,
    # comes before the polar coordinate, phi.
    Y = sph_harm(abs(m), el, phi, theta)

    # Linear combination of Y_l,m and Y_l,-m to create the real form.
    if m < 0:
        Y = np.sqrt(2) * Y.imag
    elif m > 0:
        Y = np.sqrt(2) * Y.real
    Yx, Yy, Yz = np.abs(Y) * xyz

    # Colour the plotted surface according to the sign of Y.
    cmap = plt.cm.ScalarMappable(cmap=plt.get_cmap('PRGn'))
    cmap.set_clim(-0.5, 0.5)

    ax.plot_surface(Yx, Yy, Yz,
                    facecolors=cmap.to_rgba(Y.real),
                    rstride=2, cstride=2)

    # Draw a set of x, y, z axes for reference.
    ax_lim = 0.5
    ax.plot([-ax_lim, ax_lim], [0,0], [0,0], c='0.5', lw=1, zorder=10)
    ax.plot([0,0], [-ax_lim, ax_lim], [0,0], c='0.5', lw=1, zorder=10)
    ax.plot([0,0], [0,0], [-ax_lim, ax_lim], c='0.5', lw=1, zorder=10)
    # Set the Axes limits and title, turn off the Axes frame.
    ax.set_title(r'$Y_{{{},{}}}$'.format(el, m))
    ax_lim = 0.5
    ax.set_xlim(-ax_lim, ax_lim)
    ax.set_ylim(-ax_lim, ax_lim)
    ax.set_zlim(-ax_lim, ax_lim)
    ax.axis('off')

画出单个SH基：

fig = plt.figure(figsize=plt.figaspect(1.))
ax = fig.add_subplot(projection='3d')
l, m = 3, 0
plot_Y(ax, l, m)
plt.savefig('Y{}_{}.png'.format(l, m))
plt.show()

Y3_0

如果想要画出一系列的基：

el_max = 3
figsize_px, DPI = 800, 100
figsize_in = figsize_px / DPI
fig = plt.figure(figsize=(figsize_in, figsize_in), dpi=DPI)
spec = gridspec.GridSpec(ncols=2*el_max+1, nrows=el_max+1, figure=fig)
for el in range(el_max+1):
    for m_el in range(-el, el+1):
        print(el, m_el)
        ax = fig.add_subplot(spec[el, m_el+el_max], projection='3d')
        plot_Y(ax, el, m_el)
plt.tight_layout()
plt.savefig('sph_harm.png')
plt.show()

sph_harm.png

如果想要画出 FOD ：

def plot_SH(coeffs):
    import numpy as np
    import matplotlib.pyplot as plt
    import matplotlib.gridspec as gridspec
    # The following import configures Matplotlib for 3D plotting.
    from mpl_toolkits.mplot3d import Axes3D
    from scipy.special import sph_harm

    # Grids of polar and azimuthal angles
    theta = np.linspace(0, np.pi, 100)
    phi = np.linspace(0, 2 * np.pi, 100)
    # Create a 2-D meshgrid of (theta, phi) angles.
    theta, phi = np.meshgrid(theta, phi)
    # Calculate the Cartesian coordinates of each point in the mesh.
    xyz = np.array([np.sin(theta) * np.sin(phi),
                    np.sin(theta) * np.cos(phi),
                    np.cos(theta)])

    el_max = 8
    sumX, sumY, sumZ = np.zeros((100, 100)), np.zeros((100, 100)), np.zeros((100, 100))
    for el in range(el_max + 1):
        if el%2 != 0:
            continue
        for m_el in range(-el, el + 1):
            # print(el, m_el)
            Y = sph_harm(m_el, el, phi, theta)
            if m_el < 0:
                Y = np.sqrt(2) * Y.imag
            elif m_el > 0:
                Y = np.sqrt(2) * Y.real
            else:
                Y = Y.real
            i = int(1 / 2 * el * (el + 1) + m_el)
            Yx, Yy, Yz = coeffs[i] * Y * xyz

            sumX += Yx
            sumY += Yy
            sumZ += Yz

    fig1 = plt.figure(figsize=plt.figaspect(1.))
    ax1 = fig1.add_subplot(projection='3d')
    ax1.plot_surface(sumX, sumY, sumZ)
    ax1.set_title('fod')
    ax1.axis('off')
    plt.show()

from dipy.io.image import load_nifti
fod, affine = load_nifti(fod_path)
coeffs = fod[103,89,76]
plot_SH(coeffs)

FOD

参考：
https://scipython.com/blog/visualizing-the-real-forms-of-the-spherical-harmonics/
https://mrtrix.readthedocs.io/en/latest/concepts/spherical_harmonics.html