Matlab - Watershed to extract lines - lost information

Question

I have a vein image as follow. I use watershed algorithm to extract the skeleton of the vein.

My code: (K is the original image).

level = graythresh(K);
BW = im2bw(K,level);
D = bwdist(~BW);
DL = watershed(D);
bgm = DL == 0;
imshow(bgm);

The result is:

As you can see lot of information is lost. Can anybody help me out? Thanks.

Answer 1

It looks like the lighting is somewhat uneven. This can be corrected using certain morphological operations. The basic idea is to compute an image that represents just the uneven lighting and subtract it, or to divide by it (which also enhances contrast). Because we want to find just the lighting, it is important to use a large enough structuring element, so that the operation examines more global properties rather than local ones.

%# Load image and convert to [0,1].
A = im2double(imread('http://i.stack.imgur.com/TQp1i.png'));
%# Any large (relative to objects) structuring element will do.
%# Try sizes up to about half of the image size.
se = strel('square',32);
%# Removes uneven lighting and enhances contrast.
B = imdivide(A,imclose(A,se));
%# Otsu's method works well now.
C = B > graythresh(B);
D = bwdist(~C);
DL = watershed(D);
imshow(DL==0);

Here are C (left), plus DL==0 (center) and its overlay on the original image:

Answer 2

You are losing information because when you apply im2bw, you are basically converting your uint8 image, where the pixel brightness takes a value from intmin('uint8')==0 to intmax('uint8')==255, into a binary image (where only logical values are used). This is what entails a loss of information that you observed. If you display the image BW you will see that all the elements of K that had a value greater than the threshold level turn into ones, while those ones below the threshold turn into zeros.

Answer 3

Yes, you'll need to lower your threshold likely (lower than what Otsu's method is giving you). And if the edge map is noisy when you lower the threshold, you should apply a 2-D Gaussian smoothing filter before you lower the threshold. This will move the edges slightly but will clean up noise too, so it's a tradeoff.

The 2-D Gaussian can be applied doing something like

w=gausswin(N,Alpha)  % you'll have to play with N and alpha
K = imfilter(K,w,'same','symmetric'); % something like these options

Before you apply the rest of your algorithm.