Houdiniで球面調和関数を理解してみる！

はじめに

あけましておめでとうございます。

皆さんは「球面調和関数」って聞いたことありますか？
筆者と同じ、ゲーム業界のグラフィックスに関わる方なら、リアルタイムレンダリングに関連する用語として聞いたことがあるかもしれません。
どんなものか調べてみると、大抵は下の図のような、ぷにぷにしたものがピラミッド状に並んだ画像が出てくると思います。

（https://en.wikipedia.org/wiki/Spherical_harmonics より引用）

このぷにぷに達が何に役立つかというと、具体的にはキューブマップの圧縮などに役立つのですが……。（こちらの記事が分かりやすいです↓）

今回筆者は、ちょっと違う用途で球面調和関数を扱えるようになりたいと思い、Houdiniを使って勉強してみました。
どんな用途で？という話は機会があれば別記事に書くとして、本記事ではHoudiniで球面調和関数を可視化したり、利用できるようにした過程を紹介します。

微妙に前回のフーリエ変換の記事の続編のような感じです。
（というより、球面調和関数を理解する前準備としてフーリエ変換を勉強したのでした）

最終的な成果物

球面調和関数自体の形（ぷにぷに）を可視化している様子です。

球面調和関数を用いて入力形状を「スケーリング係数」に圧縮し、さらに「スケーリング係数」と球面調和関数から「入力形状の近似」を再構築しています。
つまり、形状データの圧縮と復元の過程を表しています。

おことわり

• Houdini Indie 20.5.370 を使用しています。
• 筆者の数学知識は浅いので、不正確な部分があるかもしれません。加えて、数学的な背景や意味を理解することは必ずしも優先していません。
• VEXコードを載せていますが、あまりいいコードにはなっていないかもしれません。
• この記事を参考に何かする場合は自己責任でお願いします。
• 今回かなりMicrosoft Copilotに頼っています。

0. まずは概要をつかみたい

最初になにをしたかというと、ひたすらMicrosoft Copilotにわからないことを質問しました。
というのも、この時点では筆者も球面調和関数（による展開）について、以下のような理解しかできていなかったからです。

• フーリエ変換を球面上で行っているようなもの、ということらしい
  • つまり、球面上の任意の関数（≒グレスケ画像的なもの）を少ないデータ量に圧縮できるということ
• ぷにぷにの一つ一つが周波数を変えた正弦波に相当し、これを球面調和関数と呼んでいるらしい

具体的に以下のようなやり取りをしました。

• Q. 球面調和関数の式を教えて！
  • A. 次の式です！

   Y_{l}^{m}(\theta,\phi) = \sqrt{\frac{(2l+1)(l-m)!}{4\pi(l+m)!}} P_l^m(\cos\theta) e^{im\phi}
  

• Q. 球面調和関数の式に含まれる各項の意味を教えて！
  • A. それぞれ以下のような意味があります！
    • $\sqrt{\frac{(2l+1)(l-m)!}{4\pi(l+m)!}}$
      • 正規化定数。関数の全体的な大きさを適切にスケールする。
    • $ P_l^m(\cos\theta) $
      • ルジャンドル陪関数。球面調和関数の形状を決定する。
    • $e^{im\phi}$
      • 複素指数関数。関数の「位相」や「回転」を表す。
    • 各変数
      • $l$
        • 次数。負でない整数。
      • $m$
        • 磁気量子数。[-l,l]の範囲の整数。
      • $\theta$
        • 緯度角。[0,π]の範囲。
      • $\phi$
        • 経度角。[0,2π]の範囲。
• Q. 経度角φだけが変化するとどうなる？
  • A. 複素指数関数の部分が変化することで関数が球面上で回転します！
• Q. では、ルジャンドル陪関数の部分は緯度角θに応じて波打つ形状を決定している？
  • A. 大体そんな感じです！
• Q. θが緯度方向に沿った突起の配置を決め、Φが経度方向に沿った配置を決め、lとmが突起の数を決めているという理解で合ってる？
  • A. ほぼ合っています！ただしmは、突起の対称性や全体的な配置に影響します。
• Q. 球面調和関数におけるルジャンドル陪関数は球面上を緯度方向に進む正弦波のような性質？
  • A. ルジャンドル陪関数は周期的な関数ではなく[0,π]の範囲で特定の形状を持つ点が正弦波とは異なります！

回答について全て鵜呑みにするわけにはいきませんが、少なくとも筋は通っているように感じました。
球面調和関数の式も、確認した限り正しいようなので、信憑性はそれなりに高いと思ってよさそうです。

球面調和関数についての理解を深めることができたので、次はこれを可視化することを目標にします。
可視化の結果は、最初の画像のような、次数に応じて複雑さを増すぷにぷにになるはずです。

1. 実践編 球面調和関数の可視化

前提

先に確認した通り、球面調和関数には4つの変数が必要で、そのうち2つは球面上の緯度角θと経度角φとなっています。
残る2つである次数lと磁気量子数mはぷにぷに一つ一つで固有なので、球面上のどの地点であるかによって値が変動することはありません。
また球面調和関数は「球面上の波のようなもの」ということなので、その返り値は波の高さ情報であるはずです。

というわけで、Sphereのポイントの位置からθとφを求め球面調和関数の入力に使い、lとmは外部から与え、その計算結果を球の中心からの距離としてポイントを動かすことで、球面上の波の高さを表現する流れを考えます。

また今回は割と複雑な数式を扱うので、VEX内で関数を定義して使っていきたいです。
VEXで関数定義する場合の書き方は以下の通りです。

VEXで関数定義する場合の書き方

function int Hoge(int a)
{
    // 処理内容
}

実数球面調和関数

本来の球面調和関数では$e^{im\phi}$の部分が複素数ですが、CGでは実数のみ扱うため実数球面調和関数というものに置き換えるそうです。

実数球面調和関数の式は以下とのこと。

Y_{l}^{m}(\theta,\phi) = \left\{
    \begin{array}{}
        \sqrt{2} K_{l}^{m} \sin(-m\phi) P_{l}^{-m}(\cos\theta) &(m<0)\\
        K_{l}^{0} P_{l}^{0}(\cos\theta) &(m=0)\\
        \sqrt{2} K_{l}^{m} \cos(m\phi) P_{l}^{m}(\cos\theta) &(m>0)\\
    \end{array}
\right.

K_{l}^{m} = \sqrt{\frac{(2l+1)(l-\lvert m \rvert)!}{4\pi(l+\lvert m \rvert)!}}

実際に論文などに載っているぷにぷにも、「実数球面調和関数の絶対値」の可視化であるようです。
以降では、この実数球面調和関数を取り扱っていくことにします。

ルジャンドル陪関数

それではいよいよ、数式をVEXで実装していきます。
まず問題になるのは、やはりルジャンドル陪関数なる箇所ですね……。
Copilotによると、緯度方向に波打つ形を作っている部分ということです。

ルジャンドル陪関数について更に調べてみると、以下のように定義されるとのことでした。

P_{l}^{m}(x) = (-1)^m(1-x^2)^\frac{m}{2} \frac{d^m}{dx^m} P_{l}(x) 

P_{l}(x) = \frac{1}{2^ll!} \frac{d^l}{dx^l}(x^2-1)^l

なんだかよくわからない上に、微分演算子が出てきてしまいました……。
しかも1回微分ではなくm回、l回微分です。

複数回の微分を計算するプログラムなんて組んだことないので、非常に悩みました。
Copilotに具体的なVEXコードを提示してもらったりもしましたが、実際に動作するものは得られず……。
しかし、ルジャンドル陪関数についても詳しく書いてくださっているt-potさんの記事を改めて見てみると解決方法が書いてありました。

ルジャンドル陪関数は、次の性質を持っています(これらは、前の定義の式から計算することができます)。

P_{m}^{m}(x) = (-1)^m(2m-1)!!(1-x^2)^\frac{m}{2}

P_{m+1}^{m}(x) =(2m+1)x P_{m}^{m}(x)

(l-m)P_{l}^{m}(x) =(2l-1)xP_{l-1}^{m}(x)-(l+m-1)xP_{l-2}^{m}(x)

この関係式から、しちめんどくさい微分を計算しなくても、ルジャンドル陪関数を求めることができます。
具体的には、１つめの式から、P00(x), P11(x), P22(x)…を求めて、次に２つめの式からP10(x), P21(x), P32(x)…を求めて、その後に３番目の式からそれ以外を求めます。
この計算方法は、特にコンピュータを使って求める場合に標準的な方法になるでしょう。

とのこと。

つまり、

• $l=m$の場合は1個目の式だけで計算できる
• $l+1=m$の場合は1個目の式を計算した後に2個目の式を計算すれば求まる
• それ以外のケースでは、
  • 2個目の式$P_{m+1}^{m}$と1個目の式$P_{m}^{m}$の結果を、それぞれ$P_{l-1}^{m}(x)$と$P_{l-2}^{m}(x)$に代入して$P_{m+2}^{m}(x)$を求める（左辺は$((m+2)-m)$で割る）
  • $P_{m+2}^{m}$と$P_{m+1}^{m}$をそれぞれ$P_{l-1}^{m}(x)$と$P_{l-2}^{m}(x)$に代入して$P_{m+3}^{m}(x)$を求める（左辺は$((m+3)-m)$で割る）
  • ……というのを、$P_{l}^{m}(x)$に到達するまで繰り返すことで求める
• （$l<m$になるケースはそもそも球面調和関数の前提に反しているので無視）

ということです。
なるほど！
この3つの式を正しい順番で計算していけば、ルジャンドル陪関数はどうにかなりそうです！

!!となっている部分は二重階乗で、整数を一個飛ばしで階乗するものとのことです。
Wikipediaによると負の奇数の二重階乗はありうるということですが、実数球面調和関数の式を見ると、$m<0$のときはmのマイナスを打ち消してルジャンドル陪関数に入力するようになっています。結果、二重階乗の入力として考慮しなければならないのは$(2m-1)!!$なので -1以上の整数 となります。

実際にVEXコードにしてみたものが以下です。

ルジャンドル陪関数

// 二重階乗
function int doubleFactorial(int a) // 二重階乗はWikipediaによると「負の奇数」アリなので、a=(-1)は有効
{
    int result;
    int df;
    
    if(a==0 || a==(-1)) // 0または-1
    {
        result = 1;
    }
    else if(a%2==0 && a>0) // 正の偶数
    {
        df = 1;
        for(int i=2; i<=a; i=i+2)
        {
            df = df * i;
        }
        result = df;
    }
    else if(a%2==1 && a>0) // 正の奇数
    {
        df = 1;
        for(int i=1; i<=a; i=i+2)
        {
            df = df * i;
        }
        result = df;
    }
    else // -1未満の整数（負の奇数は定義できるが無視）。基本ないはずなので仮で0出力
    {
        result = 0;
    }
    
    return result;
}

// P_{m}^{m}(x)
function float Pmm(float x;int m)
{
    return pow(-1.0,m) * doubleFactorial(2*m-1) * pow(1.0-x*x,m/2.0);
}

// P_{m+1}^{m}(x)
function float PmmOne(float x;int m)
{
    return (2.0*m+1.0) * x * Pmm(x,m);
}

// P_{l}^{m}(x)
function float Pml(float x;int m;int l)
{
    // 実数球面調和関数ではmが負の時、その符号を正に直してルジャンドル陪関数に入力するので、この関数内ではmは常に0または正
    
    
    float pml;
    
    if( (l-m) == 0 ) // l=mならPmmの結果を使う
    {
        pml = Pmm(x,m); 
    }
    else if( (l-m) == 1 ) // l-1=mならPmmOneの結果を使う
    {
        pml = PmmOne(x,m);
    }
    else if( (l-m) < 0) // l<mの場合は未定義なので0を返す
    {
        pml = 0;
    }
    else // (l-m) > 1 の場合。つまり、lがmより2以上大きい
    {
        float pml_prev;
    
        for(int i = m; i <= l-2; i++)
        {        
            if(i==m) // ループの初回
            {
                float right = (2.0*(i+2.0)-1.0) * x * PmmOne(x,m) - ((i+2.0)+m-1.0) * Pmm(x,m); // PmmとPmmOneでPmmTwoに相当するものを計算
                    // i+2はm基準で今求めているlの値
                pml = right / ((i+2.0)-m);
            }
            else if(i==m+1) // ループの2回目
            {
                float right = (2.0*(i+2.0)-1.0) * x * pml - ((i+2.0)+m-1.0) * PmmOne(x,m); // PmmOneと初回ループの結果でPmmThreeに相当するものを計算            
                pml_prev = pml; // 次々回で使えるように退避
                pml = right / ((i+2.0)-m);           
            }   
            else // 3回目以降
            {
                float right = (2.0*(i+2.0)-1.0) * x * pml - ((i+2.0)+m-1.0) * pml_prev; // 前回と前々回のループの結果でPmmThree以降に相当するものを計算            
                pml_prev = pml; // 次々回で使えるように退避
                pml = right / ((i+2.0)-m);              
            }        
        }
    }
    return pml;
}

そしてこれを使って試しにLineのポイントを移動してみると、このような形になりました。
WolframAlphaで同じパラメータを指定してみたものと一致するので、正しいようですね。

	m=0	m=1	m=2	m=3
l=0
l=1
l=2
l=3

なお、ここまでの計算の内容は、より洗練されたコードが以下のPDFに載っているので、参考にする場合はそちらを見ることをおススメします……。

実数球面調和関数の全体

正規化定数の計算のために階乗の関数を追加する必要がありますが、それ以外は数式の通りに記述します。

実数球面調和関数

// 階乗
function int Factorial(int a)
{
    int factorial = 1;
    if(a==0)
    {
        factorial = 1;
    }
    
    for(int i=1; i<=a; i++)
    {
        factorial = factorial * i;
    }
    return factorial;
}

// Y_{l}^{m}(θ,Φ)
function float Yml(float theta; float phi; int m; int l)
{
    float normalizer = sqrt(
                            ( (2.0*l+1.0) / (4.0 * $PI) )
                            * ( float(Factorial( l-abs(m) )) / float(Factorial( l+abs(m) )) )
                            );    
    float yml;
    if(m<0)
    {
        yml =   pow(-1,m) * sqrt(2) * normalizer * sin(-m * phi) * Pml(cos(theta),-m,l); // 負のmが入ってきたら正に直す
    }
    else if(m==0)
    {
        yml =  normalizer * Pml(cos(theta),m,l);
    }
    else if(m>0)
    {
        yml =  pow(-1,m) * sqrt(2) * normalizer * cos(m * phi) * Pml(cos(theta),m,l);
    }
    
    return yml;
}

結果

ここまでの内容を元に、実際にSphereを変形させるVEXコードは以下です。

実数球面調和関数の可視化

// 球面調和関数で使う変数
int m = chi("m");
int l = chi("l"); // ルジャンドル陪関数の定義上、m<=lでなければならない

float theta = acos(v@P.y);
float phi = atan2(v@P.z,v@P.x);

// 実数球面調和関数の計算
float height =  Yml(theta,phi,m,l);

// Sphereを変形させる
v@P = v@P * abs(height); 
if(height>=0)
{
    v@Cd = set(1,0,0);
}
else
{
    v@Cd = set(0,1,0);
}

よく見るぷにぷにの図のようにするには、実数球面調和関数の結果の絶対値を取ります。
すると無事ぷにぷにの形が可視化できました！

注意点として、AttributeWrangleの「Bindings」タブから「VEX Precision」を64-bitにしておく必要があります。
これをやらないと、階乗計算が早い段階でオーバーフローしてしまい、lが高次の場合の結果が不正確になります。（64-bitにしても10次くらい？が限界のようですが……）

2. 実践編 球面調和展開

目標設定

実数球面調和関数が実装できたので、次はこれを実際に何かに利用することを考えていきたいと思います。
ここでは、入力形状をスケーリング係数（フーリエ変換の時のスペクトル的なもの）に圧縮後、スケーリング係数から元形状の近似を再構築したいと思います。
これができれば、3次元形状を表現する際のデータ量を大きく減らせるはずです。

全体的な処理の流れ

今回はぷにぷに形状を次々と変えながら、その一つ一つに対しスケーリング係数を求めていくので、二次元フーリエ変換の時のように二重ループで回していきます。

とりあえず次数lの最大は画像のように外部からいじれるようにしておきます。
このlについてのループが外側になります。

またlの値が大きくなるごとにピラミッド状にmを増やしていくので、内側のmについてのループは画像のような設定にします。

スケーリング係数を求める式

では本題に入っていきます。
入力形状を$f(\theta,\phi)$としたとき、スケーリング係数$c_{l}^{m}$を求める式は以下の通りです。

    c_{l}^{m} = \int_{0}^{2\pi}d\phi\int_{0}^{\pi}d\theta \sin\theta f(\theta,\phi) Y_{l}^{m}(\theta,\phi)

元の関数（入力形状）と実数球面調和関数の積の（球面上での）二重積分ということですね。
……はい、お恥ずかしながらどう計算すればいいかわかりません。。。

というわけで、
助けてコパえも～ん。

二重積分を代替する方法

Copilot　＜　モンテカルロ法というものがあります。

モンテカルロ法とは、
積分範囲内でランダムに点を選び、その点での関数の値の平均を取ることで積分を近似する方法とのことです。

これを使えば、二重積分の計算を代替できるというわけですね！
その場合の式は、（おそらく）こんな感じになるはずです。

    c_{l}^{m} = \frac{4\pi}{N}\sum_{i=1}^{N} f(\theta_{i},\phi_{i}) Y_{l}^{m}(\theta_{i},\phi_{i})

さて、要するに球面上にランダムにポイントをばらまいて、その合計を平均してから単位球の表面積を掛ければいいということでしょうから、これならScatterとかを駆使すればなんとかできそうです。

元形状を近似する式

そしてここで、スケーリング係数から近似形状を復元する式も見てみます。

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

スケーリング係数と球面調和関数を掛けたものをlとmについて単純に総和するだけです。
これはCosが球面調和関数になっただけで、フーリエ変換でやったのとほとんど同じですね。

結果

これで元形状を近似するネットワークは完成です！
画像のような構成になりました

shapeToSphere

vector origin = set(0,0,0);
// 各ポイントの原点からの距離
float height = length(v@P - origin);
// 形状をSphereに変える
v@P = normalize(v@P);

f@height = height;
v@Cd = height;

calcThetaPhi

float theta = acos(v@P.y);
float phi = atan2(v@P.z, v@P.x) + $PI;

f@theta = theta;
f@phi = phi;

calcSummetion

int m = detail('op:../foreach_count_m', 'value', 0);
int l = detail('op:../foreach_count_l', 'value', 0); // ルジャンドル陪関数の定義上、m<=lでなければならない

float phi = f@phi;
float theta = f@theta;

// 中略
// 実数球面調和関数の関数定義


float d = f@height // 圧縮したいデータ
            *  Yml(theta,phi,m,l) ;  
f@integral_piece = d;

calcSum

float sum = 0.0;
int lastPointNum = chi("SamplePoints"); 

// 積分をモンテカルロ法で近似
// 全ての点の座標を足す
for(int i=0;i<=lastPointNum-1;i++)
{
    float piece = point(1,"integral_piece",i);
    sum += piece;
}

// ポイントの数で割って平均をとる
sum /= lastPointNum;
// 積分範囲の面積を掛ける。ここれは単位球の表面積
sum *= 4.0 * $PI;


addpoint(geoself(),set(0,sum,0)); // スペクトル生成用に積分した値をアトリビュートに持ったポイントを作り直す

reconstructHeight

float phi = f@phi;
float theta = f@theta;

int lmax = chi("restoreLMax");


float approx_height = 0.0;

// 中略
// 実数球面調和関数の関数定義

for(int l=0; l<=lmax; l++)
{
    for(int m = -l; m<=l; m++)
    {
    
        // 対応するスケーリング係数のインデックスを求める
        int idx = 0 ;
        for(int i=1;i<=l;i++)
        {
            idx += 2*i+1;
        }
        idx -= (-m+l);
    
        // スケーリング係数を取り出す
        vector coefP = point(1,"P",idx);
        float coef = coefP.y;
    
        // 球面調和関数にスケーリング係数を掛けて足していく
        float scaledSH = Yml(theta,phi,m,l) * coef;
        approx_height += scaledSH;
    }
}

v@Cd = approx_height;
f@approx_height = approx_height;

setHeight

v@P = v@P * f@approx_height;

期待通り動作するか見てみましょう。

星

正八面体

捻じれた四角柱

テトラポット的な形

金平糖みたいな形

岩っぽいもの

今回は10次＝121通りの関数の組み合わせを最大としています。
元形状の角やシャープになっている部分は結構失われてしまっていますが、概ね近似できていそうです！
めでたしめでたし。

…と、はじめから意図通りの結果が得られたかのように書いていますが、実際はそううまいこと行きませんでした。
以下に、どのようなつまづきがあり、どのように解決したのかも書いておこうと思います。

ポイント①　正規直交

最初にネットワークを組み終わった時は、元形状とは似ても似つかない形になってしまっていました。
そこで何が原因なのか探るために行ったのは、実数球面調和関数同士が正規直交しているかを確認することです。

球面調和関数には、ある球面調和関数（とその複素共役）同士の積分は1になり、ある球面調和関数と別の球面調和関数（の複素共役）の積分は0になる、という性質があるとのこと。
これを「関数同士が正規直交している」というそうです。
複素共役が何者かについてはこの記事では触れません。

そして同様に、実数球面調和関数同士も互いに正規直交です。
（こっちの場合はそもそも複素共役とか関係ありません）
これは明言している資料を見つけられませんでしたが、WolframAlphaで実際に積分計算してみてそのようになっていることを確認しました。

つまり、実装に間違いがないのであれば、同じ実数球面調和関数同士の積を積分したら結果は1になるし、異なるもの同士なら0になるはずです。

このことを利用して、ミスがないかのチェックを行いました。
結果、複数のケアレスミスを見つけることができました。

ポイント②　変形した形状を元データにしてはいけない

当初は元形状をそのままスケーリング係数計算部分に入力していたのですが、これは間違いでした。
球面調和関数の各種の性質は、あくまで単位球面上のθとφで成り立つものなので、元形状からそれらを計算して使おうとすると意図しない結果になるようです。

なので、元形状を一度「原点からの距離」という情報に置き換えて、実際の形は単位球に変形させて、それを入力として各種の計算を行う必要があります。
具体的には、元形状の各ポイントのlength(v@P-原点)を適当な名前のアトリビュートに保存しておいて、v@Pは正規化して上書きしてしまえばよいです。

ポイント③　Scatterでバラまくポイントは充分に一様な分布でなければならない

ポイント②までの修正が済んだ時点で、結果は次の画像のようになっていました。

なんとか元形状の特徴が見えてきましたが、やけに歪です。
このことに対する主だった原因は、Scatterしたポイントが目ではわからない程度に偏って分布していたから、のようです。
モンテカルロ法で積分を近似する場合には、各サンプリング点が一様に分布していないといけない、という前提があるのですが、そこに抵触してしまっていたということですね。
ScatterのRelax Iterationsを極端に大きな値(1000くらい)にすることで、解決できました。

補足

ちなみに、球面調和関数はただのスカラー値なので、あまり複雑な形状を近似することはできないようです。
具体的には、Sphereの各ポイントを原点方向に動かしてできるような形状以外では望み通りの結果になりません。
実際にPigHeadやRubberToyのような形状を近似しようとすると……。
元形状の特徴は残っていますが、近似できているとは言えなそうな感じですね。PigHeadに至っては鼻が大きく突き出してしまっています。口内までポリゴンが作られてる影響でしょうかね？

課題

ある程度近似できていそうとはいえ、復元した結果にデコボコが残ってしまっているのが気になる点です。
また、元形状のシャープな部分がかなり消えてしまっていたのも残念ですね。各種の論文の結果画像を見る限り、10次ならもう少し維持されると思ったのですが……。
このあたり、改善方法について何か知見をお持ちの方がいらっしゃれば是非教えていただきたいです……！

終わりに

今回は球面調和関数をHoudiniで実装する過程を紹介しました。
かなり難易度の高い挑戦でしたが、なんとか形にはなっていると思います。
自信もつきました。
個人的にはMicrosoft Copilotも大活躍でした。
いい時代になりましたね……。

球面調和関数についてはWeb上の資料があまり多くないと思うので、本記事が皆さんのお役に立てたら嬉しいです。

参考

• t-pot『テイラー、フーリエ、球面調和関数』
• ルジャンドル陪関数の計算
• ☆PROJECT ASURA☆ [OpenGL] 『球面調和関数(1)』
• docs/CGのための球面調和関数.pdf at master · mebiusbox/docs · GitHub
• 球面調和関数(Spherical Harmonics)とCGへの応用
• 球面調和関数変換の理論と実装 #C++ - Qiita
• Blenderでお手軽に球面調和関数を描こう｜Melville
• ioa.s.u-tokyo.ac.jp/~kkohno/seminar2/?plugin=attach&refer=H28_KadaiKenkyu&openfile=SphericalHarmonics.pdf