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