Implement Adaptive watershed segmentation in Matlab

I will like to implement "Adaptive Watershed Segmentation" in Matlab.

There are six steps in this algorithm. Input is figure(a) and result is figure(d).

Would you please to help me check is there any mistake in my code, and I don't know how to implement the sixth step.

Thank you so much!

Load image:

 input_image = imread('test.gif');

Step 1 : Calculate D(x,y) at each (x,y), obtain the Euclidian distance map of the binary image and assign each value of M(x,y) as 0.

 DT = bwdist(input_image,'euclidean'); % Trandform distance:Euclidian distance      
 [h,w]=size(DT);
 M = zeros(h,w);

Step 2 : Smooth the distance map using Gaussian filter to merge the adjacent maxima, set M(x,y) as 1 if D(x,y) is a local maximum, and then obtain the marker map of the distance map.

 H = fspecial('gaussian');
 gfDT = imfilter(DT,H); 
 M = imregionalmax(gfDT); % maker map, M = local maximum of gfDT

Step3 : Scan the marker map pixel by pixel. If M(x0,y0) is 1, seek the spurious maxima in its neighbourhood with a radius of D(x ,y ).When M(x,y) equals 1 and sqr((x − x0)^2 + (y − y0)^2 ) ≤ D(x0, y0) , set M(x,y) as 0 if D(x,y) < D(x0,y0).

for x0 = 1:h
    for y0 = 1:w
        if M(x0,y0) == 1
            r = ceil(gfDT(x0,y0));
            % range begin:(x0-r,y0-r) end:(x0+r,y0+r)
            xb = x0-r;
            if xb <= 0
                xb =1;
            end

              yb = y0-r;
              if yb <= 0
                  yb =1;
              end

              xe = x0+r;
              if xe > h
              xe = h;
              end

              ye = y0+r;
              if ye > w
                  ye = w;
              end

              for x = xb:xe
                  for y = yb:ye
                      if M(x,y)==1 
                          Pos = [x0,y0 ;x,y];
                          Dis = pdist(Pos,'euclidean');
                          IFA = Dis<= (gfDT(x0,y0));
                          IFB = gfDT(x,y)<gfDT(x0,y0);
                          if ( IFA && IFB)
                              M(x,y) = 0;
                          end
                      end
                  end
              end
          end
      end
  end

Step 4: Calculate the inverse of the distance map,and the local maxima turn out to be the local minima.

 igfDT = -(gfDT);

STep5:Segment the distance map according to the markers by the conventional watershed algorithm and obtain the segmentation of binary image.

 I2 = imimposemin(igfDT,M);
 L = watershed(I2);
 igfDT (L==0)=0;

Step 6 : Straighten the watershed lines by linking the ends of the watershed lines with a straight line and reclassifying the pixels along the straight line.

 I don't know how to implement this step

Answer 1

[Image Analyst]

If you get the endpoints of the line, you can use imline() to burn them into a binary image then AND it with your main image. See attached demo.

[Short]

Thank you for your rely first. But in the demo you provide I should draw the line by myself, right? Is there anyway I could make it auto ? Thank you so much !

[Image Analyst]

No, the manual part was just to get the endpoints. If you have the x,y coordinates, you can pass those in to imline(). I don't know your code but you should try to identify those locations.

[Short]

I have another question is " Calculate the inverse of the distance map,and the local maxima turn out to be the local minima."

 igfDT = -(gfDT);

Is there any wrong ? Thank you!

[Image Analyst]

Yes (wrong) - that's not a correct statement. All that does is negate the distance map. And what does the author consider the local maxima and local minima of a binary image to be?

[Short]

The article said "Once the inverse distance map of the transformed image is calculated, singular markers are defined and imposed as the local minima of transformed image. It is obvious that the marker location and the distance map are dependent on the object shape."

And the article is from http://onlinelibrary.wiley.com/doi/10.1111/j.1365-2818.2009.03125.x/epdf

Thank for your reply.

[Short]

Or I should use the imcomplement()?

[Short]

There is my code without step6. And the result is something wrong in red circle.

And how can I use imline() to straighten the watershed lines?

Thank you!

[Image Analyst]

If you want, you can use bwmorph(bw, 'branchpoints') to find the crossing points and then draw lines between them, but I do not recommend it because some "lines" should stay curves.