Segmenting a grayscale image

Richard picture Richard · Nov 15, 2012 · Viewed 13.1k times · Source

I am having trouble achieving the correct segmentation of a grayscale image:

Image to be segmented

The ground truth, i.e. what I would like the segmentation to look like, is this:

Ground truth

I am most interested in the three components within the circle. Thus, as you can see, I would like to segment the top image into three components: two semi-circles, and a rectangle between them.

I have tried various combinations of dilation, erosion, and reconstruction, as well as various clustering algorithms, including k-means, isodata, and mixture of gaussians--all with varying degrees of success.

Any suggestions would be appreciated.

Edit: here is the best result I've been able to obtain. This was obtained using an active contour to segment the circular ROI, and then applying isodata clustering:

Clusters

There are two problems with this:

  • The white halo around the bottom-right cluster, belonging to the top-left cluster
  • The gray halo around both the top-right and bottom-left cluster, belonging to the center cluster.

Answer

bla picture bla · Nov 15, 2012

Here's a starter... use circular Hough transform to find the circular part. For that I initially threshold the image locally.

 im=rgb2gray(imread('Ly7C8.png'));
 imbw = thresholdLocally(im,[2 2]); % thresold localy with a 2x2 window
 % preparing to find the circle
 props = regionprops(imbw,'Area','PixelIdxList','MajorAxisLength','MinorAxisLength');
 [~,indexOfMax] = max([props.Area]);
 approximateRadius =  props(indexOfMax).MajorAxisLength/2;
 radius=round(approximateRadius);%-1:approximateRadius+1);
 %find the circle using Hough trans.
 h = circle_hough(edge(imbw), radius,'same');
 [~,maxIndex] = max(h(:));
 [i,j,k] = ind2sub(size(h), maxIndex);
 center.x = j;     center.y = i;

 figure;imagesc(im);imellipse(gca,[center.x-radius  center.y-radius 2*radius 2*radius]);
 title('Finding the circle using Hough Trans.');

enter image description here

select only what's inside the circle:

 [y,x] = meshgrid(1:size(im,2),1:size(im,1));
 z = (x-j).^2+(y-i).^2;
 f = (z<=radius^2);
 im=im.*uint8(f);

EDIT:

look for a place to start threshold the image to segment it by looking at the histogram, finding it's first local maxima, and iterating from there until 2 separate segments are found, using bwlabel:

  p=hist(im(im>0),1:255);
  p=smooth(p,5);
  [pks,locs] = findpeaks(p);

  bw=bwlabel(im>locs(1));
  i=0;
  while numel(unique(bw))<3
     bw=bwlabel(im>locs(1)+i); 
     i=i+1;
  end


 imagesc(bw);

enter image description here

The middle part can now be obtained by taking out the two labeled parts from the circle, and what is left will be the middle part (+some of the halo)

 bw2=(bw<1.*f);

but after some median filtering we get something more reasonble

 bw2= medfilt2(medfilt2(bw2));

and together we get:

 imagesc(bw+3*bw2); 

enter image description here

The last part is a real "quick and dirty", I'm sure that with the tools you already used you'll get better results...