分形图Koch

Koch 分形曲线
1.1 分形原理
这是一类复杂的平面曲线，可用算法描述。

从一条直线段开始，将线段中间三分之一部分用等边三角形的两条边代替，形成具有5个结点的图形(图1)；在新的图形中，又将图中每一直线段中间的三分之一部分都用一等边三角形的两条边代替，再次形成新的图形(图2)，这时，图形中共有17个结点。

这种迭代继续进行下去可以形成Koch 分形曲线。

在迭代过程中，图形中的点将越来越多，而曲线最终显示细节的多少将取决于迭代次数和显示系统的分辩率。

1.2 算法分析
算法分析：考虑由直线段(2个点)产生第一个图形(5个点)的过程。

设1P 和5P 分别为原始直线段的两个端点。

现在需要在直线段的中间依次插入三个点234,,P P P 产生第一次迭代的图形(图1)。

显然，2P 位于1P 点右端直线段的三分之一处， 4P 位于1P 点右端直线段的三分之二处；而3P 点的位置可以看成是由4P 点绕2P 旋转60度(逆时针方向)而得到的，故可以处理为向量24P P 经正交变换而得到向量23P P 。

算法如下：
(1) 2151()/3P P P P =+-；
(2) 41512()/3P P P P =+-；
(3) 3242()T P P P P A =+-⨯；
图2 第二次迭代
图1 第一次迭代
在(3)中， A 为正交矩阵：
c o s s i n 33sin cos 33A πππ
π⎡⎤-⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦
算法根据初始数据(1P 和5P 点的坐标)，产生图1中5个结点的坐标。

结点的坐标数组形成一个5×2矩阵，矩阵的第一行为1P 的坐标，第二行为2P 的坐标，……，第五行为5P 的坐标。

矩阵的第一列元素分别为5个结点的X 坐标，第二列元素分别为5个结点的Y 坐标。

进一步考虑Koch 曲线形成过程中结点数目的变化规律。

设第k 次迭代产生结点数为k n ，第k+1次迭代产生结点数为1k n +，则k n 和1k n +之间的递推关系式为143k k n n +=-。

1.4 MATLAB 实现
p=[0 0;10 0];n=2;
A=[cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3)];
for k=1:5
d=diff(p)/3;m=4*n-3;
q=p(1:n-1,:);p(5:4:m,:)=p(2:n,:);
p(2:4:m,:)=q+d;
p(3:4:m,:)=q+d+d*A';
p(4:4:m,:)=q+2*d;
n=m;
end
plot(p(:,1),p(:,2),'k')
axis equal
axis off
程序运行后，可得如下分形曲线图形：
图1所示的简单图形被称为是koch 曲线的生成元。

下面给出其它几种生成元及对应曲线：
图3 第五次迭代图形
图4 树形生成元及对应曲线
图5 第二种树形生成元及对应曲线
(a) 矩形生成元
(b) 第一种生成方式
(c) 第二种生成方式
图6 矩形生成元及两种对应曲线1.5 附源程序
第一种树形源程序：
p=[0 0;10 10];n=2;
line(p(:,1),p(:,2),'Color','k');
A=[cos(pi/6) -sin(pi/6);sin(pi/6) cos(pi/6)];
for k=1:4
i=1;
for j=1:2:n
p1=p(j,:);p2=p(j+1,:);
d=(p2-p1)/3;
r(i,:)=p1;i=i+1;
q1=p1+d;
r(i,:)=q1;i=i+1;
r(i,:)=q1;i=i+1;
q2=q1+d*A';
r(i,:)=q2;i=i+1;
r(i,:)=q1;i=i+1;
q3=p1+2*d;
r(i,:)=q3;i=i+1;
r(i,:)=q3;i=i+1;
q4=q3+d*A;
r(i,:)=q4;i=i+1;
r(i,:)=q3;i=i+1;
r(i,:)=p2;i=i+1;
xy=[q1;q2];
line(xy(:,1),xy(:,2),'Color','k');
xy=[q3;q4];
line(xy(:,1),xy(:,2),'Color','k');
end
p=r;
n=5*n;
end
axis equal
axis off
第二种树形源程序
clear
p=[0 0;0 15];
b=pi/8;e=1.5;
a=[cos(b) -sin(b);sin(b) cos(b)];
for k=1:4
i=1;
for j=1:2:n
p1=p(j,:);
p2=p(j+1,:);
d=(p2-p1)/3;
new(i,:)=p1;i=i+1;
q1=p1+d;
new(i,:)=q1;i=i+1;
new(i,:)=q1;i=i+1;
new(i,:)=p1+d+e*d*a';i=i+1;
new(i,:)=q1;i=i+1;
new(i,:)=p1+d+e*d*a;i=i+1;
new(i,:)=p1+d;i=i+1;
q2=p1+2*d;
new(i,:)=q2;i=i+1;
new(i,:)=q2;i=i+1;
new(i,:)=p1+2*d+d*a;i=i+1;
new(i,:)=q2;i=i+1;
new(i,:)=p1+2*d+d*a';i=i+1;
new(i,:)=p1+2*d;i=i+1;
new(i,:)=p2;i=i+1;
end
n=n*7;p=new;
end
for m=1:2:n
line([p(m,1) p(m+1,1)],[p(m,2) p(m+1,2)],'Color','k'); end
axis equal
●矩形生成元第一种曲线
p=[0 0;10 0];n=2;
A=[0 -1;1 0];
for k=1:4
d=diff(p/3);m=5*n-4;
q=p(1:n-1,:);
p(6:5:m,:)=p(2:n,:);
p(2:5:m,:)=q+d;
p(3:5:m,:)=q+d+1*d*A';
p(4:5:m,:)=q+2*d+1*d*A';
p(5:5:m,:)=q+2*d;
n=m;
end
plot(p(:,1),p(:,2),'k')
axis equal
axis off
●矩形生成元第二种曲线
p=[0 0;10 0];n=2;
A=[0 -1;1 0];
for k=1:4
d=diff(p/3);m=5*n-4;
q=p(1:n-1,:);
p(6:5:m,:)=p(2:n,:);
p(2:5:m,:)=q+d;
p(3:5:m,:)=q+d+.7*d*A';
p(4:5:m,:)=q+2*d+.7*d*A';
p(5:5:m,:)=q+2*d;
n=m;
end
plot(p(:,1),p(:,2),'k') axis equal
axis off。