分形曲线程序

*午夜流星*

下面的程序是koch分形，曲线弄了一天，总算捣鼓出来了。
function  fenxing(a1,b1,a2,b2,nn)
%fenxing(0,0,9,0,3)
%a1,b1,a2,b2直线两端点坐标，nn分形次数
[A,B]=fenx(a1,b1,a2,b2);
for  t=1:nn
    for n=1:length(A)/5;
    k=fenx2(A(1+5*(n-1):5*n),B(1+5*(n-1):5*n));
     for m=1:4  
     [AA(5*4*(n-1)+5*(m-1)+1:5*4*(n-1)+5*(m-1)+5),BB(5*4*(n-1)+5*(m-1)+1:5*4*(n-1)+5*(m-1)+5)]=fenx(k(m,1),k(m,2),k(m,3),k(m,4));
     end     
   end
   A=AA;
     B=BB;
end
plot(A,B)
hold on
axis equal

----------------------
function [A,B]=fenx(ax,ay,bx,by)
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 (ex-cx)<0
    alpha=alpha+pi;
end
dx=cx+cos(alpha+pi/3)*L;
dy=cy+sin(alpha+pi/3)*L;
A=[ax,cx,dx,ex,bx];
B=[ay,cy,dy,ey,by];
---------------------------
function k=fenx2(A,B)
a11=A(1);b11=B(1);
a12=A(2);b12=B(2);
a21=A(2);b21=B(2);
a22=A(3);b22=B(3);
a31=A(3);b31=B(3);
a32=A(4);b32=B(4);
a41=A(4);b41=B(4);
a42=A(5);b42=B(5);
k=[a11,b11,a12,b12;a21,b21,a22,b22;a31,b31,a32,b32;a41,b41,a42,b42];

*午夜流星*

下面的程序来自于“萝卜”的博客，我做了一些修改。
x1=[1 2 2.5 3 4];
y1=[0 0 0 0 0];
h1=plot(x1,y1,'linewidth',2,'erasemode','xor');
axis equal
axis off
for g=linspace(0,1,40)*sin(pi/3);
   y1(3)=g;
   set(h1,'ydata',y1);
   drawnow;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x2=x1(1);
y2=y1(1);
for k=2:length(x1);
   t=linspace(x1(k-1),x1(k),4) ;
   tt=[t(2),mean(t),t(3:4)];
   x2=[x2,tt];
   t=linspace(y1(k-1),y1(k),4);
   tt=[t(2),mean(t),t(3:4)];
   y2=[y2,tt];
end
A=angle((y2(4:4:end)-y2(2:4:end))*i+(x2(4:4:end)-x2(2:4:end)));
for g=linspace(0,1,40)*sin(pi/3)/3;
   y2(3:4:end)=(y2(4:4:end)+y2(2:4:end))/2+imag(g*exp(i*(A+pi/2)));
   x2(3:4:end)=(x2(4:4:end)+x2(2:4:end))/2+real(g*exp(i*(A+pi/2))) ;
   set(h1,'ydata',y2,'xdata',x2);
   drawnow;
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x3=x2(1);
y3=y2(1);
for k=2:length(x2);
   t=linspace(x2(k-1),x2(k),4);
   tt=[t(2),mean(t),t(3:4)];
   x3=[x3,tt];
   t=linspace(y2(k-1),y2(k),4);
   tt=[t(2),mean(t),t(3:4)];
   y3=[y3,tt];
end
A=angle((y3(4:4:end)-y3(2:4:end))*i+(x3(4:4:end)-x3(2:4:end)));
for g=linspace(0,1,40)*sin(pi/3)/9;
   y3(3:4:end)=(y3(4:4:end)+y3(2:4:end))/2+imag(g*exp(i*(A+pi/2)));
   x3(3:4:end)=(x3(4:4:end)+x3(2:4:end))/2+real(g*exp(i*(A+pi/2)));
   set(h1,'ydata',y3,'xdata',x3);
   drawnow;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x4=x3(1);
y4=y3(1);
for k=2:length(x3);
   t=linspace(x3(k-1),x3(k),4);
   tt=[t(2),mean(t),t(3:4)];
   x4=[x4,tt];
   t=linspace(y3(k-1),y3(k),4);
   tt=[t(2),mean(t),t(3:4)];
   y4=[y4,tt];
end
A=angle((y4(4:4:end)-y4(2:4:end))*i+(x4(4:4:end)-x4(2:4:end)));
for g=linspace(0,1,40)*sin(pi/3)/27;
   y4(3:4:end)=(y4(4:4:end)+y4(2:4:end))/2+imag(g*exp(i*(A+pi/2)));
   x4(3:4:end)=(x4(4:4:end)+x4(2:4:end))/2+real(g*exp(i*(A+pi/2)));
   set(h1,'ydata',y4,'xdata',x4);
   drawnow;
