Blenderスクリプトでフラクタル立体モデルを生成する

　Blenderでのモデリングの操作にもある程度慣れて、そろそろ普通にモデリングするだけじゃなくてスクリプトも活用してみたいなあ、ということで、Pythonでフラクタル図形の生成に挑戦してみようと思いました。
　バージョンについては、本記事のコードはBlender 2.79aで動作確認をしています。

# スクリプトで動かすBlender
　例えば、以下のようなPythonコードをBlenderのテキストエディタにペタッと貼って「スクリプト実行」をポチると三角形の面一つからなるオブジェクトが生成されます。

```
import bpy
import mathutils
# 初期化

verts = [mathutils.Vector([0,0,0]),
         mathutils.Vector([1,0,0]),
         mathutils.Vector([0,1,0])] #頂点のリスト
fIndexes = [[0,1,2]] #面のリスト
mesh = bpy.data.meshes.new('さんかく')
mesh.from_pydata(verts,[],fIndexes) #点と面の情報からメッシュを生成

obj = bpy.data.objects.new('さんかく', mesh) #メッシュ情報を新規オブジェクトに渡す
bpy.context.scene.objects.link(obj) #オブジェクトをシーン上にリンク

obj.select = True #作ったオブジェクトを選択状態に
```
頂点のリストvertsには、Blenderの３Ｄビュー上の座標のリストが入ります。座標はVectorオブジェクトで入れておきます。

面のリストfIndexesには、「面を構成する頂点のインデックスのリスト（またはタプル）」からなるリストを入れます。サンプルではvertsの0番目・1番目・2番目からなる三角形を作りたいので、0,1,2からなるリストが1個だけからなるリストを定義しています。
この頂点リストverts、面リストfIndexesをmesh.from_pydata()に渡してメッシュを生成してオブジェクトを作るという感じです。

実行時のイメージがこんな感じです。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272494/c0716f3f-a76f-a2a5-5510-b34b4d7040b1.jpeg" width=50%>

ということで、Pythonスクリプト中でフラクタル図形を作る時も、こんな風に面を一つ定義するたびに頂点リストに点の座標を追加していき、その点集合に対応したインデックスのリストを面リストに追加していく感じで作っていきます。

# フラクタルとは
　例えばその辺に生えている木とか、野菜のロマネスコみたいに、そいつの一部分だけに注目してみると、そいつ自身の構造がそのまま表れていて、さらにその中の一部分を見ても、また同じ構造が現れて・・・といった感じで同じ構造が繰り返し現れる性質のものをフラクタルといいます。
今回作るのは例に挙げたような複雑度の高いものではないのですが、文字だけで伝えるのは難しいので図を交えて作り方を説明していきます。

# 三角形から作るフラクタル
今回は、正三角形をなす頂点に対して以下のルールを繰り返して適用してフラクタルを作っていきます。

＜生成ルール＞
正三角形の頂点P0,P1,P2に対し、
P3,P4,P5をそれぞれP0-P1間、P1-P2間、P2-P0間の中点で、
P6,P7,P8をそれぞれP3-P4間、P4-P5間、P5-P3間の中点、
P9を、P6,P7,P8との位置関係が正四面体の頂点をなすように定義します。

図に示すとこんな感じです。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272494/6f4b6ef3-8857-9a0f-acb2-87a51f8f0d1b.jpeg" width=50%>
ついでに上からの図も。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272494/92165e17-2bbc-fba7-8906-920aa7ba3d49.jpeg" width=50%>

このルールをただ1回適用して面を貼るだけではフラクタルにはなりませんが、このP0からP9までの座標の中にも、正三角形をなす点の組があります。上から見た図の青字で示しているように、元の三角形の半分の縮尺の正三角形が3つ、その更に半分の縮尺のが6つ現れているのが見て取れると思います。

この9つの三角形それぞれに対し、上に示した＜生成ルール＞を再帰的に定義していきます。

