图算法的应用以及在Matlab中的实现

图算法的应用以及在Matlab中的实现
图算法的应用以及在Matlab中的实现
我们首先引入一个迷宫的路径求解问题。

我们用这样一个规模为N×M布尔类型的二维矩阵gaze[N][M]来描述迷宫：矩阵的每一个元素为平面上一小块方形区域，元素值为0代表该点不可通过，1代表可以通过。

下面给出一个迷宫的例子：令矩阵的规模为11×11。

maze[N][M]为：
0 0 0 0 0 0 0 0 0 0 0
1 1 0 0 1 1 1 1 1 0 0
0 1 1 0 1 0 0 0 1 1 0
0 0 1 1 1 1 1 1 0 1 0
0 0 1 0 0 0 0 1 0 1 0
0 0 1 1 1 1 0 1 0 1 0
0 0 0 0 0 0 0 1 0 1 0
0 1 1 1 1 0 0 1 1 1 0
0 0 0 0 1 1 1 1 0 1 0
0 0 0 0 0 0 1 1 0 1 1
0 0 0 0 0 0 0 0 0 0 0
其中入口为maze[2][1]，出口在maze[10][11]。

定义这个路径求解问题如下：找到一条从入口到出口的路径，只允许在值为1的区域上通过，而且允许走的过程中只有上下左右四个方向可以选择。

我们可以用一种很简单的方法找到一条路径：从入口处开始走，一直走到第一个分岔口A，假设有3条支路，分别编号为1、2、3，我们选择支路1，继续走，如果走不通，则返回A，走支路2。

如果遇到下一个分岔口按照同样的策略，直至走到出口处。

这样一个解的结构可以抽象成如下一棵树的结构。

由于入口和分岔口A之间仅有一条确定的路，可用A代替入口作为起始点。

O为出
口。

由于路径可能不止一条，故O可能不止1个。

在对这样一个结构作深入分析之前，我们对树的结构进行定义。

首先命名该树为T。

我们定义上面一棵树上的每个点为结点，任意2结点之间的连线成为一条边，起点A是树T的根结点，取任意一条边（X，Y），称X是Y的父亲，Y是X的儿子，比如上图中，B是E 的父亲，E是B的儿子。

如果2结点M、N有着相同的父结点，则称M、N互为兄弟结点。

没有子女的节点称为叶结点，上图中O、P、Q、C、H、I均为叶结点。

考虑除根外的一个结点，比如O，任何从根节点A到O的唯一路径上的结点称为O的一个祖先。

如果M是O的一个祖先，则O是M的一个后裔。

从根A到节点M的路径长度称为M 在树T中的深度。

结点在树中的高度是从结点向下导某个叶结点最长简单路径中边的条数，树的高度是其根的高度。

数的高度也等于树中结点的最大深度。

接下来我们以这棵树的角度来分析我们的问题。

可以看到，问题的解就是从树T到叶结点O的路径。

由于走的方向只有上下左右四个,故每个分岔口最多有3条支路，即树T中每个结点最多有3个儿子，最少有2个儿子。

我们采用一种称为“深度优先搜索（Depth-first Search）”的算法来搜索出所有路径。

正如“深度优先搜索（DFS）”这一名称所暗示的那样，这种搜索算法所遵循的搜索策略是尽可能“深”地搜索树T——这样的策略正符合我们问题解得特征，因为我们的解必然
是一条从根部A 到叶结点O 的一条或几条路。

该算法的思路如下：对于新发现的结点，如果它还有以此为起点而未探测到的边，就沿此边继续探测下去。

当顶点M 所有的边都被探寻过后，搜索将回溯到M 的父结点，搜索其他没有探测过的边。

这一过程一直进行到已发现从根结点可到达O 结点的所有路时为止。

可以看到，DFS 的搜索过程其实和我们上面讲的找出口的过程是
相同的。

我们通过对结点进行着色来表示结点的状态。

开始时，每个结点均为白色，被探测后置为黑色。

用数字0，1，2，3分别代表上、左、下、右四个方向，这样，对于搜索过程中的每个结点，我们均按照递增的方向进行探测，此时，该结点的父结点已被置为黑色，所以任何一个子结点不会去探测自己的父结点。

我们用一个数组parent[v]来记录每个结点v 的父亲。

下面给出Matlab 代码，从矩阵gaze 中利用DFS 算法找出所有路径。

为了减少输入量，我们先将maze 在Microsoft Office Excel 中输入，然后在Matlab 中调用。

1. >>maze(2,1)=0;
2. >>maze1=reshape(transpose(maze),1,121);
3. >>s=distance([2 2]);
4. >>f=distance([10 11]);
5. >>cl=zeros(1,50);
6. >>cl(s)=1;
7. >>DFS(s,s,f,maze1,cl,parent);
参数含义：s为起点，f为终点，cl为每个点的颜色，parent记录每个点的父结点的下标。

其中，第3行函数distance包含在distance.m文件中，函数如下：
1.function n=distance(p)
2.n=(p(1)-1)*11+p(2);
其中的数值“11”是矩阵列数，可根据实际矩阵修改。

从输入的第2行（maze1=reshape(transpose(maze),1,121);）以及distance函数可以看出，我们为简化问题，将二维矩阵上的问题转化为数组的问题。

DFS函数定义在DFS.m文件中：
1.function DFS(s,u,f,maze1,cl,parent)
2.if u==f
i.print_path(f,s,parent);
3.else
i.for j=0:1:3
1.v=direct(u,j);
2.if ((cl(v)==0)&&(maze1(v)==1))
a)parent(v)=u;
b)cl(v)=1;
c)DFS(s,v,f,maze1,cl,parent);
3.end
ii.end
4.end
函数参数之一的u代表正在被探测的结点，初始调用为DFS(s,s,f,maze1,cl,parent)。

其中第2.i行的print_path递归地按照数组parent的索引输出路径，该函数定义在print_path.m
文件中：
1.function print_path(v,s,parent)
2.if v==s
i.fprintf('\n(%d,%d) ', point(s));
3.else
i.print_path(parent(v),s,parent);
ii.fprintf('(%d,%d) ',point(v));
4.end
其中2.i及3.ii行中的point函数是distance函数的一个求逆函数，定义在point.m文件中：
1.function p=point(s)
2.if rem(s,11)==0
i.p=[fix((s-0.5)/11)+1 11];
3.else
i.p=[fix((s-0.5)/11)+1 rem(s,11)];
4.end
其中的数值“11”是矩阵列数，可根据实际矩阵修改。

函数fix(),rem()函数均是Matlab 内部数学函数，分别是取整和取余。

文件DFS.m的第3.i.1行的函数direct是根据方向数j来获得点u 某个方向的子结点v 的下标，u、v的下标均是一维形式。

direct函数定义在direct.m文件中：
function n=direct(i,j)
1.switch j
i.case 0
1.n=i-11;
ii.case 1
1.n=i+1;
iii.case 2
1.n=i+11;
iv.case 3
1.n=i-1;
2.end
运行以上程序，输出结果为：
(2,2) (3,2) (3,3) (4,3) (4,4) (4,5) (3,5) (2,5) (2,6) (2,7) (2,8) (2,9) (3,9) (3,10) (4,10) (5,10) (6,10) (7,10) (8,10) (9,10) (10,10) (10,11) (2,2) (3,2) (3,3) (4,3) (4,4) (4,5) (4,6) (4,7) (4,8) (5,8) (6,8) (7,8) (8,8) (8,9) (8,10) (9,10) (10,10) (10,11)
输入命令imagesc (maze)，可以看到由矩阵maze生成的迷宫，其中红色是通的，蓝色是障碍物。

入口在初始输入的第1行被改为maze(2,1)，maze(2,1)被封上，这是为了使矩阵的边界都不可走，使边界条件更加简单无需额外判定。

很容易验证，上述的输出结果正是所有从入口到出口的路径。

这样一来，我们的迷宫路径求解问题得到圆满解决。

接下来看一个我们更加感兴趣的问题：
假设在某即时战略游戏中（比如说时下十分流行的WarCraft III
Frozen Throne），你的部队位于点s，当你选取了一个目的地f，计算机如何自动地在这个地形复杂（有平地、山川、河流、桥梁、森林等等地形）的区域找出一条最短路径过去。

受迷宫路径解得启发，我们可以通过搜索出所有路径然后选取长度最短的一条的方式来解决这个问题。

但是我们很快就会发现，这样的计算量将非常之大。

首先，为了更接近于真实，我们可走的选择不能仅局限为上下左右4个方向，我们可以姑且认为可以沿8个方向行进：东、西、南、北、东北、东南、西南、西北。

这样一来，这地图的树结构中的每个节点的子结点的最大数目将达到7个。

又由于地形的实际情况，可能的平地将会比较多，这样一来这棵树将会较之迷宫的树丰满得多，一些值为f的节点的高度也将特别大。

如果还是从根节点到叶结点的一条一条路径找下去，计算量之大可想而知。

可以看到，最短路径的长度就是值为f的所有结点的最小深度。

所以我们可以用另一种策略，按照树的高度一层一层地进行搜索。

与深度优先搜索优先搜索子结点的策略不同，我们优先搜索结点的兄弟结点。

这种搜索称为“广度优先搜索（Breadth-First Search，BFS）”。

BFS的过程不能够递归地加以描述，所以我们需要一种合适的方法来管理结点的搜索顺序。

这实际上也很简单，我们可以仅仅用一种称为队列的数据结构就可以成功地管理结点。

队列实现了一种先进先出（First-in，First-out）的策略。

而递归实际上是使数据后进先出（Last-in，First-out）的一种方法。

为了减少函数间的数据传递数量，我们这次更多地使用全局变量。

当我们获取地图的矩
阵以及获得起点和终点坐标(x1,y1)、(x2,y2)以后，输入如下代码即可：
1.>>global cl Q head tail;
2.>>map1=reshape(transpose(map),1, SIZE);
3.>>s=distance([x1 y1]);
4.>>f=distance([x2 y2]);
5.>>inti(s);
6.>>BFS(s,f,map1);
Q是队列的数组，存储进入队列结点的一维下标。

head纪录队列首的下标，tail纪录队列尾的下标。

Q、head、tail以及颜色函数cl均为全局变量。

SIZE代表矩阵的元素个数。

第3行的init函数是对数值进行初始化的函数，在文件init.m中定义：
1.function init(s)
2.global Q tail head cl;
3.head=1;
4.tail=1;
5.cl=zeros(1,121);
6.cl(s)=1;
7.ENQ(s);
第7行的ENQ函数是队列的“入队”的操作，由文件ENQ.m定义：
1.function ENQ(u)
2.global Q head tail;
3.tail=tail+1;
4.Q(tail)=u;
可以看到，入队操作实际就是将传递来的数据放到队列尾部。

相应的出队操作是将队首的数据取出并返回，然后head值加1，表示队首向后移一个。

入队操作的函数DEQ定义在文件DEQ.m中：
1.function n=DEQ()
2.global Q head tail;
3.n=Q(head);
4.head=head+1;
下面是我们的主要函数BFS的代码，函数BFS定义在文件BFS.m 中：
1.function BFS(s,f,map)
2.global Q head tail cl;
3.flag=0;
4.u=s;
5.while u~=f
i.u=DEQ();
ii.for j=0:1:7
1.v=direct8(u,j);
2.if((v>0) && (v<=numel(map)))
a)if((cl(v)==0)&&(map(v)~=0))
i.cl(v)=1;
ii.parent(v)=u;
iii.ENQ(v);
b)end
c)if v==f
i.flag=1;
ii.break;
d)end
3.end
iii.end
iv.if v==f
1.flag=1;
2.break;
v.end
6.end
7.if flag==1
i.print_path (f,s,parent);
8.end
第5.ii.1行的direct8函数是前一个direct函数的拓展，由4个方向变为8个。

