Opencv learning notes (8)-watershed algorithm and implementation (C language functions)

I have searched a lot of articles on Watershed algorithms. Many of them are in my blog, but I have not really understood and tested the code of the watershed algorithm. Repeat this time.

A Watershed Algorithm Description

Watershed Algorithm (watershed algorithm), as its name suggests, considers Image Segmentation Based on the composition of watershed. In reality, we can or can imagine a scene where mountains and lakes exist. It must be a situation where water goes around and mountains are surrounded by water. Of course, when necessary, we need to construct a watershed manually to prevent mutual penetration between collection basins. The distinction between the plateaus and the water, and the gap between the lake and the lake, or are all connected, is our lovely watershed ). In order to get a Relatively Concentrated collection basin, it is enough to bring the water up close to the top of the highest nearby mountain, and then it will leak to the neighbors, while the neighbors, hey, the water quality is different, will confuse yourself. In this case, we can use it to obtain the Area of an object with a large border gray scale and a small gray scale in the middle. It is a collection basin.
The water immersion method is to first obtain the starting point through a small threshold, that is, the basis of the water basin. Then, it is to drown the surrounding area, that is, the water immersion process, until the watershed is obtained. Of course, if we want to flood down to the top of the hill, that is, to process the highest gray-scale images, then the dam construction will occur, and different water basins will meet here, we need to clean ourselves, so far, because we have met the border. Even if we have not met the highest gray-scale place when we met each other, we need to build a watershed to distinguish different regions. No longer going up the hill. Build your own watershed in computer graphics, you can use gray scale to characterize the Terrain Height. In the image, we can use the similarity between the gray scale height and the landform height to study the spatial variation of the image's gray scale. This is an analysis of airspace. For example, you can use various forms of gradient computing to obtain the input of algorithms and conduct leaching water treatment. The watershed has a strong edge detection capability and has a good effect on the weak edge. This is related to the setting of the watershed expansion threshold, which can determine the range of the basin expansion. However, self-building capabilities are not affected.

Binary watershed algorithm code

/* ===================================================== ================================ Function Name: watershed function: Use the Mark-watershed algorithm to separate the input image. Implementation: no input parameter description: originalimage -- input image (grayscale image, 0 ~ 255) seedimage -- Tagged Image (Binary Graph, 0-untagged, 1-tagged) labelimage -- output image (1-first split area, 2-The second split area ,...) row -- number of image rows col -- number of image columns return value description: none ============================= */void winapi CDIB:: watershed (unsigned char ** originalimage, char ** seedimage, int ** labelimage, int row, int col) {// using namespace STD; // mark the Region ID, starting from 1, int num = 0; int I, j; // saves the array vector <int *> seedcounts; // temporary seed queue <Point> quetem; // Save the array of all seed queues in the marked area, which contains the pointer vector <queue <point> *> vque; int * array; // The queue <point> * pque; point temp; for (I = 0; I <row; I ++) {for (j = 0; j <Col; j ++) labelimage [I] [J] = 0;} int m, n, k = 0; bool up, down, right, left, upleft, upright, downleft, downright; // 8 directions ctions... // preprocessing, extract and differentiate each marked area, and initialize the seed queue for each marked area. // seed refers to the points at the marked area edge. They can be drowned (or grown) when the water level rises) // Pan's words: I understand the pixel points with smaller gradient values or the points with extremely small gray values. For (I = 0; I <row; I ++) {for (j = 0; j <Col; j ++) {// if you find a marked area if (seedimage [I] [J] = 1) {// the ID of the Area plus num ++; // allocate an array and initialize it to zero, indicating that there are 256 gray-scale arrays = new int [256]; zeromemory (array, 256 * sizeof (INT )); // input the number of seed Arrays into the vector. An array is generated for each scan, and the Region ID is used as the first dimension. The gray level is the second dimension. // Indicates the number of points corresponding to a gray scale in a basin. Seedcounts. push_back (array); // assigns the priority queue of the Mark number, which contains 256 seed queues. // indicates that a queue corresponds to a gray scale, in addition, each queue can store a set of vertex information pque = new queue <point> [256]; // Add it to the queue array, corresponding to the vque of num. push_back (pque); // The current vertex is placed in the temporary seed queue of the marked area temp. X = I; temp. y = J; quetem. push (temp); // The current vertex is marked as processed labelimage [I] [J] = num; seedimage [I] [J] = 127; // indicates that the seeds in the temporary seed queue have been processed until all the seeds have been grown. // The queue information after the seeds have been grown is saved in vque, including area numbers and gray scale, the corresponding points are stored in seedcounts while (! Quetem. empty () {up = down = right = left = false; upleft = upright = downleft = downright = false; // extracts a seed from the queue temp = quetem. front (); M = temp. x; n = temp. y; quetem. pop (); // notice the influence of 127 on the scanning process, which affects the comparison below, but does not affect the scanning if (M> 0) in the while statement) {// if the top is a growth point, add the new seedimage if (seedimage [s-1] [N] = 1) {temp. X = m-1; temp. y = N; quetem. push (temp); // if so, the marked areas will be scanned again in the while loop, because the value is 127 // The new seed point is marked as a flooded area and is the current area, and the region number is recorded to labelimage [s-1] [N] = num; s Eedimage M-1] [N] = 127;} else // otherwise, it cannot grow {up = true ;}} if (M> 0 & n> 0) {If (seedimage [M-1] [n-1] = 1) // Add a new seed {temp. X = m-1; temp. y = n-1; quetem. push (temp); // The new seed point is marked as a drowned area, that is, the next cycle is marked with 127 and is not scanned, and labelimage [M-1] [n-1] = num; seedimage M-1] [n-1] = 127;} else // otherwise, the upper left corner cannot be grown {upleft = true ;}} if (M <row-1) {If (seedimage [M + 1] [N] = 1) // Add a new seed {temp. X = m + 1; temp. y = N; quetem. push (temp); // The new seed point is marked as drowned Region, which is the current region labelimage [M + 1] [N] = num; seedimage [M + 1] [N] = 127 ;} else // otherwise the lower side cannot grow {down = true;} If (M <(Row-1) & n <(Col-1 )) {If (seedimage [M + 1] [n + 1] = 1) // Add a new seed {temp. X = m + 1; temp. y = n + 1; quetem. push (temp); // The new seed point is marked as a flooded area and is labelimage [M + 1] [n + 1] = num in the current area; seedimage [M + 1] [n + 1] = 127;} else // otherwise the lower side cannot grow {downright = true ;}} if (n <col-1) {If (seedimage [m] [n + 1] = 1) // Add a new seed {temp. X = m; temp. y = n + 1; Quetem. push (temp); // The new seed point is marked as a flooded area and is labelimage [m] [n + 1] = num in the current area; seedimage [m] [n + 1] = 127;} else // otherwise the right side cannot grow {right = true ;}} if (M> 0 & n <(Col-1) {If (seedimage [M-1] [n + 1] = 1) // if it is a growth point in the upper right corner, add it as a new seed {temp. X = m-1; temp. y = n + 1; quetem. push (temp); // The new seed point is marked as a flooded area and is labelimage [M-1] [n + 1] = num in the current area; seedimage M-1] [n + 1] = 127;} else // otherwise the upper right corner cannot be grown {upright = true ;}} if (n> 0) {If (seedimage [m] [n-1] = 1) // Add a new seed if the left side is a growth point. {Temp. X = m; temp. y = n-1; quetem. push (temp); // The new seed point is marked as labelimage [m] [n-1] = num; seedimage [m] [n-1] = 127 ;} else // otherwise, the left side cannot be grown {left = true ;}} if (M <(Row-1) & n> 0) {If (seedimage [M + 1] [n-1] = 1) // Add a new seed {temp. X = m + 1; temp. y = n-1; quetem. push (temp); // The new seed point is marked as labelimage [M + 1] [n-1] = num; seedimage [M + 1] [n-1] = 127 ;} else // otherwise, the left side cannot grow {downleft = true;} // if either of the top or bottom sides cannot grow, the current point is one of the initial seed queues. // whether or not the content can grow here is determined by see Dimage value. If (up | down | right | left | upleft | downleft | upright | downright) {temp. X = m; temp. y = N; // The following Vector Array: The first dimension is the mark number; the second dimension is the gray level of the image point // m, point N corresponds to the scanned pixel in the while loop. // Num is the current region number. // this two-dimensional information indicates, the Set and number of member points corresponding to a gray level in an area are expressed by the following two quantities: vque [Num-1] [originalimage [m] [N]. push (temp); // I changed the Num-1 to num... pan's codes... seedcounts [Num-1] [originalimage [m] [N] + +;} // while ends, scanning until quetem is empty and stops. That is, all corresponding nodes cannot grow until (or the surrounding points cannot grow or have grown)} // If end // If (num = 5) // return ;}/// in the above process, if the marked point is 0, it indicates that the point has not been scanned, or it indicates that not the input seed point // This is equivalent to finding all the points of the watershed line of the Initial Region passed by seedimage; and recording each region with a number, at the same time, Edge Points of the collection basin enter the queue. // The above is the Program for collecting water basins. It is also a connected area. /*************************************// /Test the remaining non-basin points. Int seednum; for (I = 0; I <row; I ++) {for (j = 0; j <Col; j ++) {If (seedimage [I] [J] = 0) seednum ++ ;}} cstring STR; Str. format ("pre region num: % d", num); afxmessagebox (STR ); /***********************************/bool actives; // at a certain water level, int waterlevel is the sign of the growth of all marked seeds; // The water level starts from scratch and corresponds to the gray level, use the four-connection method for (waterlevel = 0; waterlevel <180; waterlevel ++) // The second-dimensional... {Actives = true; while (actives) {actives = false; // process the region corresponding to each marker in sequence, and the number of points in the region corresponding to this marker is in seedcounts for (I = 0; I <num; I ++) // The first dimension... {If (! Vque [I] [waterlevel]. empty () // the corresponding watershed is not empty set, I indicates the region number, waterlevel indicates the grayscale {actives = true; while (seedcounts [I] [waterlevel]> 0) {seedcounts [I] [waterlevel] --; // retrieve a vertex with one fewer temp = vque [I] [waterlevel]. front (); // retrieves a vertex in the region and clears the edge vertex, indicating that the pixel has been processed at the current // grayscale level. Vque [I] [waterlevel]. Pop (); M = temp. X; n = temp. Y; // If (M> 0) {If (! Labelimage M-1] [N]) // if it is not processed, it indicates that there is no label. It should have been initialized to 0 before the input. // {temp has also been initialized at the beginning of this function. X = m-1; temp. y = N; labelimage M-1] [N] = I + 1; // the top point is marked as a flooded area. // note that this mark is the same as the area number of the scan point, must it be in the region of this label? Yes // so that at least this point will be scanned in the next round to ensure that no omission is made, but will the next round of processing make it reasonable // is it categorized? The problem is that the tag does not necessarily add it to the seed queue. That is to say, it is only drowned and cannot be drowned up. Only the following growth conditions can be met. If (originalimage M-1] [N] <= waterlevel) // if the top is a growth point, the queue is added to the current queue. the queue of the current height is {vque [I] [waterlevel]. push (temp);} else // otherwise, add the originalimage [M-1] [N] To the grayscale queue. Why? {Vque [I] [originalimage [M-1] [N]. push (temp); seedcounts [I] [originalimage M-1] [N] ++ ;}} if (M <row-1) {If (! Labelimage [M + 1] [N]) // If {temp. X = m + 1; temp. y = N; labelimage [M + 1] [N] = I + 1; // mark the lower point as the drowned area if (originalimage [M + 1] [N] <= waterlevel) // if the following is a growth point, add it to the current queue {vque [I] [waterlevel]. push (temp);} else // otherwise, add originalimage [M + 1] [N]-level queue {vque [I] [originalimage [M + 1] [N]. push (temp); seedcounts [I] [originalimage [M + 1] [N] ++ ;}} if (n <col-1) {If (! Labelimage [m] [n + 1]) // if not processed on the right {temp. X = m; temp. y = n + 1; labelimage [m] [n + 1] = I + 1; // mark the right vertex as a flooded area if (originalimage [m] [n + 1] <= waterlevel) // if the right side is a growth point, add it to the current queue {vque [I] [waterlevel]. push (temp);} else // otherwise, add the originalimage [m] [n + 1] level queue {vque [I] [originalimage [m] [n + 1]. push (temp); seedcounts [I] [originalimage [m] [n + 1] ++ ;}} if (n> 0) {If (! Labelimage [m] [n-1]) // if not processed on the left {temp. X = m; temp. y = n-1; labelimage [m] [n-1] = I + 1; // The left point is marked as a flooded area if (originalimage [m] [n-1] <= waterlevel) // If the left side is a growth point, add it to the current queue {vque [I] [waterlevel]. push (temp);} else // otherwise, add the originalimage [m] [n-1] level queue {vque [I] [originalimage [m] [n-1]. push (temp); seedcounts [I] [originalimage [m] [n-1] ++ ;}}}} // while loop ends seedcounts [I] [waterlevel] = vque [I] [waterlevel]. size ();} // If end} // For Loop end} // while Ring end} // For Loop end /******************************* * ** // test whether the origional segmentation num is changed...int nonzeronum = num; for (m = 0; m <num; m ++) {for (I = 0; I <row; I ++) {for (j = 0; j <Col; j ++) {If (labelimage [I] [J] = m) break;} If (labelimage [I] [J] = m) break ;} if (j = col-1 & I = row-1 & labelimage [I] [J]! = M) nonzeronum --;} Str. format ("new region num: % d", nonzeronum); // This indicates the expandable vertex afxmessagebox (STR ); /*********************************/while (! Vque. Empty () {pque = vque. Back (); Delete [] pque; vque. pop_back () ;}while (! Seedcounts. empty () {array = seedcounts. back (); Delete [] array; seedcounts. pop_back () ;}}// the calling code in the main program is as follows (for processing 8 Bitmap) /**************************************/ // get one image... capcountmdoc * pdoc = getdocument (); hdib m_hdiboptflo = onreadimage (); // hdib is the DIB image handle lpstr lpdib2 = (lpstr): globallock (hglobal) m_hdiboptflo ); int m_imageh = viewdib. dibheight (lpdib2); int m_imagew = viewdib. dibwidth (lpdib2); // viewdib. medianfilter (lpdib2, m_imagew, m_imageh, 3,3, 1,1); // filterlpstr lpdibbits = viewdib. finddibbits (lpdib2); DWORD dwidth = widthbytes (m_imagew * 8); // watershed... lpstr lpdibg = viewdib. makegray8bmp (lpdib2); lpdibg = viewdib. getgrad (lpdibg); lpstr lpdibbitsg = viewdib. finddibbits (lpdibg); int I, K, num = 0; int m = 0, n = 0; // get oritional imageunsigned char ** ori_image; ori_image = alloc_uszchar_space (m_imageh, m_imagew); for (I = 0; I <m_imageh; I ++) {for (k = 0; k <m_imagew; k ++) ori_image [I] [k] = * (lpdibbitsg + I * dwidth + k);} // get seedimage and initialize itchar ** seed_image; seed_image = alloc_char_space (m_imageh, m_imagew ); for (I = 0; I <m_imageh; I ++) {for (k = 0; k <m_imagew; k ++) seed_image [I] [k] = 0 ;} unsigned char uszthre = 3; unsigned char usztem = 0; for (I = 2, num = 0; I <m_imageH-2; I ++) {for (k = 3; k <m_imageW-3; k ++) {usztem = (byte) * (lpdibbitsg + I * dwidth + k); If (usztem <uszthre) {seed_image [I] [k] = 1; num ++ ;}}// get labelimage, and initialize itint ** labelimage = alloc_int_space (m_imageh, m_imagew ); for (I = 0; I <m_imageh; I ++) {for (k = 0; k <m_imagew; k ++) labelimage [I] [k] = 0 ;} // call watershed .. viewdib. watershed (ori_image, seed_image, labelimage, m_imageh, m_imagew); // get the maximum number in labelimage [] [] int nregion = 0; for (I = 0; I <m_imageh; I ++) {for (k = 0; k <m_imagew; k ++) {If (nregion <labelimage [I] [k]) nregion = labelimage [I] [k];}

To better understand the algorithm, the explanation is as follows:

First, prepare mountains and initial water. This mountain is our initial image. For example, we can use the gradient of the naturally acquired image to represent the height of each point in the mountain. The initial water is recorded under the threshold value thre, all mountains below this height are filled with water until the thre height. Thus there are three initial values: unsigned char ** ori_image, char ** seed_image, and INT ** label_image. The last one is to prepare for the final result. Of course, initialization should be done well. For example, ora_image is assigned a gradient value of the original image (256-color grayscale image), and seed_image is the position with water in the initial state, which is reset without water, label_image is initialized to 0, and the region number corresponding to each point is obtained.
The next step is to record the added water into a connected area, that is, to see how many irrelevant collection basins and five separate pools, then we will go up to five lakes, and as high as possible, as long as you don't think of overflow. The algorithm records the number of connected areas into several data structures. The only difference is how to connect these connected areas into one, which is expressed by a data structure. Well, we are going to use a vector container to implement initial storage and save the arrays of all the seed queues in the marked area, it contains the pointer vector <queue <point> *> vque of the seed queue, which is composed of a series of image points in the same region, we come from a collection basin :); the storage method is as follows: queue <point>
* Pque = new queue <point> [256]; vque. push_back (pque), so that a Member is put into this area, that is, the container-set basin storage, each pointer in the container points to a set basin, that is, the connected region we want. Therefore, we can conveniently perform operations by Directly Reading the value of the container data structure, and get a region (queue pointer) with a script) and every queue is not simple, it is not as easy as a column of integer numbers, so ah, this algorithm is really a headache, a member of this queue is a point; the starting pointer of a 256 queue stored in vque is really cruel. That is, vque
[M] [N] indicates a queue, which can store a series of operation points; obviously, the capacity is 256 because there may be different gray-scale points between 0 and 256 in all the initial or final regions, so I use queues to record these Member points in one region, it is very likely that there is only one collection basin here, so there is one region for all the 256 gray-scale points. It is possible that the method for calculating the initial connected areas is, the eight-connected Neighbor Method (also has other methods: Four-Connected Domain, six-Connected Domain), that is, scanning each pixel of the input seed_image one by one, all marked initial collection basins are included in their respective areas, which is a scanning of the renovated image to form an external circulation. First, create a temporary queue quetem to process the connection of the current initial set basin, and store the growth points of the specific initial set basin area scanned one by one, and form an internal loop. The processing of the current scan point is to first determine whether the point is the point of an initial set basin, if not skipped; next, if it is the point of the initial set basin, in this case, is there any non-growth point in its eight-Connected Domain (here the non-growth point refers to the point not marked in seed_image), and the points in the scanned eight-connected Neighbor Area can grow, if yes, add it to the quetem in the temporary queue. if at least one vertex that cannot be grown exists in the scanned connected neighbor, add it to the seed queue, note that the current region is num, and the current scan point is (m, n), so that the current gray scale is ori_image [m] [N], and add the current vertex to the seed queue: vque [Num-1] [ori_image [m] [N]. push (point (m, n )). There are two loops: quetem and seed_image. After the two cycles are complete, the records of each connected initial set basin are obtained, and the storage mark is the region number num. At the same time, the initial watershed is obtained, it is placed in the storage location vque, which identifies the corresponding region number and the gray scale of the corresponding points in the region, that is, the set of points corresponding to the specific gray scale in a specific region; we can obtain the region number of the points of these watershed, and obtain all the gray-scale point information of the region. In a word, there are two functions for calculating the connected areas. One is to mark the initial area, the other is to mark the initial area of the watershed, and the other is to find the initial area of the watershed, then you can start "Shui manliang Shan. This is the drowning process. The drowning process is also implemented by an embedded loop: The External Loop increases the water level (the number of cycles must be less than 256), and The waterlevel increases, it turns out that we have already done the initial water level separation, so now we can start with thre, let the water level rise slowly, let its original lake expand slowly, and try to use its available space, instead of drowning in other neighboring lakes. The inner cycle is the watershed point (in vque [] []) that scans each initial region (the current num, and thus has num loops) and expands according to the given water level. The expansion process is as follows: the current point of the scanned watershed is checked one by one for the four connected neighborhood, if the four-connected domain is a bit unmarked (it must be a high point, and it must have been scanned before a low level ), then, first mark the vertex with num (note that it is the current num), and then judge whether it can grow under the current water level (whether the gray scale is less than or equal to waterlevel ), if it can be grown, it will be added to the vque [num] [waterlevel] seed queue and will go into the inner loop again. Otherwise, if it cannot grow under the current water level, then it is added to the vque [num] [ori_image [neighboring point] queue in the watershed set of this neighboring point. This reciprocating cycle continues until the corresponding area is complete, with a water level, scanning the watershed of all areas, so that each of them expands under one water level at the same time, this ensures that no jump occurs (that is, a water level and a region expand globally ). In the end, all the regions have been expanded in each water level, and we have obtained a segmentation chart. Our large lakes form a way to implement this watershed algorithm. It is not difficult to find that this implementation method cannot generate a new collection basin. That is to say, due to the limitations of the initial collection basin, it is likely to waste most of the collection basins with a higher starting point. In this way, we need to modify the algorithm. Implementation Method: In the process of drowning, We are drowned by the level of the threshold value thre, so we can evaluate the unlabeled points outside the initial region (from seed_image ), mark it with the condition that the internal cycle of this round is completed first, and then a new collection basin is found in the remaining unmarked areas and added to our seed queue, at the beginning of the next water level, there will be another new member who will continue to expand and grow. The members who have their own day will be separated one by one. However, the gradient image is used here. Generally, the initial segmentation of the threshold can meet our requirements, and the gray-scale changes can be smoothly obtained first. The gradient information is powerful enough; if the new basin is used for expansion, it is more suitable for the original image. The main purpose of the watershed algorithm is to find the connected area of the image. When gradient information is used as the input image, there will be a conflict. If gradient computing is performed on the original image without filtering and smoothing, it is easy to divide the object into multiple objects, this is because of the influence of noise. If filtering is performed, it is easy to synthesize several original objects into one object. Of course, the object here mainly refers to the target area with little or similar gray-scale values in the image.

 

 

I am so lucky to have a look at the algorithm today. I basically understand it, but I still don't understand some details. Let's take a look at it and come up with it first. This article is more than the source code in my blog.