【Matlab系列】之基于元胞自动机的城市规划（程序修正）

Date：2018/2/15  除夕

---

转载自：http://blog.csdn.net/sk18192449347/article/details/77312160

1、城市规模设计 

雄安新区占地总面积约为2000平方公里，涉及河北省雄县、容城、安新3个县及周边部分区域，地处北京、天津、保定腹地，通过ArcGIS地图软件搜索该区域并从中提取出来，区域图如下所示。 
 
为对雄安新区进行更好的仿真模拟，首先先在地图中截取雄安新区地图，然后进行边缘轮廓提取和白洋淀等不可开发地区的剔除，获得预处理图像。最后用MATlAB进行图像灰度化、二值化处理如下图所示。为后续元胞自动机提供演变地图。 
 
城市规划 模型总步骤： 
首先确定其组成的主要元素：元胞、元胞空间、元胞状态、元胞邻域及转变规则 
分析模拟城市空间结构； 
确定模型的参数：繁殖参数、扩散参数、传播参数及受规划约束参数 
确定模型所需元胞转换规则的定义 
进行城市发展模拟。

2、城市规模仿真模拟 

初始阶段： 

中期阶段：

远期阶段：

3、程序代码

clear;close all
I=imread('city.png');
level=graythresh(I);   % 图像灰度处理
bw=im2bw(I,level);  % 图像二值化处理
I=bw;
cells=double(I);
%定义 button
plotbutton=uicontrol('style','pushbutton',...
'string','Run',...
'fontsize',12,...
'position',[100,400,50,20],...
'callback','run=1;');
erasebutton=uicontrol('style','pushbutton',...
'string','Stop',...
'fontsize',12,...
'position',[300,400,50,20],...
'callback','freeze=1;');
number=uicontrol('style','text',...
'string','1',...
'fontsize',12,...
'position',[20,400,50,20]);
%初始化
cells(33,44)=2;
cells(88,31)=2;
cells(33,80)=2;
%参数定义
kuosan=0.2;fanzhi=0.2;chuanbo=0.2;chengshihua=0.0004;sch=[];skz=[];t=1;
[x,y]=size(cells);
run=0;
freeze=0;
stop=0;
while (stop==0)
    if(run==1)
for i=2:x-1
    for j=2:y-1
        if(cells(i,j)~=1)
        if(cells(i,j)==0) %自然增长
            if(rand<chengshihua)
                cells(i,j)=2;
            end
            if(aroundcenter(i,j,cells))
                if(rand<chuanbo)
                    cells(i,j)=2;
                end
            end
         end
        if(cells(i,j)==2 && rand<kuosan)  %边界增长
               if(exist1(i,j,cells))   
                   m=1+3*rand;
                   switch  m
                       case 1
                          ii=i-1;jj=j;
                       case 2
                          ii=i;jj=j-1;
                       case 3
                          ii=i;jj=j+1;
                       otherwise 4
                           ii=i+1;jj=j;
                   end
               if(cancity(ii,jj,cells))  
                   cells(ii,jj)=2;
               end
               end
        end
       if(cells(i,j)==2 && exist1(i,j,cells)) %新扩展中心型
           if(rand<fanzhi)
               if(cancity(i,j,cells))
                   cells(i,j)=3;
               end
           end
       end
        end
    end
end
ch=0;kzch=0;
for i=1:x
    for j=1:y
        if(cells(i,j)==2) ch=ch+1;end
        if(cells(i,j)==3) kzch=kzch+1;end
    end
end
sch(t)=ch;skz(t)=kzch;
t=t+1;
[A,B]=size(cells);
Area(1:A,1:B,1)=zeros(A,B);
Area(1:A,1:B,2)=zeros(A,B);
Area(1:A,1:B,3)=zeros(A,B);

for i=1:A
  for j=1:B
    if cells(i,j)==1
        Area(i,j,:)=[1,1,1];
    elseif cells(i,j)==0
        Area(i,j,:)= [255, 255, 255];
    elseif cells(i,j)==3
        Area(i,j,:)= [255,0,0];
    elseif cells(i,j)==2
        Area(i,j,:)= [255,177,0];
    end
  end
end
% for i=1:A
%   for j=1:B
%     if cells(i,j)==1
%         Area(i,j,:)=[255,0,0];
%     end
%   end
% end
pause(0.0005);
Area=uint8(Area);
imagesc(Area);
axis equal;
axis tight;
%计步
stepnumber=1+str2num(get(number,'string'));
set(number,'string',num2str(stepnumber));
    end
if freeze==1
run=0;
freeze=0;
end
drawnow
end

子函数aroundcenter.m：

function a=aroundcenter(i,j,cells)
  a=0;
  if(cells(i-1,j)==3)  a=1;end
  if(cells(i,j-1)==3)  a=1;end
  if(cells(i,j+1)==3)  a=1;end
  if(cells(i+1,j)==3)  a=1;end
end

子函数cancity.m：

function result=cancity(i,j,cells)
s=0
if(cells(i-1,j-1)==1)   s=s+1;end
if(cells(i-1,j)==1)     s=s+1;end
if(cells(i-1,j+1)==1)   s=s+1;end
if(cells(i,j-1)==1)     s=s+1;end
if(cells(i,j+1)==1)     s=s+1;end
if(cells(i+1,j-1)==1)   s=s+1;end
if(cells(i+1,j)==1)     s=s+1;end
if(cells(i+1,j+1)==1)   s=s+1;end
if(s>=4)
    result=0;
else 
    result=1;
end
end

子函数exist1.m：

function result= exist1(i,j,cells)
a=0;
if(cells(i-1,j)==2)  a=a+1;end
if(cells(i,j-1)==2)  a=a+1;end
if(cells(i,j+1)==2)  a=a+1;end
if(cells(i+1,j)==2)  a=a+1;end
if(a>=2)
    result=1;
else
    result=0;
end
end