How do I efficiently created a BW mask for this microscopic image?

Question

So some background. I was tasked to write a matlab program to count the number yeast cells inside visible-light microscopic images. To do that I think the first step will be cell segmentation. Before I got the real experiment image set, I developed an algorithm use a test image set utilizing watershed. Which looks like this:

The first step of watershed is generating a BW mask for the cells. Then I would generate a bwdist image with imposed local minimums generated from the BW mask. With that I can generate the watershed easily.

As you can see my algorithm rely on the successful generation of BW mask. Because I need to generate the bwdist image and markers from it. Originally, I generate the BW mask following the following steps:

1. generate the Local standard deviation of image sdImage = stdfilt(grayImage, ones(9))

2. Use BW thresholding to generate the initial BW mask binaryImage = sdImage < 8;

3. use imclearborder to clear the background. Use some other code to add the cells on the border back.

---

Background finished. Here is my problem

---

But today I received the new real data sets. The image resolution is much smaller and the light condition is different from the test image set. The color depth is also much smaller. These make my algorithm useless. Here is it:

Using stdfilt failed to generate a good clean images. Instead it generate stuff like this (Note: I have adjusted parameters for the stdfilt function and the BW threshold value, following is the best result I can get) :

As you can see there are light pixels in the center of the cells that not necessary darker than the membrane. Which lead the bw thresholding generate stuff like this:

The new bw image after bw thresholding have either incomplete membrane or segmented cell centers and make them unsuitable to the other steps.

I only start image processing recently and have no idea how should I proceed. If you have an idea please help me! Thanks!

For your convience, I have attached a link from dropbox for a subset of the images

Answer 1

I think there's a fundamental problem in your approach. Your algorithm uses stdfilt in order to binarize the image. But what that essentially means is you're assuming there is there is low "texture" in the background and within the cell. This works for your first image. However, in your second image there is a "texture" within the cell, so this assumption is broken.

I think a stronger assumption is that there is a "ring" around each cell (valid for both images you posted). So I took the approach of detecting this ring instead.

So my approach is essentially:

1. Detect these rings (I use a 'log' filter and then binarize based on positive values. However, this results in a lot of "chatter"
2. Try to remove some of the "chatter" initially by filtering out very small and very large regions
3. Now, fill in these rings. However, there is still some "chatter" and filled regions between cells left
4. Again, remove small and large regions, but since the cells are filled, increase the bounds for what is acceptable.
5. There are still some bad regions, most of the bad areas are going to be regions between cells. Regions between cells are detectable by observing the curvature around the boundary of the region. They "bend inwards" a lot, which is expressed mathematically as a large portion of the boundary having a negative curvature. Also, to remove the rest of the "chatter", these regions will have a large standard deviation in the curvature of their boundary, so remove boundaries with a large standard deviation as well.

Overall, the most difficult part will be removing regions between cells and the "chatter" without removing the actual cells.

Anyway, here's the code (note there are a lot of heuristics and also it's very rough and based on code from older projects, homeworks, and stackoverflow answers so it's definitely far from finished):

cell = im2double(imread('cell1.png'));
if (size(cell,3) == 3) 
    cell = rgb2gray(cell);
end

figure(1), subplot(3,2,1)
imshow(cell,[]);

% Detect edges
hw = 5;
cell_filt = imfilter(cell, fspecial('log',2*hw+1,1));

subplot(3,2,2)
imshow(cell_filt,[]);

% First remove hw and filter out noncell hws
mask = cell_filt > 0;
hw = 5;
mask = mask(hw:end-hw-1,hw:end-hw-1);

subplot(3,2,3)
imshow(mask,[]);

rp = regionprops(mask, 'PixelIdxList', 'Area');
rp = rp(vertcat(rp.Area) > 50 & vertcat(rp.Area) < 2000);

mask(:) = false;
mask(vertcat(rp.PixelIdxList)) = true;

subplot(3,2,4)
imshow(mask,[]);

% Now fill objects
mask1 = true(size(mask) + hw);
mask1(hw+1:end, hw+1:end) = mask;
mask1 = imfill(mask1,'holes');
mask1 = mask1(hw+1:end, hw+1:end);

mask2 = true(size(mask) + hw);
mask2(hw+1:end, 1:end-hw) = mask;
mask2 = imfill(mask2,'holes');
mask2 = mask2(hw+1:end, 1:end-hw);

mask3 = true(size(mask) + hw);
mask3(1:end-hw, 1:end-hw) = mask;
mask3 = imfill(mask3,'holes');
mask3 = mask3(1:end-hw, 1:end-hw);

mask4 = true(size(mask) + hw);
mask4(1:end-hw, hw+1:end) = mask;
mask4 = imfill(mask4,'holes');
mask4 = mask4(1:end-hw, hw+1:end);

mask = mask1 | mask2 | mask3 | mask4;

% Filter out large and small regions again
rp = regionprops(mask, 'PixelIdxList', 'Area');
rp = rp(vertcat(rp.Area) > 100 & vertcat(rp.Area) < 5000);

mask(:) = false;
mask(vertcat(rp.PixelIdxList)) = true;

subplot(3,2,5)
imshow(mask);

% Filter out regions with lots of positive concavity

% Get boundaries
[B,L] = bwboundaries(mask);

% Cycle over boundarys
for i = 1:length(B)
    b = B{i};

    % Filter boundary - use circular convolution
    b(:,1) = cconv(b(:,1),fspecial('gaussian',[1 7],1)',size(b,1));
    b(:,2) = cconv(b(:,2),fspecial('gaussian',[1 7],1)',size(b,1));

    % Find curvature
    curv_vec = zeros(size(b,1),1);
    for j = 1:size(b,1)
        p_b = b(mod(j-2,size(b,1))+1,:);  % p_b = point before
        p_m = b(mod(j,size(b,1))+1,:);    % p_m = point middle
        p_a = b(mod(j+2,size(b,1))+1,:);  % p_a = point after

        dx_ds = p_a(1)-p_m(1);              % First derivative
        dy_ds = p_a(2)-p_m(2);              % First derivative
        ddx_ds = p_a(1)-2*p_m(1)+p_b(1);    % Second derivative
        ddy_ds = p_a(2)-2*p_m(2)+p_b(2);    % Second derivative
        curv_vec(j+1) = dx_ds*ddy_ds-dy_ds*ddx_ds;
    end


    if (sum(curv_vec > 0)/length(curv_vec) > 0.4 || std(curv_vec) > 2.0)
        L(L == i) = 0;
    end
end

mask = L ~= 0;

subplot(3,2,6)
imshow(mask,[])

Output1:

Output2: