C++實現(xiàn)分水嶺算法(Watershed Algorithm)
分水嶺分割方法(Watershed Segmentation),是一種基于拓?fù)淅碚摰臄?shù)學(xué)形態(tài)學(xué)的分割方法,其基本思想是把圖像看作是測地學(xué)上的拓?fù)涞孛?,圖像中每一點像素的灰度值表示該點的海拔高度,每一個局部極小值及其影響區(qū)域稱為集水盆,而集水盆的邊界則形成分水嶺。分水嶺的概念和形成可以通過模擬浸入過程來說明。在每一個局部極小值表面,刺穿一個小孔,然后把整個模型慢慢浸入水中,隨著浸入的加深,每一個局部極小值的影響域慢慢向外擴(kuò)展,在兩個集水盆匯合處構(gòu)筑大壩,即形成分水嶺。
分水嶺的計算過程是一個迭代標(biāo)注過程。分水嶺比較經(jīng)典的計算方法是L. Vincent提出的。在該算法中,分水嶺計算分兩個步驟,一個是排序過程,一個是淹沒過程。首先對每個像素的灰度級進(jìn)行從低到高排序,然后在從低到高實現(xiàn)淹沒過程中,對每一個局部極小值在h階高度的影響域采用先進(jìn)先出(FIFO)結(jié)構(gòu)進(jìn)行判斷及標(biāo)注。
分水嶺變換得到的是輸入圖像的集水盆圖像,集水盆之間的邊界點,即為分水嶺。顯然,分水嶺表示的是輸入圖像極大值點。因此,為得到圖像的邊緣信息,通常把梯度圖像作為輸入圖像,即:
grad(f(x,y))=((f(x-1,y)-f(x+1,y))^2 + (f(x,y-1)-f(x,y+1))^2)^0.5
式中,f(x,y)表示原始圖像,grad(.)表示梯度運(yùn)算。
分水嶺算法對微弱邊緣具有良好的響應(yīng),圖像中的噪聲、物體表面細(xì)微的灰度變化,都會產(chǎn)生過度分割的現(xiàn)象。但同時應(yīng)當(dāng)看出,分水嶺算法對微弱邊緣具有良好的響應(yīng),是得到封閉連續(xù)邊緣的保證的。另外,分水嶺算法所得到的封閉的集水盆,為分析圖像的區(qū)域特征提供了可能。
為消除分水嶺算法產(chǎn)生的過度分割,通??梢圆捎脙煞N處理方法,一是利用先驗知識去除無關(guān)邊緣信息。二是修改梯度函數(shù)使得集水盆只響應(yīng)想要探測的目標(biāo)。
為降低分水嶺算法產(chǎn)生的過度分割,通常要對梯度函數(shù)進(jìn)行修改,一個簡單的方法是對梯度圖像進(jìn)行閾值處理,以消除灰度的微小變化產(chǎn)生的過度分割。即:
g(x,y)=max(grad(f(x,y)),gθ)
式中,gθ表示閾值。
程序可采用方法:用閾值限制梯度圖像以達(dá)到消除灰度值的微小變化產(chǎn)生的過度分割,獲得適量的區(qū)域,再對這些區(qū)域的邊緣點的灰度級進(jìn)行從低到高排序,然后在從低到高實現(xiàn)淹沒的過程,梯度圖像用Sobel算子計算獲得。對梯度圖像進(jìn)行閾值處理時,選取合適的閾值對最終分割的圖像有很大影響,因此閾值的選取是圖像分割效果好壞的一個關(guān)鍵。缺點:實際圖像中可能含有微弱的邊緣,灰度變化的數(shù)值差別不是特別明顯,選取閾值過大可能會消去這些微弱邊緣。
下面用C++實現(xiàn)分水嶺算法:
#define _USE_MATH_DEFINES #include <cstddef> #include <cstdlib> #include <cstring> #include <climits> #include <cfloat> #include <ctime> #include <cmath> #include <cassert> #include <vector> #include <stack> #include <queue> using namespace std; typedef void GVVoid; typedef bool GVBoolean; typedef char GVChar; typedef unsigned char GVByte; typedef short GVInt16; typedef unsigned short GVUInt16; typedef int GVInt32; typedef unsigned int GVUInt32; typedef long long GVInt64; typedef unsigned long long GVUInt64; typedef float GVFloat32; typedef double GVFloat64; const GVBoolean GV_TRUE = true; const GVBoolean GV_FALSE = false; const GVByte GV_BYTE_MAX = UCHAR_MAX; const GVInt32 GV_INT32_MAX = INT_MAX; const GVInt32 GV_INT32_MIX = INT_MIN; const GVInt64 GV_INT64_MAX = LLONG_MAX; const GVInt64 GV_INT64_MIN = LLONG_MIN; const GVFloat32 GV_FLOAT32_MAX = FLT_MAX; const GVFloat32 GV_FLOAT32_MIN = FLT_MIN; const GVFloat64 GV_FLOAT64_MAX = DBL_MAX; const GVFloat64 GV_FLOAT64_MIN = DBL_MIN; class GVPoint; class GVPoint { public: GVInt32 x; GVInt32 y; public: GVPoint() : x(0), y(0) { } GVPoint(const GVPoint &obj) : x(obj.x), y(obj.y) { } GVPoint(GVInt32 x, GVInt32 y) : x(x), y(y) { } public: GVBoolean operator ==(const GVPoint &right) const { return ((x == right.x) && (y == right.y)); } GVBoolean operator !=(const GVPoint &right) const { return (!(x == right.x) || !(y == right.y)); } }; /* * <Parameter> * <image> image data; * <width> image width; * <height> image height; * <label out> image of labeled watersheds. */ GVVoid gvWatershed( const GVByte *image, GVInt32 width, GVInt32 height, GVInt32 *label) { // Local constants const GVInt32 WSHD = 0; const GVInt32 INIT = -1; const GVInt32 MASK = -2; const GVPoint FICT_PIXEL = GVPoint(~0, ~0); // Image statistics and sorting GVInt32 size = width * height; GVInt32 *image_stat = new GVInt32[GV_BYTE_MAX + 1]; GVInt32 *image_space = new GVInt32[GV_BYTE_MAX + 1]; GVPoint *image_sort = new GVPoint[size]; ::memset(image_stat, 0, sizeof (GVInt32) * (GV_BYTE_MAX + 1)); ::memset(image_space, 0, sizeof (GVInt32) * (GV_BYTE_MAX + 1)); ::memset(image_sort, 0, sizeof (GVPoint) * size); for (GVInt32 i = 0; !(i == size); ++i) { image_stat[image[i]]++; } for (GVInt32 i = 0; !(i == GV_BYTE_MAX); ++i) { image_stat[i + 1] += image_stat[i]; } for (GVInt32 i = 0; !(i == height); ++i) { for (GVInt32 j = 0; !(j == width); ++j) { GVByte space = image[i * width + j]; GVInt32 index = image_stat[space] - (++image_space[space]); image_sort[index].x = j; image_sort[index].y = i; } } for (GVInt32 i = GV_BYTE_MAX; !(i == 0); --i) { image_stat[i] -= image_stat[i - 1]; } // Watershed algorithm GVPoint *head = image_sort; GVInt32 space = 0; GVInt32 *dist = new GVInt32[size]; GVInt32 dist_cnt = 0; GVInt32 label_cnt = 0; std::queue<GVPoint> queue; ::memset(dist, 0, sizeof (GVInt32) * size); ::memset(label, ~0, sizeof (GVInt32) * size); for (GVInt32 h = 0; !(h == (GV_BYTE_MAX + 1)); ++h) { head += space; space = image_stat[h]; for (GVInt32 i = 0; !(i == space); ++i) { GVInt32 index = head[i].y * width + head[i].x; GVInt32 index_l = ((head[i].x - 1) < 0) ? -1 : ((head[i].x - 1) + head[i].y * width); GVInt32 index_r = !((head[i].x + 1) > width) ? -1 : ((head[i].x + 1) + head[i].y * width); GVInt32 index_t = ((head[i].y - 1) < 0) ? -1 : (head[i].x + (head[i].y - 1) * width); GVInt32 index_b = !((head[i].y + 1) > height) ? -1 : (head[i].x + (head[i].y + 1) * width); label[index] = MASK; if ( (!(index_l < 0) && !(label[index_l] < WSHD)) || (!(index_r < 0) && !(label[index_r] < WSHD)) || (!(index_t < 0) && !(label[index_t] < WSHD)) || (!(index_b < 0) && !(label[index_b] < WSHD))) { dist[index] = 1; queue.push(head[i]); } } dist_cnt = 1; queue.push(FICT_PIXEL); while (GV_TRUE) { GVPoint top = queue.front(); GVInt32 index = top.y * width + top.x; GVInt32 index_l = ((top.x - 1) < 0) ? -1 : ((top.x - 1) + top.y * width); GVInt32 index_r = !((top.x + 1) > width) ? -1 : ((top.x + 1) + top.y * width); GVInt32 index_t = ((top.y - 1) < 0) ? -1 : (top.x + (top.y - 1) * width); GVInt32 index_b = !((top.y + 1) > height) ? -1 : (top.x + (top.y + 1) * width); queue.pop(); if (top == FICT_PIXEL) { if (queue.empty()) break; else { ++dist_cnt; queue.push(FICT_PIXEL); top = queue.front(); queue.pop(); } } if (!(index_l < 0)) { if ((dist[index_l] < dist_cnt) && !(label[index_l] < WSHD)) { if (label[index_l] > WSHD) { if ((label[index] == MASK) || (label[index] = WSHD)) label[index] = label[index_l]; else if (!(label[index] == label[index_l])) label[index] = WSHD; } else if (label[index] == MASK) { label[index] = WSHD; } } else if ((label[index_l] == MASK) && (dist[index_l] == 0)) { dist[index_l] = dist_cnt + 1; queue.push(GVPoint(top.x - 1, top.y)); } } if (!(index_r < 0)) { if ((dist[index_r] < dist_cnt) && !(label[index_r] < WSHD)) { if (label[index_r] > WSHD) { if ((label[index] == MASK) || (label[index] = WSHD)) label[index] = label[index_r]; else if (!(label[index] == label[index_r])) label[index] = WSHD; } else if (label[index] == MASK) { label[index] = WSHD; } } else if ((label[index_r] == MASK) && (dist[index_r] == 0)) { dist[index_r] = dist_cnt + 1; queue.push(GVPoint(top.x + 1, top.y)); } } if (!(index_t < 0)) { if ((dist[index_t] < dist_cnt) && !(label[index_t] < WSHD)) { if (label[index_t] > WSHD) { if ((label[index] == MASK) || (label[index] = WSHD)) label[index] = label[index_t]; else if (!(label[index] == label[index_t])) label[index] = WSHD; } else if (label[index] == MASK) { label[index] = WSHD; } } else if ((label[index_t] == MASK) && (dist[index_t] == 0)) { dist[index_t] = dist_cnt + 1; queue.push(GVPoint(top.x, top.y - 1)); } } if (!(index_b < 0)) { if ((dist[index_b] < dist_cnt) && !(label[index_b] < WSHD)) { if (label[index_b] > WSHD) { if ((label[index] == MASK) || (label[index] = WSHD)) label[index] = label[index_b]; else if (!(label[index] == label[index_b])) label[index] = WSHD; } else if (label[index] == MASK) { label[index] = WSHD; } } else if ((label[index_b] == MASK) && (dist[index_b] == 0)) { dist[index_b] = dist_cnt + 1; queue.push(GVPoint(top.x, top.y + 1)); } } } for (GVInt32 i = 0; !(i == space); ++i) { GVInt32 index = head[i].y * width + head[i].x; dist[index] = 0; if (label[index] == MASK) { label_cnt++; label[index] = label_cnt; queue.push(head[i]); while (!queue.empty()) { GVPoint top = queue.front(); GVInt32 index_l = ((top.x - 1) < 0) ? -1 : ((top.x - 1) + top.y * width); GVInt32 index_r = !((top.x + 1) > width) ? -1 : ((top.x + 1) + top.y * width); GVInt32 index_t = ((top.y - 1) < 0) ? -1 : (top.x + (top.y - 1) * width); GVInt32 index_b = !((top.y + 1) > height) ? -1 : (top.x + (top.y + 1) * width); queue.pop(); if (!(index_l < 0) && (label[index_l] == MASK)) { queue.push(GVPoint(top.x - 1, top.y)); label[index_l] = label_cnt; } if (!(index_r < 0) && (label[index_r] == MASK)) { queue.push(GVPoint(top.x + 1, top.y)); label[index_r] = label_cnt; } if (!(index_t < 0) && (label[index_t] == MASK)) { queue.push(GVPoint(top.x, top.y - 1)); label[index_t] = label_cnt; } if (!(index_b < 0) && (label[index_b] == MASK)) { queue.push(GVPoint(top.x, top.y + 1)); label[index_b] = label_cnt; } } } } } // Release resources delete[] image_stat; delete[] image_space; delete[] image_sort; delete[] dist; }
以上就是本文的全部內(nèi)容,希望對大家的學(xué)習(xí)有所幫助,也希望大家多多支持腳本之家。
相關(guān)文章
C語言編程技巧 關(guān)于const和#define的區(qū)別心得
盡量用const和inline而不用#define 這個條款最好稱為:“盡量用編譯器而不用預(yù)處理”,因為#define經(jīng)常被認(rèn)為好象不是語言本身的一部分。這是問題之一。再看下面的語句:2013-02-02基于malloc與free函數(shù)的實現(xiàn)代碼及分析
本篇文章介紹了malloc與free函數(shù)的實現(xiàn)代碼及分析。需要的朋友參考下2013-05-05C++基礎(chǔ)入門教程(六):為什么創(chuàng)建類的時候要用new?
這篇文章主要介紹了C++基礎(chǔ)入門教程(六):為什么創(chuàng)建類的時候要用new?本文講解了使用new創(chuàng)建動態(tài)結(jié)構(gòu)體、為什么要有new、自動存儲(自動變量、局部變量)、動態(tài)存儲、vector和array等內(nèi)容,需要的朋友可以參考下2014-11-11詳談浮點精度(float、double)運(yùn)算不精確的原因
這篇文章主要介紹了詳談浮點精度(float、double)運(yùn)算不精確的原因,具有很好的參考價值,希望對大家有所幫助。如有錯誤或未考慮完全的地方,望不吝賜教2021-12-12c++中處理相關(guān)數(shù)學(xué)函數(shù)
數(shù)學(xué)庫函數(shù)聲明在 math.h 中,主要有:2013-04-04