【影象處理筆記】影象分割之形態學分水嶺
2022-10-17 15:01:05
0 引言
迄今為止，我們討論了基於三個主要概念的分割：邊緣檢測、閾值處理和區域提取。每種方法都有優點[例如全域性閾值處理具有速度優勢]和缺點[例如在基於邊緣的分割中，需要進行後處理(如邊緣連線)]。本節討論的基於形態學分水嶺概念的方法。分水嶺分割體現了其他三種方法的許多概念，因此往往會產生更穩定的分割結果，包括連通的分割邊界。
1 原理
分水嶺方法是一種基於拓撲理論的數學形態學的分割方法，基本思想是把影象看作測地學上的拓撲地貌，將畫素點的灰度值視為海拔高度，整個影象就像一張高低起伏的地形圖。每個區域性極小值及其影響區域稱為集水盆，集水盆的邊界則形成分水嶺。演演算法的實現過程可以理解為洪水淹沒的過程：最低點首先被淹沒，然後水逐漸淹沒整個山谷；水位升高到一定高度就會溢位，於是在溢位位置修建堤壩；不斷提高水位，重複上述過程，直到所有的點全部被淹沒；所建立的一系列堤壩就成為分隔各個盆地的分水嶺。
分水嶺的計算過程是一個迭代標註過程，通過尋找集水盆和分水嶺對影象進行分割。經典的分水嶺演演算法分為排序過程和淹沒過程兩個步驟，首先對每個畫素的灰度級從低到高排序，然後在從低到高的淹沒過程中，對每一個區域性極小值在 h 階高度的影響域進行判斷及標註。
2 cv::watershed的使用
2.1 相關函數
OpenCV提供了函數watershed實現基於標記的分水嶺演演算法。輸入影象一般是原圖，或者梯度圖，由於噪聲和梯度的區域性不規則性會導致過度分割，控制過度分割的一種方法依據是標記。標記影象為CV_32S 型別，我們預先把一些區域標註好，影象中每個非零畫素代表一個標籤。對影象中部分畫素做標記，表明它的所屬區域是已知的。這些標註了的區域稱之為種子點。watershed 演演算法會把這些標記的區域慢慢膨脹填充整個影象。
void watershed( InputArray image, // 輸入8位元3通道影象。 InputOutputArray markers //輸入/輸出標記的 32 位單通道影象（地圖）。 與 image 具有相同的大小。 );
OpenCV中，函數distanceTransform()用於計算影象中每一個非零點畫素與其最近的零點畫素之間的距離，輸出的是儲存每一個非零點與最近零點的距離資訊。影象上越亮的點，代表了離零點的距離越遠。可以根據距離變換的這個性質，經過簡單的運算，用於細化字元的輪廓和查詢物體質心（中心）。
void distanceTransform( InputArray src, //輸入影象，8-bit 單通道灰度影象 OutputArray dst, //儲存了每一個點與最近的零點的距離資訊，點越亮離零點越遠。 int distanceType,//距離的型別 int maskSize, //距離變換遮罩的大小，通常取 3, 5 int dstType=CV_32F); //輸出影象型別
2.2 標記的生成
標記由3部分組成：確定前景、確定背景和未知區域。分為兩步：
（1）通過 閾值分割+形態學【例1】，距離變換【例2】，滑鼠互動等方法得到這三個區域。
（2）通過 connectedComponents【例1】或 drawContours【例2】標記區域，其中未知區域用0標記，其他區域用大於0的整數標記。
例1 閾值分割+形態學得到三區域→connectedComponents標記區域→分水嶺
#include <opencv2/opencv.hpp>
#include <iostream>
using namespace cv;
using namespace std;
int main() {
Mat src;
Mat img = imread("./16.tif");
cvtColor(img, src, COLOR_BGR2GRAY);
imshow("src", src);
// 閾值處理
Mat thresh;
threshold(src, thresh, 0, 255, THRESH_OTSU);
// 生成確定背景區域
Mat background;
Mat ele = getStructuringElement(MORPH_RECT, Size(3, 3));
dilate(thresh, background, ele, cv::Point(-1, -1), 2);
bitwise_not(background, background);
// 生成確定前景區域，並利用連通區域標記
Mat foreground;
morphologyEx(thresh, foreground, MORPH_OPEN, ele, cv::Point(-1, -1), 2);
int n = connectedComponents(foreground, foreground, 8, CV_32S);// 此時確定前景大於0，其餘為0
// 生成標記圖
Mat markers = foreground;
markers.setTo(255, background);// 將確定背景設為255，其餘為0的不動，即為unkown
Mat markers8u;
markers.convertTo(markers8u, CV_8U, 10);//灰度值*10使得差異變得明顯
imshow("Markers（輸入）", markers8u);
// 分水嶺演演算法標註目標的輪廓
watershed(img, markers);// 輪廓由-1表示
markers.convertTo(markers8u, CV_8U, 10);//灰度值*10使得差異變得明顯
imshow("Markers（輸出）", markers8u);
// 後處理（顏色填充）
Mat mark;
markers.convertTo(mark, CV_8U);//轉換後-1變成0
bitwise_not(mark, mark);
vector<Vec3b> colors;
for (size_t i = 0; i < n; i++)
{
int b = theRNG().uniform(0, 256);
int g = theRNG().uniform(0, 256);
int r = theRNG().uniform(0, 256);
colors.push_back(Vec3b((uchar)b, (uchar)g, (uchar)r));
}
Mat dst = Mat::zeros(markers.size(), CV_8UC3);
for (int i = 0; i < markers.rows; i++)
{
for (int j = 0; j < markers.cols; j++)
{
int index = markers.at<int>(i, j);
if (index > 0 && index <= static_cast<int>(n))
dst.at<Vec3b>(i, j) = colors[index - 1];
}
}
imshow("dst", dst);
waitKey(0);
return 0;
}
例2 梯度+距離變換得到三區域→drawContours標記區域→分水嶺
分水嶺分割的主要應用之一是，從背景中提取出接近一致的（團狀）目標。由變化較小的灰度表徵區域有較小的梯度之。因此，在實踐中，我們常常看到分水嶺分割應用到梯度影象，而不是應用到影象本身。在這一表述中，匯水盆地的區域極小值與對應感興趣目標的梯度的極小值密切相關。
#include <opencv2/opencv.hpp>
using namespace cv;
using namespace std;
int main() {
Mat src_bgr = imread("./18.png");
imshow("src", src_bgr);
Mat mask, src_gray;
cvtColor(src_bgr, src_gray, COLOR_BGR2GRAY);
// 1.梯度+距離變換得到前景區域
Mat blurImg;
blur(src_gray, blurImg, Size(7, 7));
Mat gx, gy, gxy;
Sobel(blurImg, gx, CV_16S, 1, 0, 3);
Sobel(blurImg, gy, CV_16S, 0, 1, 3);
convertScaleAbs(gx, gx);
convertScaleAbs(gy, gy);
gxy = gx + gy;
Mat binImg;
threshold(gxy, binImg, 110, 255, THRESH_BINARY);
Mat ele = getStructuringElement(MORPH_RECT, Size(3, 3));
morphologyEx(binImg, binImg, MORPH_OPEN, ele);
morphologyEx(binImg, binImg, MORPH_CLOSE, ele);
floodFill(binImg, Point(0, 0), Scalar(255));
bitwise_not(binImg, binImg);
Mat dist, binImg1;
distanceTransform(binImg, dist, DIST_L2, 3);
normalize(dist, dist, 0, 1.0, NORM_MINMAX);
imshow("dist", dist);
threshold(dist, binImg1, 0.4, 1.0, THRESH_BINARY);
binImg1.convertTo(binImg1, CV_8U, 255);
imshow("binImg1", binImg1);
// 2.確定背景區域
Mat markers = Mat::zeros(binImg1.size(), CV_32S);
circle(markers, Point(5, 5), 3, Scalar(255), -1);
inRange(src_gray, Scalar(0), Scalar(15), mask);
// 3.生成標記圖
vector<vector<Point>> contours;
findContours(binImg1, contours, RETR_EXTERNAL, CHAIN_APPROX_SIMPLE);
for (size_t i = 0; i < contours.size(); i++)
{
drawContours(markers, contours, static_cast<int>(i), Scalar(static_cast<int>(i) + 1), -1);
}
markers.setTo(255, mask);
Mat markers8u;
markers.convertTo(markers8u, CV_8U, 10);
imshow("Markers1", markers8u);
// 4.分水嶺
watershed(src_bgr, markers);
// 5.後處理
Mat mark;
markers.convertTo(mark, CV_8U);
bitwise_not(mark, mark);
markers.convertTo(markers8u, CV_8U, 10);
imshow("Markers2", markers8u);
vector<Vec3b> colors;
for (size_t i = 0; i < contours.size(); i++)
{
int b = theRNG().uniform(0, 256);
int g = theRNG().uniform(0, 256);
int r = theRNG().uniform(0, 256);
colors.push_back(Vec3b((uchar)b, (uchar)g, (uchar)r));
}
Mat dst = Mat::zeros(markers.size(), CV_8UC3);
for (int i = 0; i < markers.rows; i++)
{
for (int j = 0; j < markers.cols; j++)
{
int index = markers.at<int>(i, j);
if (index > 0 && index <= static_cast<int>(contours.size()))
dst.at<Vec3b>(i, j) = colors[index - 1];
}
}
imshow("dst", dst);
waitKey(0);
return 0;
}
3 watershed原始碼
3.1 大致流程
1. 對輸入的標記圖mask的邊界（1畫素）置為-1，即邊界不處理
2. 初始階段：遍歷mask中所有點，將標記為0且四鄰域中有標記大於0的點（也就是標記點外面的一圈點）進隊。有256個佇列，進哪個佇列呢？分兩種情況：（1）四鄰域中只有一個標記大於0的點，計算該點和這個鄰域點的梯度（差值），若梯度為250，就進第250個佇列。由於輸入的原圖是三通道的，算梯度要算三個通道的梯度，選最大的；（2）四鄰域中不止一個標記大於0的點，選擇梯度最小的點進隊。
3. 經過上一步，得到256個佇列，有些佇列是空的，有些不是。定位到第一個非空的佇列，遍歷佇列中的點。比較這個點的四鄰域，分三種情況：（1）四鄰域中有一個以上的鄰域標記大於0（且它們標記不同），該點為分水嶺；（2）四鄰域中只有一個鄰域標記大於0，那該點和這個鄰域點標記相同，即被這個鄰域點擴充（3）鄰域中沒有大於0的標記點，這是不可能的。因為標記點要麼是初始種子點，要麼是初始階段延伸的種子點的鄰接點，即該點一定存在一個鄰接點是標記點
4. 在上一步中，情況（1）的分水嶺點被pop出佇列，且在mask中被標記為-1；情況（2）的點為新標記點，被pop出佇列，同時檢視其四鄰域點，有沒有沒處理的點。沒處理的點被push進佇列（相當於又往外一圈），進哪個佇列呢？和上面一樣，根據梯度t。再次進入第3步，直到所有佇列為空。有所不同的是，初始階段從第0個佇列開始找非空佇列，而之後從min(t, active_queue)開始。也就是說，如果第二圈的點梯度都很大，就從第active_queue開始找非空佇列。如果第二圈的梯度都很小，那就從最小的梯隊開始找非空佇列。
3.2 原始碼註解
註解來自OpenCV 原始碼中分水嶺演演算法 watershed 函數原始碼註解，對比OpenCV 2.4.9，OpenCV 4.5中Mat取代了CvMat，cv::Mat是一個類（Class），而CvMat是一個Struct，獲取元素的寫法不同，如下：
Mat test = (Mat_<int>(3, 3) << 1, 2, 3, 4, 5, 6, 7, 8, 9 );
int* p = test.ptr<int>();//獲取元素的寫法不同，之前是test->data.i
int istep = test.step/sizeof(int);// 步長 = 一行位元組數 / sizeof(畫素資料型別)
for (int i = 0; i < test.rows; i++) {// 逐行掃描
for (int j = 0; j < test.rows; j++) {// 逐列掃描
if (p[j] == 5) {
//四鄰域
cout << p[j-1] << endl;//左
cout << p[j+1] << endl;//右
cout << p[j-istep] << endl;//上
cout << p[j+istep] << endl;//下
}
}
p += istep;
}
原始碼註解
#include "precomp.hpp"
// 結點，用於儲存原始圖img中畫素的偏移量和輸出圖mask中畫素的偏移量
typedef struct CvWSNode
{
struct CvWSNode* next;
int mask_ofs;
int img_ofs;
}
CvWSNode;
// 佇列，用於儲存結點 CvWSNode
typedef struct CvWSQueue
{
CvWSNode* first;
CvWSNode* last;
}
CvWSQueue;
// 分配空間
static CvWSNode*
icvAllocWSNodes( CvMemStorage* storage )
{
CvWSNode* n = 0;
int i, count = (storage->block_size - sizeof(CvMemBlock))/sizeof(*n) - 1;
n = (CvWSNode*)cvMemStorageAlloc( storage, count*sizeof(*n) );
for( i = 0; i < count-1; i++ )
n[i].next = n + i + 1;
n[count-1].next = 0;
return n;
}
CV_IMPL void
cvWatershed( const CvArr* srcarr, CvArr* dstarr )
{
const int IN_QUEUE = -2; // 加入到佇列q中的點定義為 -2
const int WSHED = -1; // 「分水嶺」在mask中定義為 -1
const int NQ = 256; // 佇列的數量 256，其實是對應灰度的數量
cv::Ptr<CvMemStorage> storage;
CvMat sstub, *src;
CvMat dstub, *dst;
CvSize size;
CvWSNode* free_node = 0, *node;
CvWSQueue q[NQ]; // 長度為256的CvWSQueue陣列，注意陣列中每個元素都是一個佇列，佇列中每個元素是一個節點
int active_queue; // 指明當前處理的佇列，q[active_queue]
int i, j;
int db, dg, dr;
int* mask; // 指向標記影象的指標
uchar* img; // 指向原始影象的指標
int mstep, istep; // mstep是mask對應的一行畫素數（不是位元組數），istep是img對應的一行畫素數
int subs_tab[513];
// MAX(a,b) = b + MAX(a-b,0) 取最大值
#define ws_max(a,b) ((b) + subs_tab[(a)-(b)+NQ])
// MIN(a,b) = a - MAX(a-b,0) 取最小值
#define ws_min(a,b) ((a) - subs_tab[(a)-(b)+NQ])
// 進隊操作
#define ws_push(idx,mofs,iofs) \
{ \
if( !free_node ) \
free_node = icvAllocWSNodes( storage );\
node = free_node; \
free_node = free_node->next;\
node->next = 0; \
node->mask_ofs = mofs; \
node->img_ofs = iofs; \
if( q[idx].last ) \
q[idx].last->next=node; \
else \
q[idx].first = node; \
q[idx].last = node; \
}
// 出隊操作
#define ws_pop(idx,mofs,iofs) \
{ \
node = q[idx].first; \
q[idx].first = node->next; \
if( !node->next ) \
q[idx].last = 0; \
node->next = free_node; \
free_node = node; \
mofs = node->mask_ofs; \
iofs = node->img_ofs; \
}
// 求出 ptr1 和 ptr2 指向的畫素 r,g,b 差值的最大值
#define c_diff(ptr1,ptr2,diff) \
{ \
db = abs((ptr1)[0] - (ptr2)[0]);\
dg = abs((ptr1)[1] - (ptr2)[1]);\
dr = abs((ptr1)[2] - (ptr2)[2]);\
diff = ws_max(db,dg); \
diff = ws_max(diff,dr); \
assert( 0 <= diff && diff <= 255 ); \
}
src = cvGetMat( srcarr, &sstub );
dst = cvGetMat( dstarr, &dstub );
// 對引數做檢查，要求影象src的型別是8UC3,dst的型別是32SC1，src和dst size相同
if( CV_MAT_TYPE(src->type) != CV_8UC3 )
CV_Error( CV_StsUnsupportedFormat, "Only 8-bit, 3-channel input images are supported" );
if( CV_MAT_TYPE(dst->type) != CV_32SC1 )
CV_Error( CV_StsUnsupportedFormat,
"Only 32-bit, 1-channel output images are supported" );
if( !CV_ARE_SIZES_EQ( src, dst ))
CV_Error( CV_StsUnmatchedSizes, "The input and output images must have the same size" );
size = cvGetMatSize(src); // 獲取影象的size
storage = cvCreateMemStorage();
// 步長 = 一行位元組數 / sizeof(畫素資料型別)
istep = src->step; // img是uchar型， sizeof(uchar) = 1，所以忽略除數
img = src->data.ptr; // 獲取 uchar型別指標
mstep = dst->step / sizeof(mask[0]); // mask是int（32SC1）型，sizeof(mask[0]) = 4
mask = dst->data.i; // 獲取 int型別指標
memset( q, 0, NQ*sizeof(q[0]) ); // 初始化佇列q
for( i = 0; i < 256; i++ )
subs_tab[i] = 0;
for( i = 256; i <= 512; i++ )
subs_tab[i] = i - 256;
// draw a pixel-wide border of dummy "watershed" (i.e. boundary) pixels
// 把影象四個邊的畫素畫成分水嶺
// mask的首行和末行畫成分水嶺
for( j = 0; j < size.width; j++ )
mask[j] = mask[j + mstep*(size.height-1)] = WSHED;
// initial phase: put all the neighbor pixels of each marker to the ordered queue -
// determine the initial boundaries of the basins
// 初始階段：把每個標記的所有鄰居畫素放到有序佇列中去，以確定聚水盆的初始邊界
// 即每個標記（種子，全為正值，1,2,3...）都是一個初始聚水盆，標記的周圍一圈的鄰居畫素就是聚水盆的初始邊界
// 這裡用的是一種逆向思維，不是找標記點，而是判斷每一個點是否為標記點的鄰居,若是，則該點也被擴充為與標記點同型別的標記點
// 若是多個標記點的鄰居，選擇梯度最小的標記點的型別，作為該點的標記點型別
for( i = 1; i < size.height-1; i++ )
{
img += istep; mask += mstep; // 逐行掃描
mask[0] = mask[size.width-1] = WSHED; // 每一行的首列和末列畫成分水嶺，加上前面的首行和末行，mask被分水嶺方框圍起來
for( j = 1; j < size.width-1; j++ ) // 逐列
{
int* m = mask + j; // mask的每個畫素
if( m[0] < 0 ) m[0] = 0; // 該點若為負值，先置為零（初始狀態下除了四邊是分水嶺(-1)其餘點不應該存在負值？）
if( m[0] == 0 && (m[-1] > 0 || m[1] > 0 || m[-mstep] > 0 || m[mstep] > 0) ) // 若該點為非標記點（0），且四鄰域存在標記點（>0）
{
// 求出原圖中該點到有標記點的四鄰域中，梯度值最小（idx）方向的點，將該點和對應的最小梯度值放入q[idex]佇列中
// 兩個畫素的r,g,b 三個通道中相差最大的值作為畫素間的梯度值
uchar* ptr = img + j*3;
int idx = 256, t;
if( m[-1] > 0 )
c_diff( ptr, ptr - 3, idx );
if( m[1] > 0 )
{
c_diff( ptr, ptr + 3, t );
idx = ws_min( idx, t );
}
if( m[-mstep] > 0 )
{
c_diff( ptr, ptr - istep, t );
idx = ws_min( idx, t );
}
if( m[mstep] > 0 )
{
c_diff( ptr, ptr + istep, t );
idx = ws_min( idx, t );
}
assert( 0 <= idx && idx <= 255 );
ws_push( idx, i*mstep + j, i*istep + j*3 ); // 將該點在img和mask中的座標（一維表示）儲存在q[idx]佇列中
m[0] = IN_QUEUE; // 在mask中標記該點已入隊
}
}
}
// find the first non-empty queue
// 定位到第一個非空的佇列
for( i = 0; i < NQ; i++ )
if( q[i].first )
break;
// if there is no markers, exit immediately
// 若i=256，說明陣列q中所有佇列為空
if( i == NQ )
return;
active_queue = i;
img = src->data.ptr;
mask = dst->data.i;
// recursively fill the basins
// 遞迴地填滿聚水盆
for(;;)
{
int mofs, iofs; // 將二維影象線性化後影象畫素的座標 mask_offset 和 img_offset 的縮寫
int lab = 0, t;
int* m;
uchar* ptr;
// 如果這個灰度上的佇列處理完了，就繼續找下一個非空佇列
if( q[active_queue].first == 0 )
{
for( i = active_queue+1; i < NQ; i++ )
if( q[i].first )
break;
if( i == NQ )
break;
active_queue = i;
}
ws_pop( active_queue, mofs, iofs ); // 從q[active_queue]佇列中取出一個結點資料
// 找到這個結點記錄的img和mask中的畫素點，比較該點在mask中的鄰居點
// 鄰居點中如果有標記點：該點與鄰居點的標記型別不同，則該點為分水嶺；該點與鄰居點標記型別相同，則該點不變
// 如果有非標記點：將非標記點擴充為標記點
m = mask + mofs;
ptr = img + iofs;
t = m[-1];
if( t > 0 ) lab = t;
t = m[1];
if( t > 0 )
{
if( lab == 0 ) lab = t;
else if( t != lab ) lab = WSHED; // 如果該畫素點的標記型別和鄰居畫素標記型別都 > 0 且不同，則為分水嶺
}
t = m[-mstep];
if( t > 0 )
{
if( lab == 0 ) lab = t;
else if( t != lab ) lab = WSHED;
}
t = m[mstep];
if( t > 0 )
{
if( lab == 0 ) lab = t;
else if( t != lab ) lab = WSHED;
}
// 因為標記點要麼是初始種子點，要麼是初始階段延伸的種子點的鄰接點
// 該點一定存在一個鄰接點是標記點，所以lab一定會賦值一次，不為 0
assert( lab != 0 );
// 若lab > 0 ，則該點被周圍的標記點擴充；若lab = -1（WSHED）,則該點定義為分水嶺，繼續下一個迴圈
m[0] = lab;
if( lab == WSHED )
continue;
// lab > 0 的情況，確定該點為標記點，且鄰居點中存在未標記點的情況，將未標記點擴充為標記點
if( m[-1] == 0 )
{
c_diff( ptr, ptr - 3, t ); // 計算梯度t
ws_push( t, mofs - 1, iofs - 3 ); // 將m[-1]這一未標記的點擴充為標記點，進隊
active_queue = ws_min( active_queue, t ); // 判斷，若t < 當前處理的佇列active_queue值，則下一次迴圈中處理q[t]佇列，否則繼續處理當前佇列
m[-1] = IN_QUEUE;
}
if( m[1] == 0 )
{
c_diff( ptr, ptr + 3, t );
ws_push( t, mofs + 1, iofs + 3 );
active_queue = ws_min( active_queue, t );
m[1] = IN_QUEUE;
}
if( m[-mstep] == 0 )
{
c_diff( ptr, ptr - istep, t );
ws_push( t, mofs - mstep, iofs - istep );
active_queue = ws_min( active_queue, t );
m[-mstep] = IN_QUEUE;
}
if( m[mstep] == 0 )
{
c_diff( ptr, ptr + istep, t );
ws_push( t, mofs + mstep, iofs + istep );
active_queue = ws_min( active_queue, t );
m[mstep] = IN_QUEUE;
}
}
}
void cv::watershed( InputArray _src, InputOutputArray markers )
{
Mat src = _src.getMat();
CvMat c_src = _src.getMat(), c_markers = markers.getMat();
cvWatershed( &c_src, &c_markers );
}
參考
1. 岡薩雷斯《數位影像處理（第四版）》Chapter 10（所有圖片可在連結中下載）
2.【youcans 的 OpenCV 例程200篇】180.基於距離變換的分水嶺演演算法
3.OpenCV 原始碼中分水嶺演演算法 watershed 函數原始碼註解