如何检验dmperm的结果

最初发布于：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        5      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        5        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 =  1        2        3        4        5      10        6        7        8        9      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        2        3        4        7        8        9      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 =  1        2        3        4        5      10        6        7        8        9      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        2        3        4        7        8        9      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的代码发生了一些变化，但基本方法仍然是一样的。