|
|
楼主 |
发表于 2014-7-7 18:51
|
显示全部楼层
以下开始介绍几例最简单的分形算法：
一、Cantor三分集的递归算法
选取一个欧氏长度的直线段，将该线段三等分，去掉中间一段，剩下两段。将剩下的两段分别再三等分，各去掉中间一段，剩下四段。将这样的操作继续下去，直到无穷，则可得到一个离散的点集。点数趋于无穷多，而欧氏长度趋于零。经无限操作，达到极限时所得到的离散点集称之为Cantor集。
1.给定初始直线两个端点的坐标（ax，ay）和（bx，by），按Cantor三分集的生成规则计算出个关键点的坐标如下：
cx=ax+(bx-ax)/3
cy=ay-d
dx=bx-(bx-ax)/3
dy=by-d
ay=ay-d
by=by-d
2.利用递归算法，将计算出来的新点分别对应于（ax，ay）和（bx，by），然后利用步骤1的计算关系计算出下一级新点（cx，cy）和（dx，dy），并压入堆栈。
3.给定一个小量c，当（bx，by）<c时，被压入堆栈中的值依次释放完毕，同时绘制直线段（ax，ay）-（bx，by），然后程序结束。
下面给出matlab程序：
function f=cantor(ax,ay,bx,by)
c=0.005;d=0.005;
if (bx-ax)>c
x=[ax,bx];y=[ay,by];hold on;
plot(x,y,'LineWidth',2);hold off;
cx=ax+(bx-ax)/3;
cy=ay-d;
dx=bx-(bx-ax)/3;
dy=by-d;
ay=ay-d;
by=by-d;
cantor(ax,ay,cx,cy);
cantor(dx,dy,bx,by);
end
运行cantor（0,5,5,5），出现图例如下：
二、Koch曲线的递归算法
在一单位长度的线段上对其三等分，将中间段直线换成一个去掉底边的等边三角形，再在每条直线上重复以上操作，如此进行下去直到无穷，就得到分形曲线Koch曲线。
1.给定初始直线（ax，ay）-（bx，by），按Koch曲线的构成原理计算出各关键点坐标如下：
cx=ax+(bx-ax)/3
cy=ay+(by-ay)/3
ex=bx-(bx-ax)/3
ey=by-(by-ay)/3
l=sqrt((ex-cx)^2+(ey-cy)^2)
alpha=atan((ey-cy)/(ex-cx))
dy=cy+sin(alpha+pi/3)*l
dx=cx+cos(alpha+pi/3)*l
2.利用递归算法，将计算出来的新点分别对应于（ax，ay）和（bx，by），然后利用步骤1中的计算公式计算出下一级新点（cx，cy），（dx，dy），（ex，ey），并压入堆栈。
3.给定一个小量c，当l<c时，被压入堆栈中的值依次释放完毕，同时绘制直线段（ax，ay）-（bx，by），然后结束程序。
下面给出matlab程序：
function f=Koch(ax,ay,bx,by,c)
if (bx-ax)^2+(by-ay)^2<c
x=[ax,bx];y=[ay,by];
plot(x,y);hold on;
else
cx=ax+(bx-ax)/3; cy=ay+(by-ay)/3;
ex=bx-(bx-ax)/3; ey=by-(by-ay)/3;
l=sqrt((ex-cx)^2+(ey-cy)^2);
alpha=atan((ey-cy)/(ex-cx));
if (alpha>=0&(ex-cx)<0)|(alpha<=0&(ex-cx)<0)
alpha=alpha+pi;
end
dy=cy+sin(alpha+pi/3)*l;
dx=cx+cos(alpha+pi/3)*l;
Koch(ax,ay,cx,cy,c);
Koch(ex,ey,bx,by,c);
Koch(cx,cy,dx,dy,c);
Koch(dx,dy,ex,ey,c);
end
运行Koch(0,0,100,0,10)，出现图例如下：
三、生成填充Julia集
1.设定参数a，b以及一个最大的迭代步数N。
2.设定一个限界值R，即实数R≧max（2，sqrt（a^2+b^2）。
3.对于平面上以R为半径的圆盘内的每一点进行迭代，如果对于所有的n≦N，都有|x^2+y^2|≦R，那么，在屏幕上绘制出相应的起始点，否则不绘制。
下面给出matlab程序：
a=-0.11;b=0.65;r=2;
for x0=-1:0.01:1
for y0=-1:0.01:1
x=x0;y=y0;
if x0^2+y0^2<1
for n=1:80
x1=x*x-y*y+a;
y1=2*x*y+b;
x=x1;
y=y1;
end
if (x*x+y*y)<r
plot(x0,y0);
end
hold on;
end
end
end
a=-0.11,b=0.65
a=-0.13,b=0.77
a=-0.19,b=0.6557
四、牛顿迭代
牛顿迭代是在数值求解非线性方程（组）的时候经常使用的方法。有些牛顿迭代能够绘制出漂亮的图形来，所以现在也常用于设计图形。
Matlab程序如下：
首先编写newton函数：
function y=newton(z)
if (z==0)
y=0;
return;
end
for i=1:1:2000
y=z-(z^3-1)/(3*z^2);
if (abs(y-z)<1.0e-7)
break;
end
z=y;
end
接着进入主程序：
clear all;clc;
A=1;B=0;C=1;
for a=-1:0.005:1
for b=-1:0.005:1
x0=a+b*i;
y=newton(x0);
if abs(y-A)<1.0e-6
plot(a,b,'r');hold on;
elseif abs(y-B)<1.0e-6
plot(a,b,'g');hold on;
elseif abs(y-C)<1.0e-6
plot(a,b,'y');hold on;
end
end
end
五、迭代函数系IFS
IFS是分形的重要分支。它是分形图像处理中最富生命力而且最具有广阔应用前景的领域之一。这一工作最早可以追溯到Hutchinson于1981年对自相似集的研究。美国科学家M.F.Barnsley于1985年发展了这一分形构型系统，并命名为迭代函数系统（Iterated Function System，IFS），后来又由Stephen Demko等人将其公式化，并引入到图像合成领域中。IFS将待生成的图像看做是由许多与整体相似的（自相似）或经过一定变换与整体相似的（自仿射）小块拼贴而成。
算法：
1.设定一个起始点（x0，y0）及总的迭代步数。
2.以概率P选取仿射变换W，形式为
X1=a x0+b y0 +e
Y1=c x0+d y0+f
3.以W作用点（x0，y0），得到新坐标（x1，y1）。
4.令x0=x1，y0=y1。
5.在屏幕上打出（x0，y0）。
6.重返第2步，进行下一次迭代，直到迭代次数大于总步数为止。
下面给出一些IFS植物形态的matlab程序：
a=[0.195 -0.488 0.344 0.433 0.4431 0.2452 0.25;
0.462 0.414 -0.252 0.361 0.2511 0.5692 0.25;
-0.058 -0.07 0.453 -0.111 0.5976 0.0969 0.25;
-0.035 0.07 -0.469 -0.022 0.4884 0.5069 0.2;
-0.637 0 0 0.501 0.8562 0.2513 0.05];
x0=1;y0=1;
for i=1:10000
r=rand;
if r<=0.25
x1=a(1,1)*x0+a(1,2)*y0+a(1,5);
y1=a(1,3)*x0+a(1,4)*y0+a(1,6);
end
if r>0.25 & r<=0.5
x1=a(2,1)*x0+a(2,2)*y0+a(2,5);
y1=a(2,3)*x0+a(2,4)*y0+a(2,6);
end
if r>0.5 & r<=0.75
x1=a(3,1)*x0+a(3,2)*y0+a(3,5);
y1=a(3,3)*x0+a(3,4)*y0+a(3,6);
end
if r>0.75 & r<=0.95
x1=a(4,1)*x0+a(4,2)*y0+a(4,5);
y1=a(4,3)*x0+a(4,4)*y0+a(4,6);
end
if r>0.95 & r<=1
x1=a(5,1)*x0+a(5,2)*y0+a(5,5);
y1=a(5,3)*x0+a(5,4)*y0+a(5,6);
end
x0=x1;y0=y1;
plot(x1,y1);hold on;
end
得到图例如下：
修改部分系数便可得到另一种形态：
六、三角形分形
function triangles(n);
clc;close all;
if nargin==0;
n=4;
end
rand('state',2);
C=rand(n+4,3);
figure;
axis square equal;hold on;
a=-pi/6;
p=0;
r=1;
[p,r,n,a]=tritri(p,r,n,a,C);
function [p,r,n,a]=tritri(p,r,n,a,C);
% 画一个三角形
% p 是三角形中心
% r是三角形半径
% n是递归次数
% a是三角形角度
% C是颜色矩阵
z=p+r*exp(i*([0:3]*pi*2/3+a));
zr=p+r*exp(i*([0:3]*pi*2/3+a))/2;
pf=fill(real(z),imag(z),C(n+2,:));
set(pf,'EdgeColor',C(n+2,:));
if n>0;
[p,r,n,a]=tritri(p,r/2,n-1,a+pi/3,C);
n=n+1;r=r*2;a=a-pi/3;
[zr(1),r,n,a]=tritri(zr(1),r/4,n-1,a,C);
n=n+1;r=r*4;
[zr(2),r,n,a]=tritri(zr(2),r/4,n-1,a,C);
n=n+1;r=r*4;
[zr(3),r,n,a]=tritri(zr(3),r/4,n-1,a,C);
n=n+1;r=r*4;
end
七、曼德布罗集合
Mandelbrot set是在复平面上组成分形的点的集合。Mandelbrot集合可以用复二次多项式f(z)=z^2+c来定义。其中c是一个复参数。对于每一个c，从z=0开始对f(z)进行迭代序列 (0, f(0), f(f(0)), f(f(f(0))), .......)的值或者延伸到无限大，或者只停留在有限半径的圆盘内。曼德布罗集合就是使以上序列不延伸至无限大的所有c点的集合。从数学上来讲，曼德布罗集合是一个复数的集合。一个给定的复数c或者属于曼德布罗集合M，或者不是。
1.设定参数a，b，以及一个最大的迭代步数N。
2.设定一个限界值R，不妨设实数R=2。
3.对于参数平面上的每一点c（a，b），使用以R为半径的圆盘内的每一点进行迭代，如果对于所有的n≤N，都有|x*x+y*y|≤R*R，那么，在屏幕上绘制出相应的起始点c（a，b），否则不绘制。
下面给出matlab程序：
r=4;%限界值
for a=-2:0.002:1
for b=-2:0.002:1%参数a，b取到一个范围
x=a;y=b;%初始的复数c
for n=1:20
x1=x*x-y*y+a;%复数平方加一个c的运算
y1=2*x*y+b;
x=x1;%迭代
y=y1;
end
if(x*x+y*y)<r%限界
plot(a,b);
end
hold on;
end
end
八、脑分形
作为IFS的一种应用
a=[0.03 0 0 0.45 0 0 0.05;
-0.03 0 0 -0.45 0 0.4 0.15;
0.56 -0.56 0.56 0.56 0 0.4 0.4;
0.56 0.56 -0.56 0.56 0 0.4 0.4];
x0=1;y0=1;
for i=1:100000
r=rand;
if r<=0.05
x1=a(1,1)*x0+a(1,2)*y0+a(1,5);
y1=a(1,3)*x0+a(1,4)*y0+a(1,6);
end
if r>0.05 & r<=0.2
x1=a(2,1)*x0+a(2,2)*y0+a(2,5);
y1=a(2,3)*x0+a(2,4)*y0+a(2,6);
end
if r>0.2 & r<=0.6
x1=a(3,1)*x0+a(3,2)*y0+a(3,5);
y1=a(3,3)*x0+a(3,4)*y0+a(3,6);
end
if r>0.6 & r<=1
x1=a(4,1)*x0+a(4,2)*y0+a(4,5);
y1=a(4,3)*x0+a(4,4)*y0+a(4,6);
end
x0=x1;y0=y1;
plot(x1,y1);
hold on;
end
最简单的brain fractal
|
|
发表于 2014-7-7 18:33
提升卡
变色卡
显身卡
