Connectivity in binary 3D matrix

I am using bwconncomp to find the largest connected component in a 3D binary matrix (volume). However, it doesn't work the way I want. Specifically, the code is

      conncomp = bwconncomp(volume, 6);    

      % Identify the largest component using cellfun.
      [~, maxcell] = max(cellfun(@numel, conncomp.PixelIdxList));

      % Zero the image and assign to it the largest component.
      volume = zeros(size(volume));
      volume(conncomp.PixelIdxList{1, maxcell}) = 1;

Take a look at the figure below. Imagine that these are just profiles of the "slices" in the 3rd dimension. At the top is the original volume, at the center is what I get using this code, and at the bottom is what I want. bwconncomp skips a part of the large 3D object by cutting it in the slice where it is first split in two (slice 7 from left). I want to preserve that part because it is part of the object, as shown at the bottom. Changing connectivity doesn't help.

The person who solves this is my hero! Any suggestions? Many thanks.

Edit: it seems that when I try this on small 3x3x3 matrices it works fine, but not on my large images. So I also attached an example TIFF stack and the analysed one to see what I mean. For example, slices 10 and 12 in the analysed image should but do not include some connected areas.

Answer 1

[Image Analyst]

I guess I'll be your hero. There is (now) a built-in function for this. It's called bwareafilt(). It came in the last few versions so hope you have a recent version. If not, you can try my attached function that I used for older MATLAB versions when they did not have bwareafilt(), but I have not tested it on 3D binary images, just 2D images.

[Xen]

Hey Image Analyst, that's a nice code, thanks. Actually my original code works fine as well, I was just being stupid and used it in an incorrect way and I didn't notice that immediately because it worked as expected in most cases. I didn't know about bwareafilt; this might be useful somewhere.

[Sara Salimi]

Thanks for the lines of code that you have shared, it worked for me too. I have one question, what if we want to take the 3 largest connected components in the 3d volume?

How should we refine the code

[Image Analyst]

There seems to be no bwareafilt3() yet, so you'll have to do a workaround. First find all the blob volumes with regionprops3() then sort them to find the 3 biggest ones. Then use bwareaopen() which does seem to take a 3-D binary image to filter out all but the 3 biggest.

props = regionprops(binaryImage3d, 'Volume')
sortedVolumes = sort([props.Volume], 'descend')
thirdBiggest = sortedVolumes(3);
% Pull out 3 biggest into another 3-D logical image.
binaryImage2 = bwareaopen(binaryImage3d, thirdbiggest);