How to extract rib cage for finding region of interest ?

Hello, Im trying to find lung roi from raw CT image. I cannot detect rib cage with bwarefilt() to extract because it checks just 2d, so cannot combine that info with connectivity. So I felt regionprops3 will help me but im not sure how to use regionprops3 to extract it, Im open to any alternative ways to extract that rib cage.

-no lung tissue should be lost

error in line 36:

Subscripting into a table using one subscript (as in t(i)) or three or more subscripts (as in t(i,j,k)) is not supported. Always specify a row subscript and a variable subscript, as in t(rows,vars).

my code is

source_path_raw = 'C:\Users\rawimage';
V_raw = niftiread(source_path_raw);

%to initialize
x=size(V_raw,1);
y=size(V_raw,2);
z = size(V_raw,3);
BW = zeros(x,y,z);

slice_num = size(V_raw,3);
V_raw = imbinarize(V_raw);
for i=1:slice_num
one_slice = V_raw(:,:,i);
radius =4; % radius 4,decomposition 0  veri kaybetmemek için en optimali
decomposition = 0;
se = strel('disk',radius,decomposition);
BW(:,:,i) = imerode(one_slice, se);
end
BW = imfill(BW, 'holes');
V_BW = bwareaopen(BW,100);

 volumeViewer(V_BW(:,:,:));
% Attached creenshot from here

conn=18;
n=3;

cc = bwconncomp(V_BW,conn);
if (cc.NumObjects == 0), error('%s: error - input array is empty', mfilename); end
% Getting more information about the connected-components
rp = regionprops3(cc, "Volume");

Cent = rp(1).Volume;
   crow = round(Cent(2));
   ccol = round(Cent(1));
   cplane = round(Cent(3));
   
% Sorting by size, largest first
[~, ind] = sort([rp.Volume], 'descend');
% Keeping only the n largest connected-components
rp = rp(ind(crow, ccol, cplane));
% Setting relevant indices in the output array
Xout = false(size(V_BW));
n = min(n,length(rp));
for i=1:n, 
    Xout(rp(i).PixelIdxList) = true; 
end
% Keeping only the statistics of the n largest connected-components
rp = rp(1:n);

volumeViewer(Xout(:,:,:))

Answer 1

[Matt J]

I cannot detect rib cage with bwarefilt() to extract because it checks just 2d

A 3D (and higher) version of bwareafilt is available here

https://www.mathworks.com/matlabcentral/fileexchange/90762-further-tools-for-analyzing-objects-in-n-dimensional-images

[Berkem Teksüz]

No there is no built-in function even in 2021-a version ?

help bwlareafiltn

--- bwlareafiltn not found. Showing help for bwareafilt instead.

even I looked at bwpropfiltn function that other you sent ,

bwlmaskpropfiltn(BWL,mask,attrib,varargin)

also here I do not have a mask parameter because Iam already trying to create lung tissue mask. If you want I can upload raw image with wetransfer which is anonymised.

[Matt J]

It really didn't occur to you to click on the link I gave? You weren't the least bit curious?

[Berkem Teksüz]

I'm already talking about that link's functions ,

this is just empty function?

function BWL=bwlareafiltn(BWL,varargin)

%The syntaxes below follow the same rules as in bwareafilt(), except BWL

%can be N-dimensional and can also be a label image

%

% BWL2 = bwlareafiltn(BWL,range)

% BWL2 = bwlareafiltn(BWL,n)

% BWL2 = bwlareafiltn(BWL,n,keep)

BWL=bwlpropfiltn(BWL,'Area',varargin{:});

also bwlpropfiltn not found. too

[Matt J]

No, it's not empty. As you can see from the last line, it calls another function in the package.

[Matt J]

Berken Teksuz's comment moved here:

Thank you I was not aware of it, I set your four functions in the same folder, but I was taking this error,

Unrecognized function or variable 'undefined'.
Error in bwlpropfiltn (line 67)
  if ~undefined('range')
Error in bwlareafiltn (line 8)
  BWL=bwlpropfiltn(BWL,'Area',varargin{:});
Error in lungRoi (line 24)
BWL=bwlareafiltn(V_BW,5);

so I just changed "%do sorting" part of bwlpropfiltn function to [~,idx]=maxk(propdata,n); because I just want first or second biggest blob in the case of where lung tissue bigger, however this extracts me full image, not biggest, as ı understand from bwareafilt documentaiton, this must not happen, also works for n smallest volumes(shows smallest particles true), but does not work for largest ?

Briefly, although I set biggest blob(n=1) ==> bwlareafiltn(V_BW,1);

it gave me full image.(ss2)

I hope it is clear,

THANK YOU!

Edited files are added.

[Matt J]

Yep, those were bugs. I have fixed them, however, and uploaded corrected versions to the File Exchange.

[Berkem Teksüz]

I think, you did not updated because there writes 1May also, cannot find biggest blob, and same with yesterday's code.

I still cannot solve this issue,

Unrecognized function or variable 'undefined'.

Error in bwlpropfiltn (line 71)
  if ~undefined('range')

Error in bwlareafiltn (line 9)
  BWL=bwlpropfiltn(BWL,'Area',varargin{:});

Error in bwlareaopenn (line 6)
  BWL=bwlareafiltn(BWL,[P+1,inf]);

typo,

function BWL=bwlareaopenn(BWL,P)
%The following syntax, similar to bwareaopen, is supported.
%
%BWL2=bwareaopenn(BWL,P)

  BWL=bwareafiltn(BWL,[P+1,inf]);

