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