How to efficiently create a BW mask for this microscopic image?
So, some kind of background. I was tasked with writing a matlab program to count the number of yeast cells inside microscopic images with visible light. For this, I believe the first step is to segment the cells. Before I got the real set of experimental experiments, I developed an algorithm using a set of test images using a watershed . Which looks like this:
The first step in the watershed is to create a BW mask for the cells. I would then create a bwdist image with local minima superimposed generated from the BW mask. Thanks to this, I can easily create a watershed.
As you can see, my algorithm relies on the successful creation of the BW mask. Because I need to generate the bwdist image and markers. Initially, I generate the BW mask like this:
- generate local standard deviation of the image sdImage = stdfilt (grayImage, ones (9))
- Use BW threshold to generate initial BW mask binaryImage = sdImage <8;
- use imclearborder to clear the background. Use different code to add cells on the border.
The background is complete. Here is my problem
But today I got new real datasets. The image resolution is much smaller and the lighting condition is different from the set test image. The color depth is also much shallower. This makes my algorithm useless. Here he is:
Using stdfilt failed to create nice clean images. Instead, it generates things like this (Note: I have parameters configured for the stdfilt function and BW threshold, the following is the best result I can get):
As you can see, there are light pixels in the center of the cells that are not needed darker than the membrane. Which leads to the threshold bw, generates things like this:
The new image bw after the bw threshold has either an incomplete membrane or segmented cell centers and makes them unusable for other steps.
I have just started image processing and have no idea how I should proceed. If you have an idea please help me! Thank you!
For your confidence, I have added a link from the dropbox for a subset of images
source to share
I think there is a fundamental problem with your approach. Your algorithm is used stdfilt
to binarize the image. But what this means is you are assuming there is a low "texture" in the background and inside the cell. This works for your first image. However, the second image has a "texture" in the cell, so this assumption is violated.
I think the stronger assumption is that there is a "ring" around each cell (valid for both images you posted). So I took an approach to finding this ring.
So my approach is essentially:
- Finding these rings (I use the "log" filter and then binarize based on positive values. However, this leads to a lot of "chatter"
- Try to remove some of the chatter first by filtering out very small and very large areas.
- Now fill in these rings. However, there is still some "chatter" and filled areas between the cells on the left.
- Delete the small and large areas again, but since the cells are full, increase the range.
- There are still some bad regions, most of the bad regions will be between cells. Regions between cells can be detected by observing the curvature around the border of the area. They are "bent inward" a lot, which is mathematically expressed as a large portion of the border having a negative curvature. Also, to remove the rest of the "chatter", these areas will have a large standard deviation in the curvature of their border, so remove the large standard deviation borders.
All in all, the trickiest part will be removing the areas between cells and chatting without removing the actual cells.
Anyway, here's the code (note that there are many heuristics, and also very crude and based on code from old projects, homework and stackoverflow answers, so it's definitely far from complete):
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,[])
Output 1:
Output 2:
source to share