C/C++微实践 - 细胞计数

  从一张二值化的动物组织照片中识别并统计出细胞的数量，涉及文件操作以及图的宽度优先遍历算法。

版权声明

本文可以在互联网上自由转载，但必须：注明出处(作者：海洋饼干叔叔)并包含指向本页面的链接。

本文不可以以纸质出版为目的进行改编、摘抄。

1. 细胞计数问题

  识别并统计显微镜下一幅细胞照片中的细胞数量是血液常规检查中最基本的任务，二值化处理可以帮助简化统计。 一幅细胞照片经过二值化处理后可以转化为像素值为0或1的矩阵。 下图展示了这样一个二值化矩阵，在该矩阵中，值为1的元素表示该处是细胞或细胞的一部分，该元素的上、下、左、右的相邻元素如果也是 1，则相邻元素与该元素位于同一个细胞内；矩阵中值为0的元素表示该处无细胞。

  对于上图所示的细胞照片，按上述规则，容易数出照片中包含 7个细胞。 注意，第3行第3列是一个孤立细胞（在图中已用底纹区分） ，它与第2行第5列的细胞并非同一个，因为它位于第2行第4列细胞的左下方，而不是上下左右的位置。

  在本实验对应的实验子目录中有一个名为cellpicture.txt的文本文件，其内容为细胞照片的二值化矩阵。 请编写程序，从该文件中读取矩阵内容并统计该矩阵中的细胞数量。

  该文件的内容如下：

1
2
3
4
5
6
7
8
9
10
11
12
13

12 14
10111000011100
01100110001101
00000111000011
00110000001000
00011000111000
00111100010011
10011101100111
11000100000001
00000000011000
00000000000000
10001100110000
10001000011111

  其中，第1行的12和14以空格分隔，表示该矩阵有12行14列。第2~13行则是12行元素数据，每行有14个值为0或者1的字符。

cellpicture.txt文件下载地址： http://codelearn.club/2022/11/cellcounter/cellpicture.txt

2. 从文件中读取二值化矩阵

  函数readCellPicture()负责从文件cellpicture.txt读取二值化矩阵，并将其存储在二维向量d中。

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

//Project - CellCounter
#include <iostream>
#include <fstream>
#include <vector>
#include <queue>
using namespace std;

bool readCellPicture(vector<vector<int>>& d, int& m, int& n){
    string s;

    ifstream in("D:/C2Cpp/C21_Exception/CellCounter/cellpicture.txt");
    in >> m >> n;
    getline(in,s);

    d.resize(m);
    for (int i=0;i<m;i++){
        getline(in,s);
        if (s.size()!=n)
            return false;
        vector<int>& r = d[i];
        r.resize(n);
        for (int j=0;j<n;j++)
            r[j] = s[j]-'0';
    }

    return true;
}

int main() 
{
    int m=0,n=0;
    vector<vector<int>> d;
    if (!readCellPicture(d,m,n)){
        cout << "Read cellpicture.txt failed." << endl;
        return 0;
    }

    for (int i=0;i<m;i++){
        for (int j=0;j<n;j++)
            cout << d[i][j];
        cout << endl;
    }    
    return 0;
}

第8行：形参d引用了一个类型为vector<vector<int>>的向量，该向量的元素类型为vector<int>，也是向量，子向量内保存整数。形参m,n也是引用，用于“返回”矩阵的行列数。在函数内，d、m、n预期均要被修改。

如果读取正确，函数返回true, 如果出错，返回false。

第11行：D:/C2Cpp/C21_Exception/CellCounter/cellpicture.txt是作者计算机上cellpicture.txt存储的绝对路径，读者需要根据实际情况进行修改。如果读者没有这个文件，其内容可以从本文的第1小节中复制，用记事本编辑保存即可。

第12行：从输入文件流(ifstream)对象in依次读取行数m和列数n。如本书第20章所述，ifstream也是流对象，也支持>>操作符。

第13行：在读完m和n之后，通过一个看似多余的getline()“消化”掉第1行的换行符。该行执行之后，s预期为一个空字符串。

第15行：修改d的尺寸为m，即12行。执行完成后，向量d包含12个子向量。

第16 ~ 24行：循环读入m(12)行的数据，并存储至二维向量d。

第17行：从文件中读取一行内容至s。

第18行 ~ 19行：如果单行字符串的长度不等于期望的列数n(14)，返回false表示读取失败。

第20行：获取向量d内第i个子向量的引用。

第21行：将子向量的长度改为n。

第22 ~ 23行：逐一向s内的每个字符转换成整数0或者1，存入d[i]。

第26行：返回true表示操作成功。

第33 ~ 36行：调用readCellPicture()函数读取二值化矩阵，如果失败，报错并返回退出。

第38 ~ 42行：逐行逐列打印输出二维向量d的全部元素。

3. 图的宽度优先遍历

  细胞的计数依赖于对矩阵元素的遍历。 在逐行逐列遍历矩阵元素的过程中，发现一个值为1的元素，则意味着发现了一个细胞。 此时，需要从该元素出发，向多个方向搜索，找出该元素归属细胞所包含的全部像素/元素，并将这些像素/元素标记为“已探索” ，以避免在后续遍历过程中，这些像素被错误地认为属于一个“新细胞” 。

  这个任务可以通过一个称为“图的宽度优先遍历”的算法来解决，该算法可用函数explorePixel( )来描述。其中，d为代表二值矩阵的嵌套向量，该矩阵有m行n列，搜索出发点为i行j列。

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

void explorePixel(vector<vector<int>>& d,const int m, const int n,
                  const int i, const int j)
{
    queue<pair<int,int>> q;  //先进先出队列q，元素类型为包含两个整数的值对(pair)
    q.push(std::make_pair(i,j)); //std::make_pair()函数将坐标i,j封装成一个pair

    while (!q.empty()){
        auto& e = q.front(); //q.front()返回队首元素的引用
        int x = e.first, y = e.second;  //e的类型为pair<int,int>&, e.first, e.second对应x,y座标
        q.pop();   //弹出队列q的首元素
        if (d[x][y]<0)
            continue;
        d[x][y] = 0 - d[x][y];
        if (x>0 && d[x-1][y]>0)
            q.push(make_pair(x-1,y));  //在队列末尾添加一个pair
        if (x<(m-1) && d[x+1][y]>0)
            q.push(make_pair(x+1,y));
        if (y>0 && d[x][y-1]>0)
            q.push(make_pair(x,y-1));
        if (y<(n-1) && d[x][y+1]>0)
            q.push(make_pair(x,y+1));
    }

    return;
}

【工作原理说明】

• explorePixel( )函数的任务可以描述为：从位于i行j列的像素出发，向其周边进行搜索，将与该像素有连接关系的其他像素全部找出，并标记为“已探索”。

• 当像素值=0时，该像素不属于细胞；像素值=1时，该像素属于细胞且“待探索”，像素值=-1时，该像素属于细胞且“已探索”。

• q为一个先进先出队列，其中存储那些已发现，但其自身及其相邻像素尚未被探索的像素。最初，q仅包含像素(i, j)。然后，程序将一直循环，直至队列空为止，每次循环的过程如下：

    1) 从队列中取出第0个待探索的像素(x, y)；

    2) 检查(x,y)是否为已探索，如是，continue进入下一轮循环；

    3) 将(x,y)标记为“已探索”，即将其值变更为 0 – d[x][y]；

    4) 按上、下、左、右的顺序检查(x,y)的相邻像素，如果相邻像素的值>0，说明该相邻像素也属于当前细胞，将其加入待探索像素队列q。

• 函数中的x表示行坐标，y表示列坐标。即与平面坐标系的通常习惯有所不同，这里的x表示上下方向，y表示左右方向。

  为帮助读者理解上述“图的宽度优先遍历”算法，我们手工模拟一遍explorePixel( )函数的执行过程，从像素(0,4)出发。请见下图，为了与程序一致，我们的行号、列号改为从下标0开始。

  当程序遍历矩阵元素到第0行第4列时，发现像素(0,4)值为1，即该像素状态为待探索且属于一个新细胞，执行explorePixel(d,m,n,0,4)从像素(0,4)出发搜索与该像素连接的全部细胞像素：

1. 队列q被初始化为只包含像素(0,4)，其值为：[ (0,4) ]。
2. 像素(0,4)出队列，其值等于1待探索，像素(0,4)被赋值-1。此时，q = 空列表。
3. 像素(0,4)位于第1行，其上方元素不存在。
4. 像素(0,4)的下方元素(1,4)的值为1，属于同一细胞，将(1,4)加入队列q。此时，q = [ (1,4) ]。
5. 像素(0,4)的左方和右方元素值为0，不属于细胞。
6. 像素(1,4)出队列，其值等于1待探索，赋值-1。此时，q = 空列表。
7. 像素(1,4)的上方元素(0,4)此时值为-1已探索。
8. 像素(1,4)的下方元素(2,4)值为1，属于同一细胞，加入队列q。此时，q = [ (2,4) ]。
9. 像素(1,4)的左方和右方元素值均为1，属于同一细胞，将左方元素(1,3)和右方元素(1,5)加入队列q。此时，q = [ (2,4), (1,3), (1,5) ]。
10. 像素(2,4)出队列，其值为1待探索，赋值-1。此时，q = [ (1,3), (1,5) ]。
11. 像素(2,4)的上方元素已探索，下方、左方、右方元素为0，故未发现新的待探索像素。
12. 像素(1,3)出队列，其值为1待探索，赋值-1；(1,3)的右方元素已探索 ，上方、左方、下方均为0，故未发现新的待探索像素。注意，(2,2)位于(1,3)的左下方，并不属于题目定义的相邻像素。此时，q = [ (1,5) ]。
13. 像素(1,5)出队列，其值为待探索，赋值-1；(1,5)的左方元素已探索，其它三个方向均为0，未发现新的待探索元素。
14. 此时，q = 空列表，循环结束。在循环过程中，(0,4)、(1,4)、(2,4)、(1,3)、(1,5)共5个像素被探索并标记为已探索，这5个像素构成了一个完整的细胞。

