Gray Matter Area and Volume Calculations

Hello, I am working on gray matter area and volume calculation of brain MRI image. I have some doubts with my code below:

% Area Calculation
Pixel_Numbers = sum(J(:));
DistanceinUnit = 1.796875; %Based on pixel spacing (in mm?)
DistanceinPixel = sqrt(x^2+y^2);
DistancePerPixel = DistanceinUnit/DistanceinPixel;
Area = (Pixel_Numbers)*(DistancePerPixel^2); %Is it in mm^2 or mm^2/pixel?
% Volume Calculation (Is the formula correct?)
LengthinPixel = length(find(J(:)==1)); 
Length = LengthinPixel*DistanceinUnit;
Volume = Area*Length;

The information about pixel spacing is 1.796875\1.796875. I'd be glad if you can help me to figure out the answer. Here is the picture. Thank you.

[Image Analyst]

What are the doubts? Does (the badly-named) J include that white frame around the image?

[Qvna Lhyvnav]

Sorry for the badly-named image :D No, it doesn't. I attach the new one here.

The doubts are:

1. Did I define for "DistanceinUnit = 1.796875" based on the pixel spacing information in MRI image correctly?
2. My script was based on yours (spatial calibration script), but I confused about the area unit since DistanceinPixel = pixel; DistancePerPixel = mm/pixel; Area = pixel*(mm^2/pixel^2) = mm^2/pixel. Do I define this correct or I miss something? I want to have a mm unit for area.
3. I also want to calculate the volume, but was my script correct (area x length)? I have no idea to find the height information.

Looking forward to your reply. Thank you.

Answer 1

[Ryan Comeau]

Hello, i think you're image appears to be binarized already. This means we can obtain the area well with the regionprops function. https://www.mathworks.com/help/images/ref/regionprops.html.

In terms of the volume, are you wanting to sum up the number of pixel in that area * intensity? In you case if the image is binary, the volume and area will be the same. volume=area*height but height is 1 because your image is binary.

The last thing you want now is pixels to mm? just determine the conversion for this, if a pixel is 1.796875x1.796875 we will have an area=(1.796875)^2*num_pixels.

If all the area is one block, it will be returned from regionprops, if not you'll need to sum it up. Here is my logic for this problem:

image=imread('J.jpg');
image=rescale(image,0,1); %helps for binarization
BW_im=imbinarize(image); %see documentation, it is versatile
hh=regionprops(BW_im,{'Centroid','Area',}); %centroid for locations.

total_area=0;
total_volume=0;
height=1;%your image is binary, if it's not this will change pixel by pixel.
%see the regionprops function to sum up pixel values as well. 
for i=1:length(hh)
    total_area=total_area+hh(i).Area;
    total_volume=total_volume+(hh(i).Area*height)
end

area_in_mm=(1.796875.^2)*total_area;

Hope this helps, if not, drop a comment.

RC

[Qvna Lhyvnav]

Hello, RC.

Yes, my image already binarized and your code totally works for me! But, I have some questions about the height. As I read on How to Return the Height of a Bounding Box?, the height is defined:

x = imread('coins.png')>100;
bb = regionprops(x,'BoundingBox')
bbMatrix = vertcat(bb(:).BoundingBox)
bbMatrix(:,3); %height

I applied the code to yours and I am stuck in converting the cell into structure array. I used cell2struct.

image=imread('J.jpg');
image=rescale(image,0,1); 
BW_im=imbinarize(image); 
hh=regionprops(BW_im,{'Centroid','Area','BoundingBox'}); 
total_area=0;
total_volume=0;
height = vertcat(hh(:).BoundingBox);
height = num2cell(height(:,3)); %cell2struct?
for i=1:length(hh)
    total_area=total_area+hh(i).Area;
    total_volume=total_volume+(hh(i).Area*height)
end
area_in_mm=(1.796875.^2)*total_area;

Do you have any idea about it? I also would like to know about what is your reference of the height for binary image is 1? Thank you before.

[Ryan Comeau]

Hello, in your code above, the bb is a struct and the line that has:

bbMatrix = vertcat(bb(:).BoundingBox)

Is taking a structure and transforming in into a matrix.

sidebar: in MATLAB, there are matrixes, cell matrices and structures.

This matrix is created because there is a hidden index which is not seen during the vertcat() operation. In this operation, both the structure array and the field of the structure array are indexed. In other words, the field is a vector. so for example, we can have bb(1).BoundinBox(3) or bb(3).BoundingBox(1). In this case, the height seems to be the third element of the bounding box, it's not the 3D height it's the 2D height. If you want the height of this bounding box, you don't actually need to remove it from the struct array, you could simply change the code as follows:

image=imread('J.jpg');
image=rescale(image,0,1); 
BW_im=imbinarize(image); 
hh=regionprops(BW_im,{'Centroid','Area','BoundingBox'}); 
total_area=0;
total_volume=0;
for i=1:length(hh)
    total_area=total_area+(hh(i).BoundingBox(3)*hh(i).BoundingBox(4));%change is here
    total_volume=total_volume+(hh(i).Area*1)
end
area_in_mm=(1.796875.^2)*total_area;

Please note, that the bounding box may return the height of the whole region on average(i don't know, you'll need to check the documentation) not the pixel by pixel height

A binary image is defined by have either 1 or 0. This is why the height is 1 for everything. If you had a greyscale brain scan, we could sum the height of each individual pixel to obtain a volume.

Hope this helps,

RC

[Qvna Lhyvnav]

I understand your explanation about the bounding box, but I don't understand why you put hh(i).BoundingBox(4) into the total_area. My logic is like this:

total_area=total_area+hh(i).Area;
total_volume=total_volume+(hh(i).Area*hh(i).BoundingBox(3)) %height is the 3rd element of the bounding box

Would you like to elaborate your code?

I also have another idea after read can anyone inform me on how to calculate the volume of the segmented region?. I added the slice thickness information, so the code would be like this:

image=imread('J.jpg');
image=rescale(image,0,1); 
BW_im=imbinarize(image); 
hh=regionprops(BW_im,{'Centroid','Area'}); 
voxel_numbers = 0;
for i=1:length(hh)
    voxel_numbers = voxel_numbers+hh(i).Area;
end
Pixel_Spacing = [1.796875;1.796875]; Slice_Thickness = 4;
Volume_Per_Voxel = Pixel_Spacing(1)*Pixel_Spacing(2)* Slice_Thickness;
Volume_in_mm3 = Volume_Per_Voxel*voxel_numbers;

Does my code make sense? Thank you before, RC.