问题的提出参见: 《图像工程上册图像处理和分析》(清华) P34
数字图像处理课上讲到图像和视觉基础时涉及到了离散数学中的一个概念: 传递闭包。刚好我的C语言刚看到数组，就编来看看吧。
二值矩阵是一种常用的描述关系的形式。我们可取矩阵中对应相关元素相交处的值为1，其它地方的的值为0。从R={{a, a), (a, b), (b, d), (d, b), (c, e)}}可得如下矩阵B：
a b c d e
┌ 1 1 0 0 0 ┐a
│ 0 0 0 1 0 │b
B =│ 0 0 0 0 1 │c
│ 0 1 0 0 0 │d
└ 0 0 0 0 0 ┘e
传递性指出，如果 a～R～b, b～R～c, 则a～R～c。在以上矩阵中，a与b相连且b与d相连，因为(a,b)和(b,d)在R中。但我们注意到(a,d)不在集合R中。包含这些隐含关系的集合称为R的传递闭包，并记为R+ 。所以R+ = {(a, a), (a, b), (a, d) , (b, b) , (b, d), (d, b), (d, d) , (c, e)}。这个集合包含(a, d), (b, b), (d, d)是根据传递性定义(即由a～R～b, b～R～d, 推出a～R～d; 由b～R～d, d～R～b, 推出b～R～b; 由d～R～b, b～R～d, 推出d～R～d)而来的。描述关系R+的矩阵如下：
a b c d e
┌ 1 1 0 1 0 ┐a
│ 0 1 0 1 0 │b
B+ =│ 0 0 0 0 1 │c
│ 0 1 0 1 0 │d
└ 0 0 0 0 0 ┘e
一般用B表示定义在具有n个元素的字母集A上关系R的n×n二值矩阵，则传递闭包的矩阵B+可如下计算：
B+ = B + BB + BBB + ……+ (B)n (式1)
式中(B)n=BBB……B (n次)。这里矩阵运算时所有乘法都用逻辑与代替，所有加法都用逻辑或代替。上式中的操作次序为B，B(B)，B(BB)，B(BBB)，……，所以在运算的每一步我们只需简单地把现有结果乘以B。
利用上式计算B+要O(n3)阶的与或运算。有一种只需对B中值为1的元素进行或运算的更有效的算法如下：
(1) 置j = 1;
(2) 对i = 1, 2, ……, n, 如果 b(i, j) = 1 ，
则对 k = 1, 2, ……, n, 置 b(i, k) = b(i, k) + b(j, k);
(3) 将j加1;
(4) 如果 j <= n, 转至(2); 否则转到(5);
(5) 结束。此时B转换成了B+。
实际中的二元关系常是等价关系，此时，矩阵B是对称的并且在用B或上述有效算法传递闭包前把所有主对角线的项都设为1。不同符号在矩阵B+的符号集中的等价组可以用从左到右从上到下扫描得到。在第i行第j列扫描到1时，我们将与第j列相应的符号高为与第i行相应的符号(它们等价)，将其它列高为0，并继续扫描B+矩阵。
本来想试着用(式1)所提到的方法实现的，都想了几天了，感觉还是不懂，矩阵的算法还是不太清楚。顺便复习了一下线性代数里学过的矩阵的乘法(我开始写的时候竟然是按对应项相乘算的，汗！)：
┌ a11 a12 a13 ┐
A =│ │
└ a21 a22 a23 ┘
┌ b11 b12 ┐
B =│ b21 b22 │
└ b31 b32 ┘
┌ a11b11+a12b21+a13b31 a11b12+a12b22+a13b32 ┐
AB =│ │
└ a21b11+a22b21+a23b31 a21b12+a22b22+a23b32 ┘
不能在这卡着了，用这种方法还是没能写出来，弄懂之后再补上来吧。
最后还是用书中提到的第二种算法实现了，搜了一下，这种方法应该叫做Warshell算法，好像是图论中求最短路径时提到过。
#include <stdio.h>
int main()
{
int b[5][5] = {{1, 1}, {0, 0, 0, 1}, {0, 0, 0, 0, 1}, {0, 1}};
int i, j, k;
printf("B:\n");
for(i = 0; i < 5; i++) /* 输出原矩阵 */
{
for(j = 0; j < 5; j++)
printf("%2d", b[i][j]);
printf("\n");
}
printf("\nB+:\n");
for(i = 0; i < 5; i++)
for(j = 0; j < 5; j++)
if( b[j][i] == 1 )
for(k = 0; k < 5; k++)
b[j][k] = b[j][k] || b[i][k]; /* b[j][k] = b[j][k] + b[i][k]; */
for(i = 0; i < 5; i++) /* 输出传递闭包矩阵 */
{
for(j = 0; j < 5; j++)
printf("%2d", b[i][j]);
printf("\n");
}
return 0;
}
运行结果:
==========================
B:
1 1 0 0 0
0 0 0 1 0
0 0 0 0 1
0 1 0 0 0
0 0 0 0 0
B+:
1 1 0 1 0
0 1 0 1 0
0 0 0 0 1
0 1 0 1 0
0 0 0 0 0
==========================
还搜到一位仁兄的另一种写法: 查看出处。虽然不太明白其原理，还是改写了一下，如下：
#include <stdio.h>
int main()
{
int b[5][5] = {{1, 1}, {0, 0, 0, 1}, {0, 0, 0, 0, 1}, {0, 1}};
int i, j, k;
for(i = 0; i < 5; i++)
for(j = 0; j < 5; j++)
for(k = 0; k < 5; k++)
if(b[i][j] && (b[i][k] || b[j][k]))
b[i][k] = 1;
for(i = 0; i < 5; i++)
{
for(j = 0; j < 5; j++)
printf("%2d", b[i][j]);
printf("\n");
}
return 0;
}
运行结果：(只输出了结果矩阵)
============================
1 1 0 1 0
0 1 0 1 0
0 0 0 0 1
0 1 0 1 0
0 0 0 0 0
============================
######################################补充########################################
今天(06.4.9)弄懂了矩阵相乘的方法(详见[050] 求两个矩阵的乘积矩阵),这样就可以用第一种方法(式1)求传递闭包了: B+ = B + BB + BBB + ……+ (B)n ,如下:
#include <stdio.h>
int main()
{
int b[5][5] = {{1, 1}, {0, 0, 0, 1}, {0, 0, 0, 0, 1}, {0, 1}};
int b_[5][5];
int i, j, k;
printf("B:\n");
for(i = 0; i < 5; i++) /* 输出原矩阵 */
{
for(j = 0; j < 5; j++)
printf("%2d", b[i][j]);
printf("\n");
}
printf("\n");
printf("B+:\n");
for (i = 0; i < 5; i++) /* 求传递闭包矩阵 */
{
for (j = 0; j < 5; j++)
{
b_[i][j] = 0;
for (k = 0; k < 5; k++)
b_[i][j] = b_[i][j] || b[i][k] && b[k][j]; /* 注释1 */
b_[i][j] = b[i][j] || b_[i][j]; /* 注释2 */
printf("%2d", b_[i][j]);
}
printf("\n");
}
return 0;
}
运行结果:
===============================
B:
1 1 0 0 0
0 0 0 1 0
0 0 0 0 1
0 1 0 0 0
0 0 0 0 0
B+:
1 1 0 1 0
0 1 0 1 0
0 0 0 0 1
0 1 0 1 0
0 0 0 0 0
===============================
★ 注释1 : 此句最终得到的b_[i][j]是矩阵相乘后的对应项的值, 即"式1"中的各项值(B, BB, BBB,……)。
★ 注释2 : 此句最终得到的b_[i][j]是矩阵相乘后的对应项的和, 即"式1"中的各项和(B,B+BB,B+BB+BBB,……)。
评论