Connect bwlabeled components

[EDIT: 20110623 22:06 CDT - reformat - WDR]

Hello, I need help on a very simple task...well..seemingly.

I have a series of images similarly to this: http://blogs.mathworks.com/images/steve/94/labeling_labeled_objects_03.png

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,

Answer 1

[Sean de Wolski]

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)

[D]

Awesome, I now have 729 objects. Thanks for the sanity check.

Ok, so I now have the pixellist xy position of each object as well as the centroids, I am still learning how to plot stuff

would I just use implay? What function makes a layermask over the original stack?

Thanks guys.

[Sean de Wolski]

pixellist is the xy position of each PIXEL in that object; and yes the centroids as well.

You want to plot the centroid overlaid on the original image for each slice?

[D]

Yes, I want to overlay each cetroid onto the original image and be able to scroll through the stack. I just need every centroid's xy position per slice than I am going to group the centroids together eventually , ie every centroid that has one above or below the slice within a threshold is going to be grouped or counted as one 'object'. I want to export each "object" as a list of line segments or vectors that connect centroid to centroid through the slices.

[Sean de Wolski]

Okay. Well being able to scroll through the stack is going to be slightly more difficult. I assume you want the centroid to be an RGB color and not grayscale? In order to use IMPLAY for that, you'll have to move your data into the fourth dimension; convert the 3rd to rgb for each slice and then tune the colors at the centroids. This is a more complicated problem.

If you just want to view one at a time; the task is trivial. See this thread:

http://www.mathworks.com/matlabcentral/answers/7394-find-pixel-coordinates-value-of-a-centroid

For tracing centroids through the stack, I have a bit of code that does that; however, it works in 3D. Please clarify what you want exactly with bullets.

[D]

Thanks I'll start playing around with it. 2d was no problem like you said.

For the rest:

-Get every centroid in a binary stack of n slices

-For each centroid look above and below the current slice (if possible) and see if there is a centroid 'nearby' in xy on those slices(which means i need some thresholding variable)

-for each thresholded group of centroids compute the vector connecting them slice by slice ( so I have xy coordinates of each centroid, along with a vector from each centroid in one slice to the next linearly).

-group centroids together which fall into the threshold as one continuous object (to be able to count the number of objects in the stack)

[Sean de Wolski]

So subtract compute the distance formula for each centroid in every slice adjacency (1-2),(2-3). For each centroid in the smaller number of object image, find the minimal distance - these two are a match. Insure that the distance is below your "match" criteria; subtract each component to get the displacement vector components. Record somehow (cell array probably).

[D]

Thanks for the update. I got some nice quiver graphs. I am following your suggestions with the displacements

something like:

else if dim!=threshold

displacements=nan;

I am working on exporting the cell array stuff now,...

[Sean de Wolski]

No. What you'll want to do is not throw away the minimum value in the min expression

[~,idx] = min(..)

Save the min val and if it is more than the threshold strike it from displacements:

[minval,idx] = min(..)

Do the subtraction like you were going to:

and then

displacements(minval>threshold,:) = nan;

or similar.

[D]

OK makes more sense now,

subDis = bsxfun(@minus,xyDiff(:,:,1), xyDiff(:,:,2)); %the value of the distance subtracted from (1-2). correct?

[Sean de Wolski]

yes

[Sean de Wolski]

The first slice, subDis(:,:,1), is deltaX, the second slice, subDis(:,:,2), is deltaY. Then when you do the hypot operation in the next step it does the sqrt(x.^2+y.^2)

[D]

Thanks, you are very helpful, I want to send you a cake or something haha, so doing this:

subDis = bsxfun(@minus,xyDiff(:,:,1), xyDiff(:,:,2));

[minval,idx]= min(hypot(subDis(:,:,1), subDis(:,:,2)),[],dim);

if minval<threshold

displacements=...;

else if displacements(minval>threshold,:)=nan;

Another question to anyone, when I use regionprops and I get a two column lists of centroids xy is there a way to organize the output so that it is listed by slice # instead of a long stream of values without separation?

[Sean de Wolski]

No reason for if, that would require a for-loop and those are boring.

Use logical indexing just as I did above:

