tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート

球面調和関数係数(3 x n行列)の可視化

3次元空間にひと塊の形(球と同相)があって、その表面上の点の集合を考える。何らかの方法で塊は球面に変形できて、その変形に伴って点の集合は球面上に写される、と考える。その時、元のX、Y、Z座標をそれぞれ球面上の分布($S^{2} \to \mathrm{R}$)と考えると、これは球面調和関数で分解できることになる。3つ係数の列が取れるので、それを3xnの行列とする。これを元の塊の表面に戻すのは、係数を使って色々の$\theta$、$\phi$についてX、Y、Zを計算してやるだけで良い。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.special import sph_harm as _sph_harm
import argparse

parser = argparse.ArgumentParser()
parser.add_argument('-i', '--input', required=True)
parser.add_argument('-o', '--output', required=True)
args = parser.parse_args()

# e.g. "/home/tak0kada/tmp/obj2sph/data/output/coef/t0001/cell00001_SPHARM.txt"
INPUT = args.input
# e.g. "/home/tak0kada/tmp/cell00001.obj"
OUTPUT = args.output

coef = np.loadtxt(INPUT)
coef_x = coef[:, 0]
coef_y = coef[:, 1]
coef_z = coef[:, 2]

theta = [i/50*np.pi for i in range(50)]
phi = [i/50*np.pi*2 for i in range(50)]

def sph_harm(l, m, theta, phi):
    if m > 0:
        return np.sqrt(2) * (-1)**m * _sph_harm(m, l, phi, theta).real
    if m == 0:
        return _sph_harm(m, l, phi, theta).real
    else:
        return np.sqrt(2) * (-1)**m * _sph_harm(-m, l, phi, theta).imag

def l(n):
    return int(np.sqrt(n))

def m(n):
    l_ = l(n)
    return n - l_**2 - l_

def x(theta, phi):
    return np.sum([coef_x[n] * sph_harm(l(n), m(n), theta, phi) for n in range(coef_x.size)])

def y(theta, phi):
    return np.sum([coef_y[n] * sph_harm(l(n), m(n), theta, phi) for n in range(coef_y.size)])

def z(theta, phi):
    return np.sum([coef_z[n] * sph_harm(l(n), m(n), theta, phi) for n in range(coef_z.size)])

X = []; Y = []; Z = []

for t in theta:
    for p in phi:
        X.append(x(t, p))
        Y.append(y(t, p))
        Z.append(z(t, p))

fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(X, Y, Z)
plt.show()

with open(OUTPUT, "w") as f:
    for i in range(len(X)):
        vert = "v {} {} {}\n".format(X[i], Y[i], Z[i])
        f.write(vert)

使い方

python coef2obj-v.py -i coef.txt -o sphere_v.obj