分形程序高级技巧入门教程--第五到八章

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

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

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


五.一个更有效的迭代逃逸次数插值公式
     我自己推导了一下迭代逃逸次数插值公式,根据逃逸阈值M,当前迭代次数i,逃逸值Zn和上一个没有逃逸的值Zo的关系,
   得到了一个插值公式: i-(f(|Zn|)-f(|M|))/(f(|Zn|)-f(|Zo|)) ; (该公式的适应性很强,可以用在很多不同的地方);
   这里,f(x)函数取log(abs(log(x))),并且当M值较大时(多迭代几次),就可以得到接近线性的插值效果;
   可以证明,该公式在M足够大的情况下可以约简为第三章介绍的插值公式(特例);当然f(x)函数也可以取其他很多函数,
   也可以得到平滑的插值结果,但可能插值的结果看起来不够线性(也许你就想要这个效果);
   新的插值公式函数代码如下:
inline double loglog(const double x){ return log(abs(log(x))); } inline double mandelbrot5(const double x0,const double y0,const long max_iter){ const static double M=256; const static double lnln_M=loglog(M); double x_bck=0; double y_bck=0; double x=x0; double y=y0; long i=0; for (;i<max_iter;++i){ if (x*x+y*y>=M) break; x_bck=x; y_bck=y; double tmp=x*x-y*y+x0; y=x*y*2+y0; x=tmp; } if (i!=max_iter){ const double lnln_Z=loglog(x*x+y*y); const double lnln_Zbak=loglog(x_bck*x_bck+y_bck*y_bck); return i-2-(lnln_Z-lnln_M)/(lnln_Z-lnln_Zbak); }else return i; } void draw_mandelbrot5(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=mandelbrot5(x0,y0,max_iter); dst.pixels(x,y)=coloring4(iter,max_iter); } } }
   生成的图片如下:


   新的插值效果平滑更好一些,但这种情况下可能不容易看出来,当M值较小的时候(M=4),可以对比一下两者的差异:
代码:

inline double mandelbrot4_e(const double x0,const double y0,const long max_iter){ const static double M=4; double x=x0; double y=y0; long i=0; for (;i<max_iter;++i){ if (x*x+y*y>=M) 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; } void draw_mandelbrot4_e(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=mandelbrot4_e(x0,y0,max_iter); dst.pixels(x,y)=coloring4(iter,max_iter); } } } inline double mandelbrot5_e(const double x0,const double y0,const long max_iter){ const static double M=4; const static double lnln_M=loglog(M); double x_bck=0; double y_bck=0; double x=x0; double y=y0; long i=0; for (;i<max_iter;++i){ if (x*x+y*y>=M) break; x_bck=x; y_bck=y; double tmp=x*x-y*y+x0; y=x*y*2+y0; x=tmp; } if (i!=max_iter){ const double lnln_Z=loglog(x*x+y*y); const double lnln_Zbak=loglog(x_bck*x_bck+y_bck*y_bck); return i-(lnln_Z-lnln_M)/(lnln_Z-lnln_Zbak); }else return i; } void draw_mandelbrot5_e(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=mandelbrot5_e(x0,y0,max_iter); dst.pixels(x,y)=coloring4(iter,max_iter); } } }

M=4,用第三章的插值公式,误差能够直接看出来了:
  

