はじめに
前回の記事↓
前回の改良では、WebGL2を使用したGPUによる2Dレンダリングを実装して描画処理の軽量化を行いました。今回はWebGL2をさらに活かして、Gray-Scottモデルを3D風に表示する機能を追加してみました。その他、境界条件の切り替えなどの機能も追加しています。
コチラです。→Gray-Scottモデル v3
ソースコードはGitHubに置いています。↓
mitaka1962/gray-scott-model (GitHub)
実行例
実行例はこんな感じになります。
3Dと2Dは切り替え可能です。
他にも色んなチューリングパターンの形成が3Dで確認できます。
色を変えたり、境界条件を変更したりもできるのでぜひ試してみてください。
3D化の詳細
3D化の実装には下記のサイトを参考にしました。陰影処理によって3D風に見せるというのが今回の趣旨です。
WebGL2でGray-Scott反応拡散系シミュレーション (Qiita)
まず、平行光源（Directional Light）を追加して画面に陰影を与えます。
これまでの2Dでの描画では、u(あるいはv)の濃度が1で白、0で黒としてグレースケールで表示していました。3D化する際にはこのグレースケール画像をハイトマップとして使います。ハイトマップからx方向、y方向における傾きを求めて接ベクトルとし、それらの外積をとることでその点における単位法線ベクトルを求めます。下記の記事を参考にしました。
フラグメントシェーダでの実装はこのようになっています。
vec2 onePixel = 1.0 / uTextureSize; // 1ピクセルあたりの長さを求める
// 現在の座標の上下左右のx(y)座標を取得
float dec_x = vTexCoord.x - onePixel.x;
float inc_x = vTexCoord.x + onePixel.x;
float dec_y = vTexCoord.y - onePixel.y;
float inc_y = vTexCoord.y + onePixel.y;
// 周期境界とする
float p_dec_x = (dec_x < 0.0) ? dec_x + 1.0 : dec_x;
float p_inc_x = (inc_x > 1.0) ? inc_x - 1.0 : inc_x;
float p_dec_y = (dec_y < 0.0) ? dec_y + 1.0 : dec_y;
float p_inc_y = (inc_y > 1.0) ? inc_y - 1.0 : inc_y;
// テクスチャから値を取得（Rにuの濃度, Gにvの濃度）
vec2 uv = texture(uDrawTex, vTexCoord).rg;
vec2 uv1 = texture(uDrawTex, vec2(p_dec_x, vTexCoord.y)).rg;
vec2 uv2 = texture(uDrawTex, vec2(p_inc_x, vTexCoord.y)).rg;
vec2 uv3 = texture(uDrawTex, vec2(vTexCoord.x, p_dec_y)).rg;
vec2 uv4 = texture(uDrawTex, vec2(vTexCoord.x, p_inc_y)).rg;
// 境界条件によって条件分岐
if (uBoundaryCondition == 1) {
// dirichlet boundary condition
uv1 = (dec_x < 0.0) ? vec2(2.0 - uv.r, -uv.g) : uv1;
uv2 = (inc_x > 1.0) ? vec2(2.0 - uv.r, -uv.g) : uv2;
uv3 = (dec_y < 0.0) ? vec2(2.0 - uv.r, -uv.g) : uv3;
uv4 = (inc_y > 1.0) ? vec2(2.0 - uv.r, -uv.g) : uv4;
} else if (uBoundaryCondition == 2) {
// neumann boundary condition
uv1 = (dec_x < 0.0) ? uv : uv1;
uv2 = (inc_x > 1.0) ? uv : uv2;
uv3 = (dec_y < 0.0) ? uv : uv3;
uv4 = (inc_y > 1.0) ? uv : uv4;
}
// 描画対象（u or v）の値を取得
float center = (uTarget == 1) ? uv.g : uv.r;
float left = (uTarget == 1) ? uv1.g : uv1.r;
float right = (uTarget == 1) ? uv2.g : uv2.r;
float down = (uTarget == 1) ? uv3.g : uv3.r;
float up = (uTarget == 1) ? uv4.g : uv4.r;
// 単位法線ベクトルを求める
vec3 dx = vec3(1.0, 0.0, (right - left) / (2.0 * 0.07));
vec3 dy = vec3(0.0, 1.0, (up - down) / (2.0 * 0.07));
vec3 normal = normalize(cross(dx, dy));
境界条件による分岐処理は後程説明します。
x, yについての偏微分を求める際の分母は2.0 * 0.07となっていますが、これには特に深い意味はありません。中心差分を用いて数値微分をしているので分母は2Δxのなるのですが、ここでΔxを0.07にしたらいい感じの陰影が作れたのでそうしました。
単位法線ベクトルを求める際はcrossやnormalizeなどの便利な組み込み関数を使用しています。
その後、平行光源の方向ベクトル（実際は逆向き）と単位法線ベクトルの内積をとってこれらのベクトルがなす角のcosを求めます。なす角が0だと光が真上から当たるので明るいですが、角度が開いていくと斜めから光が当たることになるので暗くなります。光の強さにcosを掛けることでこの性質を再現します。この辺の説明は下記のサイトがわかりやすいです。
WebGL 3D - Directional Lighting
コードは以下のようになります。
// 光源の方向ベクトル（逆向き）
vec3 light1 = normalize(vec3(1.0, 1.0, 1.0));
// cosが負になったときに加わる光の強さを0とするためmax関数を使用
float diffuse = max(1.1 * decay(center) * dot(normal, light1), 0.0);
// 異なる光源をもう一つ追加
vec3 light2 = normalize(vec3(-0.8, -1.2, 0.3));
diffuse += max(0.5 * decay(center) * dot(normal, light2), 0.0);
// 0から1の間に収める
diffuse = clamp(diffuse, 0.0, 1.0);
vec3 color;
// 色による分岐
// グレースケール以外は影色と表面の色をmix関数で混合
if (uColor == 0) {
// sky
color = mix(vec3(0.08, 0.35, 0.5), vec3(0.75, 0.88, 0.96), diffuse);
} else if (uColor == 1) {
// milky
color = mix(vec3(0.6, 0.4, 0.25), vec3(1.0, 0.94, 0.89), diffuse);
} else if (uColor == 2) {
// poison
color = mix(vec3(0.1, 0.0, 0.3), vec3(0.75, 0.6, 1.0), diffuse);
} else {
// greyscale
color = vec3(1.0) * diffuse;
}
ここでdecay()という関数を独自に追加しています。中身はこんな感じです。
float decay(float z) {
return 1.0 / pow(1.0 + (1.0 - z) / 2.0, 2.0);
}
u(あるいはv)の値が小さいほど、すなわち高さが低いほど、光源から遠ざかるので光が減衰すると考えて独自に追加してみました。逆2乗則を満たし、入力が1のときに1を返す関数として上のような式の形になりました。物理的に正しいかはわかりませんが、先ほどと同様、陰影がいい感じになるのでOKとします。
ここまで行って来たのは拡散反射（ランバート反射と言うそう）の再現だったので、次に鏡面反射（というより単なるハイライトかも）を追加します。
まず、先ほどの平行光源の逆方向ベクトルであるlight1と単位法線ベクトルnormalから反射光の方向ベクトルを求め、これと視点に向かうベクトルの内積をとります（今回は真上から見ているので(0, 0, 1)）。内積はこれらのベクトルがなす角のcosとなるので値が-1から1になりますが、ここで1周辺以外を0に近づけるために値を何乗かします。この方法は下記のサイトを参考にしました。
【連載】Unity時代の3D入門 – 第6回「鏡面反射ライティング」
コードは以下の通りです。
// light1を逆向きにして入射光の方向ベクトルにする
vec3 reflect = reflect(-light1, normal);
vec3 eye = vec3(0.0, 0.0, 1.0);
color += pow(max(0.0, dot(eye, reflect)), 40.0);
反射光の方向ベクトルはreflect()という組み込み関数で求めています（便利だ。。。）。内積を40乗しているのはこの位の値でちょうどいいハイライトになったからです。
最後に値を0から1の範囲に収めて終了です。
outColor = vec4(clamp(color, 0.0, 1.0), 1.0);
ちなみに2D描画の時の処理はたったこれだけです。
vec2 uv = texture(uDrawTex, vTexCoord).rg;
color = (uTarget == 1) ? vec3(uv.g) : vec3(uv.r);
outColor = vec4(clamp(color, 0.0, 1.0), 1.0);
境界条件の詳細
これまでは周期境界条件を用いていましたが、新たに他の境界条件としてディリクレ境界条件とノイマン境界条件を選択できるようにしました。ディリクレ境界条件を選択した場合は境界上の濃度値(u, v)が(1, 0)で一定になり、ノイマン境界条件を選択した場合は境界上の濃度勾配が0になります。
実装は先ほど出てきた通りです。
if (uBoundaryCondition == 1) {
// dirichlet boundary condition
uv1 = (dec_x < 0.0) ? vec2(2.0 - uv.r, -uv.g) : uv1;
uv2 = (inc_x > 1.0) ? vec2(2.0 - uv.r, -uv.g) : uv2;
uv3 = (dec_y < 0.0) ? vec2(2.0 - uv.r, -uv.g) : uv3;
uv4 = (inc_y > 1.0) ? vec2(2.0 - uv.r, -uv.g) : uv4;
} else if (uBoundaryCondition == 2) {
// neumann boundary condition
uv1 = (dec_x < 0.0) ? uv : uv1;
uv2 = (inc_x > 1.0) ? uv : uv2;
uv3 = (dec_y < 0.0) ? uv : uv3;
uv4 = (inc_y > 1.0) ? uv : uv4;
}
デフォルトを周期境界として、他の境界条件が選択されているときははみ出した座標における濃度値を上書きします。ディリクレ境界条件では、現在の座標での濃度値uvとはみ出した隣の座標での濃度値との平均がvec2(1.0, 0.0)で一定になるようにします。これによって境界線上の濃度値が(1, 0)で一定になります。ノイマン境界条件の時は、はみ出した座標における濃度値を現在の座標での濃度値uvと同じにすることで傾きを0にします。実装には以下のサイトを参考にしました。
とはいえ、境界条件の変更はパターン形成自体にはほとんど影響を与えないようなので、変更してもそれほど面白い変化は起こらないです。
おわりに
内積や外積などの比較的簡単な数学的知識を使うだけで、ここまで「リアルっぽい」描画ができるのは驚きでした。また、3D化することでパターン形成の合間に波打っている様子なんかが現実の液体（流体？）のような振る舞いをすることがよくわかりました。
