正六边形元胞自动机matlab实现
在网格结构为正方形细胞自动机的基础上，我尝试着使用matlab制作了一个网格结构为正六边形的细胞自动机。
参考链接：
参考了这位的连接，看明白了以后又自己进行了一点点改动。
基于Matlab的六边形细胞自动机
正六边形排布讲解
一、算法理论
1.绘制正六边形蜂窝网络：
要实现网格结构为正六边形的细胞自动机，首先先要绘制出正六边形蜂窝网络，这涉及到了一定的算法实现。
实现网格结构为正六边形的细胞自动机的状态转移，运算过程中需要记录储存的数据有：
1.细胞状态：0或1 ( 0代表生存，1代表死亡 )
2.细胞中心坐标：(x0, y0)
3.细胞六个边界点的坐标：(xi, yi) (i=1,2,3,4,5,6)
4.六个邻近细胞行、列号：( row, col )
考虑到以上数值都是整型或是浮点型数值，我选择用数组这一数据结构对其进行处理。
下面的过程记录会对这些数据的处理按顺序进行简单的介绍。
1 六边形网络算法简单介绍：
对于正方形来说，要构成一个网格，我们只有一种方式，但是对于正六边形来说，我们有多种方式排列，如下图所示的有边朝上的flat topped和尖朝上的pointy topped基础形式。
而在边朝上flat topped和尖朝上的pointy topped这两种排列的基础上，又可以拓展出以下偏移坐标系。
在偏移坐标系中，最常见的方法是将每一行(或列)进行偏移。每一列被命名为col(q)，每一行被称为row®。如上图所示，我们可以选择偏移奇数行或者偶数行，也可以选择偏移奇数列或者偶数列。
2 中心center点绘制：
正六边形网络的排列具有高度的规律性，在绘制网络之前，我们可以先推算出每个正六边形的中心点位置并画出，然后再根据中心点绘制完整的正六边形。
a 中心点初始化
首先我们设置正六边形的边长为1/2，如下图所示，左侧第一个正六边形的横坐标初始值可以设定为：0.5或1.25。
b 位置推演
如下图所示，正六边形的边长不变，水平方向上中心点到下一个中心点的横纵坐标变化为：
完成初始化和确定好推演规则以后，就可以编写代码先绘制出六边形网络的所有中心点了，如下图所示。
3 六个边界点坐标计算：
这里采用了flat topped类型的正六边形布局。
如下图所示，设置正六边形的边长为1/2，设正六边形的中心点坐标为(x,y),则它的六个边界点的计算如下图标注所示：
根据上一个步骤推演计算得到的中心点坐标，计算得出正六边形六个边界点坐标值，并且调用一下matlab中的patch函数进行多边形绘制：
patch(reshape(Cell(i,j,3:8,1),[1,6]),reshape(Cell(i,j,3:8,2),[1,6]),'w','EdgeColor','blue' );
至此基础的绘图完成，但是我们还需要给每一个正六边形的位置进行行列编号处理，对于边界的正六边形还需要将其设置为非细胞。
4 行列编号：
a 排除非细胞点（边界点）：
由于设定规则需要获得一个细胞周围细胞的状态，所以我们需要在整个界面边界设置一圈“非细胞”，以免计算时出现越界的情况。
对于下图网格的每个交点进行遍历，搜索该交点附近的2个网格的交点，并计算两个点之间的距离。当距离值符合一定限制时。如图所示的浅蓝色圆形为有效搜索范围，计算出大致的范围半径并判断当前正在搜索的正六边形的中心坐标是否在范围内。如果改正六边形周边六个相邻位置都有正六边形，那么它是细胞，反之不是。
部分判断代码：
%黄色区域范围计算
x_range = [1,1].* cos(PI/3* (i-1)-PI/6)* distance+ [-0.1 0.1]+ Cell(i_center, j_center,2,1);
y_range = [1,1].* sin(PI/3* (i-1)-PI/6)* distance+ [-0.1 0.1]+ Cell(i_center, j_center,2,2);
b 行列号初始化设置：
如上图所示，在确定交点为正六边形中心点，调用循环遍历的过程中，把判定为细胞的六边形的行列号一同进行标注即可。
修改正六边形的边长大小，已达到连续状态的细胞自动机可视化。
2.规则设定：
每个细胞有两种状态-存活或死亡，每个细胞与以自身为中心的周围六格细胞产生互动。
最基础的生存游戏设定的演变规则如下表所示：
1 基础生命游戏规则：
规则设定代码实现：
% 演变部分代码
if Cell( i, j, 1, 1)==1 && (N_live>3 || N_live<2)
Cell( i, j, 1, 2)= 0;
elseif Cell( i, j, 1, 1)==0 && (N_live==3 )
Cell( i, j, 1, 2)= 1;
end
按照以上规则进行随机初始化后，运行细胞自动机，最终的结果会演化呈现出几种稳定的结构：
2 规则变更-扩散效果：
规则设定：先把中间点置为1，每一时间步对每一点，如果周围六个点之和为偶数，则变为0，为奇数则变为 1。
规则设定代码实现：
% 演变部分代码
if mod(N_live,2)==0
Cell( i, j, 1, 2)= 0;
else
Cell( i, j, 1, 2)= 1;
end
按照以上规则进行初始化赋值设置中心点为1，运行细胞自动机，得到如下所示的效果：
仔细观察其实也可以发现这就是分形效果的一种演变。更改初始化点的位置，可以获得不一样的效果实现，如图蓝色为初始化四个边界点为状态1的单独效果，黄色为累积叠加效果：
3 规则变更-连续状态：
更改了一下每个细胞中存储的状态值，现在细胞中存储的不仅仅是0或1任意一个整数值，而是0-1的任意一个浮点数数值了。
下面设计了一个比较随意的规则，缩小网格数便于调试与观察。
if N_live>5
if Cell( i, j, 1, 1)-0.1>0
Cell( i, j, 1, 2)= Cell( i, j, 1, 1)-0.1;
end
elseif N_live<2
if Cell( i, j, 1, 1)+0.1<=1
Cell( i, j, 1, 2)= Cell( i, j, 1, 1)+0.1;
end
end
研究了好久Matlab官网的colormap函数，终于弄明白了这个函数的机制，自己写了一下不同数值映射到colormap的颜色上的代码：（只用了10种颜色，颜色太多网格数也太多展示效果也不明显。）
color = colormap(bone(10)); %颜色映射
x = linspace(0,1,100); %范围
color_index = ceil(x*9+1);
idx=ceil(Cell( i, j, 1, 1)/0.01);
if idx>100
idx=100;
elseif idx<1
idx=1;
end
c = color(color_index(idx),:);
set(pic{i, j},'FaceColor',c);
测试效果如下所示：
二、实现效果
挑选了一些状态值为连续的细胞自动机颜色变化的生成图出来（无聊的玩了半天= =）：
想到老师上课讲过的热传导示例程序，我决定改一下规则实现一下，尝试了以下的模板。
mark一下关键规则设定代码：
changed=Cell( i, j, 1, 1)*(-6)*delta;
for k=1:1:6
if Cell( i, j, k+8, 1)~= -1 && Cell( i, j, k+8, 2)~= -1 %非边界正六边形
changed=changed+Cell(Cell( i, j, k+8, 1), Cell( i, j, k+8, 2), 1, 1)*delta;
end
Cell( i, j, 1, 2)=Cell( i, j, 1, 1)+changed;
end
三、代码实现
备注很多，去掉备注可以自己多试试效果。
clear;
PI=3.1415926;
N_col = 10;
% 列数设置
%N_row = round(N_col/0.433*1.5/1.2);
N_row = 30;
Cell = zeros(N_row,N_col, 14,2);
% Cell - [细胞状态 细胞中点坐标 六个边点坐标 六个邻近图形序号]
% 细胞状态(1)- [当前状态 下一步状态]； 0- 死亡 ，1- 存活。
% 中方点坐标(2)- [x, y]
% 六个边点坐标(3~8)- [x，y]。
% 六个邻近图形行、列号(9~14)- [行号，列号]。
pause(1);
% 生成六边形中点坐标 %
y_point= 0;
for i_row =1:1:N_row
if(mod(i_row,2)==1)
x_point = 0.5;
else
x_point = 1.25;
end
y_point= y_point+ 0.433;
for i_col = 1:1:N_col
x_point = x_point+ 1.5;
Cell(i_row, i_col,2,:)= [x_point,y_point];
end
end
% 生成六边形角点坐标 %
Cell(:,:,3,1)= Cell(:,: ,2,1)-0.5;
Cell(:,:,[4,8],1)= repmat(Cell(:,: ,2,1)-0.25,1,1,2);
Cell(:,:,[5,7],1)= repmat(Cell(:,: ,2,1)+0.25,1,1,2);
Cell(:,:,6,1)= Cell(:,: ,2,1)+0.5;
Cell(:,:,[3,6],2)= repmat(Cell(:,: ,2,2),1,1,2);
Cell(:,:,[4,5],2)= repmat(Cell(:,: ,2,2)+0.433,1,1,2);
Cell(:,:,[7,8],2)= repmat(Cell(:,: ,2,2)-0.433,1,1,2);
% 设置背景颜色 %
figure('color',[0 0 0.2])
% 绘制中点 %
plot(Cell(:, :,2,1),Cell(:, :,2,2),'LineStyle','none','Marker','.','color','r' );
hold on;
if(Cell(N_row,N_col,2,1)>1.2*Cell(N_row,N_col,2,2))
xlim([0, Cell(N_row,N_col,2,1)]+1.5);
ylim([0, Cell(N_row,N_col,2,1)]./1.2);
else
xlim([0, Cell(N_row,N_col,2,2).*1.2+1]);
ylim([0, Cell(N_row,N_col,2,2)]);
end
% 绘制六边形 %
pic= cell(N_row,N_col);
for i=1:1:N_row
for j=1:1:N_col
pic{i,j}= patch(reshape(Cell(i,j,3:8,1),[1,6]),reshape(Cell(i,j,3:8,2),[1,6]),'w','EdgeColor',[0.9 0.9 0.9],'LineWidth',1 );
end
end
% 寻找邻近图形的行、列号（没有细胞的位置为[-1,-1]）
for i_center=1:1:N_row
for j_center=1:1:N_col
for i=1:1:6
found = 0;
% 在目标细胞附近 士2 格的范围内寻找，减少运算量
for i_find= max(1, i_center-2):min(N_row, i_center+2)
for j_find= max(1, j_center-2):min(N_col, j_center+2)
distance = sqrt((Cell(i_center, j_center,2,1)- Cell(i_find, j_find,2,1))^2 ...
+ (Cell(i_center, j_center,2,2)- Cell(i_find, j_find,2,2))^2);
% 初步根据距离判断
if distance> 0.85 && distance< 0.9 %计算范围限制
x_range = [1,1].* cos(PI/3* (i-1)-PI/6)* distance+ [-0.1 0.1]+ Cell(i_center, j_center,2,1);
y_range = [1,1].* sin(PI/3* (i-1)-PI/6)* distance+ [-0.1 0.1]+ Cell(i_center, j_center,2,2);
% 再根据夹角判断
if Cell(i_find, j_find,2,1)> x_range(1) && Cell(i_find, j_find,2,1)< x_range(2) ...
&& Cell(i_find, j_find,2,2)> y_range(1) && Cell(i_find, j_find,2,2)< y_range(2)
Cell(i_center, j_center,i+8,:) = [i_find, j_find];
found= 1;
end
end
end
end
if found ==0
Cell(i_center, j_center,i+8,:)= -1;
end
end
end
end
% 赋予细胞随机初试状态0或1
%Cell( :, :, 1, 1)=round(rand(N_row, N_col));
%对角初始化为1--------------------------------------
% Cell( 2, 2, 1, 1)=1;
% Cell( N_row-1, N_col-1, 1, 1)=1;
% Cell( 2, N_col-1, 1, 1)=1;
% Cell( N_row-1, 2, 1, 1)=1;
%中心初始化为1--------------------------------------
%Cell( N_row/2, N_col/2, 1, 1)=1;
%非连续状态--------------------------------------
%Cell( :, :, 1, 1)=rand(N_row, N_col);
Cell( N_row/2, N_col/2, 1, 1)=20;
delta=0.05;
changed=delta;
% 开始迭代
while 2>1
for i=1:1:N_row
for j=1:1:N_col
% 计算细胞周围存活细胞数量
% N_live = 0;
% for k=1:1:6
% if Cell( i, j, k+8, 1)~= -1 && Cell( i, j, k+8, 2)~= -1 %非边界正六边形
% N_live = N_live+ Cell(Cell( i, j, k+8, 1), Cell( i, j, k+8, 2), 1, 1);
% end
% end
% mod1 判断生死
% if Cell( i, j, 1, 1)==1 && (N_live>3 || N_live<2)
% Cell( i, j, 1, 2)= 0;
% elseif Cell( i, j, 1, 1)==0 && (N_live==3 )
% Cell( i, j, 1, 2)= 1;
% end
% mod2 扩散
% if Cell( i, j, 1, 1) == 0&&mod(N_live,2)==0
% Cell( i, j, 1, 2)= 0;
% elseif mod(N_live,2)==1
% Cell( i, j, 1, 2)= 1;
% elseif Cell( i, j, 1, 1) == 0&&mod(N_live,2)==1
% Cell( i, j, 1, 2)= 1;
% elseif Cell( i, j, 1, 1) == 1&&mod(N_live,2)==0
% Cell( i, j, 1, 2)= 0;
% elseif Cell( i, j, 1, 1) == 1&&mod(N_live,2)==1
% Cell( i, j, 1, 2)= 0;
% end
%mod3 连续测试
% if N_live>5
% if Cell( i, j, 1, 1)-0.1>0
% Cell( i, j, 1, 2)= Cell( i, j, 1, 1)-0.05;
% end
% elseif N_live<3
% if Cell( i, j, 1, 1)+0.1<=1
% Cell( i, j, 1, 2)= Cell( i, j, 1, 1)+0.05;
% end
% end
%mod4 扩散
disp(Cell( N_row/2, N_col/2, 1, 2));
changed=Cell( i, j, 1, 1)*(-6)*delta;
for k=1:1:6
if Cell( i, j, k+8, 1)~= -1 && Cell( i, j, k+8, 2)~= -1 %非边界正六边形
changed=changed+Cell(Cell( i, j, k+8, 1), Cell( i, j, k+8, 2), 1, 1)*delta;
end
Cell( i, j, 1, 2)=Cell( i, j, 1, 1)+changed;
end
% 更新显示状态
color = colormap(hot(20)); %颜色映射
x = linspace(0,1,1000); %范围
color_index = ceil(x*19+1);
idx=ceil(Cell( i, j, 1, 1)/0.001);
if idx>1000
idx=1000;
elseif idx<1
idx=1;
end
c = color(color_index(idx),:);
set(pic{i, j},'FaceColor',c);
% if Cell( i, j, 1, 1) >0.5
% set(pic{i, j},'FaceColor','flat');
% else
% set(pic{i, j},'FaceColor',[1 1 1]);
% end
end
end
drawnow;
% 更新生死状态
Cell( :, :, 1, 1)= Cell( :, :, 1, 2);
%pause(0.1);
end
1018
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
