C++實現(xiàn)分水嶺算法（Watershed Algorithm）
更新時間：2018年01月22日 09:34:37 作者：nagao_kagetora
這篇文章主要為大家詳細(xì)介紹了C++實現(xiàn)分水嶺算法Watershed Algorithm，具有一定的參考價值，感興趣的小伙伴們可以參考一
下
分水嶺分割方法（Watershed Segmentation），是一種基于拓?fù)淅碚摰臄?shù)學(xué)形態(tài)學(xué)的分割方法，其基本思想是把圖像看作是測地學(xué)上的拓?fù)涞孛玻瑘D像中每一點(diǎn)像素的灰度值表示該點(diǎn)的海拔高度，每一個局部極小值及其影響區(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)注。
分水嶺變換得到的是輸入圖像的集水盆圖像，集水盆之間的邊界點(diǎn)，即為分水嶺。顯然，分水嶺表示的是輸入圖像極大值點(diǎn)。因此，為得到圖像的邊緣信息，通常把梯度圖像作為輸入圖像，即：
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ū)域的邊緣點(diǎn)的灰度級進(jìn)行從低到高排序，然后在從低到高實現(xiàn)淹沒的過程，梯度圖像用Sobel算子計算獲得。對梯度圖像進(jìn)行閾值處理時，選取合適的閾值對最終分割的圖像有很大影響，因此閾值的選取是圖像分割效果好壞的一個關(guān)鍵。缺點(diǎ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語言 基本語法,這里提供簡單的示例和代碼來詳細(xì)講解C語言的基本語法，開始學(xué)習(xí)C語言的朋友可以看一下2016-08-08
這篇文章主要用代碼介紹了C語言線性表的順序線性表，對于學(xué)習(xí)數(shù)據(jù)結(jié)構(gòu)與算法的朋友很有參考借鑒價值,需要的朋友可以參考下2021-09-09
這篇文章主要介紹了C++ 中二分查找遞歸非遞歸實現(xiàn)并分析的相關(guān)資料,需要的朋友可以參考下2017-06-06
本文主要介紹C 語言數(shù)據(jù)類型轉(zhuǎn)換,這里通過代碼實例進(jìn)行詳解,這是C語言基礎(chǔ)部分,需要的朋友可以參考下2016-07-07
這篇文章主要給大家介紹了關(guān)于c語言/c++溢出問題的相關(guān)資料，文中通過圖文介紹的非常詳細(xì)，對大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價值，需要的朋友們下面隨著小編來一起學(xué)習(xí)學(xué)習(xí)吧2021-02-02
這篇文章主要介紹了C++11 變參模板的相關(guān)資料，幫助大家更好的理解和學(xué)習(xí)c++11，感興趣的朋友可以了解下2020-08-08
簡單了解C語言中直接插入排序與直接選擇排序?qū)崿F(xiàn)
這篇文章主要介紹了C語言中直接插入排序與直接選擇排序?qū)崿F(xiàn),插入排序的基本操作就是將一個數(shù)據(jù)插入到已經(jīng)排好序的有序數(shù)據(jù)中,從而得到一個新的、個數(shù)加一的有序數(shù)據(jù),需要的朋友可以參考下2016-03-03
C++語言的模板技術(shù)包括函數(shù)模板和類模板，模板技術(shù)是一種代碼重用技術(shù)，函數(shù)和類是C++語言中兩種主要的重用代碼形式，這篇文章主要介紹了C++函數(shù)模板和類模板,需要的朋友可以參考下2022-09-09
這篇文章主要介紹了C++ CopyFile,MoveFile用法案例詳解,本篇文章通過簡要的案例,講解了該項技術(shù)的了解與使用,以下就是詳細(xì)內(nèi)容,需要的朋友可以參考下2021-09-09
這篇文章主要為大家詳細(xì)介紹了C語言實現(xiàn)關(guān)機(jī)小程序，文中示例代碼介紹的非常詳細(xì)，具有一定的參考價值，感興趣的小伙伴們可以參考一下2020-02-02