11

ImagePy解析：24 -- 骨架图转图论sknw解析

 3 years ago
source link: https://qixinbo.info/2020/11/20/imagepy_24/
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
neoserver,ios ssh client

ImagePy解析：24 -- 骨架图转图论sknw解析

sknw是一个从骨架图中创建图网络的库，代码在这里。

它不仅可以实现将单线转变成图graph的效果，而且里面的trace函数还可以实现像素追踪，将图像中的单线的坐标序列依次提取出来，从而将图像转变为矢量图。（sknw可以对闭合曲线进行坐标提取，对于闭合曲线，也可以使用find_contour来提取这些坐标序列）

输入图像必须是一个二值的骨架图。
比如，这里的示例图像矩阵为：

1
2
3
4
5
6
7
8
9
10
img = np.array([
    [0,0,0,1,0,0,0,0,0],
    [0,0,0,1,0,0,0,1,0],
    [0,0,0,1,0,0,0,0,0],
    [1,1,1,1,0,0,0,0,0],
    [0,0,0,0,1,0,0,0,0],
    [0,1,0,0,0,1,0,0,0],
    [1,0,1,0,0,1,1,1,1],
    [0,1,0,0,1,0,0,0,0],
    [0,0,0,1,0,0,0,0,0]])


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
node_img = mark_node(img)

def mark_node(ske):
    buf = np.pad(ske, (1,1), mode='constant')
    nbs = neighbors(buf.shape)
    acc = np.cumprod((1,)+buf.shape[::-1][:-1])[::-1]
    mark(buf, nbs)
    return buf

def mark(img, nbs): # mark the array use (0, 1, 2)
    img = img.ravel()
    for p in range(len(img)):
        if img[p]==0:continue
        s = 0
        for dp in nbs:
            if img[p+dp]!=0:s+=1
        if s==2:img[p]=1
        else:img[p]=2

这一步是将上面的骨架图中的特有节点标识出来（注意：新形成的图是在原图周围附加了一圈0作为缓冲）：
（1）像素值原来为0的地方，仍然为0；
（2）如果某非0像素，其八邻域有2个非0像素，那么标识该像素为1；这种像素代表了骨架图中的中间段的像素（其中有一部分1代表了环形闭合结构，在后面会特殊处理）；
（3）如果某非0像素，其八邻域中的非0像素个数不是2（比如是1、3等），那么标识该像素为2，这种像素代表了骨架图中端点和交点部分的像素。

经过标识后，得到的图像矩阵为：

1
2
3
4
5
6
7
8
9
10
11
[[0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 2 0 0 0 0 0 0]
[0 0 0 0 1 0 0 0 2 0 0]
[0 0 0 0 2 0 0 0 0 0 0]
[0 2 1 2 2 0 0 0 0 0 0]
[0 0 0 0 0 1 0 0 0 0 0]
[0 0 1 0 0 0 2 0 0 0 0]
[0 1 0 1 0 0 2 2 1 2 0]
[0 0 1 0 0 1 0 0 0 0 0]
[0 0 0 0 2 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0]]


1
2
3
4
5
6
7
8
9
10
img = img.ravel()
buf = np.zeros(131072, dtype=np.int64)
num = 10
nodes = []
for p in range(len(img)):
    if img[p] == 2:
        isiso, nds = fill(img, p, num, nbs, acc, buf)
        if isiso and not iso: continue
        num += 1
        nodes.append(nds)

依次观察值为2的像素，对其进行等值填充（主要函数就是fill函数），并记录这些节点在原图中的坐标位置。
fill函数为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def fill(img, p, num, nbs, acc, buf):
    img[p] = num
    buf[0] = p
    cur = 0; s = 1; iso = True;

while True:
        p = buf[cur]
        for dp in nbs:
            cp = p+dp
            if img[cp]==2:
                img[cp] = num
                buf[s] = cp
                s+=1
            if img[cp]==1: iso=False
        cur += 1
        if cur==s:break
    return iso, idx2rc(buf[:s], acc)

原理就是探究这些值为2的像素的邻居是否仍然是2，若是，则将其亦纳入探究范围，这样就标识出了这些节点。
经过上述代码后，得到的nodes数值为：

1
2
3
4
5
[array([[0, 3]], dtype=int16), array([[1, 7]], dtype=int16), array([[2, 3],
       [3, 2],
       [3, 3]], dtype=int16), array([[3, 0]], dtype=int16), array([[5, 5],
       [6, 5],
       [6, 6]], dtype=int16), array([[6, 8]], dtype=int16), array([[8, 3]], dtype=int16)]

如何理解呢？很好理解，比如第一个坐标[0,3]就是第一个值为2的像素在原图中的位置，而坐标组合([[2, 3], [3, 2], [3, 3]])代表那三个相邻的值为2的像素。
同时img中的像素值也发生了变化，比如第一个值为2的像素，它的值由2变成了10（如上代码硬编码），而第一个值为2的像素则变成了11，同理，那三个相邻的值为2的像素变成了12，依次类推，最后一个值为2的像素变成了16。

1
2
3
4
5
6
7
edges = []
for p in range(len(img)):
    if img[p] <10: continue
    for dp in nbs:
        if img[p+dp]==1:
            edge = trace(img, p+dp, nbs, acc, buf)
            edges.append(edge)