代码在direct8.m文件中：
1.function n=direct8(i,j)
2.switch j
i.case 0
1.n=i-11;
ii.case 1
1.n=i+1;
iii.case 2
1.n=i+11;
iv.case 3
1.n=i-1;
v.case 4
1.n=i-11+1;
vi.case 5
1.n=i+11+1;
vii.case 6
1.n=i+11-1;
viii.case 7
1.n=i-11-1;
3.end
其中每个case中的数值“11”是矩阵的列数，可根据实际地图情况加以修改。

我们对8个方向的搜索顺序是北、东、南、西、东北、东南、西南、西北。

由于我们对网格的划分以及方向的选取都很粗糙，所以按照这个顺序可以有效减小实际的长度，这是因为非东南西北方向的步长为2，而我们把它们与东南西北（步长为1）一样对待处理。

所以要优先搜索东南西北。

BFS.m文件的第7.i行的函数print_path沿用之前的print_path。

如果我们依然选取之前的迷宫作为输入，则需要将BFS函数的第5.ii行的循环改为：
for j=0:1:3
v=direct(u,j);
……
然后我们输入
1.>>global cl Q head tail;
2.>>maze1=reshape(transpose(maze),1,numel(maze));
3.>>s=distance([2 2]);
4.>>f=distance([10 11]);
5.>>inti(s);
6.>>BFS(s,f,maze1);
得到的结果为：
(2,2) (3,2) (3,3) (4,3) (4,4) (4,5) (4,6) (4,7) (4,8) (5,8) (6,8) (7,8) (8,8) (8,9) (8,10) (9,10) (10,10) (10,11)
即为DFS算法得出的最短路径。

接下来我们尝试使数据可视化。

我们接下来要做的是，调用BFS后，我们得到的结果不再是路径上点的坐标集合，而是在迷宫上画出的一条路径。

这实际上很简单，只要改变一下函数print_path函数中对数据处理的语句即可。

改写后的函数为print_path_plot，保存在print_path_plot.m文件中：
1.function print_path_plot(v,s,parent)
2.if v~=s
i.print_path_plot(parent(v),s,parent);
ii.p1=point(v);
iii.p2=point(parent(v));
iv.hold on;
v.plot([p1(2) p2(2)],[p1(1) p2(1)]);
3.end
然后在BFS.m文件的第7.i行将print_path (f,s,parent);改为print_path_plot(f,s,parent);然后在主窗口中输入：
1.>> imagesc (maze); figure(gcf);
2.>> init(s)
3.>> BFS(s,f,maze1)
就得到了显示了路径的迷宫。

我们对于这个迷宫已经进行了十分多的分析，可以看到，最后的图形结果是最理想的。

但是有个缺点，需要用户输入的代码比较多。

下面我们写一个入口函数Find_Shortest_Path，通过调用这一函数就可以完成所有工作。

我们首先把矩阵maze保存到当前目录下的maze.mat数据文件中。

下面是Find_Shortest_Path的代码：
1.function Find_Shortest_Path(s,f)
2.load maze.mat;
3.global Q head tail;
4.s1=distance(s);
5.f1=distance(f);
6.imagesc (maze);
7.r=input('Do you want to see the shortest path?Y/N[Y]:','s');
8.if isempty(r)
i.r='Y';
9.end
10.if r=='Y'
i.maze1=reshape(transpose(maze),1,numel(maze));
ii.init(s1);
iii.BFS(s1,f1,maze1);
11.end
我们只要在窗口中输入Find_Shortest_Path([2 2],[10 11])就可以看到迷宫，然后得到询问是否显示路径，我们直接按回车（默认为Y）就可以看到地图上的路径了。

起点和终点可以自由设定，只要在红色区域就可以了。

