基于递归法的Sierpinski“垫片”Matlab编程（添加上Matlab源码）

 

Sierpinski“垫片”的构造原理如下：首先取一个正三角形，将其等分为四个小的正三角形，并舍去中间的一个，然后将余下3个正三角形中的每一个再等分为4个更小的正三角形，并舍去各自中间的一个，进一步再将余下的9个正三角形分别按同样的方法操作取舍，如此反复操作下去，直至无穷。最后所得图形就是谢尔宾斯基三角形，亦即Sierpinski“垫片”。

由上述原理可知，该图形具有严格的相似性，可通过递归法实现计算机绘制。算法原理如下：一个正三角形的顶点位置可以由其中心O和半径R确定。为方便起见，用复数来表示正三角形的各个特征值，设正三角形的中心位置为p = p，半径为r = r，三个顶点的位置可用向量z=[x1+i*y1,x2+i*y2,x3+i*y3]表示，则有：

z = p + r * exp( i * ( [0:3] * pi * 2/3 + a ) )

其中，a 是正三角形的偏转角度，当 a = 0 时，正三角形有一条底边平行于坐标轴纵轴；当 a = -pi/6 时，正三角形有一条底边平行于坐标轴横轴。

由上式即可确定一个正三角形A，当半径缩小为 r/2 时，所得的正三角形B的面积恰为A的 1/4 ，而且易见，B的中心及其顶点恰好是将A四等分的四个小的正三角形C1、C2、C3和C4的中心，由此，通过正三角形B又可确定正三角形C1、C2、C3和C4。根据这个原理，通过递归法就能将正三角形A进行分形，得到Sierpinski“垫片”。

这种方法也可以推广到其它多边形的分形中，只要确定多边形的中心和半径，根据一定的偏转角度，就可以画出多边形分形和“雪花”，下面的程序将给出Sierpinski“垫片”、“地毯”以及六边形“雪花”的分形生成过程。

程序源码及注释：

l         Sierpinski“垫片”

function triangles(n);

%  三角形分形;

%  n 为递归次数，亦即分形的层数;

 

clc;close all;clear;

if nargin==0;        % 设定递归次数n

    n=5;

end

rand('state',2);     % 将当前状态重设到状态2，

% 使每次程序运行时生成的各个分形三角形的颜色不变

C=rand(n+4,3);    % 生成一个(n+4)*3的随机矩阵C，对不同阶段的分形填色

figure;

axis square equal;       % 保持纵、横坐标的刻度单位一致

axis([-1 11 -1 1.5]);     % 设定标度范围

grid;hold on;

a=-pi/6;          % 设定三角形的偏转角度。a=-pi/6时，三角形有一条底边平行于x轴

p=-2;            % 设定三角形中心位置。p=j时，中心位置为(j,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 是颜色矩阵;

 

s=2;                                                 % 设定三角形的中心偏移量，将不同阶段的分形图形分开

z=p+s+r*exp(i*([0:3]*pi*2/3+a));     % z 是表示主三角形各个顶点位置的矩阵

zr=p+s+r*exp(i*([0:3]*pi*2/3+a))/2;  % zr 是表示主三角形内各个分形的顶点位置的矩阵

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/2,n-1,a,C);           % 绘制主三角形右下角的分形

    n=n+1;r=r*2;

    [zr(2),r,n,a]=tritri(zr(2),r/2,n-1,a,C);           % 绘制主三角形上部的分形

    n=n+1;r=r*2;

    [zr(3),r,n,a]=tritri(zr(3),r/2,n-1,a,C);           % 绘制主三角形左下角的分形

    n=n+1;r=r*2;

end