M=4,新的插值公式(还保持着平滑,但不够线性):

 
六.使用误差扩散来杜绝色差感
     可以看到很多分形图像的颜色过渡非常的流畅亮丽,这是很多分形图像引人入胜的重要因素;
   我们的显示器每个颜色分量有256级亮度,大部分情况下眼睛已经很难看出颜色梯度; 但在某些
   特殊情况下,生成的图像还是可能会产生肉眼可见的色差;而程序内部实际上可以生成更精确的
   颜色,只是在转化到整数颜色时丢失了,这时我们可以利用误差扩散来处理转化过程;
   (误差扩散的原理可以参看我的相关3篇文章: <图形图像处理－之－误差扩散 上 中 下 >
     http://blog.csdn.net/housisong/archive/2008/04/23/2316924.aspx )

   为了简单,我用的最简单的*1扩散模板,代码如下:
   //浮点RGB颜色 struct Colorf{ double R; double G; double B; inline Colorf(){} inline Colorf(const Colorf& c):R(c.R),G(c.G),B(c.B){} inline Colorf(const Color32& c):R(c.r*(1.0/255)),G(c.g*(1.0/255)),B(c.b*(1.0/255)){} inline Colorf(const double _R,const double _G,const double _B):R(_R),G(_G),B(_B){} inline void addColor(const Colorf& c){ R+=c.R; G+=c.G; B+=c.B; } inline void subColor(const Colorf& c){ R-=c.R; G-=c.G; B-=c.B; } inline void mulColor(const double p){ R*=p; G*=p; B*=p; } inline static long fToIColor(const double fc){ if (fc<=0) return 0; else if (fc>=1) return 255; else return (long)(fc*255+0.49999); } inline Color32 toColor32()const{ return Color32(fToIColor(R),fToIColor(G),fToIColor(B)); } }; inline double sinColorf(double iter){ return (sin(iter*2*3.1415926/510-3.1415926*0.5)+1)*0.5; } inline Color32 coloring6(const double iter,const long max_iter,const Colorf& errorColorIn,Colorf& errorColorOut){ Colorf color=errorColorIn; if (iter==max_iter){ color.addColor(Color32(255,0,0)); }else{ color.addColor(Colorf(sinColorf(iter*20),sinColorf(iter*15+85),sinColorf(iter*30+171))); } Color32 resultColor=color.toColor32(); errorColorOut=color; errorColorOut.subColor(resultColor); return resultColor; } void draw_mandelbrot6(const TPixels32Ref& dst,const TViewRect& rect,const long max_iter){ for (long y=0;y<dst.height;++y){ Colorf errorColor(0,0,0); 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=mandelbrot5(x0,y0,max_iter); dst.pixels(x,y)=coloring6(iter,max_iter,errorColor,errorColor); } } }
   生成的图片如下:
   (很多人眼睛可能看不出差异,但这一幅的颜色过度更好,我也很难随手构造出一幅特殊情况的图像让大家能直接看出差异)

七.集合内部的颜色
     前面我们给集合外的区域利用迭代逃逸次数赋予了比较丰富的颜色,集合内部只是简单的设置为红色;其实集合内部也是可以发挥想象
力赋予漂亮的颜色的,比如根据迭代最后的复数值来生成颜色: (任何迭代过程中生成的数据都可以作为颜色函数的参数)
inline double mandelbrot7(double* xList,double* yList,const long max_iter){ const static double M=256; const static double lnln_M=loglog(M); const double& x0=xList[0]; const double& y0=yList[0]; double x_bck=0; double y_bck=0; double x=x0; double y=y0; long i=0; for (;i<max_iter;++i){ if (x*x+y*y>=M) break; x_bck=x; y_bck=y; double tmp=x*x-y*y+x0; y=x*y*2+y0; x=tmp; xList[i+1]=x; yList[i+1]=y; } if (i!=max_iter){ const double lnln_Z=loglog(x*x+y*y); const double lnln_Zbak=loglog(x_bck*x_bck+y_bck*y_bck); return i-2-(lnln_Z-lnln_M)/(lnln_Z-lnln_Zbak); }else return i; } inline Color32 coloring7(const double iter,const long max_iter,double* xList,double* yList,const Colorf& errorColorIn,Colorf& errorColorOut,double k=1){ Colorf color=errorColorIn; if (iter==max_iter){ const double x=xList[max_iter]; const double y=yList[max_iter]; double z=sqrt(x*x+y*y); double zd=z-sqrt(xList[max_iter-1]*xList[max_iter-1]+yList[max_iter-1]*yList[max_iter-1]); color.addColor(Colorf(sinColorf(z*2000*k),sinColorf(y*x*1000*k),sinColorf(zd*1000*k))); }else{ color.addColor(Colorf(sinColorf(iter*20*k),sinColorf(iter*15*k+85),sinColorf(iter*30*k+171))); } Color32 resultColor=color.toColor32(); errorColorOut=color; errorColorOut.subColor(resultColor); return resultColor; } void draw_mandelbrot7(const TPixels32Ref& dst,const TViewRect& rect,const long max_iter){ std::vector<double> xList; std::vector<double> yList; xList.resize(max_iter+1); yList.resize(max_iter+1); for (long y=0;y<dst.height;++y){ Colorf errorColor(0,0,0); 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; xList[0]=x0; yList[0]=y0; double iter=mandelbrot7(&xList[0],&yList[0],max_iter); dst.pixels(x,y)=coloring7(iter,max_iter,&xList[0],&yList[0],errorColor,errorColor); } } }

