Watershed by immersion (Vincent-Soille watershed algorithm)

procedure: Watershed-by-Immersion
Input: digital grey scale image G = (D;E; im).
Output: labelled watershed image lab on D.

#define INIT   1 // initial value of lab image
#define MASK   2 // initial value at each level
#define WSHED 0 // label of the watershed pixels
#define FICTITIOUS (-1, -1)  // fictitious pixel 2 D

curlab <= 0; // curlab is the current label
fifo_init(queue);

for all p in D do
   lab[p] <= INIT;
   dist[p] <= 0; // dist is a work image of distances
end for

SORT pixels in increasing order of grey values (minimum hmin, maximum hmax)

// Start Flooding
for h = hmin to hmax do // Geodesic SKIZ of level h   1 inside level h

   for all p in D with im[p] = h do // mask all pixels at level h
             // these are directly accessible because of the sorting step
      lab[p] <= MASK;
      if p has a neighbour q with (lab[q] > 0 or lab[q] == WSHED) then
         // Initialize queue with neighbours at level h of current basins or watersheds
         dist[p] <= 1;
         fifo_add(p, queue);
      end if
   end for

   curdist <= 1;
   fifo_add(FICTITIOUS, queue);

   loop // extend basins
      p <= fifo_remove(queue);
      if p == FICTITIOUS then
         if fifo_empty(queue) then
            break;
         else
            fifo_add(FICTITIOUS, queue);
            curdist <= curdist + 1;
            p <= fifo_remove(queue);
         end if
      end if
      for all q in neighbours of p do // labelling p by inspecting neighbours
         if dist[q] < curdist and (lab[q] > 0 or lab[q] == WSHED) then
            // q belongs to an existing basin or to watersheds
            if lab[q] > 0 then
               if lab[p] == MASK or lab[p] == WSHED then
                  lab[p] <= lab[q];
               else if lab[p] != lab[q] then
                  lab[p] <= WSHED;
               end if
            else if lab[p] == MASK then
               lab[p] <= WSHED;
            end if
         else if lab[q] == MASK and dist[q] == 0 then // q is plateau pixel
            dist[q] <= curdist + 1;
            fifo_add(q, queue);
         end if
      end for
   end loop

   // detect and process new minima at level h
   for all p in D with im[p] = h do
      dist[p] <= 0; // reset distance to zero
      if lab[p] == MASK then // p is inside a new minimum
         curlab <= curlab + 1;  // create new label
         fifo_add(p, queue);
         lab[p] <= curlab;
         while not fifo_empty(queue) do
            q <= fifo_remove(queue);
            for all r in neighbours of q do // inspect neighbours of q
               if lab[r] == MASK then
                  fifo_add(r, queue);
                  lab[r] <= curlab;
               end if
            end for
         end while
      end if
   end for
end for
// End Flooding

Reference

The Watershed Transform: De nitions, Algorithms and
Parallelization Strategies
Jos B.T.M. Roerdink and Arnold Meijster
Institute for Mathematics and Computing Science
University of Groningen
P.O. Box 800, 9700 AV Groningen, The Netherlands
Email: roe@cs.rug.nl,a.meijster@rc.rug.nl