C++实现分水岭算法（Watershed Algorithm）

分水岭分割方法（Watershed Segmentation），是一种基于拓扑理论的数学形态学的分割方法，其基本思想是把图像看作是测地学上的拓扑地貌，图像中每一点像素的灰度值表示该点的海拔高度，每一个局部极小值及其影响区域称为集水盆，而集水盆的边界则形成分水岭。分水岭的概念和形成可以通过模拟浸入过程来说明。在每一个局部极小值表面，刺穿一个小孔，然后把整个模型慢慢浸入水中，随着浸入的加深，每一个局部极小值的影响域慢慢向外扩展，在两个集水盆汇合处构筑大坝，即形成分水岭。
　　分水岭的计算过程是一个迭代标注过程。分水岭比较经典的计算方法是L. Vincent提出的。在该算法中，分水岭计算分两个步骤，一个是排序过程，一个是淹没过程。首先对每个像素的灰度级进行从低到高排序，然后在从低到高实现淹没过程中，对每一个局部极小值在h阶高度的影响域采用先进先出（FIFO）结构进行判断及标注。
　　分水岭变换得到的是输入图像的集水盆图像，集水盆之间的边界点，即为分水岭。显然，分水岭表示的是输入图像极大值点。因此，为得到图像的边缘信息，通常把梯度图像作为输入图像，即：
　　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(.)表示梯度运算。
　　分水岭算法对微弱边缘具有良好的响应，图像中的噪声、物体表面细微的灰度变化，都会产生过度分割的现象。但同时应当看出，分水岭算法对微弱边缘具有良好的响应，是得到封闭连续边缘的保证的。另外，分水岭算法所得到的封闭的集水盆，为分析图像的区域特征提供了可能。
　　为消除分水岭算法产生的过度分割，通常可以采用两种处理方法，一是利用先验知识去除无关边缘信息。二是修改梯度函数使得集水盆只响应想要探测的目标。
　　为降低分水岭算法产生的过度分割，通常要对梯度函数进行修改，一个简单的方法是对梯度图像进行阈值处理，以消除灰度的微小变化产生的过度分割。即：
　　g(x,y)=max(grad(f(x,y)),gθ)
　　式中，gθ表示阈值。
　　程序可采用方法：用阈值限制梯度图像以达到消除灰度值的微小变化产生的过度分割，获得适量的区域，再对这些区域的边缘点的灰度级进行从低到高排序，然后在从低到高实现淹没的过程，梯度图像用Sobel算子计算获得。对梯度图像进行阈值处理时，选取合适的阈值对最终分割的图像有很大影响，因此阈值的选取是图像分割效果好坏的一个关键。缺点：实际图像中可能含有微弱的边缘，灰度变化的数值差别不是特别明显，选取阈值过大可能会消去这些微弱边缘。

　　下面用C++实现分水岭算法：

#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;
} 

以上就是本文的全部内容，希望对大家的学习有所帮助，也希望大家多多支持我们。