生成的图片如下:


八.julia集合
  复数函数 f(z)=z^2+zM, 其中zM(jx0,jy0)为常量，对平面上的一点z0进行迭代，经足够多次迭代
后函数值不扩散，这类z0点组成的集合为Julia集，对每一个特定的zM值都构成一个相应的Julia集；
　可以证明:mandelbrot集是使Julia集为连通的参数(zM)的集合。
  把复数方程改写为对应的实数方程:
    x(i+1)=x(i)*x(i)-y(i)*y(i)+jx0;  //实部
    y(i+1)=x(i)*y(i)*2+jy0;          //虚部
  julia迭代函数:  (和 mandelbrot 迭代函数对比一下 )
long julia(const double x0,const double y0,const long max_iter,double jx0,double jy0){ 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+jx0; y=x*y*2+jy0; x=tmp; } return i; }
  融合前几章的内容,实际代码如下:
inline double julia8(double* xList,double* yList,const long max_iter,double jx0,double jy0){ const static double M=256; const static double lnln_M=loglog(M); const double& x0=xList[0]; const double& y0=yList[0]; double x_bck=0; double y_bck=0; double x=x0; double y=y0; long i=0; for (;i<max_iter;++i){ if (x*x+y*y>=M) break; x_bck=x; y_bck=y; double tmp=x*x-y*y+jx0; y=x*y*2+jy0; x=tmp; xList[i+1]=x; yList[i+1]=y; } if (i!=max_iter){ const double lnln_Z=loglog(x*x+y*y); const double lnln_Zbak=loglog(x_bck*x_bck+y_bck*y_bck); return i-2-(lnln_Z-lnln_M)/(lnln_Z-lnln_Zbak); }else return i; } void draw_julia8(const TPixels32Ref& dst,const TViewRect& rect,const long max_iter,double jx0,double jy0){ std::vector<double> xList; std::vector<double> yList; xList.resize(max_iter+1); yList.resize(max_iter+1); for (long y=0;y<dst.height;++y){ Colorf errorColor(0,0,0); 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; xList[0]=x0; yList[0]=y0; //double iter=julia(xList[0],yList[0],max_iter,jx0,jy0); double iter=julia8(&xList[0],&yList[0],max_iter,jx0,jy0); dst.pixels(x,y)=coloring7(iter,max_iter,&xList[0],&yList[0],errorColor,errorColor,0.5); } } }
  调用该函数生成一幅图片,代码如下: (尺寸640x480,中心点(0,0),r=1.6;最大迭代次数1000)
    const long max_iter=1000; TViewRect rect; rect.x0=0; rect.y0=0; rect.r=1.6; TPixels32 dstPic; dstPic.resizeFast(640,480); draw_julia8(dstPic.getRef(),rect,max_iter,jx0,jy0); { //save pic TFileOutputStream bmpOutStream(dstFileName); TBmpFile::save(dstPic.getRef(),&bmpOutStream);//保存结果图 片 } 
生成的图片如下:
   
    (Julia集 zM=(-0.74543,0.11301) )
   
    (Julia集 zM=(0.28888,0.012325) )
   
    (Julia集 zM=(-0.81442,0.19633) )


ps:文章没有过多的关注例子图片是否好看(而偏向于达到该目的一些手段和技巧),都用了统一的着色方案,实际

上 生成好看的分形图像 才是本系列文章的真正重点,这也是编写分形程序的最重要的动力!
   要生成好看的分形图像,我觉得要点在于图形和颜色;图形(也就是构图和形状)还稍好办一些,而赋予合适的颜
色 更有挑战性!有时候也需要一些运气成分;
   比如,我对第3副Julia集图片的颜色生成函数进行了一些调整(后面的文章会涉及一些),得到了一幅新的图片: