image - この顕微鏡画像の BW マスクを効率的に作成するにはどうすればよいですか?

Question

だからいくつかの背景。私は、可視光顕微鏡画像内の酵母細胞の数をカウントするための matlab プログラムを作成することを任されました。そのためには、最初のステップが細胞のセグメンテーションになると思います。実際の実験画像セットを入手する前に、流域を利用したテスト画像セットを使用するアルゴリズムを開発しました。次のようになります。

流域の最初のステップは、セルの BW マスクを生成することです。次に、BW マスクから生成された局所最小値を適用して bwdist 画像を生成します。これで、流域を簡単に生成できます。

ご覧のとおり、私のアルゴリズムは BW マスクの生成の成功に依存しています。そこから bwdist イメージとマーカーを生成する必要があるためです。もともと、次の手順に従って BW マスクを生成します。

1. 画像のローカル標準偏差を生成しますsdImage = stdfilt(grayImage, ones(9))

2. BW しきい値を使用して、初期 BW マスクを生成します。 binaryImage = sdImage < 8;

3. 背景をクリアするにはimclearborderを使用します。他のコードを使用して、境界線にセルを追加します。

---

背景完成。これが私の問題です

---

しかし、今日、新しい実際のデータ セットを受け取りました。画像の解像度ははるかに小さく、光の状態はテスト画像セットとは異なります。色深度もはるかに小さいです。これらは私のアルゴリズムを役に立たなくします。はい、これ：

stdfiltを使用しても、きれいなイメージを生成できませんでした。代わりに、次のようなものを生成します (注: stdfilt関数と BW しきい値のパラメーターを調整しました。以下が得られる最良の結果です)。

ご覧のとおり、細胞の中心に明るいピクセルがあり、膜よりも暗い必要はありません。bw のしきい値処理が次のようなものを生成する原因は次のとおりです。

bw しきい値処理後の新しい bw 画像には、不完全な膜またはセグメント化された細胞中心があり、他のステップには適していません。

画像処理を最近始めたばかりで、どのように進めたらよいかわかりません。アイデアがあれば、私を助けてください！ありがとう！

ご参考までに、画像のサブセットへのドロップボックスからのリンクを添付しました

Answer 1

あなたのアプローチには根本的な問題があると思います。アルゴリズムはstdfilt、画像を二値化するために使用します。しかし、それが本質的に意味することは、背景とセル内に低い「テクスチャ」があると仮定していることです。これは最初の画像で機能します。ただし、2 番目の画像では、セル内に「テクスチャ」があるため、この仮定は破られています。

より強い仮定は、各セルの周りに「リング」があることだと思います（投稿した両方の画像に有効です）。そこで、代わりにこのリングを検出するというアプローチを取りました。

したがって、私のアプローチは基本的に次のとおりです。

1. これらのリングを検出します (「ログ」フィルターを使用し、正の値に基づいて 2 値化します。ただし、これにより多くの「チャター」が発生します)。
2. 非常に小さい領域と非常に大きい領域を除外することにより、最初に「チャター」の一部を削除してみてください
3. 次に、これらのリングを埋めます。ただし、セル間の「チャター」と塗りつぶされた領域がまだ残っています。
4. ここでも、小さな領域と大きな領域を削除しますが、セルが塗りつぶされているため、許容できる範囲を広げます。
5. まだいくつかの悪い領域がありますが、ほとんどの悪い領域はセル間の領域になります。セル間の領域は、領域の境界付近の曲率を観察することで検出できます。それらは大きく「内側に曲がり」ます。これは、負の曲率を持つ境界の大部分として数学的に表現されます。また、残りの「チャター」を除去するために、これらの領域は境界の曲率に大きな標準偏差を持つため、標準偏差が大きい境界も除去します。

全体として、最も難しい部分は、実際のセルを削除せずに、セルと「チャター」の間の領域を削除することです。

とにかく、ここにコードがあります（多くのヒューリスティックがあり、非常に大まかで、古いプロジェクト、宿題、およびスタックオーバーフローの回答のコードに基づいているため、完成にはほど遠いことに注意してください）：

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: