バーンズリーのシダ(バリエーション)

　「バーンズリーのシダ」では、「さっぱりわかりません」で終わらせていた以下の生成規則について掘り下げてみたいと思います。



　各々、大雑把には以下の様な部分を形成するための数式となっています。

• (x’，y’)＝f1　：　茎、確率＝1%
• (x’，y’)＝f2　：　左の葉、確率＝7%
• (x’，y’)＝f3　：　右の葉、確率＝7%
• (x’，y’)＝f4　：　細かいひだ、確率＝85%

　f1〜f4の各数式は、以下の様な形をしています。



　2×2行列(a〜d)は回転と拡大・縮小変換、加算項(e，f)は平行移動の、各アフィン変換を実施するパラメータです。

アフィン変換

　加算項(平行移動、e，f)の影響は自明なので、解析すべきは変換行列(a〜d)の方です。回転や相似変換などの変換行列は以下の様な形をしています。

• 回転変換行列：Mr(θ)
• 原点の周りに反時計回りに角度θで回転する変換行列です。
• 相似変換行列：Mk(Kx，Ky)
• 原点を中心に、x軸方向へ×Kx倍、y軸方向へ×Ky倍の拡大・縮小を実施する行列です。
• シアー変換行列：Ms(Sx，Sy)
• 斜めに歪ませるような変換を実施する行列です。

　バーンズリーのシダを生成するための数式f1〜f4で使われている行列(a〜d)は、上の様なアフィン変換の合成だと考えられます。アフィン変換を合成するには、各行列の積を取ります。
　例えば、回転(Mr)して→拡大・縮小(Mk)して→シアー変換(Ms)する場合、その変換行列Mt＝Ms(Sx，Sy)×Mk(Kx，Ky)×Mr(θ)となります。行列は掛け算の順番で結果が異なるので、演算順番は重要です。
　変換順番が、拡大・縮小(Mk)して→シアー変換(Ms)して→回転(Mr)する、だった場合、その変換行列は、Mt’＝Mr(θ)×Ms(Sx，Sy)×Mk(Kx，Ky)と、なって、一般にその計算結果は異なります(Mt≠Mt’)。

変換行列の解析

　本ページでは、バーンズリーのシダの変換行列を上例のMtの順番だと仮定して、そこからθ，Kx，Ky，Sx，Syの各パラメータを逆算してみたいと思います。

　変換行列は2行×2列なので、変換行列(a〜d)＝Mtとおくと、4個の方程式を立てることが出来ます。これに対して未知数はθ，Kx，Ky，Sx，Syの5個と、これでは方程式の数が足りず、解が求まりません。
　そこで、シアー変換のSy＝0に固定してしまいます。Sy＝0に固定しても、恐らく必要な自由度を失わないので、問題ないでしょう。これによって未知数が4個となり、4元4連立方程式が立ちますので、バーンズリーのシダの生成規則からθ，Kx，Ky，Sxのパラメータを逆算することが出来ます。

　さてさて、言うは易しですが、まず、4×4行列Mtの計算がしんどいです。その結果立てた4元4連立方程式も、三角関数を含む方程式なので、解くのはそこそこ大変です。

Mathematicaを使って解く

　そこで、前節の代数計算をMathematicaと言うソフトウェアにやらせようと思います。Mathematicaや、それを実行するための超小型コンピュータRaspberry Piについては、以下の記事をご参照ください。

• Raspberry Piの導入に関して：「デザートにラズベリーパイはいかが？」
• Mathematicaの使い方について：「Mathematicaを使ってみよう」

　Mathematicaで変換行列Mtを計算させた様子は以下の通りです。

• 回転変換行列Mrと、相似変換行列Mk、シアー変換行列Msを関数として定義した後、
• バーンズリーのシダ生成規則に登場する2×2変換行列をmf1〜mf4として宣言、
• 最後に合成変換行列Mt＝Ms(ss)×Mk(kx，ky)×Mr(th)を計算させました。
• 薄い字で「Out[8]＝」となっているところがMtを計算した結果です。
• ※都合により、θ⇒th、Sx⇒ssと打ち換えています。

　これで準備が完了しましたので、MathematicaのSolve[ ]関数を使って、4元4連立方程式を4回解きます(↓)。



※上のSolve[ ]一回で、1つの4元4連列方程式を解いています。「4元4連列方程式を解くためにSolve[ ]を4回実行している」のではありません！Solve[ ]が4回なので、4元4連列方程式を4回解いています。

　それぞれ、複数個の解が示されました。どれを使っても同じですが、θ(＝th)の絶対値が最小となる解を選んで使用しました。……ちなみに、上の4元4連立方程式×4回を手計算で解くのは御免蒙りますｍ（１＿＿）ｍ。

Xi-TRACEデータの作成

　意味不明だったバーンズのシダの生成規則(変換行列)を、回転変換、相似変換、シアー変換に分解することが出来ましたので、各々のパラメータがシダの葉の生成にどのように影響を与えるのか？各パラメータを少しづつ変更して見れば、知ることが出来るはずです。

　そこで、th，kx，ky，ssのパラメータから、一次変換行列を生成する手続きを作りました(↓)。

　　procedure　bsmat　th, kx,ky, ss, mx,my
　　{
　　　　para　mx,my,0;
　　　　matrix　1, 　　0,　　0,　　0,
　　　　　　　　ss,　　1,　　0,　　0,
　　　　　　　　0, 　　0,　　1,　　0,
　　　　　　　　0, 　　0,　　0,　　1;
　　　　sim　kx,ky,1;
　　　　rotz　th*180/PI;
　　}

　回転変換(rotz命令)と相似変換(sim命令)は、Xi-TRACEの命令を使って実装しました。シアー変換はXi-TRACEの命令ではサポートが無いので、matrix命令を使って等価な変換を実施しています。(mx，my)は、平行移動(para命令)のパラメータです。

　上の手続きを使うと、例えば、f2の変換処理は以下の様になります。

　　＝＝＝＝旧(↓)
　　　　：
　　elif (rn<pr1 + pr2)
　　{
　　　　// f2: 左葉
　　　　matrix　0.20,　-.26,　0,　0.00,
　　　　　　　　0.23,　0.22,　0,　1.60,
　　　　　　　　0, 　　0, 　　1,　0,
　　　　　　　　0, 　　0, 　　0,　1;
　　　　linear　nx,ny,nz;
　　}
　　　　：

　　＝＝＝＝新(↓)
　　　　：
　　elif (rn<pr1 + pr2)
　　{
　　　　// f2: 左葉
　　　　static : bsmat　0.915101,　0.328024,0.31644　-0.104089,　0,1.60;
　　　　linear　nx,ny,nz;
　　}
　　　　：

　手続きbsmatに、Mathematicaで求めたth＝0.915101，kx＝0.328024，ky＝0.31644，ss＝−0.104089及び、平行移動項(0，1.60)を渡して変換行列を生成するように変更しました。

パラメータを変更して確かめる

　意味不明だった変換行列に、一応理屈が通りました。
　ここから先は、地道にパラメータを変更して見て、バーンズリーのシダがどのような形に変形するのか？を解析して行きます。

• 右葉を形成するf3の発生確率pr3に関して：
• pr3＝0%に設定して描画してみたところ、右葉が描画されなくなりました。
• 自己相似図形の右半分もまた、描画されなくなります。
• 0%＜pr3≦7%の範囲では、pr3→大となるにつれ、葉の形が段々濃くなってくるように変化しました。
• また、葉の形は、先端部が先に濃くなって、根元に近い部分は後から濃くなるイメージです。
• 他の部分の確率pr1，pr2に関しても、概ね同じです。
• pr1＝0だと茎の部分が消滅します。
• pr2＝0だと左葉が消滅します。
• pr4→0とすると、小葉が生成されないことで、形が判然としなくなってしまいます。
• 左右の葉を形成するf2，f3に関して：
• 平行移動項f2(e，f)＝(0，1.6)とf3(e，f)＝(0，0.44)は、y軸方向の平行移動で、1.6および0.44と不均衡に設定されています(オリジナル)。
• これをy軸方向に1.0に揃えてみると、左右の葉が同じ位置から発生しました。
• つまり、y軸方向の平行移動は、左・右葉の出る位置を決めているようです。
• 相似変換Kxの符号が左・右を決めていることも判りました。
• また、Kxが葉の長さ、Kyが葉の幅を決めるパラメータになっています。
• 回転変換θは、葉の延びる角度(方向)を決定します。
• 小葉を形成するf4に関して：
• 回転角度θがシダの葉の傾きになっています。
• バーンズリーのシダのオリジナルパラメータではθ≒−0.047ラジアン(約−5.4°)ですが、この状態では全体が右側に傾いています。
• θ＞0にすると、逆に左側に傾きます。ただし、自己相似図形の葉も逆側(下側)に反ってしまうため、「シダの葉」とはイメージが変わってしまいます。
• θ＝0にすると、真っすぐな葉ができます。同じく自己相似図形の葉もすべて真っすぐになります。
• 平行移動成分(y軸方向)を変化させると、全体の大きさが変化します。
• 右葉、左葉を個別に操作しようとしても、全ての自己相似図形中の右葉と左葉が一斉に変化してしまいます。
• 従って、画像の右半分と左半分を個別に操作することはできません。

　解析したところ、ざっとこんな所(↑)です。言葉ばかりでわかり難いと思いますので、下に示す動画でシダの葉が変化する様子を確かめてみてください。

完成データ

　完成したデータをリンクしておきます。

　データを実行すると、以下の様なアニメーションデータを作成することが出来ます。生成式の変換行列をダイナミックに変化させながら、シダの形がどうなるかを動画にしてあります。



　2Dのグラフィックですが、計算には結構な時間がかかりました。Xi-TRACEのデータ(スクリプト)はインタプリタ方式で遅いですし、高解像度でそれなりにシダらしき形状にするためには繰り返し回数NNを大きくとる必要があったためです。さすがに、これでレイトレーシングする気にはなれませんでした（１＾＾；）。