Detecting and isolating part of an image

Question

I'm working with cross-sectional CT scans of the leg using matlab's image processing toolbox. Part of what I'm attempting to do is to automatically isolate the leg in any image, which requires removing any extraneous objects that are present in the image. So far my code is able to take out most objects, which works for the majority of my images. However, some of the images have an object that is flush with the leg, which my code isn't able to remove.

Here's an example of the problem: The first image shows the original image, and the second one shows the "filtered" image, after it has gotten rid of extraneous objects. But as you can see, there is still something outside of the leg cross-section that I would like to get rid of.

I've thought of using something to detect ellipses, but unfortunately some of my images are not very good homogenous shapes (like the one below).

And here's the current code I have:

% Read input file
inputfile = dlmread(filename,';');

% Normalise input file (0 to 1) & show original image
I_original = mat2gray(inputfile,[min(inputfile(:)) max(inputfile(:))]);
figure; imshow(I_original)

% Detect leg
BW1 = im2bw(I_original,.2);
BW2 = imfill(BW1,'holes');

% Remove misc. objects around leg (i.e., stray objects <5000 pixels)
BW_final = bwareaopen(BW2,5000);

% Generate filtered image
I_filtered = I_original;
I_filtered(imcomplement(BW_final)) = NaN;
figure; imshow(I_filtered)

Any suggestions on how to fix this would be appreciated! Thanks!

Answer 1

Since the Extraneous part is connected with leg section via a small thread like thing, Try converting it to binary Image, and then eroding the image to get rid of smaller objects. That should disconnect your both objects. Then choose the object with biggest area.

Helpful functions - im2bw, imerode, rogionprops.

Answer 2

Behold the power of shortest path.

Prepare image

I'm using the image in the question for this so there are some extra steps

I = imread('leg.jpg');
I = rgb2gray(I);
I = double(I);
I = I(56:682,411:1027);

Bandpass filter

We are interested in finding the edges of the oval sections. The first step is to make a bandpass filter, in this case a difference of Gaussians, to isolate the smaller features.

H1 = fspecial('gaussian',[60,60],5);
H2 = fspecial('gaussian',[60,60],2);
H = H2-H1;

Edge filter Use the Hilbert transform to make two edge filters, one in the x and one in the y direction.

Hedge = imag(hilbert(H));

Filtered images

Filter the image with the edge filters. Square each of the edge filter responses, add them together, then take the square root. The result is the magnitude of the edges.

Iodd = sqrt(imfilter(I,Hedge).^2 + imfilter(I,Hedge.').^2);

Shortest Path

Shortest phase looks for the path of lowest energy (but not negative), so invert the edge values. We use the log of the edge values here as that makes strong and weak edges a little closer together in value

G = log(Iodd);
G = max(G(:)) - G;

Perform shortest path. I won't explain this fully as it would require an entire paper.

[ path, energy, lastIndex, pathImage ] = ...
    circularShortestPath(G, 0.1, [300,300], [100,400], [400,720]);

imagesc(G);
hold on;
plot(path(:,1),path(:,2),'r','LineWidth',2);
hold off;
pause;

imagesc(G);
hold on;
fill(path(:,1),path(:,2),'r');
hold off;
pause;

Segment

Use the path to segment the image

J = zeros(size(I));
mask = roipoly(J,path(:,1),path(:,2));
segmentedI = I .* mask;
imagesc(segmentedI); title('segmented image'); colormap gray(256); 

More

Changing the centre point of the shortest path you can segment other things as well:

[ path, energy, lastIndex, pathImage ] = ...
    circularShortestPath(G, 0.1, [114,394], [10,400], [400,720]);

Code available on my neglected website:

http://imageprocessing.com.au/research/code/polarTransform.m http://imageprocessing.com.au/research/code/circularShortestPath.m http://imageprocessing.com.au/research/code/linearShortestPath.m