end

*午夜流星*

以上2种程序，编程思路不一样。
我的程序采用递归法，可以任意分形，“萝卜”的是机械的限定的，分一次形，程序重复一次，但是图形生成比较漂亮，有兴趣的可以修改一下，使程序简化。

*午夜流星*

function levy(n)
%  levy(9),n表示levy曲线分形次数。
x1=0;y1=0;
x2=1;y2=0;
W=zeros(2^n,4);
plot([x1,x2],[y1,y2])
hold on
[X,Y]=levy2(x1,y1,x2,y2);
hold on
W=levy1(X,Y);
w1=W(:,1:4);
w2=W(:,5:8);
w=[w1;w2];
for k=1:n
         for i=1:2^k
         [X,Y]=levy2(w(i,1),w(i,2),w(i,3),w(i,4));
         W(i,:)=levy1(X,Y);
          end         
     w1=W(:,1:4);
   w2=W(:,5:8);
   w=[w1;w2];
end
------------------------
function w=levy1(X,Y)
a1=X(1);b1=Y(1);
a2=X(3);b2=Y(3);
a3=X(3);b3=Y(3);
a4=X(2);b4=Y(2);
w=[a1,b1,a2,b2,a3,b3,a4,b4];
-----------------------
function [X,Y]=levy2(x1,y1,x2,y2)
x3=1/2*(x1+x2-y1+y2);
y3=1/2*(x1-x2+y1+y2);
a=[x1,x3,x2];
b=[y1,y3,y2];
plot(a,b)
axis equal
hold on
X=[x1,x2,x3];
Y=[y1,y2,y3];

levy曲线有很强的自适应性，比较美观，具有装饰性。这个程序仅仅只是把算法表达出来了，由于用了2层循环，n比较大的时候运行就不理想了，没有办法表现levy曲线的美了。很可惜。

[ 本帖最后由 bainhome 于 2008-6-3 20:12 编辑 ]

yycgy

厉害，流星兄可以搞个关于分形的讲座了。

*午夜流星*

function tree(n,a,b)
%  tree(8,pi/8,pi/8)
%  tree(10,pi/4,pi/4)
x1=0;y1=0;
x2=0;y2=1;
plot([x1,x2],[y1,y2])
hold on
[X,Y]=tree2(x1,y1,x2,y2,a,b);
hold on
W=tree1(X,Y);
w1=W(:,1:4);
w2=W(:,5:8);
w=[w1;w2];
for k=1:n
         for i=1:2^k
         [X,Y]=tree2(w(i,1),w(i,2),w(i,3),w(i,4),a,b);
         W(i,:)=tree1(X,Y);
          end   
   w1=W(:,1:4);
   w2=W(:,5:8);
   w=[w1;w2];
end
-------------------
function w=tree1(X,Y)
a1=X(1);b1=Y(1);
a2=X(2);b2=Y(2);
a3=X(1);b3=Y(1);
a4=X(3);b4=Y(3);
w=[a1,b1,a2,b2,a3,b3,a4,b4];
----------------------------
function [X,Y]=tree2(x1,y1,x2,y2,a,b)
L=sqrt((x2-x1)^2+(y2-y1)^2);
if (x2-x1)==0
    a=pi/2;
   else if (x2-x1)<0
   a=pi+atan((y2-y1)/(x2-x1));
   else
   a=atan((y2-y1)/(x2-x1));  
  end
end
x3=x2+L*2/3*cos(a+b);
y3=y2+L*2/3*sin(a+b);
x4=x2+L*2/3*cos(a-b);
y4=y2+L*2/3*sin(a-b);
a=[x3,x2,x4];
b=[y3,y2,y4];
plot(a,b)
axis equal
hold on
X=[x2,x3,x4];
Y=[y2,y3,y4];

[ 本帖最后由 bainhome 于 2008-6-3 20:13 编辑 ]

*午夜流星*

改变夹角参数的树型图。

WaitingForMe

试了一下那个tree函数,发现实在是非常慢
优化了一下
并提供一个比较实用的旋转数据的函数:

rotatedata

function [newx,newy,newz]=rotatedata(xdata,ydata,zdata,azel,alpha,origin)
%在三维空间种将数据进行旋转，修改自matlab自带的rotate函数
%
%   Modified by WaitingForMe
%   heroaq_2002@163.com
%   2005/3/4
%   

%ROTATE Rotate objects about specified origin and direction.
%   ROTATEDATA(Xdata,Ydata,Zdata,[THETA PHI],ALPHA,ORIGIN) rotates the objects with handles H
%   through angle ALPHA about an axis described by the 2-element
%   direction vector [THETA PHI] (spherical coordinates).  
%   All the angles are in degrees.  The handles in H must be children
%   of the same axes.
%
%   THETA is the angle in the xy plane counterclockwise from the
%   positive x axis.  PHI is the elevation of the direction vector
%   from the xy plane (see also SPH2CART).  Positive ALPHA is defined
%   as the righthand-rule angle about the direction vector as it
%   extends from the origin.
%
%   ROTATEDATA(Xdata,Ydata,Zdata,[X Y Z],ALPHA,ORIGIN) rotates the objects about the direction
%   vector [X Y Z] (Cartesian coordinates). The direction vector
%   is the vector from the center of the plot box to (X,Y,Z).
%
%   See also SPH2CART, CART2SPH.
%

%判断输入变量的个数
if nargin<6
    error('Not enough input arguments! Type ''help rotatedata'' to get some help!')
end
%找到旋转的单位轴向量
if prod(size(azel)) == 2 % theta, phi
    theta = pi*azel(1)/180;
    phi = pi*azel(2)/180;
    u = [cos(phi)*cos(theta); cos(phi)*sin(theta); sin(phi)];
elseif prod(size(azel)) == 3 % direction vector
    u = azel(:)/norm(azel);
end

alph = alpha*pi/180;
cosa = cos(alph);
sina = sin(alph);
vera = 1 - cosa;
x = u(1);
y = u(2);
z = u(3);
%旋转矩阵
rot = [cosa+x^2*vera x*y*vera-z*sina x*z*vera+y*sina; ...
        x*y*vera+z*sina cosa+y^2*vera y*z*vera-x*sina; ...
        x*z*vera-y*sina y*z*vera+x*sina cosa+z^2*vera]';
[m,n] = size(xdata);
if isempty(zdata)
    zdata=zeros(size(xdata));%在对二维数据进行旋转时，z参数输入可以是空矩阵
end
newxyz = [xdata(:)-origin(1), ydata(:)-origin(2), zdata(:)-origin(3)];
newxyz = newxyz*rot;
newx = origin(1) + reshape(newxyz(:,1),m,n);
newy = origin(2) + reshape(newxyz(:,2),m,n);
newz = origin(3) + reshape(newxyz(:,3),m,n);

这个函数修改自matlab内部函数rotate,用法是一样的,只不过返回输入输入旋转以后的数据

下面是同样层数计算的时间
这个是我的程序的运行时间:
tic;mytree(10,45,0.618);toc
elapsed_time =
0.2810
这个是上面的原程序运行的时间:
tic;fsimwe(10,pi/4,pi/4);toc
elapsed_time =
51.4680
不过程序中生成树枝的方式有一点不一样,上面程序中的树枝好象不往下面长,呵呵
========================================
function mytree(n,ang,p,isrand)
%生成一个数
%n,层数
%ang,左右旋转的角度,[ang1,ang2],如果输入的是一个标量数据,则自动修改为[-ang,ang]
%p,每次长出新芽的生长比例
rt=[0;0];
br=[0;1];
plot([0 0],[0 1]);
if nargin<4 isrand=0; elseif isrand>1 isrand=0.99; end
for i=1:n
    [rt,br]=grow(rt,br,ang,p,isrand);
    treeplot(rt,br);
end

function [rt1,br1]=grow(rt,br,ang,p,isrand)
l=br-rt;
brext=p.*l;
[rt1,br1]=rot(br,brext,ang);
%加入一个随机成活的控制
if size(rt1,2)>1000 & isrand
    alive=rand(1,size(rt1,2));
    rt1(:,alive<isrand)=[];
    br1(:,alive<isrand)=[];
end

function [rt,brg]=rot(br,brext,ang)
%分别向左右偏离一个角度
if length(ang)==1
    ang=[ang -ang];
end
[x,y]=rotatedata(brext(1,:),brext(2,:),[],[0 0 1],ang(1),[0 0 0]);
l=br+[x;y];
[x,y]=rotatedata(brext(1,:),brext(2,:),[],[0 0 1],ang(2),[0 0 0]);
r=br+[x;y];
brg=[l,r];
rt=[br,br];

function treeplot(r,b)
%画出树枝
line([r(1,:);b(1,:)],[r(2,:);b(2,:)]);

==========================================

后来加入了一个随机树控制生长树枝的概率,给一个图:

*午夜流星*

感谢WaitingForMe 的优化，我所做的只是表现了递归算法，没有深层次探索。

*午夜流星*

下面是台湾的單維彰先生写的koch曲线，是我发现的最简单的程序，函数内部自调用，很简洁。
http://libai.math.ncu.edu.tw/bcc16/B/matlab/d15.shtml
function koch(p,q,n)
%koch([0;0],[1;0],5)
if (n==0)
    plot([p(1);q(1)], [p(2);q(2)], 'LineWidth',4,'Color','red');
    hold on;
else
    c = q-p;
    c = [-c(2); c(1)];
    c = (p+q)/2 + c/sqrt(12);   % 求出「向左側翹起 1/3」的頂點座標向量 c
    a = (2*p+q)/3;              % 求出從 p 到 q 的 1/3 處端點座標向量 a
    b = (p+2*q)/3;              % 求出從 p 到 q 的 2/3 處端點座標向量 b
    koch(p, a, n-1);            % 對 pa 線段做下一回合
    koch(a, c, n-1);            % 對 ac 線段做下一回合
    koch(c, b, n-1);            % 對 cb 線段做下一回合
    koch(b, q, n-1);            % 對 bq 線段做下一回合
end
axis equal

*午夜流星*

http://www.math.ncu.edu.tw/~liyl/
http://www.math.ncu.edu.tw/~liyl/fractal/MatlabProgram.html
台湾李易霖先生的博客，关于分形的一些有趣的东西，做的很不错。

*午夜流星*

基于LS算法的曲线程序。
function kochLS(alpha0,n)
syms F S T  
%  A=[F S S F S S F];
%  C=[F S F T T F S F];
% A=[F S F S F S F S];  
%C=[F F S F S F S F S F S F T F];alpha0=pi/2
x1=0; y1=0;
x2=0.1; y2=0;  %起始点
L=0.1;         %段长

A=[F T T F T T F ];  % 初始公理
C=[F S F T T F S F]; %生成规则
%alpha0旋转角度
   for j=1:n
      k=1;
    for i=1:length(A)
     if A(i)==F
         B(k+length(C):length(A)+length(C)-i+k-1)=(A(i+1:end));
         B(k:k+length(C)-1)=C;
         k=k+length(C);         
     else         
        B(k)=A(i);
         k=k+1;         
     end
   end
A=B;
end

  for kk=1:length(A)
      if A(kk)==F
          plot([x1,x2],[y1,y2]),hold on,axis equal
              if (x2-x1)==0
                 alpha=pi/2;
                else if (x2-x1)<0
                  alpha=pi+atan((y2-y1)/(x2-x1));
                 else
                  alpha=atan((y2-y1)/(x2-x1));
                 end
              end
        x1=x2;
        y1=y2;
         else if A(kk)==S
             alpha=alpha+alpha0;
              else if A(kk)==T
               alpha=alpha-alpha0;
              end
          end
       end
       x2=x1+L*cos(alpha);
       y2=y1+L*sin(alpha);
end

*午夜流星*

关于IFS玛的分形
function sierpinski_ifs(n,w1,w2,w3)
%sierpinski_ifs(10000,1/3,1/3,1/3)
%w1,w2,w3出现频率
M1=[0.5 0  0   0 0.5  0];
M2=[0.5 0 0.5  0 0.5  0];
M3=[0.5 0 0.25 0 0.5  0.5];
x=0;y=0;
r=rand(1,n);
B=zeros(2,n);
k=1;
for i=1:n   
    if r(i)<w1
        a=M1(1);b=M1(2);c=M1(3);d=M1(4);e=M1(5);f=M1(6);
    else if  r(i)<w1+w2
        a=M2(1);b=M2(2);c=M2(3);d=M2(4);e=M2(5);f=M2(6);  
        else if  r(i)<w1+w2+w3
        a=M3(1);b=M3(2);c=M3(3);d=M3(4);e=M3(5);f=M3(6);  
        end
    end
    end
x=a*x+b*y+c;
y=d*x+e*y+f;
B(1,k)=x;
  B(2,k)=y;
  k=k+1;
end   
plot(B(1,:),B(2,:),'.','markersize',0.1)

*午夜流星*

