关闭

如何检验dmperm的结果

1110人阅读 评论(0) 收藏 举报
最初发布于:2011年7月12日,http://blog.sina.com.cn/xialulee
    实际上在整个《读Matlab7.7 IPT的bwlabel函数》系列中,我们还没有给出用Python和csparse实现的bwlabel的完整代码。虽然在《读Matlab7.7 IPT的bwlabel函数(四)——bwlabel大结局:dmperm现身》中用csparse的dmperm消除错误标签,最终得到了标记图像,但是那一切都是以交互的方式运行的,如果想要一个完整的python版的bwlabel,还得再写写代码。
    由于一直没有给出完整的代码,所以很多观众对于bwlabel的运作方式还存在一些误解。这里就来说明一些主要的疑问。
    1. dmperm的作用是什么?
   《读Matlab7.7 IPT的bwlabel函数(一)——bwlabel1的Python描述》中描述的bwlabel1的算法可知,在bwlabel运行的第一个阶段,bwlabel1这个函数是这样对图像进行标记的。对于如下的图像(x代表非零像素):
x 0
x 0
x x
0 x
x x
x 0
x 0
bwlabel1会从左往右按列扫描(因为Matlab是column-major的系统),每次发现非零像素簇时,就会进行标记,比如扫描完第一列后,会得到:
1 0
1 0
1 x
0 x
2 x
2 0
2 0
然后对第二列进行扫描,一开始发现第二列的那一簇与第一列的标记为1的簇是相邻的,因此将其标记为1:
1 0
1 0
1 1
0 1
2 1
2 0
2 0
但是接着就发现,第二列的这一簇与第一列中标记为2的簇也是相邻的。因此,标记2打错了!标记为2的簇与标记为1的簇是同一簇,不应该被打上不同的标签。
    为了能让后续的代码消除这种错误,bwlabel1的返回值中包含了两个列表rowEquivalences和colEquivalences,用于记录错误的标记。比如上面的例子,标签1和标签2属于同一个非零像素簇,那么此时rowEquivalences = [1]
colEquivalences = [2]
在执行完bwlabel1之后,rowEquivalences和colEquivalences中就记录了所有的错误标签。rowEquivalences[k]与colEquivalences[k]就是属于同一簇,但是打上的不同的标签。
    rowEquivalences和colEquivalences就是消除错误标签的关键数据了。看看下面的这个例子:
rowEquivalences = [1 3 5]
colEquivalences = [2 4 2]
可以看出1和2属于同一簇,3和4属于同一簇,5和2属于同一簇。再仔细观察,可以发现1,2,5是同一簇。于是合并错误标签的问题就变成了一个典型的无向非连通图的connected components的问题。在《读Matlab7.7 IPT的bwlabel函数(三)——如果是xialulee,会怎样》中,我们是用图的遍历来解决的。而Matlab则使用Dulmage-Mendelsohn decomposition来实现相同的功能。
    要使用Dulmage-Mendelsohn decomposition来获取图的connected components,首先需要用矩阵的形式来表达图。表达的方法在图论里好像叫做邻接矩阵。例如:
A[m, n] = 1
表示元素m到n有通路;
A[m, n] = 0
表示元素m到n没有通路。
    对于无向图,如果m到n有通路,则n到m也有通路,所以:
A[n, m] = 1 if A[m, n] = 1
因此无向图的矩阵是一个对称矩阵。又由于每个元素都可以到达自身,因此有:
A[n, n] = 1
因此无向图的矩阵是一个对角线元素为1的对称矩阵。
    Dulmage-Mendelsohn decomposition通过对这个对称矩阵进行行和列的交换,使其变成分块对角阵。则每一个分块对应图的一个connected component。在Matlab中,直接使用Matlab的dmperm函数即可实现Dulmage-Mendelsohn decomposition,如果使用其它的语言和平台,比如C/C++,Python,则可以借助csparse库(参见http://people.sc.fsu.edu/~burkardt/c_src/csparse/csparse.html)的cs_dmperm函数来实现。
    有很多观众在此提出了一个疑问:有时候csparse的cs_dmperm与Matlab的dmperm结果不同!这样还能给出最后的结果吗?要回答这个问题,还是得回顾Dulmage-Mendelsohn decomposition在bwlabel中的作用:找出错误标签的connected components。比如,如果元素1 2 5构成一个connected component,则可以写为:
1 2 5
1 5 2
2 1 5
2 5 1
都是一样的,在这里,元素出现的顺序并不重要,重要的是将属于一个connected component的元素放在一起。所以不要直接比较Matlab的dmperm和csparse的cs_dmperm的结果,而是要比较它们的分组结果
    这里选择一个实际的例子来说明如何检验dmperm的输出。数据来自《VC和matlab 测试dmperm函数的图像》
rowEquivalences = [28    28    13    13    28    13    13    13    28    13    28    28    13    13    13    28    13    28    28    13    28    13    28    28    13    28    13    13    13    20    20    20    20    26    20    13    13    26    20    13    20    20    26    13    20    20    20    20    13    20    13    20    26    20    13    20    20    23    20    23    20    15    13    23    20    20    18    19    20    20    20    19    20    13    10       20
]
colEquivalences = [20    20    28    28    20    28    28    28    20    28    20    20    28    28    28    20    28    20    20    28    20    28    20    20    28    13    20    20    20    26    26    26    23    20    26    20    20    20    26    20    23    23    20    20    23    23    23    26    20    23    20    23    20    23    20    23    24    18    24    18    25    22    20    18    19    19    12    18    19    19    13    18    13    16        4
]
这相当于是Matlab的bwlabel1的输出中包含的错误标签,比如(Matlab下标从1开始)rowEquivalences(1)=28, colEquivalences=20,说明标签为28的像素簇和标签为20的像素簇实际上属于同一个像素簇;rowEquivalences(3)=13, colEquivalences=28,说明标签为13的像素簇和标签为28的像素簇实际上属于同一个像素簇(于是13 28 20都属于同一个像素簇)。我用人工的方式找出了其中包含的所有的connected components:
28 20 13 20 23 24 18 25 19 12 16属于同一个像素簇;
22 15属于同一个像素簇;
5 10 4属于同一个像素簇。
dmperm的作用就是,用计算机来获取这三个connected components。
    但是,在使用dmperm之前,必须将rowEquivalences和colEquivalences表达的错误标签图构造一个对称的,对角线元素为1的稀疏矩阵,这些之前就已经讨论过了。
    Matlab的dmperm的输出结果是:
p =                10                11    12    13    16    18    19    20    23    24    25    26    28    14    15    22    17    21    27    29    30    31    32    33

r = 1                      10    11    12    23    24    26    27    28    29    30    31    32    33    34
    这里p和r就是一串数字,其中蕴含了什么关于connected components的玄机呢?r的作用是用来分隔p(在Dulmage-Mendelsohn decomposition中,p代表了行交换和列交换的次序,r代表了变换之后,每个分块的起始行)。r(k)代表了第k个分组第一个对应元素在p中的下标,r(k+1)-r(k)则等于第k个分组包含的元素个数,比如,看r(4)=4, r(5)=7,说明第4个分组从p的第4个元素开始,r(5)-r(4)=3,说明第4个分组包含了三个元素,因此,第四个分组即是由p(4) p(5) p(6)构成,即:4 5 10。
    再看r(10)=12, r(11)=23,说明第10个分组从p(12)开始,一个包含r(11)-r(10)=11个元素,即p(12)..p(22): 12 13 16 18 19 20 23 24 25 26 28。
    r(12)=24, r(13)=26,说明第12个分组从p(24)开始,到p(25)结束,即:15 22。
    综上,dmperm输出的r和p告诉了我们3个connected components包含的元素:
12 13 16 18 19 20 23 24 25 26 28
4 5 10
15 22
再看看之前我们人工做的分组:
28 20 13 20 23 24 18 25 19 12 16
22 15
5 10 4
虽然元素排列的顺序不一样,但分组是完全一样的。只要分组正确,就可以认为dmperm成功完成了任务。
    接下来再试试csparse的cs_dmperm,这里我使用了之前在《让Python也有dmperm(续)》中给出的代码来做的。但是这里不能直接使用rowEquivalences和colEquivalences构造图的邻接矩阵,因为基于C的系统,下标都是从0开始的,而bwlabel做的标签,是从1开始计数的,所以要将rowEquivalences和colEquivalences的所有元素减一才行:
In [13]: r = array([28,   28,   13,   13,   28,   13,   13,   13,   28,   13,
28,   28,   13,   13,   13,   28,   13,   28,   28,   13,   28,   13,   28,   28,   13,   28,   13,   13,   13,   20,   20,   20,   20,   26,   20,   13,   13,
  26,   20,   13,   20,   20,   26,   13,   20,   20,   20,   20,   13,   20,
13,   20,   26,   20,   13,   20,   20,   23,   20,   23,   20,   15,   13,   23,   20,   20,   18,   19,   20,   20,   20,   19,   20,   13,   10,    5])
In [14]: c = array([20,   20,   28,   28,   20,   28,   28,   28,   20,   28,
20,   20,   28,   28,   28,   20,   28,   20,   20,   28,   20,   28,   20,   20,   28,   13,   20,   20,   20,   26,   26,   26,   23,   20,   26,   20,   20,
  20,   26,   20,   23,   23,   20,   20,   23,   23,   23,   26,   20,   23,
20,   23,   20,   23,   20,   23,   24,   18,   24,   18,   25,   22,   20,   18,   19,   19,   12,   18,   19,   19,   13,   18,   13,   16,    5,    4])

In [15]: r -= 1

In [16]: c -= 1
减一之后,才能构造矩阵(以下代码用于构造对称性):
In [18]: r1 = hstack((r, c))

In [19]: c1 = hstack((c, r))
以下代码用于将对角线元素置1:
In [22]: r1 = hstack((r1, range(33)))

In [23]: c1 = hstack((c1, range(33)))
构造矩阵:
In [29]: A = csc_matrix(([1]*len(r1), (r1, c1)))
使用csparse的cs_dmperm执行分解:
In [30]: from csparse import dmperm

In [31]: p, q, r, s, cc, rr = dmperm(A)
看看结果:
In [48]: print p
[0, 1, 2, 3, 4, 9, 5, 6, 7, 8, 10, 11, 17, 18, 19, 25, 24, 23, 22, 12, 27, 15, 13, 14, 21, 16, 20, 26, 28, 29, 30, 31, 32]

In [49]: print r
[0, 1, 2, 3, 6, 7, 8, 9, 10, 11, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33]
比较Matlab的dmperm输出:
p =                10                11    12    13    16    18    19    20    23    24    25    26    28    14    15    22    17    21    27    29    30    31    32    33

r = 1                      10    11    12    23    24    26    27    28    29    30    31    32    33    34
除了两者在数值上相差1之外,Matlab输出的p和csparse的cs_dmperm输出的p竟然如此不同!但是不必惊慌,我已经不止一次提到,Dulmage-Mendelsohn decomposition的结果不唯一,它的作用就是找出connected components,而一个connected components内部各个元素出现的顺序,则是无关紧要的。那么我们如何检验cs_dmperm的结果是否正确呢?当然是用r对p进行分组,看看各个connected components包含的元素是否一致。对于
p = [0, 1, 2, 3, 4, 9, 5, 6, 7, 8, 10, 11, 17, 18, 19, 25, 24, 23, 22, 12, 27, 15, 13, 14, 21, 16, 20, 26, 28, 29, 30, 31, 32]
r = [0, 1, 2, 3, 6, 7, 8, 9, 10, 11, 22, 23, 25, 26, 27, 28, 29, 30, 31, 32, 33]
还是和前面一样的方法进行操作,只不过这里p的下标是从0开始的:
r[3]=3, r[4]=6, 说明第3个分组对应p[3]..p[5],即3 4 9;
r[9]=11, r[10]=22,说明第9个分组对应p[11]..p[21],即[11, 17, 18, 19, 25, 24, 23, 22, 12, 27, 15];
r[11]=23, r[12]=25,说明第11个分组对应p[23]..p[24],即[14, 21]
即三个connected components对应的元素分别为:
3, 4, 9
11, 17, 18, 19, 25, 24, 23, 22, 12, 27, 15
14, 21
为了方便和Matlab比较,对所有元素加一,得到:
4, 5, 10
12, 18, 19, 20, 26, 25, 24, 23, 13, 28, 16
15, 22
再看看Matlab找出的三个connected components:
12 13 16 18 19 20 23 24 25 26 28
4 5 10
15 22
以及我们人工做的结果:
28 20 13 20 23 24 18 25 19 12 16
22 15
5 10 4
可以看出,虽然各个元素排列的顺序不同,但是每组结果中各个connected components包含的元素是一样的。
    所以,想要检验Dulmage-Mendelsohn decomposition的结果,不要直接比较输出的p和r,而是应该使用r对p进行分组,看看各个connected components包含的元素是否与期望一致,这才是正确的做法。
    从读者们的反馈来看,似乎有很多场合需要用C/C++来实现Matlab的bwlabel的功能。在自己实现一个bwlabel之前,必须将bwlabel算法的各个步骤以及意义完全弄明白才行。xialulee的这个《读bwlabel》系列就是想帮助大家弄明白bwlabel所使用的方法。可惜,虽然bwlabel从本质上说,并不复杂,但是xialulee的表达能力有限(而且有很多饶舌的地方,比如Matlab使用column-major,C和Python使用row-major,C和Python下标从0开始,Matlab下标从1开始,以及在消除错误标签时相关的图的理论,写着写着就一塌糊涂了),所以有很多细节仍然含混不清。因此各位需要实现bwlabel功能的同学,仍然需要狠狠地啃一啃代码了(以及一些数学和算法书籍,因为涉及到图论中常见的connected components的问题)。
    另外,本系列基于Matlab7.7进行讨论的,在后续的版本中,bwlabel的代码发生了一些变化,但基本方法仍然是一样的。
0
0

查看评论
* 以上用户言论只代表其个人观点,不代表CSDN网站的观点或立场
    个人资料
    • 访问:256473次
    • 积分:3407
    • 等级:
    • 排名:第9868名
    • 原创:83篇
    • 转载:43篇
    • 译文:0篇
    • 评论:36条
    最新评论