分形程序高级技巧入门教程--第一到四章

                          分形程序高级技巧入门教程--第一到四章
                                  housisong@gmail.com

摘要:
    本系列文章是写给分形编程爱好者的一个入门教程;文章章节包括(规划中的,可能增删):
    一.复数迭代的mandelbrot集合; 二.颜色平滑的简单周期算法; 三.迭代逃逸次数插值的颜色平滑;
    四.使用sin函数做颜色平滑; 五.一个更有效的迭代逃逸次数插值公式; 六.使用误差扩散来杜绝色差感;
    七.集合内部的颜色; 八.julia集合; 九.迭代生成的复数值的插值;
    十.迭代生成的复数值的高阶插值; 十一.图形的放大和旋转; 十二.复数初始值的变换;
    十三.固定颜色表; 十四.设计图案的映射; 十五.纹理映射; 十六.并行计算;
    十七.更多的迭代方程和颜色方案; 十八.高次mandelbrot集合; 十九.其他分形的介绍、分形动画

正文:

   (1-4章完整项目源代码: http://cid-10fa89dec380323f.office.live.com/self.aspx/.Public/hssFractal1%5E_4.zip )
   (我的其他分形相关文章: http://blog.csdn.net/housisong/category/382024.aspx )

一.复数迭代的mandelbrot集合

对于复数函数 f(z)=z^2+z0; 给定任意的z0，迭代可以得到序列：z0^2 + z0,
z0^4 + 2*z0^3 + z0^2 + z0,...... 我们把经过任意次迭代以后而不发散的z0定义为
Mandelbrot集合; 把复数方程改写为对应的实数方程:
    x(i+1)=x(i)*x(i)-y(i)*y(i)+x(0);  //实部
    y(i+1)=x(i)*y(i)*2+y(0);          //虚部

Mandelbrot迭代函数:我们认为迭代一定次数max_iter以后,没有逃逸的点就属于集合,返回max_iter; 如果发生逃逸,返回当前迭代次数,代码如下:  (当|Z(i)|^2>=4时,点必然逃逸)

inline long mandelbrot1(const double x0,const double y0,const long max_iter){ double x=x0; double y=y0; long i=0; for (;i<max_iter;++i){ if (x*x+y*y>=4) break; double tmp=x*x-y*y+x0; y=x*y*2+y0; x=tmp; } return i; }

 

迭代逃逸次数映射成颜色值:我们根据不同的逃逸次数,简单生成一个对应的RGB颜色值,代码如下:

//color=i%256 // / / / / / // / / / / / // / / / / / inline Color32 coloring1(const long iter,const long max_iter){ if (iter==max_iter){ return Color32(255,0,0); }else{ return Color32((iter*20)%256,(iter*15+85)%256,(iter*30+171)%256); } }  

我们把复平面上以点(x0,y0)为中心横轴2*r宽的复数区域均匀采样with(宽)*height(高)个点, 获取其Mandelbrot迭代逃逸次数,并映射成相应的颜色填充到一个with(宽)*height(高)的图片中,代码如下:
struct TViewRect{ double x0,y0,r; }; void draw_mandelbrot1(const TPixels32Ref& dst,const TViewRect& rect,const long max_iter){ for (long y=0;y<dst.height;++y){ for (long x=0;x<dst.width;++x) { double x0=(2*rect.r)*x/dst.width+rect.x0-rect.r; double yr=rect.r*dst.height/dst.width; //y轴半径 double y0=(2*yr)*y/dst.height+rect.y0-yr; long iter=mandelbrot1(x0,y0,max_iter); dst.pixels(x,y)=coloring1(iter,max_iter); } } }

调用该函数生成一幅图片,代码如下: (尺寸640x480,中心点(-0.5,0),r=2;最大迭代次数1000)

const long max_iter=1000; TViewRect rect; rect.x0=-0.5; rect.y0=0; rect.r=2; TPixels32 dstPic; dstPic.resizeFast(640,480); draw_mandelbrot1(dstPic.getRef(),rect,max_iter); { //save pic TFileOutputStream bmpOutStream(dstFileName); TBmpFile::save(dstPic.getRef(),&bmpOutStream);//保存结果图片 }  

生成的图片如下:


二.颜色平滑的简单周期算法
  直接的%求余容易造成较强的颜色阶梯,我们用一个平滑一点的函数来处理迭代值到颜色的转化: color=(i%510)-255;   当i连续变化的时候,得到的color值也是连续变化的: (其他代码相同)
//color=(i%510)-255 // // // / // / / / / / // / // // inline long modColor2(long iter){ return abs((iter+255)%510-255); } inline Color32 coloring2(const long iter,const long max_iter){ if (iter==max_iter){ return Color32(255,0,0); }else{ return Color32(modColor2(iter*20),modColor2(iter*15+85),modColor2(iter*30+171)); } }
生成的图片如下:


对于相邻的逃逸次数,现在的颜色还不够平滑,可以把那些*乘法系数取消掉,这样就能保证得到 的颜色的平滑(但颜色的对比就比较小了):
inline Color32 coloring3_s(const long iter,const long max_iter){ if (iter==max_iter){ return Color32(255,0,0); }else{ return Color32(modColor2(iter),modColor2(iter+85),modColor2(iter+171)); } }
生成的图片如下:


三.迭代逃逸次数插值的颜色平滑
   相邻点的逃逸次数差为1,对应到颜色最大亮度的1/255时才能看起来平滑(如:coloring1_s);
如果想增大他们之间的颜色对比,又有可能产生颜色梯度(如:coloring1); 我们如何同时满足
这两个要求呢?答案是对逃逸次数进行插值! 
  逃逸次数 插值公式: i +1-log(log(|Z|))/log(P); (P为复数方程的指数,这里为2)
  这是一个近似插值公式,有一定的误差;如果多迭代几次的话,可以减小该公式的误差;
  (公式来源参见: http://linas.org/art-gallery/escape/ray.html ) 代码如下: (增加了迭代次数)

static const double _divLog2=1.0/log(2.0); inline double log2(double x){ return log(x)*_divLog2; } inline double mandelbrot3(const double x0,const double y0,const long max_iter){ double x=x0; double y=y0; long i=0; for (;i<max_iter;++i){ if (x*x+y*y>=256) break; double tmp=x*x-y*y+x0; y=x*y*2+y0; x=tmp; } if (i!=max_iter){ return i +1-log2(log2(x*x+y*y)); }else return i; } //color=(i%510)-255 inline long modColor3(double iter){ return (long)(abs(fmod(iter+255,510)-255)); } inline Color32 coloring3(const double iter,const long max_iter){ if (iter==max_iter){ return Color32(255,0,0); }else{ return Color32(modColor3(iter*20),modColor3(iter*15+85),modColor3(iter*30+171)); } } void draw_mandelbrot3(const TPixels32Ref& dst,const TViewRect& rect,const long max_iter){ for (long y=0;y<dst.height;++y){ for (long x=0;x<dst.width;++x) { double x0=(2*rect.r)*x/dst.width+rect.x0-rect.r; double yr=rect.r*dst.height/dst.width; double y0=(2*yr)*y/dst.height+rect.y0-yr; double iter=mandelbrot3(x0,y0,max_iter); dst.pixels(x,y)=coloring3(iter,max_iter); } } }  

生成的图片如下:


四.使用sin函数做颜色平滑
  sin函数也是个不错的颜色平滑调节函数,函数值连续,处处可导... 产生的颜色更柔和一些,代码如下:

//color=sin(i) inline long sinColor(double iter){ return (long)( (sin(iter*2*3.1415926/510-3.1415926*0.5)+1)*0.5*255 ); } inline Color32 coloring4(const double iter,const long max_iter){ if (iter==max_iter){ return Color32(255,0,0); }else{ return Color32(sinColor(iter*20),sinColor(iter*15+85),sinColor(iter*30+171)); } } void draw_mandelbrot4(const TPixels32Ref& dst,const TViewRect& rect,const long max_iter){ for (long y=0;y<dst.height;++y){ for (long x=0;x<dst.width;++x) { double x0=(2*rect.r)*x/dst.width+rect.x0-rect.r; double yr=rect.r*dst.height/dst.width; double y0=(2*yr)*y/dst.height+rect.y0-yr; double iter=mandelbrot3(x0,y0,max_iter); dst.pixels(x,y)=coloring4(iter,max_iter); } } }  
生成的图片如下: