離散型物体データの球面調和関数展開

球面調和関数を使って.obj データを級数展開してくれるやつを後輩が書いてくれたのでやってみる。
例えばこれ(J Med Syst. 2014 May;38(5):20.)では、病気の肝臓と健康な肝臓をCT画像から構成して、球面調和関数による比較分類を行っている。
肝硬変が進んだ病気の肝臓は、線維化といって見た目がカチカチのゴツゴツになる。臨床的には辺縁(エッジ)が鈍と言ったり、ミクロには肝表面に再生結節というブツブツができたりする。
CT画像から肝臓の形を再構成したはいいが、物体同士の比較をするのに球面調和関数で級数展開しておく。degree とorder によって球面調和関数が決まる。ただし、 が小さすぎると大雑把にしかならない。次元の高い は物体表面の細かい正確さを決めているとされる。

Fig.5 からパクった。 が大きくなるとより肝臓っぽくなる。この論文では以前の報告から でよさそうと決めている。他に、医学方面で球面調和関数を使っている論文は例えば脳MRI とかだが、12 から25 くらいでやっていることが多そう。
ただし、 をひたすら大きく取れば物体細部の分解能が上がるかというとそうではないっぽくて、例えば球面調和関数は で計算されるから、 ひとつにつき計算されるorder が となる(球面調和関数の展開図がピラミッド状に表現されていることからわかる)。
　
球面調和関数の係数はspintransformation の最後に吐き出される仕様になっている。球面調和関数自体はR ではやるのが面倒くさそうなのでこれまたrPython でゴリ押しする。
くらいでそれ以上は細かいところがグダグダになるようなのでここが限界のようだが、かといって元のうさぎが表現できているかというとあやしい。
　

元リンクのサンプル。
　

元のスタンフォードうさぎ

library(rgl)
library(geometry)
library(SphericalCubature)
d <- read.table("bunny.obj") # とりあえずスタンフォードうさぎ
d0 <- split(d[,-1], d[,1])
m0 <- tmesh3d(t(d0$v), t(d0$f), homogeneous=FALSE)

# 球面調和関数の係数データ
sh_coef <- as.numeric(strsplit(readLines("coef.txt"), " ")[[1]])
library(rPython)

# Python で球面調和関数を計算できるようにおまじない
python.exec("import numpy as np")
python.exec("from scipy.special import sph_harm as _sph_harm")
python.exec("
def sph_harm(l, m, theta, phi):\
  if m >= 0:\
    return np.sqrt(2) * _sph_harm(m, l, phi, theta).real\
  else:\
    return np.sqrt(2) * _sph_harm(-m, l, phi, theta).imag")

# 球面調和関数の係数を与えて元の物体の級数展開をゴリ押しする
python.exec("
def SH(Lmax, c, u, v):\
  r = np.zeros([len(u), len(v)])\
  for l in range(Lmax + 1):\
    for m in range(-l, l + 1):\
      r += c[l**2 + l + m] * sph_harm(l, m, np.array(u, ndmin=2).T, v).real\
  \
  return map(list, r)")

# 球面座標用の theta, phi
u <- seq(0, 180, length=80)*pi/180 # u = np.linspace(0, np.pi, 80)
v <- seq(0, 360, length=60)*pi/180 # v = np.linspace(0, 2 * np.pi, 60)

# 元リンクのデータを再現するなら
cvec <- c(10, 1, 2, 0, 1, 5, 3, 4, 0)

# rPython でPython に投げてゴリ押しする
r <- c(do.call(rbind, python.call("SH", 2, sh_coef, u, v)))

x <- r*c(outer(sin(u), cos(v)))
y <- r*c(outer(sin(u), sin(v)))
z <- r*c(outer(cos(u), rep(1, length(v))))
xyz <- cbind(x, y, z)
Z <- convhulln(xyz/r)
m1 <- tmesh3d(t(xyz), tmesh3d(t(xyz/r), t(Z), homogeneous=F)$it, homogeneous=FALSE)
open3d()
shade3d(m1, col="grey")
rgl.viewpoint(30, 20, zoom=0.8)
movie3d(spin3d(axis=c(0, 1, 0)), duration=12, fps=5)