Detecting and isolating part of an image

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

% 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!

• The images look the same? Which is the the bit that must be removed? The long bit thin part with the blob on the end? Can you post an example of your code? If you are just looking for the large oval section I would use circular shortest path (CSP) with the centre at the centre of the image. CSP is also useful for segmenting the white parts (bones?) if you can also find a point near the centre. – geometrikal Jun 10 '14 at 14:13
• I edited my post to show my code for isolating the leg so far. Yes, I am looking to get rid of the long thin part with the blob in that example. Unfortunately my images are fairly different from each other, and some of them don't have nice oval shapes. Would circular shortest path still work for those? I have found a nice way to segment out the bones using edge detection. Thank you for your help! – user10167 Jun 11 '14 at 8:41
• while implementing this code , its showing this error: Undefined function or variable 'polarTransform'. Error in circularShortestPath (line 17) P = polarTransform(M, centrePointV, radiusV, stepsV); Error in isolation (line 12) [ path, energy, lastIndex, pathImage ] = ... what to do? – user21250 May 30 '16 at 7:48

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.

• This. Although I would probably do opening operator i.e ersion then dilation to maintain the bit you want to keep the same size – nivag Jun 11 '14 at 13:18
• This is what I was going to post as well. Also, solidity may be a good metric to use if that connection point is sometimes larger. – Ed S. Jun 11 '14 at 21:19

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));
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:

• +1 for your awesome explanation for the new concept. I myself will try it for several other problems. :) – shreelock Jun 13 '14 at 6:45