如何有效地为此显微图像创建 BW 蒙版?
How do I efficiently created a BW mask for this microscopic image?
所以一些背景。我的任务是编写一个 matlab 程序来计算可见光显微图像中酵母细胞的数量。为此,我认为第一步是细胞分割。在我得到真正的实验图像集之前,我开发了一种算法,使用利用 watershed 的测试图像集。看起来像这样:
watershed 的第一步是为细胞生成 BW mask。然后我将生成一个 bwdist 图像,其中强加了从 BW 掩码生成的局部最小值。这样我就可以轻松生成分水岭了。
如您所见,我的算法依赖于成功生成 BW 掩码。因为我需要从中生成 bwdist 图像和标记。最初,我按照以下步骤生成 BW 掩码:
- 生成图像的局部标准差sdImage = stdfilt(grayImage, ones(9))
- 使用 BW 阈值生成初始 BW mask binaryImage = sdImage < 8;
- 使用imclearborder清除背景。使用其他一些代码将边框上的单元格添加回来。
后台完成。这是我的问题
但是今天我收到了新的真实数据集。图像分辨率小得多,光照条件与测试图像集不同。颜色深度也小得多。这些使我的算法无用。在这里:
使用 stdfilt 未能生成清晰的图像。相反,它会生成这样的东西(注意:我已经调整了 stdfilt 函数的参数和 BW 阈值,以下是我可以获得的最佳结果):
如您所见,细胞中心有一些亮像素,不一定比细胞膜暗。哪个导致 bw 阈值生成这样的东西:
bw thresholding 后的新bw 图像有不完整的膜或分割的细胞中心,使它们不适合其他步骤。
我最近才开始图像处理,不知道该如何进行。如果您有想法,请帮助我!谢谢!
为了您的方便,我附上了来自保管箱的 link 用于 subset of the images
我认为您的方法存在根本问题。您的算法使用 stdfilt
来对图像进行二值化。但这实质上意味着您假设单元格中背景 和 中存在低 "texture" 和 。这适用于您的第一张图片。但是,在您的第二张图片中,单元格内有一个 "texture",因此这个假设被打破了。
我认为更有力的假设是每个单元格周围都有一个 "ring"(对您发布的两张图片都有效)。所以我采用了检测这个环的方法。
所以我的方法本质上是:
- 检测这些环(我使用'log'过滤器,然后根据正值进行二值化。但是,这会导致很多"chatter"
- 最初尝试通过过滤掉非常小和非常大的区域
来删除一些 "chatter"
- 现在,填写这些圆环。但是,仍然有一些 "chatter" 和单元格之间的填充区域 left
- 同样,删除小区域和大区域,但由于单元格已填充,因此增加可接受的范围。
- 还有一些坏区,大部分坏区是单元格之间的区域。细胞之间的区域可以通过观察区域边界周围的曲率来检测。他们 "bend inwards" 很多,这在数学上表示为边界的大部分具有负曲率。另外,要删除其余的 "chatter",这些区域的边界曲率标准差较大,因此也应删除标准差较大的边界。
总的来说,最困难的部分是在不移除实际单元格的情况下移除单元格和 "chatter" 之间的区域。
无论如何,这是代码(注意有很多启发式方法,而且它非常粗糙,并且基于旧项目、家庭作业和 Whosebug 答案的代码,所以它肯定远未完成):
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,[])
输出 1:
输出2:
所以一些背景。我的任务是编写一个 matlab 程序来计算可见光显微图像中酵母细胞的数量。为此,我认为第一步是细胞分割。在我得到真正的实验图像集之前,我开发了一种算法,使用利用 watershed 的测试图像集。看起来像这样:
watershed 的第一步是为细胞生成 BW mask。然后我将生成一个 bwdist 图像,其中强加了从 BW 掩码生成的局部最小值。这样我就可以轻松生成分水岭了。
如您所见,我的算法依赖于成功生成 BW 掩码。因为我需要从中生成 bwdist 图像和标记。最初,我按照以下步骤生成 BW 掩码:
- 生成图像的局部标准差sdImage = stdfilt(grayImage, ones(9))
- 使用 BW 阈值生成初始 BW mask binaryImage = sdImage < 8;
- 使用imclearborder清除背景。使用其他一些代码将边框上的单元格添加回来。
后台完成。这是我的问题
但是今天我收到了新的真实数据集。图像分辨率小得多,光照条件与测试图像集不同。颜色深度也小得多。这些使我的算法无用。在这里:
使用 stdfilt 未能生成清晰的图像。相反,它会生成这样的东西(注意:我已经调整了 stdfilt 函数的参数和 BW 阈值,以下是我可以获得的最佳结果):
如您所见,细胞中心有一些亮像素,不一定比细胞膜暗。哪个导致 bw 阈值生成这样的东西:
bw thresholding 后的新bw 图像有不完整的膜或分割的细胞中心,使它们不适合其他步骤。
我最近才开始图像处理,不知道该如何进行。如果您有想法,请帮助我!谢谢!
为了您的方便,我附上了来自保管箱的 link 用于 subset of the images
我认为您的方法存在根本问题。您的算法使用 stdfilt
来对图像进行二值化。但这实质上意味着您假设单元格中背景 和 中存在低 "texture" 和 。这适用于您的第一张图片。但是,在您的第二张图片中,单元格内有一个 "texture",因此这个假设被打破了。
我认为更有力的假设是每个单元格周围都有一个 "ring"(对您发布的两张图片都有效)。所以我采用了检测这个环的方法。
所以我的方法本质上是:
- 检测这些环(我使用'log'过滤器,然后根据正值进行二值化。但是,这会导致很多"chatter"
- 最初尝试通过过滤掉非常小和非常大的区域 来删除一些 "chatter"
- 现在,填写这些圆环。但是,仍然有一些 "chatter" 和单元格之间的填充区域 left
- 同样,删除小区域和大区域,但由于单元格已填充,因此增加可接受的范围。
- 还有一些坏区,大部分坏区是单元格之间的区域。细胞之间的区域可以通过观察区域边界周围的曲率来检测。他们 "bend inwards" 很多,这在数学上表示为边界的大部分具有负曲率。另外,要删除其余的 "chatter",这些区域的边界曲率标准差较大,因此也应删除标准差较大的边界。
总的来说,最困难的部分是在不移除实际单元格的情况下移除单元格和 "chatter" 之间的区域。
无论如何,这是代码(注意有很多启发式方法,而且它非常粗糙,并且基于旧项目、家庭作业和 Whosebug 答案的代码,所以它肯定远未完成):
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,[])
输出 1:
输出2: