Connect bwlabeled components

8 views (last 30 days)
D
D on 24 Jun 2011
[EDIT: 20110623 22:06 CDT - reformat - WDR]
Hello, I need help on a very simple task...well..seemingly.
except I have a tiff stack with numerous layers in series. In this image the blobs move around as one goes deeper into the tiff stack series (it is a spatial not temporal series). I am trying to do the following:
Find all unique objects per layer, connect objects from one plane to another that 'overlap'(thus indicating that both object x1 in layer a and x2 in layer b are the same object just at different spatial depths) as a line segment, generate a output with orientation in z of each line segment.
I have done the first instant. I am looking for advice on how to go about the object grouping, ie. connecting the same 'real object' in the serialized sections.
Any functions I should look at? Particle tracking functions perhaps?
Cheers,

Accepted Answer

Sean de Wolski
Sean de Wolski on 27 Jun 2011
Using bwconncomp on a three dimensional image will group them in 3d. You'll have to call it on each individual slice to see how many objects per slice there are:
sliceObjs = zeros(size(stack,3),1);
for ii = 1:size(stack,3)
CC = bwconncomp(stack(:,:,ii));
sliceObjs(ii) = CC.NumObjects;
end
EDIT This is a script, with the engine completely generalized and ready to go. I saved you above image; cropped it and translated it for visualization error checking purposes. The translation was done using my function FEX:imtranslate
%%Data and stuff
I = imread('ans627.png'); %your image
I = I(35:276,86:381,1)==255; %keep relevant part; convert to binary
I(:,:,2) = imtranslate(I,[10 -14]); %artificially translate for testing
I(5:11,270:275,2) = true; %add an eleventh object with no friends
szI = size(I);
Ds = cell(szI(3)-1); %preallocate place to store stuff
CCold = bwconncomp(I(:,:,1)); %cc of first slice
RPold = regionprops(CCold,'centroid');
centsA = vertcat(RPold(:).Centroid); %extract centroids
for ii = 1:(szI(3)-1)
CCnew = bwconncomp(I(:,:,ii+1)); %cc of next slice
RPnew = regionprops(CCnew,'centroid');
centsB = vertcat(RPnew(:).Centroid);
dim = 2; %which dimension will control in the bsxfun expression?
if CCold.NumObjects>CCnew.NumObjects
dim = 1;
end
%Engine:
xyDiff = bsxfun(@(x,y)abs(x-y),reshape(centsA,[],1,2),reshape(centsB,1,[],2)); %Centroids distances in each dimension
[~,idx] = min(hypot(xyDiff(:,:,1),xyDiff(:,:,2)),[],dim);
if dim==2;
displacements = centsB(idx,:)-centsA;
H = quiver(centsA(:,1),centsA(:,2),displacements(:,1),displacements(:,2));
else
displacements = centsB-centsA(idx,:);
H = quiver(centsB(:,1),centsB(:,2),displacements(:,1),displacements(:,2));
end
RPold = RPnew; %get ready to move on
CCold = CCnew;
centsA = centsB;
Ds{ii} = displacements; %can store other stuff too.
end
%SCd 06/27/2011
If you now look at displacements, you'll see that the artificial translation was recovered. I didn't place a threshold on distance - just used the minimum one. This could be added easily enough. (replace objects without a match to nan to maintain order)
  20 Comments
D
D on 1 Jul 2011
In addtion, when I run your code above, I occasionally get the error on this line:
[~,idx] = min(hypot(xyDiff(:,:,1),xyDiff(:,:,2)),[],dim);
which generally says Expression or statement is
incorrect--possibly unbalanced (, {, or [.
Also I sometimes get :Subscript indices must either be real
positive integers or logicals on the displacement lines:
displacements = centsB(idx,:)-centsA;
My centsA and centsB are 258x2 and 272x2 respectively. I am working on padding the smaller of the two with nans (then convert to 0s when needed for some manipulation), is there slick way of doing this?

Sign in to comment.

More Answers (3)

Wolfgang Schwanghart
Wolfgang Schwanghart on 24 Jun 2011
I'd start with using the function bwlabeln to label the connected components in the 3d. This, however, requires them to never overlap or touch. If this occurs, you might want to use imerode prior to bwlabeln to let each object shrink a little.
  2 Comments
Wolfgang Schwanghart
Wolfgang Schwanghart on 27 Jun 2011
Oh right. Somehow I thought bwconncomp works in 2D only, but I was wrong.

Sign in to comment.


D
D on 24 Jun 2011
I would like to use bwconncomp, but I am running Matlab 7.5(R2007B).
Thus far I have this :
a=zeros([500 500 90], 'uint8'); % 90 frames deep
for i=1:length(a)
[Stack(:,:,i),map]=imread('file.tif',i);
BW=bwlabeln(a);
s=regionprops(BW,'Centroid');
end
That should find all the centroids in the volume right?
  1 Comment
Sean de Wolski
Sean de Wolski on 24 Jun 2011
Yes - in 3D.
If you want the output structure identical to bwconncomp (for storage) purposes and you have a label matrix, you can use my label2CC function.
http://www.mathworks.com/matlabcentral/fileexchange/30702-label2cc

Sign in to comment.


D
D on 24 Jun 2011
Thanks, I think I almost got it working...I guess I will use this space for some more questions on this code(I am still getting the basics in Matlab).
It seems to be the case that the 'orientation' function in regionprops only seems to work in 2D. Is there a 3D variant? Could I take an orientation of each labeled object as an orthogonal vector to the plane (in polar coordinates?)?
I have some more interesting questions coming up, more about morphology, which is still beyond me, but perhaps will interest the gurus. I will update with code as I figure at everything.
  7 Comments
Sean de Wolski
Sean de Wolski on 27 Jun 2011
length(..) returns the size of the _largest_ dimension. Therefore, don't use it, since it's unreliable and can change!
force SIZE to output the size of the dimension you want, in this case 3, size(L,3);

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by