こんにちは、NUMです。
何だか肌寒い日が続くようになって、やっと夏が終わりそうな感じがしますね。
僕の最近のトピックは、剱岳という日本で一番難しいとされている山に登ったことです。
登山というより崖を登っている感覚に近くて、命の危機を感じるような場面がいくつかありました、、、
頂上まで登ってみると今まで見たこともないような景色でした。
9月中旬でも剱岳付近は既に植物が紅葉が始まっていて、頂上まで登らずとも景色が最高です。自然は一番の芸術だとつくづく思いました。
今回は山と植物に関連して、数理生物学の分野でもある縞枯れ現象という森林の更新過程を表したモデルで遊んでいこうと思います。
縞枯れ現象とは
縞枯れ現象とは亜高山帯の針葉樹であるシラビソ、オオシラビソなどの優先林で見られる森林更新現象です。
（優先林：洪水時に流木発生源とならない木が優先して植えられている森林なので、同じ種類の木が植えられています。）
縞枯れと言われる理由は、木々が立ち枯れたり倒れたりして、生きている木と枯れている木が交互になっている状態が縞模様に見えるからです。
この現象は恒常風が主な要因と言われています。（恒常風：年間を通して決まった方向に吹く風）
縞枯れ現象が発生するのは悪い事ではなく、木々の世代交代を促す良い役割があると言われています。
モデル化
それでは縞枯れ現象の再現方法を考えていきます。
縞枯れ現象は恒常風によって発生するものと仮定し、2次元のセルオートマトンを使って再現していきます。今回はセルを樹木と見なし、近傍の木々の状態によって、樹木が枯れるか成長するかを判定していくのでセルオートマトンを使っています。
判定対象の樹木が成長するかどうかは、風上の樹木（風上サイト）と判定対象の樹木（サイト）の樹高差が関係していると考えられています。
風上サイトよりサイトの方が樹高が高い場合、恒常風を受けやすくなり、樹木が枯れる可能性があるからです。
以下イメージ図です。
サイト樹高：h(x,y)
風上サイト樹高：h(x-1,y-1)（恒常風が左上から吹くと仮定）
近隣の樹高差の制限値：d
判定条件：
h(x,y)とh(x-1,y-1)の樹高差がdより大きい場合、h(x,y)は枯れる。
それ以外の場合、h(x,y)は成長する。
ただし、実際にはh(x-1,y)とh(x,y-1)も風上サイトであると考えられます。
そのため、h(x-1,y-1)を最も影響度の高い風上サイト、h(x-1,y)とh(x,y-1)を影響度の低い樹木と考えます。
判定では樹高の値を使用するため、3つの樹高の平均を影響度定数を使用して計算します。
Aを影響度の定数（Aの値は0~1、1の場合h(x-1,y-1)と同じ影響度）とすると
以下数式となります。
ah(x,y)= h(x-1,y-1) + A*(h(x-1,y)+h(x,y-1))/1+2*A
今までの説明を擬似プログラムで表現してみます。（補足：t+1というのは次の世代のサイトの樹高、0は枯れる、1は1世代あたりの成長値）
h(x,y)(t+1) = 0：if(h(x,y)-ah(x,y) > d)
h(x,y)(t+1) = h(x,y)(t) + 1：else
それでは待ちに待ったProcessingで実装していきましょう！
実装
/*
縞枯モデル
立ち枯れ帯付近の樹高が高い樹は卓越風にされされることで枯れる
風上の樹が枯れた場合、付近の樹が風にさらされて枯れる
風にさらされない樹は一定の速度で成長していく
①風上の樹の樹高
付近の樹が①より高い：枯れる
付近の樹が①より低い：成長する
*/
class WaveRegeneration {
float[][] trees; // 現世代樹高
float[][] nextTrees; // 次世代樹高
int [][] windAngles; // 樹木群と風源の角度
int windWardX; // 風源x座標
int windWardY; // 風源y座標
int scl; // 拡大率
int len; // 配列要素-1
float limit = 20; // 樹高の閾値
float grow = 5; // 1世代あたりの成長値
float alfa = random(0, 1); // 風上サイトの影響度
//コンストラクタ
WaveRegeneration(int treesSize, int x, int y) {
this.windWardX = x;
this.windWardY = y;
this.scl = width/treesSize;
this.len = treesSize-1;
setTrees(treesSize);
}
/**
* 樹木群二次元配列（現世代と次世代）と各樹木と風源の角度の二次元配列を設定
* @param1 int 配列の要素数
*/
void setTrees(int treesSize) {
this.trees = new float[treesSize][treesSize];
this.nextTrees = new float[treesSize][treesSize];
this.windAngles = new int[treesSize][treesSize];
for (int col = 0; col<treesSize; col++) {
for (int row = 0; row<treesSize; row++) {
this.trees[col][row] = int(random(20, 255)); // 樹高
this.windAngles[col][row] = calcWindDirectionAngle(col, row);
}
}
}
/**
* 指定の値まで世代交代をさせる
* @param1 int 世代交代数
*/
void passageTime(int time) {
for (int t = 0; t<time; t++) {
updateTrees();
swap();
}
}
/**
* 樹木群の状態の更新
*/
void updateTrees() {
for (int col = 1; col<this.len; col++) {
for (int row = 1; row < this.len; row++) {
float siteH = this.trees[col][row]; // サイトの樹高
int windAngle = this.windAngles[col][row]; // 風源からの角度
float windWardSiteAvgH = getAvgWindWardTreeHeight(col, row, windAngle); //3つの風上サイトから樹高を計算して取得
// サイトが枯れるか成長するかを判定
if (siteH - windWardSiteAvgH > this.limit) {
this.nextTrees[col][row] = 0; // 枯れる
} else {
this.nextTrees[col][row] = siteH + this.grow; // サイトの現在の樹高に成長値を加算
}
}
}
}
/**
* 可視化
*/
void drawSites() {
for (int col = 1; col<len; col++) {
for (int row = 1; row < len; row++) {
float mapTree = map(this.trees[col][row], 0, limit, 0, 1);
color c = lerpColor(color(0), color(255, 255, 255), mapTree);
fill(c);
rect(col*this.scl, row*this.scl, this.scl, this.scl);
}
}
}
/**
* 風上サイトの樹高を計算する
* 風は南東から吹いていると仮定すると、風向に垂直なtrees[col+1][row+1]が最も風避けの効果が大きい
* trees[col[row+1]とtrees[col+1[row]は風避けの効果が小さいため、alfaで影響度を調整している
* @param1 int 参照列
* @param2 int 参照行
* @param3 int 風向の角度
*/
float getAvgWindWardTreeHeight(int col, int row, int windAngle) {
float eastTree = this.trees[col][row+1];
float southEastTree = this.trees[col+1][row+1];
float southTree = this.trees[col+1][row];
float southWestTree = this.trees[col+1][row-1];
float westTree = this.trees[col][row-1];
float northWestTree = this.trees[col-1][row-1];
float northTree = this.trees[col-1][row];
float northEastTree = this.trees[col-1][row+1];
float mostProtectionTree = 0; // 最も風避けの影響が高い樹
float lowProtectionTree1 = 0; // 風避けの影響が低い樹①
float lowProtectionTree2 = 0; // 風避けの影響が低い樹②
if (windAngle >= 0 & windAngle <90) {
//風向：東、南東
mostProtectionTree = southEastTree;
lowProtectionTree1 = eastTree;
lowProtectionTree2 = southTree;
} else if (windAngle >= 90 & windAngle <180) {
//風向：南、南西
mostProtectionTree = southWestTree;
lowProtectionTree1 = southTree;
lowProtectionTree2 = westTree;
} else if (windAngle >= 180 & windAngle <270) {
//風向：西、北西
mostProtectionTree = northWestTree;
lowProtectionTree1 = westTree;
lowProtectionTree2 = northTree;
} else if (windAngle >= 270 & windAngle <360) {
//風向：北、北東
mostProtectionTree = northEastTree;
lowProtectionTree1 = northTree;
lowProtectionTree2 = eastTree;
}
return (mostProtectionTree + alfa*(lowProtectionTree1 + lowProtectionTree2))/(1+2*alfa);
}
/**
* ①現世代配列
* ②次世代配列（計算済）
* ②を①に代入して次の数値計算に使用する
* 次世代配列は計算結果が代入されて更新されるためtempを設定しても問題なし
*/
void swap() {
float[][] temp = this.trees;
this.trees = this.nextTrees;
this.nextTrees = temp;
}
/**
* 対象座標から風源の角度を取得
* @param1 int 配列のx座標
* @param2 int 配列のy座標
*/
int calcWindDirectionAngle(int col, int row) {
int windAngle = int(degrees(atan2((row*this.scl)-(this.windWardY), (col*this.scl)-(this.windWardX))));
//181°以降は-180~0を180~360に変換する必要があるため
//lightAngle変数に360°を加算して調整している
if (windAngle<0) {
windAngle = 360+windAngle;
}
return windAngle;
}
}
上記が縞枯れモデルのプログラムです。
基本的にプログラムのコメントに概要は書いているので、重要なポイントだけを説明していきます。
WaveRegenerationクラス
縞枯れ現象は2次元配列のセルオートマトンの仕組みを使って再現しています。
WaveRegenerationクラスのメンバtreesとnextTreesが操作対象の2次元配列です。
二つの2次元配列には樹高の値が格納されています。
コンストラクタ
WaveRegeneration(配列の要素数、風源のx座標、風源のy座標)を指定します。
setTreesメソッド
2つの現世代、次世代の2次元配列を指定の要素数で作成します。
また、windAnglesの2次元配列には要素の位置から見た風源との角度を計算して設定します。
これは風源の位置は1度設定したら変更されないので、事前に計算しておくことで処理の効率化をしています。
（ただし、風源の位置座標を動的させる場合、この配列は使えないので注意です。）
updateTreesメソッド
2次元配列処理のためにfor文をネストしています。
for文のカウンタ変数の初期値が1で終了条件の値も要素数-1としているのには理由があります。
風上サイトの樹高を計算する際にサイトの近傍にアクセスします。
2次元配列の最初と最後の列と行から処理を行うとどうなるでしょうか？
配列外の要素(赤塗り)にアクセスしてしまうためエラーになります。
そのため、青塗りの部分を処理していく必要があります。
getAvgWindWardTreeHeightメソッド
各セルからの風源の角度を指定して、風上サイトの近傍がどの範囲になるかを選択し、最終的に風上サイトの樹高を計算しています。
角度の計算方法は以前の[球に影をつける](https://www.num-kt.com/post/shadowboll)記事で書いているので気になった方は参照して下さい。
swapメソッド
次世代の樹高はtreesの樹高を参照、計算してnextTreesに入れていました。
図の3世代目を計算するには2世代目の値を使用したいですが、nextTreesに入っていて計算出来ません。
そのため、treesにnextTreesを代入してあげると、この問題が解決できます。
複数世代計算が必要な場合でも、2つの配列さえ作成すれば良いので効率がいいですよね。
あとはお決まりのsetup関数にWaveRegenerationクラスのインスタンスを用意してあげて
passageTimeメソッドで世代交代数を指定、drawSitesメソッドで可視化すればOKです。
WaveRegenerationクラスのメンバalfaで影響度の値をランダムにしているので、描画するたびに見え方が違って面白いです。
是非好みの模様を生成してジェネラティブアートの面白さを体感してみて下さい！
まとめ
今回は縞枯れモデルという自然現象のモデル化と実装に挑戦してみました。
僕は元々、大学では生物を学んでいましたが、生物を数学やシミュレーションという観点で学んだことはありませんでした。
正直こっちの方がメチャクチャ面白いです。（大学で数理生物学を専攻すればよかった、、、、、泣）
まあ、趣味程度でぼちぼちやっていこうかと思います笑
縞枯れモデルを使用した作品の紹介になります。
この作品は剱岳に行った時の思い出を表現しました。
頂上までの苦しい道のり、登り切った後の達成感の二つを表しています。
縞枯れモデルは山のイメージを表現するのにとても使えると思いました。
以上が縞枯れモデルの実装になります。
最後まで読んでいただきありがとうございました！
参考文献
コメント