creating a sphere (as a mask) on a 3D image

Question

I am using a 3D raw image (size of PET image(float) is [84*71*103] with spacing size of 4.07*4.07*3 mm).

I want to:

• Use the voxel coordinates [116,177,86] that correspond to center of the tumor (defined by me), and create a sphere that is concentric with the center of the tumor.
• Change the diameter of the sphere to be a little bigger than the tumor itself.

So the sphere is going to serve as a mask which is overlaid with the image, and I can get the statistics of the values of the image within that spherical shape.

I hope I have made it clear enough. I would really appreciate your help. My code is not working correctly, and I don't know how to proceed.

fid_1 = fopen('I:\PatientData\patient3\case1\DIR_CT1toCT2\New-     
crop\PET1FEM_PET2_DIFF.img','r','l');
pet1fem_pet2_diff = fread(fid_1,'float32');
pet1fem_pet2_diff = reshape(pet1fem_pet2_diff, [84 71 103]);
% interpolation is nearest neighbour, and 'crop' Makes output the same size as the input image
pet1fem_pet2_diff = imrotate(pet1fem_pet2_diff,90,'nearest','crop');

% create the image 
imageSizeX = 84;
imageSizeY = 71;
imageSizeZ = 103;
% columns are in the x direction and, rows are in the y direction
[columnsInImage rowsInImage pagesInImage] = meshgrid(1:imageSizeX, 1:imageSizeY,1:imageSizeZ);

% Next create the sphere in the image.
centerX = 29;
centerY = 26;
centerZ = 74;
radius = 5;

nonSphereVoxels = (rowsInImage - centerY).^2 ...
    + (columnsInImage - centerX).^2 + (pagesInImage - centerZ).^2 > radius.^2;
    pet1fem_pet2_diff(nonSphereVoxels) = 0;

figure(1);
imagesc(pet1fem_pet2_diff(:,:,30)); 
title('PET1FEM-PET2-DIFF');

Answer 1

I have solved a simplified 2D version of your problem, I hope that this helps to answer your question. For starters, let's create an input

A = rand(128,128); % this is your pet1_diff data
% you can plot it like this:
clf reset
imagesc(A)
daspect([1 1 1])

Next, we apply the cutoff:

[x,y] = meshgrid(1:128,1:128);
r2 = ((x-40).^2 + (y-40).^2) > 10^2; % radius 10 cutoff from position (40,40)
A(r2) = 0; % delete pixels outside the cutoff radius

% plot the filtered data
clf reset
imagesc(A)
daspect([1 1 1])

Is this the kind of result you were aiming at?