How to measure the thickness of an object having inner & outer irregular circle like shapes

I have segmented left ventricle in a MRI image using Random-Walker as shown in below image. Now I want to measure the thickness of Myocardium i.e. the gap between the inner & outer shape as shown in the following figure. Can someone please help me with this. Thanks in advance for your help! http://img690.imageshack.us/img690/9152/output1.jpg

[Matt J]

How would you define "thickness" seeing as the gap between the inner and outer shape does not have a uniform cross-section?

[Prashant Kashid]

Actually this is the point of concern. I have pointed out the area between the right ventricle & left ventricle (outlined in red shapes) which is somewhat uniform compare to rest of the thickness. but at least is there a way out to measure this gap between the inner & outer shapes at the pointed location in an image. Thanks in advance!

[Walter Roberson]

Look at the very top of the inner area. Should the "thickness" there be the distance to the outer area directly above it ("north" of it), or should the "thickness" there be the distance between the point and the nearest point in the outer circle ?

[Prashant Kashid]

Yes you got it right. As the inner area is not intersected with the outer one; the "Thickness" is considered as the distance between any point "A" of inner area & immediate point "B"(nearest point to "A") on outer area.

Answer 1

[Matt J]

You might also define an average thickness as follows

thickness = mean(min(interdists(A,B),[],2));

where the columns of the matrix A are the coordinates of the points on the outer boundary and B are the points on the inner boundary. This requires my interdists() utility function below.

function Graph=interdists(A,B)
%Finds the graph of distances between point coordinates
%
% (1) Graph=interdists(A,B)
%
% in:
% 
% A:  matrix whose columns are coordinates of points, for example 
%          [[x1;y1;z1], [x2;y2;z2] ,..., [xM;yM;zM]]
%     but the columns may be points in a space of any dimension, not just 3D.
% 
% B:  A second matrix whose columns are coordinates of points in the same
%     Euclidean space. Default B=A.
% 
% 
% out:
% 
% Graph: The MxN  matrix of separation distances in l2 norm between the coordinates. 
%        Namely, Graph(i,j) will be the distance between A(:,i) and B(:,j).
%
%
% (2) interdists(A,'noself') is the same as interdists(A), except the output 
%     diagonals will be NaN instead of zero. Hence, for example, operations
%     like min(interdists(A,'noself')) will ignore self-distances.
%
% See also getgraph

noself=false;
if nargin<2  
    B=A; 
elseif ischar(B)&&strcmpi(B,'noself')
    noself=true;
    B=A;
end

N=size(A,1);
B=reshape(B,N,1,[]);

Graph=sqrt(sum(bsxfun(@minus, A, B).^2,1));

Graph=squeeze(Graph);

if noself
    n=length(Graph);
   Graph(linspace(1,n^2,n))=nan; 
end

[Image Analyst]

I think though that this requires that the number of coordinates in B is the same as in A, right (it's not really obvious)? He may not have that. The inner boundary may have fewer points than the outer boundary. Would this still work in that case?

[Matt J]

No, it is not required that B and A have the same number of coordinates. The expression

min(interdists(A,B),[],2);

will be a vector of min. distances of length size(A,2), the number of coordinates in the outer boundary. This will be true regardless of the size of B.

Obviously, the definition is not really symmetric - you would get a different average thickness if you reversed the roles of A and B, but hopefully both boundaries are well-sampled enough to make the difference negligible.

Or, you could make the definition symmetric by doing

    G=interdists(A,B);
    tBA=min(G,[],1);
    tAB=min(G,[],2).';
    thickness = mean([tBA,tAB]) ;

[Image Analyst]

I see. It would be interesting to compare the values you get from this with the values from my two methods.

Say, could interdists be used to generate flexible envelopes, like you'd get from alpha shapes, concave hull, restricted convex hull, balloons, snakes, active contours, that kind of thing? You could have one contour be the actual shape, and the other be the convex hull, and generate a curve that is x% of the way in between, so you could tailor (by varying x) how much the envelope "hugs" the original shape.

[Matt J]

Sounds like you could do that something like the following

 [~,idx]=min(interdists(A,B)[],2);
 envelope=x*A+(1-x)*B(:,idx);

However, for large numbers of points, you would probably want to use John D'Errico's IPDM which is probably better optimized than interdists(). E.g., you can get the nearest points without computing the entire distance matrix.

Answer 2

[Image Analyst]

What I would do, since you already have the coordinates of the inner and outer boundary (which is the really tough part), would be to get the area of the zone and divide by the length of the zone. You can do this:

1. Use polyarea() on both boundaries to get both areas, or use regionprops if you have a binary image instead of boundary coordinates
2. Subtract the areas to get the area of the "in between" zone
3. Determine the perimeter of both boundaries. There are a variety of ways such as regionprops()
4. Take the average of those two perimeters. This can be considered as the average length of the midpoint of the zone.
5. Consider the zone as a rectangle with a length of the perimeter you just calculated and an area you just calculates. So the average width is simply the (zone area) / (zone perimeter).

Would that figure of merit work for you?

An alternative way is to take double the mean value of the skeleton of the Euclidean distance map. "Huh?" you say? What's that? Well let's take it one step at a time. The Euclidean distance map is the distance from some part of a binary shape to the nearest background pixel. So it's zero at the edge, and it is maximum along the "spine" or "skeleton" of the shape because that's where it's farthest from the background. So right there in the middle (on the spine or skeleton) you have basically the radius or half the distance to the edge, which is half the thickness. So to get the thickness you double that number. So now you need to do this along every point in the skeleton and you will have an array of thicknesses. Take the mean of those thicknesses and that is your average thickness of the zone. So, basically it's this:

1. Get a binary image of only the annualar zone
2. Call bwdist() on the zone to get the EDM
3. Call bwmorph(, 'skel') to get the skeleton
4. Multiply the skeleton by the EDM image (mask it)
5. Extract the non-zero pixels of that masked image
6. Take their mean and multiply by 2 to get mean thickness.

They're both fairly simple. You might do it both ways and see which you like best. They should both be fairly similar in values.