Mandelbrot集合及其渲染

什么是Mandelbrot集合？

Mandelbrot集合是在復數平面上組成分形的點的集合，它正是以數學家Mandelbrot命名。

Mandelbrot集合可以用復二次多項式

\[f_c(z)=z^2+c \]

來定義
其中c是一個復數。對於每一個c，從\(z = 0\),開始對\(f_c(z)\)進行迭代。

序列\((0, f_ c(0), f_c(f_ c(0)), f_ c(f_ c(f_ c(0))), \ldots)\)的元素的模（復數具有模的概念）或者延伸到無窮大，或者只停留在有限半徑的圓盤內。Mandelbrot集合就是使以上序列不延伸至無限大的所有c點的集合。
從數學上來講，Mandelbrot集合是一個復數的集合。一個給定的復數c或者屬於Mandelbrot集合M，或者不屬於。比如，取c = 1，那么這個序列就是(0, 1, 2, 5, 26, ...)，顯然它的值會趨於無窮大；而如果取c = i，那么序列就是(0, i, -1+i, -i, -1+i, -i,...)，它的值會一直停留在有限半徑的圓盤內。

事實上，一個點屬於Mandelbrot集合當且僅當它對應的序列（由上面的二項式定義）中的任何元素的模都不大於2。這里的2就是上面提到的“有限半徑”。

繪制Mandelbrot集合

可以將屏幕上的一個像素映射為坐標系中的一點，如果該點屬於Mandelbrot集合，就將該像素着為黑色，這樣逐一對每個像素進行判斷和着色，就可以模擬繪制Mandelbrot集合了。

完成映射后來考慮如何判斷一個點是否屬於該集合。其根據就是上面的結論：一個點屬於Mandelbrot集合當且僅當它對應的序列（由上面的二項式定義）中的任何元素的模都不大於2，由於序列的的元素有無窮多個，我們只能取有限的迭代次數來模擬了，比如取100或1000次。

下面的代碼shader代碼完成了上面的思想。
其中迭代次數為200.
fragCoord.xy傳入的當前要計算顏色的像素點的坐標。iResolution.xy是顯示區域的寬高。iTime是時間變量。
代碼的前面一部分將像素點坐標轉換到了邏輯坐標系的坐標，並且隨着時間變化對顯示區域不斷縮放，使得我們可以動態的觀察不同大小尺度的Mandellbrot集合。
后面一部分計算了當前點是否位於Mandelbrot集合。

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 p = -1.0 + 2.0 * fragCoord.xy / iResolution.xy;//[0, 1]
    p.x *= iResolution.x/iResolution.y;//keep with h/w radio

    // animation	
	float tz = 0.5 - 0.5*cos(0.225*iTime);//[0,1]
    float zoo = pow( 0.5, 13.0*tz );//[1/2^13,1]
	vec2 c = vec2(-0.05,.6805) + p*zoo;

    // iterate
    vec2 z  = vec2(0.0);
    float m2 = 0.0;
    int isIn = 1;
    for( int i=0; i<200; i++ )
    {
        if( m2 > 2.0 ) { isIn = 0; break; }
			
        // Z -> Z² + c			
        z = vec2( z.x*z.x - z.y*z.y, 2.0*z.x*z.y ) + c;
			
        m2 = dot(z,z);
    }

    vec3 col = vec3(isIn==1?0:1 );
    
    fragColor = vec4( col, 1.0 );
}

可以進入這個地址觀看效果：https://www.shadertoy.com/view/3lsXWH

我們已可以從可視化的效果中觀察到數學的美妙，從簡單的遞歸定義中可以得的這種在大小尺度上不斷重復的圖案。
然而渲染效果不夠好，我們觀察到很多黑色的散點，和主體並不是連接的，而且在縮放的過程中，這些散點不斷忽然消失，又忽然出現。

改進渲染效果

前面一個假設“將屏幕上的一個像素映射為坐標系中的一點”其實是不對，屏幕中的像素是有限的，一個像素點對應的應該是坐標系中一塊區域，而不是一個點，然而計算中將像素點映射到了對應區域的中心點，如果這塊區域位於Mendelbrot集合分形圖形的邊緣分支上，它包含了無窮的分形的細節，區域中心點可能剛好位於Mandelbrot集合中也可能剛好位於集合外。所以在分形圖形的邊緣分支處出現了很多的散點。隨着縮放，在細節豐富地方的黑色點出來類似噪聲的效果，其實是由於前后兩幀像素點對應的區域的中心位置是否在集合中發生了變化。