＜生成ルール＞を2段階目まで適用するとこうなります。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272494/484f602f-c6be-27b9-aafd-07ef43f2ae53.jpeg" width=50%>

9つの三角形それぞれに対し、また新たな三角形が9つずつ定義されるので、9×9=81個の三角形で面が作られました。

この生成ルールを何度も繰り返していくことで、細かいトゲトゲがびっしり生えたフラクタルが出来上がるというわけです。



これをコードで実現するとこんな感じになります。

```
import bpy
import math
import mathutils
# 初期化

def normal(face): #法線ベクトル
    v1 = face[1] - face[0]
    v2 = face[2] - face[0]
    normal = v1.cross(v2)
    normal.normalize()#正規化
    return normal

def tetraFractal(level, face): #level: フラクタルの複雑度 face: フラクタル生成のターゲットとなる面
    if level <= 0:# 頂点と面の情報をリストに追加
        global idxV
        global verts
        global fIndexes

        verts.append(face[0])
        verts.append(face[1])
        verts.append(face[2])
        fIndexes.append((idxV,idxV+1,idxV+2))
        idxV += 3
    else: # 一つ下のレベルに分解
        pts = face[:] #p0~p2
        pts.append((pts[0] + pts[1])/2) #p3
        pts.append((pts[1] + pts[2])/2) #p4
        pts.append((pts[2] + pts[0])/2) #p5
        pts.append((pts[3] + pts[4])/2) #p6
        pts.append((pts[4] + pts[5])/2) #p7
        pts.append((pts[5] + pts[3])/2) #p8
        nrm = 0.5*math.sqrt(1/6.0)*(pts[0]-pts[1]).length * normal(face)
        pts.append((pts[0] + pts[1] + pts[2])/3 + nrm) #p9

        level -= 1 # 
        tetraFractal(level,[pts[0],pts[3],pts[5]])
        tetraFractal(level,[pts[1],pts[4],pts[3]])
        tetraFractal(level,[pts[2],pts[5],pts[4]])
        tetraFractal(level,[pts[3],pts[6],pts[8]])
        tetraFractal(level,[pts[4],pts[7],pts[6]])
        tetraFractal(level,[pts[5],pts[8],pts[7]])
        tetraFractal(level,[pts[6],pts[9],pts[8]])
        tetraFractal(level,[pts[6],pts[7],pts[9]])
        tetraFractal(level,[pts[7],pts[8],pts[9]])
#
# グローバル変数定義
#
verts = [] #頂点のリスト
idxV = 0 # 頂点のインデックス番号
fIndexes = [] #面のリスト
mesh = bpy.data.meshes.new('tetraFractal') #メッシュ情報

# フラクタルを作成する起点となる三角形の頂点座標
tetraV  = [mathutils.Vector([math.cos(0) , math.sin(0) , 0.0]),
     mathutils.Vector([math.cos(2*math.pi/3) , math.sin(2*math.pi/3), 0.0]),
     mathutils.Vector([math.cos(4*math.pi/3) , math.sin(4*math.pi/3), 0.0])
     ]
#
# メインルーチン
#
tetraFractal(4,tetraV)
mesh.from_pydata(verts,[],fIndexes) #　追加した点と面情報からメッシュを定義

obj = bpy.data.objects.new('tetraFractal', mesh) #メッシュ情報からオブジェクトを生成
bpy.context.scene.objects.link(obj)

obj.select = True
```

実行結果：
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272494/4e8884c6-7386-8e27-04cf-37ec514b1b35.jpeg" width=50%>

このコードで肝となっている関数tetraFractal()を呼び出す事で、頂点リストと面リストにフラクタル図形を構成している頂点と面の情報がどんどん追加されていきます。一つ目の引数levelに、上で示した生成ルールを適用する回数を入れ、faceには3角形の座標をしめすVectorオブジェクトからなるリストを入れます。

tetraFractal()のelse文のブロック内のpts[0]~pts[9]には、上の説明で言うP0~P9にあたる座標を計算して入れます。そして、このptsから上の説明で述べた9つの三角形にあたる座標の組み合わせを選び出し、その三角形それぞれに対して更にtetraFractal自身を呼び出しています。呼び出す前にlevelを1つ下げているので、最終的にはlevelが0になって頂点と面の情報をリストに追加する処理が呼び出されて終了するわけです。

上のコードで「# フラクタルを作成する起点となる三角形の頂点座標」の部分で、三角形の座標tetraVを定義してtetraFractalを呼び出してるわけですが、その部分を以下の様に書き換えて、四面体をなす4つの3角形に対して実行してみます。

```
# フラクタルを作成する起点となる四面体の頂点座標
tetraV  = [mathutils.Vector([math.cos(0) , math.sin(0) , 0.0]),
     mathutils.Vector([math.cos(2*math.pi/3) , math.sin(2*math.pi/3), 0.0]),
     mathutils.Vector([math.cos(4*math.pi/3) , math.sin(4*math.pi/3), 0.0]),
     mathutils.Vector([0, 0, -1*math.sqrt(2)])
     ]

# tetraVからなる四面体の面情報
tetraF = [[tetraV[0],tetraV[1],tetraV[2]],
        [tetraV[1],tetraV[0],tetraV[3]],
        [tetraV[2],tetraV[1],tetraV[3]],
        [tetraV[0],tetraV[2],tetraV[3]]]

#
# メインルーチン
#
for i in range(4):# 四面体をなす4つの三角形に対して実行
    tetraFractal(4,tetraF[i])
```
実行結果：
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272494/85b4d387-7f49-9d74-3e16-deb84affd46f.jpeg" width=50%>
手に握ったらものすごく痛そうな物体ができましたね。

この関数は、levelの値によって計算量が爆発的に増えるので注意が必要です。今はlevelは4に設定していますが、7くらいでもう自分の環境では実行に2～3分くらいかかってしまいますし、作られる面の数も半端ない事になります。スペックやメモリの容量と相談した上でlevelの値は設定してください。

## 補足説明
例に挙げたコード中の関数tetraFractal()の引数faceに入れるVectorオブジェクトは、こんな感じで普通の数値型の変数を扱う感覚で加減算や定数倍などが計算できます。

![コンソール.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272494/35a430a8-a741-5c9b-4ec1-eab71d5e21cf.jpeg)


2つの座標間の差を計算してlengthメソッドを呼び出せば距離も計算できます。
![コンソール2.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272494/a30e48e0-193a-2911-b4f7-2fa22b59ea9a.jpeg)

また、＜生成ルール＞中のP9の座標を計算するために、p0~p2からなる三角形の面に垂直な、法線ベクトルの計算が必要になるので、外積を使って定義してます。下記の部分です。

```
def normal(face): #法線ベクトル
    v1 = face[1] - face[0]
    v2 = face[2] - face[0]
    normal = v1.cross(v2)
    normal.normalize()#正規化
    return normal
```
4行目のv1.cross(v2)は外積v1×v2を計算してます。
5行目のnormalizeは正規化、つまり長さを1に揃えています。
図に示すとこんな感じ。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/272494/b7d5f5e9-ddc0-0ab2-399e-08a495cd1d94.jpeg" width=50%>

これで法線ベクトルが求められるので、これを適当な長さで定数倍して、あとは三角形の中央の座標に足し合わせる事でP9を求める事ができます。この部分ですね。

```
nrm = 0.5*math.sqrt(1/6.0)*(pts[0]-pts[1]).length * normal(face)
pts.append((pts[0] + pts[1] + pts[2])/3 + nrm) #p9
```

# 参考資料
**公式リファレンス**
Math Types & Utilities (mathutils)
https://docs.blender.org/api/current/mathutils.html

**フラクタル**
雑科学ノート　－ フラクタルの話 －
https://hr-inoue.net/zscience/topics/fractal/fractal.html
解説 -- フラクタルImaginary Cube --
https://www.i.h.kyoto-u.ac.jp/users/tsuiki/cfractal/kaisetu.html