How to get the intensity values (data stored) that is stored in the voxels (3D image) at the certain distances from the origin of the sphere?

Hello, I have a 3D medical image with data values stored in its each voxel. I would like to pull out the data values of image at certain radial distances. The origin of my image is (0, 0, 0) and radius of the sphere is 20 mm. I would like to generate 19 concentric shells of equal thickness 1 mm (so total is 20 mm) in between center of sphere and its surface. At each concentric shell, I am interested to extract the value of the data stored in 3D image. How to write simple Matlab code for this situation? Has anyone done this before? Please provide me some ideas and suggestions, may be a little piece of code as well.

[Walter Roberson]

Are these DICOM images? If so then the first thing to do is to pull information about the distance between centers out of the headers. It should not generally be assumed to be the same in all directions, as it is not unusual in some modalities for the Z spacing to be different than the X or Y spacing.

[blues]

Thank you Walter for your reply. Image is in MetaImage (.mhd and .raw) format. The header files looks as follows:

ObjectType = Image
NDims = 3
BinaryData = True
BinaryDataByteOrderMSB = False
CompressedData = False
TransformMatrix = 1 0 0 0 1 0 0 0 1
Offset = -0.156024 -0.156024 -0.156024
CenterOfRotation = 0 0 0
ElementSpacing = 0.034672 0.034672 0.034672
DimSize = 10 10 10
AnatomicalOrientation = ???
ElementType = MET_FLOAT
ElementDataFile = dose1e6events-Dose.raw

[Walter Roberson]

It looks like the ElementSpacing has the relevant information about pixel sizes.

Answer 1

[Walter Roberson]

https://www.mathworks.com/matlabcentral/fileexchange/29344-read-medical-data-3d

filename = ...;
info = mha_readheader(filename);
V = mha_read_volume(info);

spx = info.ElementSpacing(1);
spy = info.ElementSpacing(2);
spz = info.ElementSpacing(3);
%careful, first dimension is y not x.
sx = size(V,2);
sy = size(V,1);
sz = size(V,3);
cx = sx/2; cy = sy/2; cz = sz/2;
xv = ((1:sx) - cx) * spx;
yv = ((1:sy) - cy) * spy;
zv = ((1:sz) - cz) * spz;
[X, Y, Z] = ndgrid(xv, yv, zv);
D = sqrt(X.^2 + Y.^2 + Z.^2);

rvals = 0 : 19;
rthick = 1;
Nr = length(rvals);
shells = cell(Nr, 2);
for ridx = 1 : Nr
  r = rvals(ridx);
  mask = (D >= r & D < r+thick);
  shells{ridx, 2} = mask;
  shells{ridx, 1} = V(mask);
end

Now the first column of the cell array shells contain the extracted data values, and the second column contains the mask of locations that were extracted (in case you need it later.)

[blues]

I put my file name: acbd_efh.mhd in place of filename and when I ran it, it says undefined variable abcd_efh or class abcd_efh.mhd. Am I missing something? And how can I store the output array in a .txt file? Sorry for my stupid question!

[Walter Roberson]

filename = 'acbd_efh.mhd';

You will need to define the output format for the text file.

When you do, keep in mind that you are extracting at non-rectangular coordinates, so the result of extracting the pixels will just be a list of values, ordered by linear index but (unless you specifically command) not identified as to where they came from. Just a whole bunch of numbers, with a different number of them for each of the different shells, with the larger radii tending to have more associated values.

[blues]

Thank you. I got data values in shells first column, second column is the mask of locations. Now I would like to extract the average pixel data values covered by the radius of say 1 mm spherical shell, 2 mm spherical shell, ...up to 20 mm (i.e the radius of sphere). I want to pull out [(average pixel data value, radial distance 1 mm), ( average pixel data value, radial distance 2 mm ), ( , )........ up to 20 sets, so that I could make a plot and visualize it. Is there a simplistic way to do it?

[Walter Roberson]

avgpixelvals = cellfun(@mean, shells(:,1));
plot(rvals, avgpixelvals);

[blues]

This avgpixelvals is just the mean values of pixel data stored in whole sphere, not at each concentric shells, right?

But I want to extract average pixel values at each concentric shells. In brief; average pixel values of data in a spherical shell of radius 1 mm, average pixel (it should be cumulated value) values of data in a spherical shell of radius 2 mm and so on, up to 20 mm. Hope I describe my question clearly.

Thank you very much for your great help.

[Walter Roberson]

No, avgpixelvals is the mean for each concentric shell. shells is a cell array in which the first column contains only values that lie within the particular concentric sphere for that one radius.

Answer 2

[blues]

One more baby query: my avgpixelvals min and max values are 5*10^-11 and 8*10^-8 (in y-axis), so I tried to set axis([0 120 5.^(-11) 8.^(-8)]); but this is not working. How to set this y-axis?

[Walter Roberson]

Your graph appears to show a y value peaking about 8.4*10^-8

Anyhow, you used the wrong numeric syntax

axis([0 120 5E-11 8E-8])

or

axis([0 120 5*10.^(-11) 8*10.^(-8)]);

With 5E-11 being only about 1/1600 of 8E-8, you should expect that MATLAB will probably round the lower bound to 0. After all, you would need your y resolution to be at least 1600 for it to make a 1 pixel difference.