Fig.2 输入为Find_Shortest_Path([2 2],[10 11])的输出图形
以上所有M文件以及数据文件maze.mat均保存在MAZE目录下。

下面回到我们最感兴趣的在3维地形下的最短路径求解问题。

首先要做的是，利用Matlab的3维绘图功能为我们绘制出一块3维地形图。

我们首先想到利用函数的方法，z代表高度且z=f(x,y)，3维空间里的每个点坐标为(x,y,z)。

为了实现平地、高山、峡谷等地形，我们想到用分段函数来描述不同区域的高度。

但是这样不仅使我们的大量精力消耗在寻找合适的函数上，而且做出的图形将十分呆板。

我们想通过一个规模较小的、每个元素值代表高度的矩阵来描述一块地形。

但是同时又要要求地图比较精细，比较光滑，于是我们想到了利用二维插值来获得规模更大的矩阵从而很好地模拟地形。

有了以上所有的算法和思路，我们现在已经很容易来解决这个问题。

而二维插值和3维绘图在Matlab中的易于实现也为我们绘制3维地形创造了很好的条件。

我们首先在Excel中输入一个大小为的矩阵map，元素的值为高度。

在Matlab中导入数据，并保存在map.mat文件中。

和函数Find_Shortest_Path类似，我们写一个Find_Shortest_Path3函数，用来在3维平面上找出最短路径并以图形化的方式输出。

1.function Find_Shortest_Path3(s,f)
2.load map.mat;
3.global Q head tail cl;
4.s=s.*20;
5.f=f.*20;
6.s1=distance(s);
7.f1=distance(f);
8.[x y]=meshgrid(0:0.2:10);
9.[xi yi]=meshgrid(0:0.05:10);
10.zi=interp2(x,y,map,xi,yi,'spline');
11.mesh(xi,yi,zi);
12.xlabel(‘Y’);
13.ylabel(‘X’);
14.r=input('Do you want to see the shortest path?Y/N[Y]:','s');
15.if isempty(r)
i.r='Y';
16.end
17.if r=='Y'
i.set(0,'RecursionLimit',5000);
ii.map1=reshape(transpose(zi),1,numel(zi));
iii.init(s1);
iv.BFS3(s1,f1,map1);
18.end
除了6~8行对我们的矩阵进行二维插值外，Find_Shortest_Path3与Find_Shortest_Path 差别很小。

其中第13.i 行是设定最大递归数量。

我们首先分析改变后的广度优先搜索，即函数BFS3。

比起二维情况，三维地形下还要考虑两点之间的高度差，如果高度差大于某一设定值，就认为同样也是不能通过的。

这样，我们只要将原来BFS第 5.ii.2.a行的if((cl(v)==0)&&(map(v)~=0))改为if((cl(v)==0)&&(map(v)~=0)&&(abs(map(u)-map(v))<0.05))即可。

同时还要修改的是输出路径的函数print_path_plot。

因为我们的输出路径图形也要有高度并随着地形起伏而起伏。

我们设定路径高度的是比地面高0.2。

修改后的函数如下：
1.function print_path_plot3(v,s,parent,map)
2.if v~=s
i.print_path_plot3(parent(v),s,parent,map);
ii.p1=point(v);
iii.p2=point(parent(v));
iv.hold on;
v.plot3([p1(2) p2(2)]./20,[p1(1) p2(1)]./20,[map(v)+0.2 map(parent(v))+0.2],'*');
3.end
第2.v行的./100是因为我们矩阵的下标并不就是坐标，而是差一个常数因子。

以上在3维环境下的操作的所有M文件以及数据文件map.mat
均保存在3D_MAP目录下。

我们的矩阵是201×201，所以相关涉及到行、列操作的，都要改动。

直接调用Find_Shortest_Path3([x1 y1],[x2 y2])，即可，其中，（x1,y1）,（x2,y2）分别是起点、终点的实际坐标值。

