connected component analysis error

Question

I'm trying to do connected component analysis.but I'm getting error. I need the vertebral body ;but I'm getting some other objects. Image is:

Result is:

    im= imread('im.bmp');
    figure,imshow(im);
    K1=imadjust(im);
    figure, imshow(K1), title('After Adjustment Image')
    threshold = graythresh(K1);
    originalImage = im2bw(K1, threshold);

    originalImage = bwareaopen(originalImage,100);
    se = strel('disk', 2); %# structuring element
    closeBW = imclose(originalImage,se);

    figure,imshow(closeBW);

    CC = bwconncomp(closeBW);
    L = labelmatrix(CC);
    L2 = bwlabel(K1);
    figure, imshow(label2rgb(L));

Answer 1

Segmentation isn't my area, so I'm not sure what the best approach is. Here are a couple heuristic ideas I came up with:

1. Discard regions that are too big or too small.

It looks like you can expect a certain size from the vertebra.

regionIdxs = unique(L(:));
regionSizes = accumarray(L(:)+1,1);

If we look at regionSizes, we see the region sizes in pixels:

  213360
     919
     887
     810
     601
     695
   14551
     684
    1515
     414
     749
     128
     173
   26658

The regions you want (rows 2-6) are on the range 500-1000. We can probably safely discard regions that are <200 or >2000 in size.

goodRegionIdx = (regionSizes>200) & (regionSizes<2000);
regionIdxs = regionIdxs(goodRegionIdx);
regionSizes = regionSizes(goodRegionIdx);

2. Look at the image moments of the desired regions.

The eigenvalues of the covariance matrix of a distribution characterize its size in its widest direction and its size perpendicular to that direction. We are looking for fat disk-shapes, so we can expect a big eigenvalue and a medium-sized eigenvalue.

[X,Y] = meshgrid(1:size(L,2),1:size(L,1));
for i = 1:length(regionIdxs)
    idx = regionIdxs(i);
    region = L==idx;
    totalmass = sum(region(:));

    Ex(i)  = sum( X(1,:).*sum(region,1) )  / totalmass;
    Ey(i)  = sum( Y(:,1).*sum(region,2))   / totalmass;
    Exy(i) = sum(sum( X.*Y.*region ))      / totalmass;
    Exx(i) = sum(sum( X.*X.*region ))      / totalmass;
    Eyy(i) = sum(sum( Y.*Y.*region ))      / totalmass;

    Varx(i)  = Exx(i) - Ex(i)^2;
    Vary(i)  = Eyy(i) - Ey(i)^2;
    Varxy(i) = Exy(i) - Ex(i)*Ey(i);

    Cov = [Varx(i) Varxy(i); Varxy(i) Vary(i)];

    eig(i,:) = eigs(Cov);
end

If we look at the eigenvalues eig:

  177.6943   30.8029
  142.4484   35.9089
  164.6374   26.2081
  112.6501   22.7570
  138.1674   24.1569
   89.8082   58.8964
  284.2280   96.9304
   83.3226   15.9994
  113.3122   33.7410

We are only interested in rows 1-5, which have eigenvalues on the range 100-200 for the largest and below 50 the the second. If we discard these, get the following regions:

goodRegionIdx = (eig(:,1)>100) & (eig(:,1)<200) & (eig(:,2)<50);
regionIdxs = regionIdxs(goodRegionIdx);

We can plot the regions by using logical OR |.

finalImage = false(size(L));
for i = 1:length(regionIdxs)
    finalImage = finalImage | (L==regionIdxs(i) );
end

We seem to get one false positive. Looking at the ratio of the eigenvalues eig(:,1)./eig(:,2) is one idea but that seem to be a little problematic too.

You could try some sort of outlier detection like RANSAC to try and eliminate the region you don't want, since true vertebra tend to be spatially aligned along a line or curve.

I'm not sure what else to suggest. You may have to look into more advanced segmentation methods like machine learning if you can't find another way to discriminate the good from the bad. Having a stricter preprocessing method might be one thing to try.

Hope that helps.