很是感謝☆Ronny丶博主在其博文《圖像分析:二值圖像連通域標記》中對二值圖像連通域的介紹和算法闡述,讓我這個毫無數據結構算法底子的小白可以理解和復現代碼。本文的目的是基於我本身的理解,對該博文中Two-Pass算法的一些優化和補充,同時也但願幫助更多像我同樣的人較快地掌握連通域標記。html
連通域標記是圖像分割計數的一個重要環節,在工業上應用很是地多。例如像硬幣的計件,在二值化處理後,爲了可以感知數量,就得對硬幣區域進行標記(固然標記前可能還要通過一系列的形態學處理)。另外,還有一個我想到的,更有趣、也更具備挑戰性的例子——二維碼連通域標記,這用來檢驗算法的性能是再合適不過了。言歸正題——本文介紹了兩大流行算法,一個是利用DFS的Seed-Filling算法,另外一個是Two-Pass算法。後者由於處理等價對的方法不一樣,又細分爲DFS Two-Pass(使用DFS處理等價對)和Union-Find Two-Pass(使用並查集處理等價對)。若是硬要給這三種算法排序的話,大概是Union-Find Two-Pass > Seed-Filling > DFS Two-Pass,反正我寫的程序是這樣的速度排序。算法
這個算法其實實質就是DFS,筆者曾經有幸作過一個「水窪連通」的算法題,當時就是用DFS或者BFS來作的,顯然,「水窪連通」也是屬於連通域標記問題的。DFS在這個問題上的思路是:優先地尋找一個完整連通域,在找的同時把他們都標記一下,找完一個完整連通域, 再去找下一個連通域。按照這個想法,程序無非就是維護一個堆棧或者隊列罷了,寫起來相對簡潔易懂。要說缺點的話,就是頻繁的堆棧操做可能會拉低程序的性能。數據結構
簡要地說明一下這部分代碼含義,故事就定義成小明踩水坑吧,雖然小明對我表示本身很文靜,只喜歡作數學題。首先,定義了一個二維矩陣labels,大小跟二值圖同樣。一開始labels都是標籤0,這是一個無效標籤,能夠理解爲是充滿迷霧的未知區域或者是已肯定的非水坑區域。小明每到達一個新的單位域(也就是一個新像素),首先要先看看這個域是否是不曾踩過的水坑(不曾踩過的水坑其標籤爲0且灰度值爲255),若是是的話,那麼小明就原地開心地踩水坑了,踩過之後還不忘給它畫上一個大於0的標記(以標籤1爲例)。接下來,小明回顧四周,又發現了接壤的另外一個水坑, 因而又在該水坑上留下了標記1······這樣看似單調的循環,在小明眼裏倒是一次次奇妙的冒險。愉快的時光很短暫,小明不一下子就發現身邊已經沒有「新鮮」的水坑了,傷心的同時回到最初的那個水坑,繼續朝遠方走去。漸漸地,眼前依稀出現了陌生又熟悉的水坑,重現微笑的小明決定要開啓新的旅途,所以標記1.0進化至2.0。函數
故事的結束,要額外補充一點,程序裏要不停地將新的單位域加入隊列, 所以隊列遍歷其上限是動態的。性能
1 vector<vector<int>> seedFilling(Mat src) 2 { 3 4 // 標籤容器,初始化爲標記0 5 vector<vector<int>> labels(src.rows, vector<int>(src.cols, 0)); 6 // 當前的種子標籤 7 int curLabel = 1; 8 // 四連通位置偏移 9 pair<int, int> offset[4] = {make_pair(0, 1), make_pair(1, 0), make_pair(-1, 0), make_pair(0, -1)}; 10 // 當前連通域中的單位域隊列 11 vector<pair<int, int>> tempList; 12 13 for (int i = 0; i < src.rows; i++) 14 { 15 for (int j = 0; j < src.cols; j++) 16 { 17 // 當前單位域已被標記或者屬於背景區域, 則跳過 18 if (labels[i][j] != 0 || src.at<uchar>(i, j) == 0) 19 { 20 continue; 21 } 22 // 當前單位域未標記而且屬於前景區域, 用種子爲其標記 23 labels[i][j] = curLabel; 24 // 加入單位域隊列 25 tempList.push_back(make_pair(i, j)); 26 27 // 遍歷單位域隊列 28 for (int k = 0; k < tempList.size(); k++) 29 { 30 // 四連通範圍內檢查未標記的前景單位域 31 for (int m = 0; m < 4; m++) 32 { 33 int row = offset[m].first + tempList[k].first; 34 int col = offset[m].second + tempList[k].second; 35 // 防止座標溢出圖像邊界 36 row = (row < 0) ? 0: ((row >= src.rows) ? (src.rows - 1): row); 37 col = (col < 0) ? 0: ((col >= src.cols) ? (src.cols - 1): col); 38 39 // 鄰近單位域未標記而且屬於前景區域, 標記並加入隊列 40 if (labels[row][col] == 0 && src.at<uchar>(row, col) == 255) 41 { 42 labels[row][col] = curLabel; 43 tempList.push_back(make_pair(row, col)); 44 } 45 } 46 } 47 // 一個完整連通域查找完畢,標籤更新 48 curLabel++; 49 // 清空隊列 50 tempList.clear(); 51 } 52 } 53 54 return labels; 55 }
關於Two-Pass的算法原理能夠參考上面提到的博文,原文仍是很詳細的,惟一的遺憾就是後面程序的註釋有點少,看起來會吃力些,說白了就是本身菜。要找一張二維圖像中的連通域,很容易想到能夠一行一行先把子區域找出來,而後再拼合成一個完整的連通域,由於從每一行找連通域是一件很簡單的事。這個過程當中須要記錄每個子區域,爲了知足定位要求,而且節省內存,咱們須要記錄子區域所在的行號、區域開始的位置、結束的位置,固然還有一個表徵子區域總數的變量。須要注意的就是子區域開始位置和結束位置在行首和行末的狀況要單獨拿出來考慮。優化
1 // 查找每一行的子區域 2 // numberOfArea:子區域總數 stArea:子區域開始位置 enArea:子區域結束位置 rowArea:子區域所在行號 3 void searchArea(const Mat src, int &numberOfArea, vector<int> &stArea, vector<int> &enArea, vector<int> &rowArea) 4 { 5 for (int row = 0; row < src.rows; row++) 6 { 7 // 行指針 8 const uchar *rowData = src.ptr<uchar>(row); 9 10 // 判斷行首是不是子區域的開始點 11 if (rowData[0] == 255) 12 { 13 numberOfArea++; 14 stArea.push_back(0); 15 } 16 17 for (int col = 1; col < src.cols; col++) 18 { 19 // 子區域開始位置的判斷:前像素爲背景,當前像素是前景 20 if (rowData[col - 1] == 0 && rowData[col] == 255) 21 { 22 // 在開始位置更新區域總數、開始位置vector 23 numberOfArea++; 24 stArea.push_back(col); 25 // 子區域結束位置的判斷:前像素是前景,當前像素是背景 26 }else if (rowData[col - 1] == 255 && rowData[col] == 0) 27 { 28 // 更新結束位置vector、行號vector 29 enArea.push_back(col - 1); 30 rowArea.push_back(row); 31 } 32 } 33 // 結束位置在行末 34 if (rowData[src.cols - 1] == 255) 35 { 36 enArea.push_back(src.cols - 1); 37 rowArea.push_back(row); 38 } 39 } 40 }
另一個比較棘手的問題,如何給這些子區域標號,使得同一個連通域有相同的標籤值。咱們給單獨每一行的子區域標號區分是很容易的事, 關鍵是處理相鄰行間的子區域關係(怎麼判別兩個子區域是連通的)。spa
主要思路:以四連通爲例,在上圖咱們能夠看出BE是屬於同一個連通域,判斷的依據是E的開始位置小於B的結束位置,而且E的結束地址大於B的開始地址;同理能夠判斷出EC屬於同一個連通域,CF屬於同一個連通域,所以能夠推知BECF都屬於同一個連通域。設計
迭代策略:尋找E的相連區域時,對前一行的ABCD進行迭代,找到相連的有B和C,而D的開始地址已經大於了E的結束地址,此時就能夠提早break掉,避免沒必要要的迭代操做;接下來迭代F的時候,因爲有E留下來的基礎,所以對上一行的迭代能夠直接從C開始。另外,當前行以前的一行若是不存在子區域的話,那麼當前行的全部子區域均可以直接賦新的標籤,而不須要迭代上一行。指針
標籤策略:以上圖爲例,遍歷第一行,A、B、C、D會分別獲得標籤一、二、三、4。到了第二行,檢測到E與B相連,以前E的標籤仍是初始值0,所以給E賦上B的標籤2;以後再次檢測到C和E相連,因爲E已經有了標籤2,而C的標籤爲3,則保持E和C標籤不變,將(2,3)做爲等價對進行保存。同理,檢測到F和C相連,且F標籤仍是初始值0,則爲F標上3。如此對全部的子區域進行標號,最終能夠獲得一個等價對的列表。code
下面的代碼實現了上述的過程。子區域用一維vector保存,沒辦法直接定位到某一行號的子區域,所以須要用curRow來記錄當前的行,用firstAreaPrev記錄前一行的第一個子區域在vector中的位置,用lastAreaPrev記錄前一行的最後一個子區域在vector中的位置。在換行的時候,就去更新剛剛說的3個變量,其中firstAreaPrev的更新依賴於當前行的第一個子區域位置,因此還得用firstAreaCur記錄當前行的第一個子區域。
1 // 初步標籤,獲取等價對 2 // labelOfArea:子區域標籤值, equalLabels:等價標籤對 offset:0爲四連通,1爲8連通 3 void markArea(int numberOfArea, vector<int> stArea, vector<int> enArea, vector<int> rowArea, vector<int> &labelOfArea, vector<pair<int, int>> &equalLabels, int offset) 4 { 5 int label = 1; 6 // 當前所在行 7 int curRow = 0; 8 // 當前行的第一個子區域位置索引 9 int firstAreaCur = 0; 10 // 前一行的第一個子區域位置索引 11 int firstAreaPrev = 0; 12 // 前一行的最後一個子區域位置索引 13 int lastAreaPrev = 0; 14 15 // 初始化標籤都爲0 16 labelOfArea.assign(numberOfArea, 0); 17 18 // 遍歷全部子區域並標記 19 for (int i = 0; i < numberOfArea; i++) 20 { 21 // 行切換時更新狀態變量 22 if (curRow != rowArea[i]) 23 { 24 curRow = rowArea[i]; 25 firstAreaPrev = firstAreaCur; 26 lastAreaPrev = i - 1; 27 firstAreaCur = i; 28 } 29 30 // 相鄰行不存在子區域 31 if (curRow != rowArea[firstAreaPrev] + 1) 32 { 33 labelOfArea[i] = label++; 34 continue; 35 } 36 // 對前一行進行迭代 37 for (int j = firstAreaPrev; j <= lastAreaPrev; j++) 38 { 39 // 判斷是否相連 40 if (stArea[i] <= enArea[j] + offset && enArea[i] >= stArea[j] - offset) 41 { 42 if (labelOfArea[i] == 0) 43 // 以前沒有標記過 44 labelOfArea[i] = labelOfArea[j]; 45 else if (labelOfArea[i] != labelOfArea[j]) 46 // 以前已經被標記,保存等價對 47 equalLabels.push_back(make_pair(labelOfArea[i], labelOfArea[j])); 48 }else if (enArea[i] < stArea[j] - offset) 49 { 50 // 爲當前行下一個子區域縮小上一行的迭代範圍 51 firstAreaPrev = max(firstAreaPrev, j - 1); 52 break; 53 } 54 } 55 // 與上一行不存在相連 56 if (labelOfArea[i] == 0) 57 { 58 labelOfArea[i] = label++; 59 } 60 } 61 }
經過上面的努力,標記任務並無作完,最核心的部分正是如何處理等價對。這裏簡單貼上原博主說的DSF方法,又是深搜啊。相比於直接DFS標記連通域,先找等價對再深搜減小了大量的堆棧操做,由於前者深度取決於連通域的大小,然後者是連通域數量,顯然這兩個數量級的差距挺大的。
原博主的想法是創建一個Bool型等價對矩陣,用做深搜環境。具體作法是先獲取最大的標籤值maxLabel,而後生成一個$maxLabel*maxLabel$大小的二維矩陣,初始值爲false;對於例如(1,3)這樣的等價對,在矩陣的(0,2)和(2,0)處賦值true——要注意索引和標籤值是相差1的。就這樣把全部等價對都反映到矩陣上。
深搜的目的在於創建一個標籤的重映射。例如四、五、8是等價的標籤,都重映射到標籤2。最後重映射的效果就是標籤最小爲1,且依次遞增,沒有缺失和等價。深搜在這裏就是優先地尋找一列等價的標籤,例如一口氣把四、五、8都找出來,而後給他們映射到標籤2。程序也維護了一個隊列,當標籤在矩陣上值爲true,並且沒有被映射過,就加入到隊列。
固然不必定要創建一個二維等價矩陣,通常狀況,等價對數量要比maxLabel來的小,因此也能夠直接對等價對列表進行深搜,但不管採用怎樣的深搜,其等價對處理的性能都不可能提升不少。
1 // 等價對處理,標籤重映射 2 void replaceEqualMark(vector<int> &labelOfArea, vector<pair<int, int>> equalLabels) 3 { 4 int maxLabel = *max_element(labelOfArea.begin(), labelOfArea.end()); 5 // 等價標籤矩陣,值爲true表示這兩個標籤等價 6 vector<vector<bool>> eqTab(maxLabel, vector<bool>(maxLabel, false)); 7 // 將等價對信息轉移到矩陣上 8 vector<pair<int, int>>::iterator labPair; 9 for (labPair = equalLabels.begin(); labPair != equalLabels.end(); labPair++) 10 { 11 eqTab[labPair->first -1][labPair->second -1] = true; 12 eqTab[labPair->second -1][labPair->first -1] = true; 13 } 14 // 標籤映射 15 vector<int> labelMap(maxLabel + 1, 0); 16 // 等價標籤隊列 17 vector<int> tempList; 18 // 當前使用的標籤 19 int curLabel = 1; 20 21 for (int i = 1; i <= maxLabel; i++) 22 { 23 // 若是該標籤已被映射,直接跳過 24 if (labelMap[i] != 0) 25 { 26 continue; 27 } 28 29 30 labelMap[i] = curLabel; 31 tempList.push_back(i); 32 33 for (int j = 0; j < tempList.size(); j++) 34 { 35 // 在全部標籤中尋找與當前標籤等價的標籤 36 for (int k = 1; k <= maxLabel; k++) 37 { 38 // 等價且未訪問 39 if (eqTab[tempList[j] - 1][k - 1] && labelMap[k] == 0) 40 { 41 labelMap[k] = curLabel; 42 tempList.push_back(k); 43 } 44 } 45 } 46 47 curLabel++; 48 tempList.clear(); 49 } 50 // 根據映射修改標籤 51 vector<int>::iterator label; 52 for (label = labelOfArea.begin(); label != labelOfArea.end(); label++) 53 { 54 *label = labelMap[*label]; 55 } 56 57 }
若是讀者看到了這裏,真的要感謝一下您的耐心。Two-Pass算法的代碼要比直接深搜來得多,用很差甚至性能還遠不如深搜。原博主在文中說起了能夠用稀疏矩陣來處理等價對,奈何我較爲愚鈍,讀者能夠自研之。
講到等價對,實質是一種關係分類,於是聯想到並查集。並查集方法在這個問題上顯得很是合適,首先將等價對進行綜合就是合併操做,標籤重映射就是查詢操做(並查集能夠看作一種多對一映射)。並查集具體算法我就不嘮叨了,畢竟不怎麼打程序設計競賽。不過,採用並查集的話,函數定義估計就滿天飛了,這裏我包裝了一下,作成了類——實在是有點強迫症,其中等價對生成的函數方法跟上面的是同樣的。
網上有一些代碼也實現了這個算法,可是有的犧牲了秩優化,合併時讓樹指向較小的根,我的認爲這樣作太不值了。因此爲了解決這個,我在並查集映射後,又用labelReMap來進行第二次的映射,主要的步驟跟前面的差很少。
而後,本身跑了一下這代碼,不算畫圖標記的時間,效率要比上面的快四五倍左右,實時性上確定是綽綽有餘了。
1 class AreaMark 2 { 3 public: 4 AreaMark(const Mat src,int offset); 5 int getMarkedArea(vector<vector<int>> &area); 6 void getMarkedImage(Mat &dst); 7 8 private: 9 Mat src; 10 int offset; 11 int numberOfArea=0; 12 vector<int> labelMap; 13 vector<int> labelRank; 14 vector<int> stArea; 15 vector<int> enArea; 16 vector<int> rowArea; 17 vector<int> labelOfArea; 18 vector<pair<int, int>> equalLabels; 19 20 void markArea(); 21 void searchArea(); 22 void setInit(int n); 23 int findRoot(int label); 24 void unite(int labelA, int labelB); 25 void replaceEqualMark(); 26 }; 27 28 // 構造函數 29 // imageInput:輸入待標記二值圖像 offsetInput:0爲四連通,1爲八連通 30 AreaMark::AreaMark(Mat imageInput,int offsetInput) 31 { 32 src = imageInput; 33 offset = offsetInput; 34 } 35 36 // 獲取顏色標記圖片 37 void AreaMark::getMarkedImage(Mat &dst) 38 { 39 Mat img(src.rows, src.cols, CV_8UC3, CV_RGB(0, 0, 0)); 40 cvtColor(img, dst, CV_RGB2HSV); 41 42 int maxLabel = *max_element(labelOfArea.begin(), labelOfArea.end()); 43 vector<uchar> hue; 44 for (int i = 1; i<= maxLabel; i++) 45 { 46 // 使用HSV模型生成可區分顏色 47 hue.push_back(uchar(180.0 * (i - 1) / (maxLabel + 1))); 48 } 49 50 for (int i = 0; i < numberOfArea; i++) 51 { 52 for (int j = stArea[i]; j <= enArea[i]; j++) 53 { 54 dst.at<Vec3b>(rowArea[i], j)[0] = hue[labelOfArea[i]]; 55 dst.at<Vec3b>(rowArea[i], j)[1] = 255; 56 dst.at<Vec3b>(rowArea[i], j)[2] = 255; 57 } 58 } 59 60 cvtColor(dst, dst, CV_HSV2BGR); 61 } 62 63 // 獲取標記過的各行子區域 64 int AreaMark::getMarkedArea(vector<vector<int>> &area) 65 { 66 searchArea(); 67 markArea(); 68 replaceEqualMark(); 69 area.push_back(rowArea); 70 area.push_back(stArea); 71 area.push_back(enArea); 72 area.push_back(labelOfArea); 73 return numberOfArea; 74 } 75 76 void AreaMark::searchArea() 77 { 78 for (int row = 0; row < src.rows; row++) 79 { 80 // 行指針 81 const uchar *rowData = src.ptr<uchar>(row); 82 83 // 判斷行首是不是子區域的開始點 84 if (rowData[0] == 255) 85 { 86 numberOfArea++; 87 stArea.push_back(0); 88 } 89 90 for (int col = 1; col < src.cols; col++) 91 { 92 // 子區域開始位置的判斷:前像素爲背景,當前像素是前景 93 if (rowData[col - 1] == 0 && rowData[col] == 255) 94 { 95 // 在開始位置更新區域總數、開始位置vector 96 numberOfArea++; 97 stArea.push_back(col); 98 // 子區域結束位置的判斷:前像素是前景,當前像素是背景 99 }else if (rowData[col - 1] == 255 && rowData[col] == 0) 100 { 101 // 更新結束位置vector、行號vector 102 enArea.push_back(col - 1); 103 rowArea.push_back(row); 104 } 105 } 106 // 結束位置在行末 107 if (rowData[src.cols - 1] == 255) 108 { 109 enArea.push_back(src.cols - 1); 110 rowArea.push_back(row); 111 } 112 } 113 } 114 115 116 117 void AreaMark::markArea() 118 { 119 int label = 1; 120 // 當前所在行 121 int curRow = 0; 122 // 當前行的第一個子區域位置索引 123 int firstAreaCur = 0; 124 // 前一行的第一個子區域位置索引 125 int firstAreaPrev = 0; 126 // 前一行的最後一個子區域位置索引 127 int lastAreaPrev = 0; 128 129 // 初始化標籤都爲0 130 labelOfArea.assign(numberOfArea, 0); 131 132 // 遍歷全部子區域並標記 133 for (int i = 0; i < numberOfArea; i++) 134 { 135 // 行切換時更新狀態變量 136 if (curRow != rowArea[i]) 137 { 138 curRow = rowArea[i]; 139 firstAreaPrev = firstAreaCur; 140 lastAreaPrev = i - 1; 141 firstAreaCur = i; 142 } 143 144 // 相鄰行不存在子區域 145 if (curRow != rowArea[firstAreaPrev] + 1) 146 { 147 labelOfArea[i] = label++; 148 continue; 149 } 150 // 對前一行進行迭代 151 for (int j = firstAreaPrev; j <= lastAreaPrev; j++) 152 { 153 // 判斷是否相連 154 if (stArea[i] <= enArea[j] + offset && enArea[i] >= stArea[j] - offset) 155 { 156 if (labelOfArea[i] == 0) 157 // 以前沒有標記過 158 labelOfArea[i] = labelOfArea[j]; 159 else if (labelOfArea[i] != labelOfArea[j]) 160 // 以前已經被標記,保存等價對 161 equalLabels.push_back(make_pair(labelOfArea[i], labelOfArea[j])); 162 }else if (enArea[i] < stArea[j] - offset) 163 { 164 // 爲當前行下一個子區域縮小上一行的迭代範圍 165 firstAreaPrev = max(firstAreaPrev, j - 1); 166 break; 167 } 168 } 169 // 與上一行不存在相連 170 if (labelOfArea[i] == 0) 171 { 172 labelOfArea[i] = label++; 173 } 174 } 175 } 176 177 178 // 並查集初始化 179 void AreaMark::setInit(int n) 180 { 181 for (int i = 0; i <= n; i++) 182 { 183 labelMap.push_back(i); 184 labelRank.push_back(0); 185 } 186 } 187 188 // 查根 189 int AreaMark::findRoot(int label) 190 { 191 if (labelMap[label] == label) 192 { 193 return label; 194 } 195 else 196 { 197 //路徑壓縮優化 198 return labelMap[label] = findRoot(labelMap[label]); 199 } 200 } 201 202 // 合併 203 void AreaMark::unite(int labelA, int labelB) 204 { 205 labelA = findRoot(labelA); 206 labelB = findRoot(labelB); 207 208 if (labelA == labelB) 209 { 210 return; 211 } 212 // 秩優化,秩大的樹合併秩小的樹 213 if (labelRank[labelA] < labelRank[labelB]) 214 { 215 labelMap[labelA] = labelB; 216 } 217 else 218 { 219 labelMap[labelB] = labelA; 220 if (labelRank[labelA] == labelRank[labelB]) 221 { 222 labelRank[labelA]++; 223 } 224 } 225 226 } 227 228 // 等價對處理,標籤重映射 229 void AreaMark::replaceEqualMark() 230 { 231 int maxLabel = *max_element(labelOfArea.begin(), labelOfArea.end()); 232 233 setInit(maxLabel); 234 235 // 合併等價對,標籤初映射 236 vector<pair<int, int>>::iterator labPair; 237 for (labPair = equalLabels.begin(); labPair != equalLabels.end(); labPair++) 238 { 239 unite(labPair->first, labPair->second); 240 } 241 242 // 標籤重映射,填補缺失標籤 243 int newLabel=0; 244 vector<int> labelReMap(maxLabel + 1, 0); 245 vector<int>::iterator old; 246 for (old = labelMap.begin(); old != labelMap.end(); old++) 247 { 248 if (labelReMap[findRoot(*old)] == 0) 249 { 250 labelReMap[findRoot(*old)] = newLabel++; 251 } 252 } 253 // 根據重映射結果修改標籤 254 vector<int>::iterator label; 255 for (label = labelOfArea.begin(); label != labelOfArea.end(); label++) 256 { 257 *label = labelReMap[findRoot(*label)]; 258 } 259 260 }
最後的最後,這些代碼都沒有經歷過「歲月的歷練」,若是存在不合理之處,請讀者指正!