Fig.2 输入为Find_Shortest_Path3([9 2],[9 8])的输出图形
Fig.3 输入为Find_Shortest_Path3([1.5 8.6],[1 5])的输出图形
正如我们整篇文章的内容一样，我们对各种方法的选取、对Matlab的应用是逐步深入的。

我们的处理问题的方式是随着我们的想法、随着我们对Matlab特性的不断掌握而逐渐深入的。

由于时间所限，我们现在所呈现的结果依然是十分粗糙的。

目前我们处理问题还停留在完全依赖由初始矩阵通过二维插值生成的新矩阵上。

这样不仅影响了我们的步长在不同方向上的不同，更严重的是限制了我们的方向选择。

我们想到了用圆来选取方向，并且相邻每步之间不再严格是矩阵上的点，而是该为实际坐标，高度值用周围矩阵的点来拟合。

这样的话可以根据要求的路径精度来选取方向个数以及圆的半径大小（圆半径大小即为步长，且此时步长严格相等）。

此时再对网格涂色的时候，不仅要考虑每步两点所在网格的颜色，还要考虑两点连线通过的网格的颜色。

这个过程可能更像波的传播，像惠更斯原理描述的那样一个过程。

不再限制落脚点为矩阵的上的点，可能对print_path造成一些困难，集中体现在存储空间的利用，我们的路径高度还要单独存储。

当然也可能直接用矩阵上的点来近似的方法来节约空间，而如何取近似又将是一个小的问题。

而且我们没有对数据溢出做处理。

我们不恰当地假定路径始终是存在的，而且地图四周均
是不可达区域。

虽然这些小问题可以简单地通过一两个函数或者if语句进行处理，但是我们最终还是因为嫌繁琐没有做这项工作。

不过我们意在演示如何找到路径，这并不影响。

当路径不存在，我们的算法将不作回应，这也是不恰当的，应该找到一条离目的地最近的路。

而且通过更为细致的观察，游戏中真实情况是，计算机找到的路径并不总是最短的。

这使我们相信，程序采用的并非BFS算法，而是某种贪心算法。

而且带权最短路径的搜索用的也正是贪心算法，我们又可以想到一种策略是动态地设置一定的关键点和相应权值，然后用Dijkstra算法来找到最短路径。

而我们已经没有时间来做进一步探索，包括对Dijkstra算法的研究以及用Matlab实现。

我们写了一共大概十几个主要的M文件，不包括各个M文件的变种，因而对Matlab 已有相当的认识。

作为一种描述性的语言，Matlab在实现算法上确实不方便，在我们从最初写好的C语言程序转换为Matlab语言时遇到了相当的困难，甚至一度产生了在Matlab上写个接口文件然后直接用我们的C函数完成工作的想法。

即使已经将所有的M文件写好，仍然十分令人不满的是Matlab对这些语句的执行速度。

BFS的时间复杂度是仅仅是O(V+E)*，但用Matlab运行起来实在使人郁闷。

Matlab的优势是对数据处理上，比如可以方便的对矩阵进行各种变换，比如可以随意的使用变量而不必拘泥于定义它的类型、尺寸大小，比如可以很简单地根据数据进行插值、绘图等等。

这都在我们处理过程中以及呈现的结果中得到了很好的体现。

而且通过Matlab在自动控制上的一些简单应用，我们深切体会到了Matlab在科学计算上的全面而强大功能。

而且Matlab提供了简单的创建GUI的途径，但是我们并没有时间再去学习GUI的相关内容，没有将我们的程序以一种最友好的方式呈现出来，很遗憾……
参考文献:
[1] Thomas H. Cormen , Charles E. Leiserson , Ronald L. Rivest , Clifford Stein : Introduction
to Algorithms(Second Edition). The MIT press,2001.
[2] 蔡子经,施伯乐: 数据结构教程. 复旦大学出版社,1994.
[3] 楼顺天,姚若玉,沈俊霞: MATLAB 7.X程序设计语言. 西安电子科技大学出版社,2008. *有关时间复杂度的概念可以参照Introduction to Algorithms。

其中V为树的结点数，E为边数。