1
2

if (d[x][y]<0)
    continue;

  读者可能会上述代码感到疑惑：既然一个待探索元素在加入队列前其值确定为1即待探索，那么在该像素出队后，其值可能变为负数（已探索）吗？

  我们考虑上图所展示的情况。当我们从(1,1)出发搜索该细胞的全部像素时，(2,2)作为像素(2,1)的相邻元素，在探索(2,1)时会被加入队列。同时，作为(1,2)的下方元素，在探索(1,2)时也会被加入队列。这样，队列中就会存在两个(2,2)，当第1个(2,2)被取出时，其像素值为1未探索，而当第2个(2,2)被取出时，其像素值已经是-1已探索。此时，再去考虑(2,2)的相邻元素已不具实践意义，故略过。

4. 循环遍历与搜索

  在实现了函数explorePixel( )之后，下述程序逐行逐列地遍历全部矩阵元素，如果发现值>0的像素，说明遇到了一个属于细胞且“待探索”的元素，将细胞计数变量iCellCount加1，然后再调用explorePixel( )函数从该像素出发探索整个细胞。由于explorePixel( )函数会将所有探索过的细胞像素置为已探索，所以当下述程序第11行、第12行的主循环遍历到已探索的细胞像素时，该像素值为-1，不会将其视为一个新细胞。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

int main()
{
    int m=0,n=0;
    vector<vector<int>> d;
    if (!readCellPicture(d,m,n)){
        cout << "Read cellpicture.txt failed." << endl;
        return 0;
    }

    int iCellCount = 0;
    for (int i=0;i<m;i++)
        for (int j=0;j<n;j++){
            if (d[i][j]<=0)
                continue;
            iCellCount+=1;
            explorePixel(d,m,n,i,j);
        }

    cout << "Found " << iCellCount << " cells.";
    return 0;
}

5. 完整源代码

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
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

//Project - CellCounter
#include <iostream>
#include <fstream>
#include <vector>
#include <queue>
using namespace std;

bool readCellPicture(vector<vector<int>>& d, int& m, int& n){
    string s;

    ifstream in("D:/C2Cpp/C21_Exception/CellCounter/cellpicture.txt");
    in >> m >> n;
    getline(in,s);

    d.resize(m);
    for (int i=0;i<m;i++){
        getline(in,s);
        if (s.size()!=n)
            return false;
        vector<int>& r = d[i];
        r.resize(n);
        for (int j=0;j<n;j++)
            r[j] = s[j]-'0';
    }

    return true;
}

void explorePixel(vector<vector<int>>& d,const int m, const int n,const int i, const int j)
{
    queue<pair<int,int>> q;
    q.push(std::make_pair(i,j));

    while (!q.empty()){
        auto& e = q.front();
        int x = e.first, y = e.second;
        q.pop();
        if (d[x][y]<0)
            continue;
        d[x][y] = 0 - d[x][y];
        if (x>0 && d[x-1][y]>0)
            q.push(make_pair(x-1,y));
        if (x<(m-1) && d[x+1][y]>0)
            q.push(make_pair(x+1,y));
        if (y>0 && d[x][y-1]>0)
            q.push(make_pair(x,y-1));
        if (y<(n-1) && d[x][y+1]>0)
            q.push(make_pair(x,y+1));
    }

    return;
}

int main()
{
    int m=0,n=0;
    vector<vector<int>> d;
    if (!readCellPicture(d,m,n)){
        cout << "Read cellpicture.txt failed." << endl;
        return 0;
    }

    int iCellCount = 0;
    for (int i=0;i<m;i++)
        for (int j=0;j<n;j++){
            if (d[i][j]<=0)
                continue;
            iCellCount+=1;
            explorePixel(d,m,n,i,j);
        }

    cout << "Found " << iCellCount << " cells.";
    return 0;
}