単位球面に対してランダムに点を生成する

概要

単位球面に対してランダムに点を生成する方法はたくさんありますが、Marsagliaの方法が大変エレガントで感動したので、それをまとめておこうと思います。

単純な方法

一つの単純な方法としては、単位球内の乱数を生成してそれを正規化することです。単位球内の乱数は、単に単位球外を棄却する方法が使えます。例えば以下のような実装が考えられます。

inline glm::dvec3 sample_on_unit_sphere(PeseudoRandom *random) {
    glm::dvec3 d;
    double sq = 0.0;
    do {
        d.x = random->uniform(-1.0, 1.0);
        d.y = random->uniform(-1.0, 1.0);
        d.z = random->uniform(-1.0, 1.0);

        sq = glm::length2(d);
    } while (sq < 1.0e-9 || 1.0 < sq);
    d /= glm::sqrt(sq);
    return d;
}

非常に直感的で、めちゃくちゃパフォーマンスが悪いかといえば、そうでもないでしょう。 ただ実際成功する確率はそこそこ低く、2x2x2の立方体の体積部の単位球の体積なので、

 
\displaystyle {
P_{ inside }=\frac { \frac { 4 }{ 3 } \pi  }{ { 2 }^{ 3 } } =\frac { 1 }{ 6 } \pi \simeq 0.523598776
}

となり、わりと棄却されて乱数生成が無駄になっていることがわかります。 実際のところ何個程度の乱数生成が必要なのか?というのの期待値を考えてみると、 単位球内の乱数生成に失敗し、乱数生成が継続される確率を、

 
\displaystyle {
1-\frac { 1 }{ 6 } \pi =r
}

とおくと、失敗する場合乱数が3個消費されるので、1回失敗する場合、2回失敗する場合、・・・・無限回失敗する場合全部を考えれば良いので乱数生成回数の期待値は、

 
\displaystyle {
N_{ avg }=3{ r }^{ 0 }+3{ r }^{ 1 }+3{ r }^{ 2 }+...3{ r }^{ \infty  }\\ =3\left( { r }^{ 0 }+{ r }^{ 1 }+{ r }^{ 2 }+...{ r }^{ \infty  } \right) \\ =3\sum _{ i=1 }^{ \infty  }{ r } \\ =3\left( \frac { 1 }{ 1-r }  \right) \\ =3\left( \frac { 1 }{ 1-\left( 1-\frac { 1 }{ 6 } \pi  \right)  }  \right) \\ =3\left( \frac { 1 }{ \frac { 1 }{ 6 } \pi  }  \right) \\ =\frac { 18 }{ \pi  } \simeq 5.72957795
}

となり、実際理想的な乱数生成の個数が最低2回であることを考えると、ちょっと無駄が多い気もします。乱数生成の負荷が高い場合は、少し気になるかもしれません。(ただ乱数生成に計算負荷の高いものを使う場合、そもそも棄却するような手法を使わない気がしますが・・・)

Marsagliaの方法

Marsagliaの方法は、2つの事実を利用して、目的を達成します。

A. 単位球面に対してランダム生成された点 は、それぞれ単一の軸にだけ着目すると-1~+1で密度は均一である

これはいろいろな証明方法があるようです。なれている方法で確認するとして、積分を使った方法で考えてみます。 下の図のような  \theta_{end} で切り取られた単位球の表面積を考えてみます。

f:id:ushiostarfish:20180712232315p:plain

これは半球積分の範囲をちょっと調整すると、

 
\displaystyle {
\int _{ 0 }^{ 2\pi  }{ \int _{ 0 }^{ \theta _{ end } }{ \sin { \theta  } d\theta  } d\phi  } \\ =2\pi \int _{ 0 }^{ \theta _{ end } }{ \sin { \theta  } d\theta  } \\ =2\pi \int _{ 0 }^{ \theta _{ end } }{ \sin { \theta  } d\theta  } 
}

のように考えることができます。 ここで、

 
\displaystyle {
x=\cos { \theta  } \\ \frac { dx }{ d\theta  } =-\sin { \theta  } \\ -dx=\sin { \theta  } d\theta 
}

の変数変換を施して、

 
\displaystyle {
=2\pi \int _{ \cos { 0 }  }^{ \cos { \theta _{ end } }  }{ -x } \\ =2\pi \int _{ 1 }^{ \cos { \theta _{ end } }  }{ -x } \\ =2\pi \left( \cos { \theta _{ end } } -1 \right) \left( -1 \right) \\ =2\pi \left( 1-\cos { \theta _{ end } }  \right) 
}

もちろん  \theta_{end} = 0 のとき 0であり、 \theta_{end} = \pi のときは  4 \pi であることはすぐにわかります。