displacements(minval>threshold,:) = nan;

all done!

For your regionprops question, look at SORTROWS.

[D]

Thanks,

I think I did something wrong I get this return:

the intial file is [200 200 50] I find the centroids, and send it to these functions yet I get this return:

Index exceeds matrix dimensions.

Error in ==>

[minval,idx]= min(hypot(subDis(:,:,1),

subDis(:,:,2)),[],dim);

the inputs are:

xyDiff = bsxfun(@(x,y)abs(x-y),reshape(centsA,...[],1,2),reshape(centsB,1,[],2));

subDis = bsxfun(@minus,xyDiff(:,:,1), xyDiff(:,:,2));

As for the SORTROWS function:

[A, index] =SORTROWS(centroidA,i);& where i is the slice number?

[Sean de Wolski]

you don't need subDis, and you wouldn't need bsxfun anyway. just use diff(xyDiff,1,3) if you needed that matrix. If you name your 200x200x50 binary file, 'I', copy and paste the above with the added minval and displacements(minval>threshold,:) = nan lines this should run in its entirety with no issues. You can store the centroids of each slice in a second cell array if you want.

[D]

Thanks, seems to work now. I just need to organize the output data for some post-processing analysis.

[D]

Quick question, (can I post more here now that the main question is vaguely answered?)

Can I append existing arrays from other arrays with a header put into the updated array? Basically I am try [C,ii]=sortrows(centsA,centsB);

I want the array C's header on the first row to look like:

X1 Y1 X2 Y2 ... XN YN?

with X1 Y1 from centsA and X2 Y2 from centB and to continue on this until the loop ends.

basically in Sean's code I have two constantly updating columns CentA and CentB, which then iterate over the computation. I want to pull out this information before it is lost in the loop.

[Sean de Wolski]

Cwith_header = vertcat({['x' num2str(ii)],['y' num2str(ii)]},C);

And yes, it's okay to continue old threads but if they've been accepted few people will probably look.

[D]

So would this work:

A=sortrows(centsA);

B=sortrows(centsB);

C=vertcat(A,B);

Cwith_header = vertcat({['x' num2str(ii)],['y' num2str(ii)]},C);

[D]

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?

Answer 2

[Wolfgang Schwanghart]

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.

[Sean de Wolski]

Personally I'd recommend using bwconncomp and if necessary labelmatrix to get the label. bwlabeln (and label matrices in general) take up a lot of memory.

[Wolfgang Schwanghart]

Oh right. Somehow I thought bwconncomp works in 2D only, but I was wrong.

Answer 3

[D]

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?

[Sean de Wolski]

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

Answer 4

[D]

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.

[Sean de Wolski]

You want the orientation of a circle?

[D]

I would like the orientation of each labeled object relative to the plane (normal to it). Or if I can figure out how to get a 500x500 pixel slice through z and then take the orientation by flattening all objects onto the plane that is a possibility.

Also, quick question, can I export line segments (x and y coordinates) of each 3d labeled object?

[Sean de Wolski]

I don't understand your goal. The plane normal to a 2d image is axis into the third dimension.

xy coordinates are the 'pixellist' option in regionprops or you could use ind2sub on CC.PixelIdxList and CC.ImageSize from a connected components structure.

[Wolfgang Schwanghart]

If you have the three dimensional label matrix (L), I think it would be easiest to calculate the centroids of each object using regionprops for L(:,:,1) and L(:,:,end). It will then be easy to obtain vectors in 3d.

[D]

I'll give it a try. I'll post some code in a bit. Thanks everyone for the help.

[D]

Ok, so I got a 3d labelmatrix, L,

BW=bwconncomp(Stack);

L=labelmatrix(BW);

I took a look at the objects :

BW.NumObjects

it returned

2614

So I have 2614 objects which may or may not be connected together precisely. Is there a way to project the centroids onto a tiff stack as an overlay?

anyway to try Wolfgang's method

I then want to compute the centroids:

s=regionprops(BW,centroid) %in a for loop?

or

using labelmatrix ->

x1=L(:,:,1)...

xn=L(:,:,length(L))

I want to see if there is a difference in num of objects

[Sean de Wolski]

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