如何改進渲染效果呢？應該讓像素反映出改點對應區域的密度，而不是僅僅取中心點計算，即如果該區域內有更多比例的點位於集合，使它表現出來的顏色更趨於黑色。
但是如何計算出對應區域的位於集合內的點的比例？想要得到這個精確的數值貌似是個不可能的任務，因為在我們的認知中分形圖形並不是一個確定性的圖形，它具有無窮的細節。但是我們可以采用某種方法估計這個比例。
假設我們可以得到某一點距離Mandelbrot集合圖形邊界的最短距離Dist,如果點位於集合內，Dist等於0，如果位於集合外，Dist是一個正值。Dist的計算方法在Inigo Quilez的文章distance rendering for fractals中給出了，在本文中就不具體闡述來了，我們把它當做一個黑盒，假設它的結果是正確的。
如果最短距離大於像素區域的大小，說明沒有像素區域內不包含集合中的點，該區域的密度為零；
如果最短距離小於像素區域的大小，可以說明以像素區域中心為中心以最小距離為半徑的圓內一定不包含集合點，我們假設剩余的面積都是被集合點充滿，然后我們可以根據面積的比例來估計集合點的比例。如下圖所示：

像素區域的寬高分別為w,h,中心點距離集合邊緣的距離為D，則區域內密度估計為\(Density = (w*h - Pi*D^2)/(w*h)\)。
這個估計值顯然會比真實值更大。
當圖像縮小時，像素對應的區域會變大，如圖中的大的紅色框所示，這時的密度估計會變大，在渲染中表現就是隨着視口放大，原先有很多細節的地方逐漸縮小，並且變得更黑。
下面是根據上面討論實現的效果：
觀看地址：https://www.shadertoy.com/view/wtlSW8、

從shaderToy的實時演示和截圖中可以看出顯示效果得到了很大提高，如果有更好的密度估計方法，顯示效率和顯示質量還能得到進一步提升。

void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    vec2 p = -1.0 + 2.0 * fragCoord.xy / iResolution.xy;
    p.x *= iResolution.x/iResolution.y;

    // animation	
	float tz = 0.5 - 0.5*cos(0.225*iTime);
    float zoo = pow( 0.5, 13.0*tz );//[1/2^12,1]
	vec2 c = vec2(-0.05,.6805) + p*zoo;

    // iterate
    bool isIn = true;
    vec2 z  = vec2(0.0);
    float m2 = 0.0;
    vec2 dz = vec2(0.0);
    for( int i=0; i<200; i++ )
    {
        if( m2>1024.0 ) { isIn=false; break; }

		// Z' -> 2·Z·Z' + 1
        dz = 2.0*vec2(z.x*dz.x-z.y*dz.y, z.x*dz.y + z.y*dz.x) + vec2(1.0,0.0);
			
        // Z -> Z² + c			
        z = vec2( z.x*z.x - z.y*z.y, 2.0*z.x*z.y ) + c;
			
        m2 = dot(z,z);
    }

    // distance	
	// d(c) = |Z|·log|Z|/|Z'|
    float d = 0.5*sqrt(dot(z,z)/dot(dz,dz))*log(dot(z,z));
    if (isIn)
        d = 0.0;
	
    // estimate density in the pixel area
    vec2 pixelScale = vec2(zoo/iResolution.y*iResolution.x/iResolution.y,zoo/iResolution.y);
    float pixelArea = pixelScale.x * pixelScale.y;
    
    float k = 1.0;
    float density = clamp((pixelArea - 3.1415926*d*d*k)/pixelArea, 0.0, 1.0);
    
    density = pow(density,3.0);             
    vec3 col = vec3( mix(1.0,0.0,density) );
    
    fragColor = vec4( col, 1.0 );
}

其中變量d就是當前像素點代表的區域中心點距離集合邊界的距離，其計算原理本文暫且不談。
vec2 pixelScale = vec2(zoo/iResolution.y*iResolution.x/iResolution.y,zoo/iResolution.y);
這行代碼根據當前的縮放比率以及顯示分辨率計算了像素代表的區域的大小。
float density = clamp((pixelArea - 3.1415926*d*d*k)/pixelArea, 0.0, 1.0);
這行代碼根據上文討論計算了估計密度。
然后剩余代碼根據估計密度計算出了像素的灰度值，這部分代碼可以更加感性的進行調整，只要能獲得更好的顯示效果。

着色

上面實現的效果是一張灰度圖，如果我們根據某種規則為像素賦予顏色，也許可以獲得更加具有美感的顯示效果。
todo