实现高阶分形的算法

Question

我正在尝试为n阶Koch分形建立代码,定义为: 为此,我在Matlab中实现了一个代码.我已经能够开发出能够绘制二阶分形的代码(零和第一阶很简单).为此,我在这篇文章中基于代码.到目前为止,我已经有了这个,它有效,但显然只是为了第二阶段(或者一个零,只需编辑一下代码)

clear all;
clc;
close all;
hold on;
% Counter-clockwise rotation matrix
R = [cosd(60) -sind(60);sind(60) cosd(60)];
angles = [0 60 -60 0 0]; % main angles
n = 2; % order
% Length of each straight segment
seglength = 3^-(n-1);
% Initialization of variables
a=0;
b=0;
c=1/3;
d=0;
for i=1:4^(n-2)
    for j=1:4
    p1 = [a;b];
    p5 = [c;d];
    p2 = (2.*p1+p5)/3;
    p4 = (2.*p5+p1)/3;
    p3 = p2 + R*(p4-p2);
    x = [p1(1) p2(1) p3(1) p4(1) p5(1)];
    y = [p1(2) p2(2) p3(2) p4(2) p5(2)];
    line(x,y);

    a=a+seglength*cosd(angles(j));
    c=c+seglength*cosd(angles(j+1));
    b=b+seglength*sind(angles(j));
    d=d+seglength*sind(angles(j+1));
    end
end
axis equal;

我知道它不是一个干净和最佳的代码.我只是对任何一种能够帮助我开发n阶Koch分形的技巧感兴趣,只需要改变n的值.

任何帮助表示赞赏!

Answer 1

您当然可以改进算法,特别是利用向量化将计算减少到单个for循环(并使其易于扩展到任何顺序).你甚至不需要递归(尽管这是另一种选择).这是一个矢量化函数koch:

function [x, y] = koch(N, x, y)
  R = [cosd(60) -sind(60); sind(60) cosd(60)];
  x = x(:).';
  y = y(:).';
  for iOrder = 1:N  % Loop N times
    x1 = x(1:(end-1));
    x5 = x(2:end);
    x2 = (2.*x1+x5)./3;
    x4 = (x1+2.*x5)./3;
    y1 = y(1:(end-1));
    y5 = y(2:end);
    y2 = (2.*y1+y5)./3;
    y4 = (y1+2.*y5)./3;
    temp = R*[x4-x2; y4-y2];
    x = [x1; x2; x2+temp(1, :); x4];
    y = [y1; y2; y2+temp(2, :); y4];
    x = [x(:).' x5(end)];
    y = [y(:).' y5(end)];
  end
end

对于每次迭代的一组M点(两者x和y每个),每条线的起点由下式给出,x(1:(M-1))终点由下式给出x(2:M).您可以通过构建一个4 M-1×矩阵来交错这3个新点,其中顶行是您的起点,每个连续行是您的3个新点之一(按链接计算).使用(:).'(并添加最后一个终点)重新生成结果矩阵将为该行提供1-by- 4*M-3set的点,可以在下一次迭代中使用.

以下是行动中的代码:

[x, y] = koch(0, [0 1], [0 0]);
plot(x, y);
axis equal;
hold on;
[x, y] = koch(1, x, y);  % or [x, y] = koch(1, [0 1], [0 0]);
plot(x, y);
[x, y] = koch(1, x, y);  % or [x, y] = koch(2, [0 1], [0 0]);
plot(x, y);
[x, y] = koch(1, x, y);  % or [x, y] = koch(3, [0 1], [0 0]);
plot(x, y);
legend({'0' '1' '2' '3'});

请注意,你可以得到一个n^个由要么通过顺序分n和你的出发点或传球1,为分n-1^次分形.