ここで、 1-\cos { \theta _{ end } } はなにかと考えてみると、縦の長さと捉えることができます。

f:id:ushiostarfish:20180712233531p:plain

一方  2 \pi はなにかといえば、円周であり、 つまり、単位球と、単位球をピッタリとすっぽりと覆う円柱の側面積とが-1~+1の 任意の  \cos \theta_{ end } で一致することになります。

f:id:ushiostarfish:20180712234410p:plain

積分値が任意の範囲で一致するということはやはり微分値でも一致するので、 言い換えれば、上の図のような配置を考える場合は、水平にどこで輪切りにしても、球面の輪切りと対応する円柱の側面積が一致するということであり、 つまり面積の密度は-1~+1で均一であり、これはつまり確率も均一であるということがわかります。

※ここでは単位球で考えましたが、この性質は任意の半径で成り立ちます。

B. 単位円に対してランダム生成された点 (a, b) の  a^2+b^2 は0~1で密度は均一である

あるrの円に、ある点(x, y) が含まれる条件は、

 
\displaystyle {
{ x }^{ 2 }+{ y }^{ 2 }\le { r }^{ 2 }
}

となります。 ここで、単位円に対してランダム生成された点 (x, y) が、ある半径r (0~1)の円に含まれている確率は、

 
\displaystyle {
P\left( x^{ 2 }+y^{ 2 }\le r^{ 2 } \right) =\frac { \pi r^{ 2 } }{ \pi  } \\ P\left( x^{ 2 }+y^{ 2 }\le r^{ 2 } \right) =r^{ 2 }
}

と考えることができます。ここで、

 
\displaystyle {
x^{ 2 }+y^{ 2 }=X\\ r^{ 2 }=U
}

とでもおいてみると、

 
\displaystyle {
P\left( X\le U \right) =U
}

ここで、Xの範囲は0~1で、Uの範囲も0~1です。 そう、つまりこれはXが 0~1に均一に分布していることを示しているわけです。

組み合わせる

単位円に均一に発生させた点を \left( u_{ 1 },u_{ 2 } \right) とすると、 まず簡単にわかることは、zについては、

 
\displaystyle {
z=1-2\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right) 
}

とすれば、zについては均一に-1~+1の範囲にできるので、これでOKということになります。 ここで、 u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } の成分、つまり大きさ成分を使ってしまったので、正規化して向きだけのものを考えると、

 
\displaystyle {
\vec { u } =\frac { \left( u_{ 1 },u_{ 2 } \right)  }{ \sqrt { u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } }  } 
}

ここで、球面上に乱数を発生させる場合の、x, y平面上のでの半径  r_{ xy } は、

f:id:ushiostarfish:20180713002206p:plain

 
\displaystyle {
r_{ xy }=\sqrt { 1-{ z }^{ 2 } } 
}

であり、これは具体的には、

 
\displaystyle {
r_{ xy }=\sqrt { 1-{ z }^{ 2 } } \\ =\sqrt { 1-{ \left( 1-2\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  \right)  }^{ 2 } } \\ =\sqrt { 1-{ \left( 1-4\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right) +4{ \left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  }^{ 2 } \right)  } } \\ =\sqrt { 1-1+4\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right) -4{ { \left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  }^{ 2 } } } \\ =\sqrt { 4\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right) -4{ { \left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  }^{ 2 } } } \\ =2\sqrt { \left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right) -{ { \left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  }^{ 2 } } } \\ =2\sqrt { \left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right) \left( 1-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  \right)  } \\ =2\sqrt { \left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  } \sqrt { 1-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  } 
}

となるので、実際にほしい x, y は結局、

 
\displaystyle {
\vec { u_{ xy } } =\frac { \left( u_{ 1 },u_{ 2 } \right)  }{ \sqrt { u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } }  } r_{ xy }\\ =\frac { \left( u_{ 1 },u_{ 2 } \right)  }{ \sqrt { u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } }  } 2\sqrt { \left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  } \sqrt { 1-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  } \\ =2\left( u_{ 1 },u_{ 2 } \right) \sqrt { 1-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  } 
}

となるわけです。 個別にかけば、

 
\displaystyle {
x=2u_{ 1 }\sqrt { 1-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  } \\ y=2u_{ 2 }\sqrt { 1-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  } \\ z=1-2\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right) 
}

なんてシンプルでしょうか!

プログラミング

以上をふまえ、コードを考えてみます。まずは棄却法をつかって、単位円に均一な乱数を生成します。

inline glm::dvec3 sample_on_unit_sphere(PeseudoRandom *random) {
    double x1;
    double x2;
    double S;
    do {
        x1 = random->uniform(-1.0, 1.0);
        x2 = random->uniform(-1.0, 1.0);
        S = x1 * x1 + x2 * x2;
    } while (S >= 1.0);

するとこの時点で、 u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } がすでに計算されていることに気づきます! なのでこれをそのまま使って、

inline glm::dvec3 sample_on_unit_sphere(PeseudoRandom *random) {
    double x1;
    double x2;
    double S;
    do {
        x1 = random->uniform(-1.0, 1.0);
        x2 = random->uniform(-1.0, 1.0);
        S = x1 * x1 + x2 * x2;
    } while (S >= 1.0);

    double two_sqrt_one_minus_s = 2.0 * std::sqrt(std::max(1.0 - S, 0.0));
    return glm::dvec3(
        x1 * two_sqrt_one_minus_s,
        x2 * two_sqrt_one_minus_s,
        1.0 - 2.0 * S);
}

エレガント!

f:id:ushiostarfish:20180713005833p:plain

乱数生成回数の期待値

単位円内の乱数を棄却法で生成する場合、成功する確率は、

 
\displaystyle {
P_{ inside }=\frac { \pi  }{ { 2 }^{ 2 } } =\frac { \pi  }{ { 4 } } \simeq 0.785398163
}

であり、乱数生成が継続される確率を、

 
\displaystyle {
1-\frac { \pi  }{ { 4 } } =r
}

とおくと、乱数生成回数の期待値は、

 
\displaystyle {
N_{ avg }=2{ r }^{ 0 }+2{ r }^{ 1 }+2{ r }^{ 2 }+...2{ r }^{ \infty  }\\ =2\left( { r }^{ 0 }+{ r }^{ 1 }+{ r }^{ 2 }+...{ r }^{ \infty  } \right) \\ =2\sum _{ i=1 }^{ \infty  }{ r } \\ =2\left( \frac { 1 }{ 1-r }  \right) \\ =2\left( \frac { 1 }{ 1-\left( 1-\frac { \pi  }{ { 4 } }  \right)  }  \right) \\ =2\left( \frac { 1 }{ \frac { \pi  }{ { 4 } }  }  \right) \\ =2\left( \frac { 4 }{ \pi  }  \right) \\ =\frac { 8 }{ \pi  } \simeq 2.54647909
}

であり、最初の方法に比べると乱数生成もとても少なくて済みます。

半球の場合

ついでに、上のアイディアをz+の半球の場合考えてみます。 一旦zを

 
\displaystyle {
z=1-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right) 
}

として使うことにします。(わざわざ逆転させているのは、逆転させないと式が単純にならないためです) とすると、 r_{ xy } は、

 
\displaystyle {
r_{ xy }=\sqrt { 1-{ z }^{ 2 } } \\ =\sqrt { 1-{ \left( 1-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  \right)  }^{ 2 } } \\ =\sqrt { 1-\left( 1-2\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right) +{ \left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  }^{ 2 } \right)  } \\ =\sqrt { 2\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right) -{ \left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  }^{ 2 } } \\ =\sqrt { \left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right) \left( 2-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  \right)  } \\ =\sqrt { u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } } \sqrt { 2-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  } 
}

ということはつまり、

 
\displaystyle {
\vec { u_{ xy } } =\frac { \left( u_{ 1 },u_{ 2 } \right)  }{ \sqrt { u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } }  } r_{ xy }\\ =\frac { \left( u_{ 1 },u_{ 2 } \right)  }{ \sqrt { u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } }  } \sqrt { u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } } \sqrt { 2-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  } \\ =\left( u_{ 1 },u_{ 2 } \right) \sqrt { 2-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  } 
}

となり、とってもシンプルになりました! まとめると、

 
\displaystyle {
x=u_{ 1 }\sqrt { 2-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  } \\ y=u_{ 2 }\sqrt { 2-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right)  } \\ z=1-\left( u_{ 1 }^{ 2 }+u_{ 2 }^{ 2 } \right) 
}

これをコードにすると、

inline glm::dvec3 sample_on_unit_hemisphere(PeseudoRandom *random) {
    double x1;
    double x2;
    double S;
    do {
        x1 = random->uniform(-1.0, 1.0);
        x2 = random->uniform(-1.0, 1.0);
        S = x1 * x1 + x2 * x2;
    } while (S >= 1.0);

    double c = std::sqrt(2.0 - S);
    return glm::dvec3(
        x1 * c,
        x2 * c,
        1.0 - S);
}

となります。さらに計算が減って効率が良くなり、いい感じです!

f:id:ushiostarfish:20180714144957p:plain

参考文献

GEORGE Marsaglia, "CHOOSING A POINT FROM SURFACE OF A SPHERE"

geometry - 3 random numbers to describe point on a sphere - Mathematics Stack Exchange