you missed bwlareafiltn (l)

[Matt J]

It still is May 1 at the MathWorks. In any case, I've now uploaded again and downloaded to double-check that my changes were checked in.

[Berkem Teksüz]

Im still trying to mask biggest volumed object(lung tissue extracted rib cage) so I made boundary test and find that after the value of 15 takes full image, untill 15 finds smallest elements whats wrong?

BWL=bwlareafiltn(BW,15,'smallest'); (look at screenshots)

source_path_raw = 'C:\Users\90507\Desktop\calisilacak_dosyalar\S19_002\4 Unnamed Series';
V_raw = niftiread(source_path_raw);

%to initialize
x=size(V_raw,1);
y=size(V_raw,2);
z = size(V_raw,3);
BW = zeros(x,y,z);

slice_num = size(V_raw,3);
V_raw = imbinarize(V_raw);
for i=1:slice_num
one_slice = V_raw(:,:,i);
radius =4; % radius 4,decomposition 0  veri kaybetmemek için en optimali
decomposition = 0;
se = strel('disk',radius,decomposition);
BW(:,:,i) = imerode(one_slice, se);
end
BW = imfill(BW, 'holes');


BW = bwareaopen(BW,50);
%BW = bwlareaopenn(BW,100); % Remove objects containing fewer than 100 pixels using bwareaopen function.

%BWL2=bwlmaskpropfiltn(BWL,mask,attrib,varargin); % göğüs kafesini
%bulduğunda maskeyi ham resimden kesişimle çıkarman için

BWL=bwlareafiltn(BW,15,'smallest');

volumeViewer(BWL(:,:,:));

Thank you very much Matt!

[Matt J]

I'm not going to be able to understand what's going on without access to V_raw. Regardless, how are you using bwareaopen on a 3D BW?

[Matt J]

Also, why are you reoding in a loop, slice by slice? imerode() can take a 3D image as input.

[Matt J]

Berkem Teksüz commented:

You are right about imerode() but,

At the stackoverflow, a highly voted comment says that,

The matlab function bwareaopen allows to remove the groups of %pixels that are smaller than the number you specify, and it works in 3D.

an example below,

% Remove white regions with less than 10 pixel
 my3DimageNEW = bwareaopen(my3Dimage, 10);

Works in 3D, what do you mean by "how are you using?" it works?

My code,

source_path_raw = 'PATH';
V_raw = niftiread(source_path_raw);

%to initialize
x=size(V_raw,1);
y=size(V_raw,2);
z = size(V_raw,3);
BW = zeros(x,y,z);

slice_num = size(V_raw,3);
V_raw = imbinarize(V_raw);

radius =4; % radius 4,decomposition 0, optimum for not to lost any data from lung tissue!
decomposition = 0;
se = strel('disk',radius,decomposition);
BW = imerode( V_raw, se);

%BW = imfill(BW, 'holes');
BW = bwlareaopenn(BW,16); % Remove objects containing fewer than 100 pixels using bwareaopen function.
BWL=bwlareafiltn(BW,2);
volumeViewer(BWL(:,:,:));

V_raw's link, C is Covid, S is healthy, both anonymised, you can download freely:

https://we.tl/t-z7bopdZ2K1

Thank you, I really need this, sorry for taking your time.

[Matt J]

It would be easiest if you would just attach a .mat file containing the binarized V-raw. Even a 512^3 logical array wouldn't create a very big .mat file.

[Berkem Teksüz]

Wisely :)

Thank you for suggestion.

Attached.

[Matt J]

This seems to do a decent job of isolating the lungs. I assume the ribcage won't be very far away...

load(websave('t','https://www.mathworks.com/matlabcentral/answers/uploaded_files/606075/V_raw_binarized.mat'));

E=entropyfilt(V_raw)>0.5;

E=imclose(E,strel('sphere',5));

BW=bwlareafiltn(E,1);

V_proc=(E>0).*BW.*V_raw;

imshow(V_proc(:,:,128))

[Berkem Teksüz]

WORKS! THANK YOU!

[Berkem Teksüz]

Does imerode in a loop provides me an airway with around 6-7 generations, is it possible, or do you have any file for this?

I saw that answer but did not work for me.

[Matt J]

I don't know what "6-7 generations" means...

[Berkem Teksüz]

Each three branch untill to final point, increases generation+1

Same as the attached blue airway branches.

[Matt J]

So you want to keep only the airways that branch at least 18 times?

[Berkem Teksüz]

No just I'm trying to find as possible as I can go to deeper levels of tree, do not care about number. Generation means just how detailed an airway's branching. Lung segmentaiton literature has a big problem about detecting deeper levels of airway tree because airway wall is getting thin for each generation so no one can detect it further 10th generation, for years. Indeed they exist up to 26th generation!. To solve this, my project's first stage is coming until , max of 7-8 generation or whatever it is close to it. After that my algorithm will start to find them by my anatomical patterns samples. In conclusion, that blue image's a bit more complex airway tree works fine for me.

Thanks wise man!

[Matt J]

I'm afraid I don't understand the task. You mean you are trying to count the number of branches in different airway paths?

In any case, it doesn't sound like a simple application of erosion. think it would be a good idea to start a new thread where you can more fully explain the new task and where the answers offered by people can focus on that.

[Berkem Teksüz]

Okay I will open a new question, We do not need to count anything. its for just a sign for developer to stop, its done. I just need an airway segmentation as complex as that blue one.Thats all.Thank you!