前陣子看到 codeParade 的 marble marcher，這幾天又翻到 Optie 製作的 Aleph-0。想說來玩看看碎形。

Julia fractal

Julia 碎形由一個複數遞迴式產生

z_i + 1 = z_i² + c

其中 Julia set 也就是一般所見的碎形圖片中黑色區域為  ∣ z_n ∣  ≤ 2 的集合，外面彩色部分則根據遞迴次數來上色。 詳細資料可以看維基百科。

因為GLSL不支援複數型態，所以使用vec2來實作。將x部分作為實部；y部分作為虛部。所以公式中 z² 需要另外寫函數來實作。

vec2 complexMultiply(vec2 a, vec2 b){
    return vec2(a.x * b.x - a.y * b.y, a.y * b.x + a.x * b.y);  
}
vec2 complexSquare(vec2 z){
    return complexMultiply(z,z);    
}

這邊將螢幕座標 ([−0.5,0.5],[−0.5,0.5]) 作為z，c則是從主程式中傳入的時間來計算。

vec4 julia(void){
    float f = mod(float(iTime)/100.0, towPi) * sign(iTime);
    vec2 z = translate + texCoord / scale;
    vec2 c = vec2(sin(f)*(1+cos(f*2.371)),sin(f));
    int i;
    for(i = 0; i < iter; i++) {
        z = complexSquare(z) + c;
        if((z.x * z.x + z.y * z.y) > 4.0) 
            break;
    }
    vec3 color = texture(colorBar, float(i) / 80.0).rgb;
    return (i == iter)? vec4(0,0,0,1) : vec4(color ,1.0);
}

將圖片輸出後做成 gif 的結果

Mandelbrot fractal

Mandelbrot set 的公式跟 Julia 相同，不過此時的參數來源不一樣：

• z: (0, 0)
• c: 螢幕座標

vec4 mandelbrot(void){
    vec2 c = translate + texCoord / scale;
    vec2 z = vec2(0);
    int i;
    for(i = 0; i < iter; i++) {
        z = complexSquare(z) + c;
        if((z.x * z.x + z.y * z.y) > 4.0) 
            break;
    }
    vec3 color = texture(colorBar, float(i) / 60.0).rgb;
    return (i == iter)? vec4(0,0,0,1) : vec4(color ,1.0);
}

產生的結果

Colormap and color smoothing

接著用 colormap 圖片來讓顏色更好看，直接從 matplotlib 的 document 抓 jet 跟 hot 來試試看。

另外圖中有很明顯的 contour ，參考 wiki 來修正看看

Distance function

經過嘗試後發現上面的結果無法同時適用於 Julia set 跟 Mandelbrot set，查了之後才發現除了用遞迴次數計算顏色以外，有公式可以直接算出複數平面上各點與集合的最短距離。 參考了好幾個 references (1, 2, …) 後，使用 Inigo Quilez 的方法加上些許修正來計算顏色。

$$dist=\lim_{n\to\infty}\frac{|z_n|log(|z_n|)}{|z_n'|}$$

為了讓結構在放大後仍然清晰，不使用絕對的距離，而是使用

float zlen = length(z);
float dist = zlen / length(dz) * log(zlen);
float dist_scaled = log(min(dist, 0.5) / iScale) / log(Bailout);

上式中的 iScale 為放大倍率、 Bailout 為 Mandelbrot set 的計算閥值。因為這種方法不是用離散的遞迴次數去內插，所以整體的顏色很平滑。

跟之前一樣，輸出成 gif

自己玩得很開心，截了一堆圖～

小記

這次新碰到的東西

• 用 stb_image_write.h 存圖片
  • 要記得加上 #define STB_IMAGE_WRITE_IMPLEMENTATION 再 include
• 把 framebuffer 中的資料抓到主程式的 buffer 中： glReadPixels

Fragment Shader Code

#version 410 core

#define JULIA       0
#define MANDELBROT  1

uniform int         iMaxIter;
uniform float       iScale;
uniform sampler1D   iColorbar;
uniform int         iTime;
uniform vec2        iTranslate;
uniform int         iSelector;

const float Bailout = 100.0f;
const float Speed   = 90.0f;
const float Pi      = 3.1415926f;
const float TwoPi   = 6.2831852f;

in  vec2 _texCoord;
out vec4 _fragColor;

vec2 complexMultiply(vec2 a, vec2 b){
    return vec2(a.x * b.x - a.y * b.y, a.y * b.x + a.x * b.y);  
}
vec2 complexSquare(vec2 z){
    return complexMultiply(z, z);   
}

vec4 colorCalculation(float i, vec2 z, vec2 dz){
    if(i >= iMaxIter)
        return vec4(0);
    float zlen = length(z);
    float dist = zlen / length(dz) * log(zlen);
    float dist_scaled = log(min(dist, 0.5) / iScale) / log(Bailout);
    vec3 color = texture(iColorbar, dist_scaled).rgb;
    return vec4(color, 1.0);
}

vec4 mandelbrot(vec2 coord){
    float f = mod(float(iTime)/Speed, TwoPi) * sign(iTime);
    vec2 c = iTranslate + coord / iScale;
    vec2 z = vec2(0.0f);
    vec2 dz = vec2(0.0f);
    float i;
    for(i = 0; i < iMaxIter; i++) {
        dz = 2.0f * complexMultiply(z, dz) + vec2(1.0f, 0.0f);
        z = complexSquare(z) + c;
        if((z.x * z.x + z.y * z.y) > Bailout * Bailout) 
            break;
    }
    return colorCalculation(i, z, dz);
} 

vec4 julia(vec2 coord){
    float f = mod(float(iTime)/Speed, TwoPi) * sign(iTime);
    vec2 z = iTranslate + coord / iScale;
    vec2 c = vec2(sin(f) * (1 + cos(f)), 2 * cos(f));
    vec2 dz = vec2(1.0f, 0.0f);
    float i;
    for(i = 0; i < iMaxIter; i++) {
        dz = 2.0f * complexMultiply(z, dz);
        z = complexSquare(z) + sin(2.5*f) * c;
        if((z.x * z.x + z.y * z.y) > Bailout * Bailout) 
            break;
    }
    return colorCalculation(i, z, dz);
}

vec4 fractal(vec2 coord){
    switch(iSelector){
        case JULIA:
            return julia(coord);
                break;
        case MANDELBROT:
            return mandelbrot(coord);
                break;
        default:
            return vec4(0);
                break;
    }
}

void main() {
    _fragColor = fractal(_texCoord);
    return; 
}