フラクタル図形の拡大

解像度を明示せずにフラクタル図形を描画する

解像度を明示せずにフラクタル図形を無限に拡大するにはどうすれば良いか考えていた。例えば、再帰呼び出しを使って描くときに、再帰の回数を人間の方で指定せずに、拡大するとその都度座標が再計算されるようにできないか。さらに、自然な書式でフラクタル図形を定義する方法として、関数型言語の遅延評価のようなものを考えていた。

各ピクセルがその集合に載っていることが分かれば良いのでシェーダーを使えば良さそうだと思った。シェルピンスキーのギャスケットの場合、面積を持たないので、線分に幅を持たせて一番近い線分までの距離を求める必要がある。そこで次のようなフラグメントシェーダーを書いた。

precision mediump float;

uniform vec2 uOrigin;
uniform float uScale;

float cross2(vec2 a, vec2 b){
    return a.x * b.y - a.y * b.x;
}

float distLine(vec2 p, vec2 a, vec2 b){
    float x = distance(p, a);
    float y = distance(p, b);
    float z = distance(a, b);
    float s = (x + y + z) / 2.0;
    float area = sqrt(s * (s - x) * (s - y) * (s - z));
    return 2.0 * area / z;
}

float distTriangle(vec2 p, vec2 vert[3]){
    float dist = 1e9;
    dist = min(dist, distLine(p, vert[0], vert[1]));
    dist = min(dist, distLine(p, vert[1], vert[2]));
    dist = min(dist, distLine(p, vert[2], vert[0]));
    return dist;
}

bool inTriangle(vec2 p, vec2 v0, vec2 v1, vec2 v2){
    float sign0 = cross2(p - v0, v1 - v0);
    float sign1 = cross2(p - v1, v2 - v1);
    float sign2 = cross2(p - v2, v0 - v2);
    if(sign0 > 0.0 && sign1 > 0.0 && sign2 > 0.0) return true;
    if(sign0 < 0.0 && sign1 < 0.0 && sign2 < 0.0) return true;
}

bool inGasket(vec2 p, vec2 v0, vec2 v1, vec2 v2){
    float width = 1.0;
    vec2 mid[3];
    for(int i = 0; i < 100; i++){
        mid[0] = (v1 + v2) / 2.0;
        mid[1] = (v2 + v0) / 2.0;
        mid[2] = (v0 + v1) / 2.0;
        if(inTriangle(p, mid[0], mid[1], mid[2])) return distTriangle(p, mid) < width;
        if(inTriangle(p, v0, mid[1], mid[2])){
            v1 = mid[1];
            v2 = mid[2];
        }
        else if(inTriangle(p, mid[0], v1, mid[2])){
            v0 = mid[0];
            v2 = mid[2];
        }
        else if(inTriangle(p, mid[0], mid[1], v2)){
            v0 = mid[0];
            v1 = mid[1];
        }
    }
}

void main(){
    vec2 a = uOrigin + vec2(0.0, uScale);
    vec2 b = uOrigin + vec2(-uScale * sqrt(3.0) / 2.0, -uScale / 2.0);
    vec2 c = uOrigin + vec2(uScale * sqrt(3.0) / 2.0, -uScale / 2.0);
    gl_FragColor = (inGasket(gl_FragCoord.xy, a, b, c) ? vec4(0.0, 0.0, 0.0, 1.0) : vec4(1.0, 1.0, 1.0, 1.0));
}

glslで再帰呼び出しやwhile(true)が使えれば、有限時間で計算が終了し、適切にフラクタルを離散化できる。実際にはそのどちらも使えないようなので、結局forループで最も大きい三角形から分割する回数を指定することになった。この方法では拡大していくといつか限界が来る。

スクリーンと交わる部分をセグメント木のように求めていく方法も考えたが、同じ理由でできないと気付いた。

解像度を明示してフラクタル図形を描画する

解像度を明示する場合でも、最も大きい三角形からどれだけ拡大したかを記録しておかなければ、縮小しても戻ってくることはできない。拡大率が可変長のデータとなってしまうので、これも難しい。

一方で両方向に無限に拡縮できる場合は絶対的なサイズを持っておく必要がないので描画が可能。常に三角形がスクリーンを覆うように頂点を持っておく。

if(isIncluded(tri)){
    if(isIncluded([tri[0], inMid(tri[0], tri[1]), inMid(tri[0], tri[2])])){
        tri[1] = inMid(tri[0], tri[1]);
        tri[2] = inMid(tri[0], tri[2]);
        return;
    }
    if(isIncluded([inMid(tri[0], tri[1]), tri[1], inMid(tri[1], tri[2])])){
        tri[0] = inMid(tri[0], tri[1]);
        tri[2] = inMid(tri[1], tri[2]);
        return;
    }
    if(isIncluded([inMid(tri[0], tri[2]), inMid(tri[1], tri[2]), tri[2]])){
        tri[0] = inMid(tri[0], tri[2]);
        tri[1] = inMid(tri[1], tri[2]);
        return;
    }
}else{
    if(isIncluded([tri[0], exMid(tri[0], tri[1]), exMid(tri[0], tri[2])])){
        tri[1] = exMid(tri[0], tri[1]);
        tri[2] = exMid(tri[0], tri[2]);
        return;
    }
    if(isIncluded([exMid(tri[1], tri[0]), tri[1], exMid(tri[1], tri[2])])){
        tri[0] = exMid(tri[1], tri[0]);
        tri[2] = exMid(tri[1], tri[2]);
        return;
    }
    if(isIncluded([exMid(tri[2], tri[0]), exMid(tri[2], tri[1]), tri[2]])){
        tri[0] = exMid(tri[2], tri[0]);
        tri[1] = exMid(tri[2], tri[1]);
        return;
    }
}

function inMid(a, b){
    return { x: (a.x + b.x) / 2, y: (a.y + b.y) / 2 };
}

function exMid(a, b){
    return { x: 2 * b.x - a.x, y: 2 * b.y - a.y };
}

逆三角形に接している部分は自己相似となっていない（反対側に頂点が存在するかどうかが絶対的な基準となり、拡大率を覚えておく必要がある）ので、拡大し続けるには適宜位置をずらす。

自然界のフラクタル

良く考えると自然界にも理想的な形でフラクタルが存在しているわけではないから、コンピュータで計算できなくても仕方ないのかもしれない。視線に平行な直線を考えると、これらは全て視線の先の消失点で交わるが、光速は有限なのでいつまで経ってもレンダリングが終了しない。