function helix_ifs(n,w1,w2,w3)
%helix_ifs(20000,0.9,0.05,0.05)
%w1,w2,w3出现频率
M1=[0.787879  -0.424242  1.758647 0.242424  0.859848  1.408065];
M2=[-0.121212 0.257576  -6.721654 0.05303   0.05303   1.377236];
M3=[0.181818  -0.136364  6.086107 0.090909  0.181818  1.568035];
x=0;y=0;
r=rand(1,n);
B=zeros(2,n);
k=1;
for i=1:n   
    if r(i)<w1
        a=M1(1);b=M1(2);c=M1(3);d=M1(4);e=M1(5);f=M1(6);
    else if  r(i)<w1+w2
        a=M2(1);b=M2(2);c=M2(3);d=M2(4);e=M2(5);f=M2(6);  
        else if  r(i)<w1+w2+w3
        a=M3(1);b=M3(2);c=M3(3);d=M3(4);e=M3(5);f=M3(6);  
        end
    end
    end
x=a*x+b*y+c;
y=d*x+e*y+f;
B(1,k)=x;
  B(2,k)=y;
  k=k+1;  
end      
plot(B(1,:),B(2,:),'.','markersize',0.1)

[ 本帖最后由 bainhome 于 2008-6-3 20:12 编辑 ]

*午夜流星*

julia集合IFS算法，算法看起来好像没有问题，但是内部图像没有显示完全
function julia_ifs(n,cx,cy)
%  julia_ifs(100000,-0.77,0.08)
%%%  f(z)=z^2+c,cx=real(c);cy=image(c);
x=1;y=1;   %z^2+c=z0,x=real(z0);y=image(z0);
B=zeros(2,n);
k=1;
A=randn(1,n);
for i=1:n
    wx=x-cx;
    wy=y-cy;
    if wx>0
        alpha=atan(wy/wx);
    end
    if wx<0
        alpha=pi+atan(wy/wx);
    end
    if wx==0
        alpha=pi/2;
    end
    alpha=alpha/2;
    r=sqrt(wx^2+wy^2);
      if A(i)<0
        r=-sqrt(r);
      else
        r=sqrt(r);
      end
    x=r*cos(alpha);
    y=r*sin(alpha);
    B(1,k)=x;
    B(2,k)=y;
    k=k+1;   
end
plot(B(1,:),B(2,:),'.','markersize',0.1)

*午夜流星*

逃逸时间算法，sierpinski垫片
function sierpinski(a,b,c,d,n,m,r)
%sierpinski(0,0,1,1,12,200,200)
%(a,b),(c,d)收敛区域左上角和右下角坐标，m为分辨率
% n为逃逸时间，需要反复试探，r逃逸半径
B=zeros(2,m*m);
w=1;
for i=1:m
    x0=a+(c-a)*(i-1)/m;
    for j=1:m
        y0=b+(d-b)*(j-1)/m;
        x=x0;
        y=y0;
      for k=1:n
          if y>0.5
              x=2*x;
              y=2*y-1;
           elseif x>=0.5
               x=2*x-1;
               y=2*y;  
           else
                x=2*x;
                y=2*y;
            end
                 if x^2+y^2>r
                  break;
                  end
         end
         if k==n
              B(1,w)=i;
              B(2,w)=j;
              w=w+1;  
          end
    end
end
plot(B(1,:),B(2,:),'.','markersize',0.1)

bainhome

再鼓励一下,虽然我已经被这些算法搞得晕乎了
:D:D:D

*午夜流星*

元胞自动机算法
function sierpinski_ca(m,n)
%sierpinski_ca(1000,3000)
%一维元胞自动机
x=1;y=1;
t=1;w=zeros(2,m*n);
s=zeros(m,n);
s(1,fix(n/3))=1;
for i=1:m-1
    for j=2:n-1
         if (s(i,j-1)==1&s(i,j)==0&s(i,j+1)==0)|...
                     (s(i,j-1)==0&s(i,j)==0&s(i,j+1)==1)
            s(i+1,j)=1;
            w(1,t)=x+3+3*j;
            w(2,t)=y+5*i;
            t=t+1;
        end   
   end
end
plot(w(1,:),w(2,:),'.','markersize',1)

*午夜流星*

function sierpinski_ca2(m,n)
%sierpinski_ca2(400,400)
%二维元胞自动机
t=1;w=zeros(2,m*n);
s=zeros(m,n);
s(m/2,n/2)=1;
for i=[m/2:-1:2,m/2:m-1]
     for j=[n/2:-1:2,n/2:n-1]
           if mod(s(i-1,j-1)+s(i,j-1)+s(i+1,j-1)+s(i-1,j)+s(i+1,j)+...
                s(i-1,j+1)+s(i,j+1)+s(i+1,j+1),2)==1
            s(i,j)=1;
            w(1,t)=i;
            w(2,t)=j;
           t=t+1;
       end
    end
end
plot(w(1,:),w(2,:),'.','markersize',0.1)

alexqxp

学习中。