Создание сферы (в виде маски) на трехмерном изображении

Question

Я использую 3D необработанное изображение (размер изображения PET (float) составляет [84*71*103] с размером интервала 4,07*4,07*3 мм).

Я бы хотел:

• Используйте координаты вокселей [116,177,86] которые соответствуют центру опухоли (определено мной), и создают сферу, которая является концентрической с центром опухоли.
• Измените диаметр сферы, чтобы она была немного больше, чем сама опухоль.

Таким образом, сфера будет служить маской, которая накладывается на изображение, и я могу получить статистику значений изображения в пределах этой сферической формы.

Надеюсь, я прояснил это достаточно. Буду очень признателен за вашу помощь. Мой код не работает правильно, и я не знаю, как поступить.

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

Я решил упрощенную 2D версию вашей проблемы, надеюсь, это поможет ответить на ваш вопрос. Для начала давайте создадим вход

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

Далее применяем обрезку:

[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])

Это тот результат, к которому вы стремились?