C++ watershed algorithm (Watershed Algorithm)
- 2020-06-01 10:46:26
- OfStack
Watershed segmentation method (Watershed Segmentation), is a kind of mathematical morphology segmentation method based on topology theory, the basic idea is to put the image as the topology of landform on geodesy, each 1 point in the image pixel gray value indicates that point elevation, each one local minimum values and effect area known as the reception basin, and set the boundary of the basin form a watershed. The concept and formation of the watershed can be illustrated by simulating the immersion process. On the surface of each local minimum, a small hole is punctured, and then the whole model is slowly immersed in the water. With the deepening of immersion, the influence domain of each local minimum slowly expands outwards, and a dam is built at the confluence of two water basins to form a watershed.
The calculation process of watershed is an iterative labeling process. The classical calculation method of watershed is proposed by L. Vincent. In this algorithm, the watershed calculation is divided into two steps, one is the sorting process, the other is the inundation process. First, the gray level of each pixel is sorted from low to high. Then, in the process of submersion from low to high, the first in first out (FIFO) structure is used to judge and label the influence domain of each local minimum at the height of order h.
Watershed transformation results in the sink image of the input image, and the boundary point between sinks is the watershed. Clearly, the watershed represents the maximum point of the input image. Therefore, in order to obtain the edge information of the image, the gradient image is usually used as the input image, that is:
grad (f (x, y)) = ((f (x - 1, y) - f (y x + 1)) ^ 2 + (f (x y - 1) - f (x, y + 1)) ^ 2) ^ 0.5
Where, f(x,y) represents the original image, and grad(.) represents the gradient operation.
Watershed algorithm has a good response to the weak edge, the noise in the image, the object surface of the subtle gray change, will produce the phenomenon of excessive segmentation. But at the same time, it should be seen that the watershed algorithm has a good response to the weak edge and is guaranteed to get the closed continuous edge. In addition, the closed basin obtained by watershed algorithm provides the possibility to analyze the regional features of the image.
In order to eliminate the excessive segmentation generated by watershed algorithm, two processing methods can be adopted. One is to remove irrelevant edge information by using prior knowledge. 2 is to modify the gradient function so that the sink only responds to the target to be probed.
In order to reduce the excessive segmentation generated by watershed algorithm, the gradient function is usually modified. A simple method is to conduct threshold processing on the gradient image to eliminate the excessive segmentation caused by the small change in gray level. That is:
g (x y) = max (grad (f (x y)), g theta)
In this formula, g represents the threshold value.
The program can adopt the following methods: the threshold value is used to limit the gradient image so as to eliminate the excessive segmentation caused by the small change in the gray value, to obtain an appropriate amount of regions, and then the gray level of the edge points of these regions is sorted from low to high, and then the gradient image is calculated by Sobel operator in the process of realizing the inundation from low to high. When threshold processing is carried out on the gradient image, the selection of the appropriate threshold has a great influence on the final image segmentation, so the selection of the threshold is a key to the image segmentation effect. Disadvantages: the actual image may contain weak edges, and the difference in the value of the grayscale change is not particularly obvious. If the selection threshold is too large, these weak edges may be eliminated.
The watershed algorithm is implemented by C++ below:
#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;
}