像素追溯部分的观察点就变成了与上述节点相邻且值为1的那些像素，即骨架图中的中间部分的像素（原理就是通过是否大于10而选择过滤出上面那些节点，然而通过其邻居是否是1来过滤得到这些值为1的像素），然后通过trace函数寻找其两端所连接的具体节点标识。

trace函数为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
def trace(img, p, nbs, acc, buf):
    c1 = 0; c2 = 0;
    newp = 0
    cur = 1
    while True:
        buf[cur] = p
        img[p] = 0
        cur += 1
        for dp in nbs:
            cp = p + dp
            if img[cp] >= 10:
                if c1==0:
                    c1 = img[cp]
                    buf[0] = cp
                else:
                    c2 = img[cp]
                    buf[cur] = cp
            if img[cp] == 1:
                newp = cp
        p = newp
        if c2!=0:break
    return (c1-10, c2-10, idx2rc(buf[:cur+1], acc))

trace的原理就是：观察这些值为1的像素的邻居，若为大于10，即找到了那些节点nodes，分别通过c1和c2来存储两侧的节点；若为1，则也将其设为进一步的观察点。
经过上述代码后，得到的edges的数值为：

1
2
3
4
5
6
7
8
9
10
11
[(0, 2, array([[0, 3],
       [1, 3],
       [2, 3]], dtype=int16)), (3, 2, array([[3, 0],
       [3, 1],
       [3, 2]], dtype=int16)), (2, 4, array([[3, 3],
       [4, 4],
       [5, 5]], dtype=int16)), (4, 6, array([[6, 5],
       [7, 4],
       [8, 3]], dtype=int16)), (4, 5, array([[6, 6],
       [6, 7],
       [6, 8]], dtype=int16))]

如何理解呢？也很好理解。比如第一个元组((0, 2, array([[0, 3], [1, 3], [2, 3]])，前两个元素0和2分别是img中的值为10和12的像素减去10所得，第三个元素就是第一个值为1的元素的坐标，以及它所连接的两个值为2的节点的坐标。其他元组也是这个意思。通过这样的元组表示，就很自然地为后面的图graph中的首尾节点及中间连接做了准备。

同时img中的像素值又发生了变化：与上面值为2的节点相连接的值为1的像素的值都变为了0（见trace函数中的硬编码），这是为了接下来的闭合曲线的处理。

闭合曲线的处理

上面的代码可以处理非闭合的曲线，因为很容易根据八邻域中的节点数目来获得交点部分所在。而对于闭合曲线，比如原图中左下方的四个1所形成的闭合曲线，其中无法找到值为2的这种像素，且无法对应到图graph中的节点的概念，因此需要特殊处理一下（这里是否处理这种闭合曲线，是用ring这个参数来指定）：

1
2
3
4
5
6
7
8
for p in range(len(img)):
    if img[p]!=1: continue
    img[p] = num; num += 1
    nodes.append(idx2rc([p], acc))
    for dp in nbs:
        if img[p+dp]==1:
            edge = trace(img, p+dp, nbs, acc, buf)
            edges.append(edge)

注意，因为上面已经将值为1且与2相连的像素置为0，所以这里寻找的是剩下的值为1的像素。若为1，然后会将它继续添加到节点nodes中。
然后再对其邻居点进行追溯trace，不断地将邻居的为1的像素加入到edge中，最终得到的edge结果为：

1
2
3
4
5
(7, 7, array([[5, 1],
       [6, 0],
       [7, 1],
       [6, 2],
       [5, 1]], dtype=int16))

代表由7号节点到7号节点的一个循环。

创建图graph

经过上述节点标识和像素追溯，可以得到节点及其边为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
nodes =
[array([[0, 3]], dtype=int16), array([[1, 7]], dtype=int16), array([[2, 3],
       [3, 2],
       [3, 3]], dtype=int16), array([[3, 0]], dtype=int16), array([[5, 5],
       [6, 5],
       [6, 6]], dtype=int16), array([[6, 8]], dtype=int16), array([[8, 3]], dtype=int16), array([[5, 1]], dtype=int16)]

edges =
[(0, 2, array([[0, 3],
       [1, 3],
       [2, 3]], dtype=int16)), (3, 2, array([[3, 0],
       [3, 1],
       [3, 2]], dtype=int16)), (2, 4, array([[3, 3],
       [4, 4],
       [5, 5]], dtype=int16)), (4, 6, array([[6, 5],
       [7, 4],
       [8, 3]], dtype=int16)), (4, 5, array([[6, 6],
       [6, 7],
       [6, 8]], dtype=int16)), (7, 7, array([[5, 1],
       [6, 0],
       [7, 1],
       [6, 2],
       [5, 1]], dtype=int16))]

使用networkx库来构建graph：

1
2
3
4
5
6
7
8
9
import networkx as nx
def build_graph(nodes, edges, multi=False):
    graph = nx.MultiGraph() if multi else nx.Graph()
    for i in range(len(nodes)):
        graph.add_node(i, pts=nodes[i], o=nodes[i].mean(axis=0))
    for s,e,pts in edges:
        l = np.linalg.norm(pts[1:]-pts[:-1], axis=1).sum()
        graph.add_edge(s,e, pts=pts, weight=l)
    return graph



report

Recommend

About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK