I wish to process around 2000 CT images. these are CT image of Material. I wish to find the size of 3D inclusion.

Hello,

I am new to MATLAB ...

I wish to do following things ....

1. Read around 2000 CT images of a metal specimen.
2. defect in metal are dark in color. I wish to calculate the 3d size of the defect.

Please guide me for the following steps:

1. How to read 2000 CT images ..
2. how to set the threshold intensity for identification of defect in the CT image
3. how to calculate the 3D size of the defect and the CT image in which the 2d area of the inclusion is maximum.

[Rik]

1. One image at a time. If your computer memory permits, load them all to a single 3D variable.
2. You can either determine a threshold manually, or you need to implement some method to automatically determine the threshold. Depending on your data you may need to do some filtering if a threshold depends on regional values (eg, in this array [1 2 3 1 3 4 5 6 7 5 7 8], both 1 and 5 could be considered low values, because of the surrounding values).
3. Once you have a logical array that is true for every voxel with the defect and false otherwise, sum(L(:)) will count the number of voxels. You can do a similar thing to create a vector of defect-voxels per slice and use the second output of max.

But the more pressing question: did you solve the issues you described in your previous question? Because as I mentioned there: without solving those, any questions about processing your data is moot.

[Vikas Verma]

Thanks alot Rik, for your prompt reply.

I am sorry for bothering you with my trivial questions. But I am trying to learn. I was able to resolve this by copying files to a different directory..... i have attached one image which i have to process ... the central circular part is the specimen.... and the black dots are the defects....As you can see, there are multiple defects in the specimen.

This image is just one slice of the volume.... these defects spread across multiple slices... i need to find the volume of each defect, the maximum area projected by the defect in Z direction.... and the location of the centroid of each defect...

clc;
clear all;
image_folder='E:\spring 8 2021 January\testing';
filenames=dir(fullfile(image_folder,'*.tif'))
total_images=numel(filenames);
for n=1:total_images
    f=fullfile(image_folder, filenames(n).name);
    our_images=imread(f);
    figure(n);
    imshow(our_images)
end

[Rik]

You're currently replacing the previous image by the next, instead of storing them in a 3D array. You're also opening 2000 figures. I don't see any indication that you've read and understood the demo Image Analyst provided you in the other thread. I suggest you have another look at that one.

[Vikas Verma]

1, I saw trhe lkink provided by image analyst in the previous thread. But i could not understand how to store images in a 3D array.

[Vikas Verma]

clc;

clear all;

image_folder='E:\testing';

filenames=dir(fullfile(image_folder,'*.tif'))

total_images=numel(filenames);

for n=1:total_images

f=fullfile(image_folder, filenames(n).name);

our_images(n,:,:)=imread(f);

figure(n);

end

[Vikas Verma]

I hope this correct ? @Rik

[Rik]

Generally you would store the slices as the third dimension, but that is indeed the general idea. You might also want to check whether the files aren't stored as color images (despite being grayscale). In that case you should read the images to a temporary variable and only store one color channel in your 3D array.

Also, there is no need to post 3 separate comments. You can edit your comments if you need to. Please also use the editor tools to format your code as code.

[Vikas Verma]

@Rik The following code i have prepared:

I am facing the following problems:

1. Should i do the contrast adjustmend of each image before importing them to a 3d Array. or after the formation of 3D array.
2. The regionprops gives the dimension like volume and other details of the defects. but is there anyother way to identify them so that it is easier to identify the value corresponds to whivch defect.
3. How can i find the Z-projection area of the defect. U explained in previous comment but i could not understand.
4. I need to ignore the defects below a certain size-projection area how can i modify the program to see only the defects above that size.

%To read the tif files and write the images to a 3D array 
image_folder='E:\vikas\testing';
filenames=dir(fullfile(image_folder,'*.tif'))
total_images=numel(filenames);
for i=1:total_images
    f=fullfile(image_folder, filenames(i).name);
    Image=imread(f);
    Volume(:,:,i)=Image; % storing images in to a volume(3D array)
end

%Thresholding and Binarizing the image
binarizeVolume=imbinarize(AdjVolume,"global");
%to get the data of the image segmentation
stats=regionprops3(binarizeVolume,"all")

[Rik]

1. Doesn't matter. As long as you're consistent.
2. The centroid property gives you the centroid of the region (so it contains an approximation of the center point). That should work well enough for the location.
3. If you use a labeled image as input (by using bwlabel on binarizeVolume), you can use the centroid to determine the label of a particular defect. Then you can use that on the result of bwlabel to make a logical array with only that defect. Using sum or max will give you a projection.
4. Loop through all defects to check the volume. If it is too small, mark it with false on the binarized image.

[Vikas Verma]

@Rik

The regionprops3 is giving me only one row table with very large volume which means only one defect which is incorrect. Can you please suggest me what is wrong.

1. For preprocessing of Image : I have following questions:

a. I have shown one of the slice in my previous comment. will "iamadjust" be suffiicient for contrast adjustment of each slice or 3D image.

b. How to define the connectivity for thresholding.

c. Can you please tell me "

1. If you use a labeled image as input (by using bwlabel on binarizeVolume), you can use the centroid to determine the label of a particular defect. Then you can use that on the result of bwlabel to make a logical array with only that defect. Using sum or max will give you a projection.
2. Loop through all defects to check the volume. If it is too small, mark it with false on the binarized image."

How to do these two steps ? please guide.

clc;
clearvars;
cimage_folder='E:\spring 8 2021 January\testing';
filenames=dir(fullfile(image_folder,'*.tif'))
total_images=numel(filenames);

for i=1:total_images
    f=fullfile(image_folder, filenames(i).name);
    Image=imread(f);
    
    Volume(:,:,i)=Image; % storing images in to a volume(3D array)
   % AdjVolume(:,:,i)=AdjImage; %storing adjusted Images in to a volume(3D array)
end

AdjVolume= imadjustn(Volume); %adjusting the contrast of the image

%Thresholding and Binarizing the image
binarizeVolume=imbinarize(AdjVolume,"global");


%Analysing the defects -finding size
stats=regionprops3(binarizeVolume,"all")

[Benjamin Thompson]

Can you post more data on this? It is hard for anyone else to follow this long thread. Are you looking at the binary 3D image output and comparing to the input image to see if it makes sense, or maybe your threshold needs adjusting in order to better segment the image?

[Vikas Verma]

i understand the thread has become long. can you please tell me what dat should i post ?

[Rik]

I sugested using bwlabel, which you haven't done yet. Actually, I should have said bwlabeln, as that is the multi-dimensional version. imadjust should not be doing anything, as you should control the binarization yourself. You want to mark certain values with true and others with false. What are those conditions exactly? What value range do you want to be 1?

You need to do this step by step. It is good to have the end goal in mind, but you need to focus on each small problem at a time.

You write regionprops3 only finds a single object. That means your binarization is not correct. Only after that does it make sense to explore how you can remove small defects.

[Vikas Verma]

@Rik

1. i agree with you there is some problem in binarization, Can you tell me how do i make pixel intensity below 40 as white and above 40 as black. what synatx should i use.
2. 2. the other problem which i am facing is that at the circumference of the material there is dark line, which is also turned black during binarization.
3. how do i remove this dark line at the circumference ?
4. How do i define connectivity of pixel in region props ?

I understand that i will have to use bwlabeln

[Rik]

BW=data<40;

Although you will have to adapt this, because you want to remove some values (that dark line you describe).

Did you read the documentation for the regionprops3 function? It will describe the options and give examples.

You also don't need to send me emails. I'm aware of this thread and I will respond when I have time and I have something to say.

[Vikas Verma]

@Rik

i have tried to find the Z projection area as shown in the code below. How do i modify the label for only the objects whose projection area is above certain limit. can you please tell me.

filenames=dir(fullfile(image_folder,'*.tif'));
total_images=numel(filenames)

Volume=zeros(2048,2048,total_images, "int8"); %predeifining the volume array to reduce processing time
for i=1:total_images
    f=fullfile(image_folder, filenames(i).name);
    Image=imread(f);
   %AdjImage=imadjust(Image); %Adjusting the contrast
    ImageBW=Image<65;
    Volume(:,:,i)=ImageBW; % storing images in to a volume(3D array)
    
    
end
%labeling the objects in volume
[L,n]=bwlabeln(Volume);
stats=regionprops3(L,"Volume","Centroid")
stats = sortrows(stats,'Volume','descend')


% Z projection calculator
zproj=zeros(n,total_images,"int8");
for j =1:n
    for k=1:total_images
        for m= 1:2048
            for n= 1:2048
                if(L(m,n,k))==j
                    zproj(j,k)=zproj(j,k)+1;
                    
                end
            end
        end
    end
end

for s=1:n
    Zprojectionarea(s) = max(zproj(s,:));
end
toc

[Rik]

The z projection can probably be replaced with a call to sum. (probably zproj=sum(L>0))

To remove certain objects, find the label of the background (probably 0), and then loop through all objects. If the volume is too small, mark it as background in your labeled volume.

[Vikas Verma]

@Rik

Thezproj=sum(L>0) is not working to calculate zprojection of each object in the label individually.

can you please guide

[Rik]

Of course that doesn't calculate a projection for each object separately, as it uses all labels. If you want a specific label, you need to select that label (e.g. sum(L==3);)

[Vikas Verma]

@Rik

sum(L==3) returns a 1x2048X7 array, when the total number of images is 7 and size of each image is 2048x2048.

how can i get numerical value from this. sum(L==3)

I tried using the "Image" property of regionprops3. It gives me the bounding box binary volume of the each object. from that it can be calculated.But i am not able to understand how to use the data from "Image" property column given by regionprops3.

1. Can u guide How can i write more data to stats table generated by regionprops3.
2. can i specify limits in regionprops3 to get the properties of object above a certain limit(such as volume above 500 or equivalentdiameter above 5 ?
3. Can i specify limits in bwlabeln to label only objects which are above certain size ?

[Rik]

Ah, then you need to specify the dimension sum must operate on:

L=randi(20,[2048,2048,7]);
projection=sum(L==5,3);
size(projection)

ans = 1×2

2048 2048

The easiest way to remove objects from regionprops is to determine your criterion for each of them, and mark as background.

L=randi([0 8],[30,20,7]);%small synthetic data
stats=regionprops3(L,{'EquivDiameter','Volume'})

stats = 8×2 table
    Volume    EquivDiameter
    ______    _____________

     531         10.047    
     464         9.6052    
     473         9.6669    
     474         9.6737    
     464         9.6052    
     454         9.5357    
     433         9.3863    
     479         9.7076    

for n=1:max(L,[],'all')
    if ... % remove if any of the conditions are true
            stats.Volume(n)<500 || ...
            stats.EquivDiameter(n)<5
        L(L==n)=0;%mark as background
    end
end
stats=regionprops3(L,{'EquivDiameter','Volume'})

stats = 1×2 table
    Volume    EquivDiameter
    ______    _____________

     531         10.047    

[Vikas Verma]

@Rik

thanks for your response. I understand this logic. it would have worked well but i had to remove the L=bwlabel(Volume) command Because L was 2048x2048x2048 and my computer ran out of memory on running the code.

I directly applied regionprops to Volume which was 3D binary Image.

1. is there a method of labeling before regionprops without running out of memory.

What i was thinking of first modifying the stats table(obtained by regionprops3) ( by deleteing rows based on table properties) and then labeling as per the modified stats table. I have alread yadded the zprojection to each row of the table(zprojection area of each object)

1. Can you please guide me how can i modify stats table(obtaines by regionprops3) based on the table properties. (i wish to delete certain rows based on the one or two attributes of the table).

[Rik]

The output is a table, which is one of the basic data types in Matlab. What basic Matlab tutorial did you try?

Also, there might be some othere error. Do you indeed have 2048 slices, 2048*2048 pixels each? That is a huge array. The data I'm used to work with (medical data) tends to be a factor 1000 smaller in total voxel count.

[Vikas Verma]

@Rik

i tried the MATLAB onramp and Image processing onramp tutorial.

The problem with table operation is following:

The stats table generated by regionprops3 contain one column "Image" in which each cell is a 3D array of the bounding box of the object. like(200X312x81) , (120x234x102) ........

I wish to delete the row in which the Image has 3d array above certain dimensions

i wish to delete row in which the image has bounding box with x>150 or y>200.

How to label the objects in 3D array which have been shortlisted in the stats after deleting rows from stats

[Rik]

You need to choose:

1. Delete the rows in the table (without modifying the label array)
2. Modify the label array

The last request ("How to label the objects in 3D array which have been shortlisted") is mutually exclusive with what you posted before ("i had to remove the L=bwlabel(Volume) command Because L was 2048x2048x2048 and my computer ran out of memory").

The code below shows you how to modify the label array by modifying the conditional statement. You can edit the table regionprops3 returns just like every other table you have encountered in the Onramp course.

rng(1)%set the seed to get repeatable results for random functions
L=randi([0 8],[30,20,7]);%small synthetic data
stats=regionprops3(L,{'EquivDiameter','Volume','Image'})

stats = 8×3 table
    Volume          Image          EquivDiameter
    ______    _________________    _____________

     451      {30×20×7 logical}       9.5146    
     510      {30×20×7 logical}       9.9127    
     436      {30×20×7 logical}        9.408    
     444      {30×20×7 logical}       9.4652    
     491      {30×20×7 logical}        9.788    
     485      {30×20×7 logical}        9.748    
     474      {30×20×7 logical}       9.6737    
     447      {30×20×7 logical}       9.4864    

for n=1:max(L,[],'all')
    if ... % remove if any of the conditions are true
            stats.Volume(n)<500 || ...
            size(stats.Image,1) > 150 || ...
            size(stats.Image,2) > 200 || ...
            stats.EquivDiameter(n)<5
        L(L==n)=0;%mark as background
    end
end
stats=regionprops3(L,{'EquivDiameter','Volume','Image'})

stats = 2×3 table
    Volume          Image          EquivDiameter
    ______    _________________    _____________

       0      { 0×0    logical}            0    
     510      {30×20×7 logical}       9.9127