三维种子点生长算法（以及Python递归深度问题）

前言

种子点生长算法是从一个种子点向周围的遍历图像的每个像素的图遍历算法。
通常在二维中有8邻域方法，三维中有6邻域与26邻域方法。
本文实现了三维种子点生长的6邻域算法。
种子点生长算法本质上是对图像的连通部分进行遍历，因此可以分别利用深度优先搜索与广度优先搜索的图遍历算法实现。
其中，深度优先搜索算法应用了递归的技术，当三维图形较大的时候，随着递归层数的加深，计算机堆栈区内存不足会导致程序无法继续，因此对于大的三维图形建议采用广度优先搜索算法。
两种算法的都是高效的图遍历算法，都可以实现O(n)时间的遍历。

Python的递归问题

Python默认的递归深度是1000层，当三维矩阵较大时，1000层肯定是不够的，使用如下代码可以设置递归层数限制：

import sys  #引入系统模块
sys.setrecursionlimit(3000) #设置为3000层
1
2

如果程序运行到限制的递归层数时，会出现下面情况：

RecursionError: maximum recursion depth exceeded in comparison
1

出现上面情况的时候，我们可以通过将限制递归层数的值调大。
但是递归层数也不是可以无限大的，当递归使用的堆栈区内存达到最大时，会出现以下情况。
使用Spyder运行程序时，如果递归层数太深导致计算机堆栈区内存不足，会出现以下情况：

Restarting kernel... 
[SpyderKernelApp] WARNING | No such comm: 27435501f4ff11ebaedb1208b1daf499
1
2

如果堆栈区的内存缺的不多的时候，可以不使用Spyder等IDE执行程序，建议使用CMD直接执行程序。

三维种子点生长算法

本实现的输入是一个给定的三维二值图像的矩阵，即只有0，1值的三维矩阵（0代表黑色，1代表白色）。输出是从一个白色的生长点周围寻找到的同为白色的连通区域中的所有点坐标。
定义的生长规则为，每次从当前生长点的上下左右前后6个相邻点中查看，如果相邻点中有白色（1）那就也搜索相邻白色点的邻域，以此类推。

深度优先搜索算法

算法实现

伪代码：

FloodFill(P,bitmap,replaceColor,targetColor)
　　foreach point T in Adj(P)
　　    if color at T is equal to targetColor
　　    set color at T to replaceColor
　　    FloodFill(T,bitmap,replaceColor,targetColor)
　　end
1
2
3
4
5
6

Python代码：

def FloodFill(P,figure,target_pos,targetcolor): 
    # function: the algorithm of grow
    # input: P is growing point (type is tuple)
    #        figure is 3D picture matrix (type is numpy.ndarray)
    #        target_pos is a set of position which we need
    #        targetcolor is the value of color which we want (type is int)
    # output: None
    x,y,z = P
    adj = [] #邻域列表
    # 如果在边界上，则不能越界添加邻域
    if x>0:
        adj.append((x-1,y,z))
    if x<L-1:
        adj.append((x+1,y,z))
    if y>0:
        adj.append((x,y-1,z))
    if y<W-1:
        adj.append((x,y+1,z))
    if z>0:
        adj.append((x,y,z-1))
    if z<H-1:
        adj.append((x,y,z+1))
        
    for xn,yn,zn in adj:# 六邻域
        if (xn,yn,zn) not in target_pos:# 属于邻域且未遍历过
            if figure[xn][yn][zn] == targetcolor: #如果该点为目标值，将该邻近点加入集合中，并且递归调用
                target_pos.add((xn,yn,zn))
                FloodFill((xn,yn,zn),figure,target_pos,targetcolor)
        
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

输出所有目标点到.txt文件:

with open('target_pos.txt','w') as f:
    for x in target_set:
        for s in x:
            f.write(str(s))
            f.write(" ")
        f.write('\n')

1
2
3
4
5
6
7

广度优先搜索算法

算法实现

伪代码：

FloodFill(seed,bitmap,replaceColor)
    get targetColor from bitmap using seed position
    set color at seed point to replaceColor
    set queue Q to empty
    push seed into Q
    while Q is not empty
        pop a point P from Q
        foreach point T in Adj(P)
            if color at T is equal to targetColor
                set color at T to replaceColor
                push T into Q
    end 
1
2
3
4
5
6
7
8
9
10
11
12

下面的代码稍微做了一些改动，主要是添加了一个edge与感知集sense，目的是当一个图像中有一个非封闭的连通区域，但是这个区域与外界的连通部分（开口）很小时，我们希望程序不遍历该连通区域外部的区域。
那么这个开口大小如何定义呢，如果我们认为这个开口很大，那么edge应该设置的较大，如果我们认为这个开口很小，那么edge应该设置的小一些。

def FloodFill(P,figure,target_pos,targetcolor):
    # function: the algorithm of grow
    # input: P is growing point (type is tuple) 
    #        figure is 3D picture matrix (type is numpy.ndarray)
    #        target_pos is a set of position which we need
    #        targetcolor is the value of color which we want (type is int)
    # output: None
    
    Q = collections.deque()  # 建立一个双边队列
    Q.append(P) # 添加队列元素
    
    while Q: # 当Q队列不为空时
        seed = Q.popleft() #弹出种子
        x,y,z = seed  #种子
        sense = set() # 感知集
        
        # 搜索到距离边界的edge范围就放弃搜索
        edge = 3
        if x<edge: 
            continue
        if x>L-edge-1: # 要注意到python中索引从0开始，那么长度为N的矩阵，末端的长度索引是N-1
            continue
        if y<edge:
            continue
        if y>W-edge-1:
            continue
        if z<edge:
            continue
        if z>H-edge-1:
            continue
        
        for i in range(-edge,edge+1):  #遍历以xyz为原点，半径为edge的范围  (edge*2+1)^3的方块
            for j in range(-edge,edge+1): # 要注意到range(0,2)会生成0，1两个数，没有2，python中的0:2也是这样，口诀：顾头不顾尾
                for p in range(-edge,edge+1):
                        sense.add(figure[x+i,y+j,z+p])  #将遍历得到的值加入集合，利用集合的去重特性，可以减少内存占用
        if 0 in sense:  #如果感知到黑色（0） 
            continue  #当前点周围不再查找，直接开始下一轮循环

        # 六邻域
        adj = [(x-1,y,z),(x+1,y,z),(x,y-1,z),(x,y+1,z),(x,y,z-1),(x,y,z+1)] #邻域列表
            
        for xn,yn,zn in adj:  #对于所有的邻域点
            if (xn,yn,zn) not in target_pos: # 该点已经在目标点集合
                if figure[xn][yn][zn] == targetcolor: #如果该点为目标值，将该邻近点加入集合中，并且递归调用
                    target_pos.add((xn,yn,zn)) # 该点添加到目标点集合
                    Q.append((xn,yn,zn)) # 将该点作为生长点加入队列
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

输出.txt文件：

with open('target_pos.txt','w') as f:  #打开文件  下次跑这段代码会覆盖上次跑这段代码生成的txt文件（因为没用调为追加写入形式，我们需要的正是覆盖所以不用动）
    for x in target_set:  # 遍历目标点集合所有点
        for s in x: #遍历点元组（x,y,z）中的每个值
            f.write(str(s)) # 将其转化为字符形式，写入文件，只能写字符数据类型
            f.write(" ") #用空格隔开
        f.write('\n') #不同点用回车隔开
1
2
3
4
5
6

参考资料

https://www.cnblogs.com/chnhideyoshi/p/SeedGrow2d.html
https://www.itdaan.com/blog/2016/05/17/84e0a189821d.html