How do I resolve this issue with 3D image visualization?

Question

I have this 3D image data that I need to visualize. I have been able to visualize it with 2D slices using imshow3D, but I would like to see the image data in 3D space.

The code I used is as follows (courtesy: How do i create a rectangular mask at known angles?), but I can't tell why it isn't displaying:

% create input image
imageSizeX = 120;
imageSizeY = 200;
imageSizeZ = 50

% generate 3D grid
[columnsInImage, rowsInImage, pagesInImage] = meshgrid(1:imageSizeX, 1:imageSizeY, 1:imageSizeZ);

% create the sphere in the image.
centerY  = imageSizeY/2;
centerX  = imageSizeX/2;
centerZ  = imageSizeZ/2;
diameter = 56;
radius   = diameter/2;

sphereVoxels = (rowsInImage - centerY).^2 ...
    + (columnsInImage - centerX).^2 + (pagesInImage - centerZ).^2 <= radius.^2;


% change image from logical to numeric labels.
Img   = double(sphereVoxels);
for ii = 1:numel(Img)
    if Img(ii) == 0
        Img(ii) = 2;  % intermediate phase voxels
    end 
 end



% specify the desired angle
angle = 60;                     

% specify desired pixel height and width of solid
width  = imageSizeX;    
height = imageSizeY;
page   = imageSizeZ;

% Find the row point at which theta will be created
y = centerY - ( radius*cos(angle * pi/180) ) 

% determine top of the solid bar
y0 = max(1, y-height); 

% label everything from y0 to y to be = 3 (solid)
Img(y0:y, 1:width, 1:page)=3;   
% figure, imshow3D(Img);
% axis on;
% grid on;


% display it using an isosurface 
fv = isosurface(Img, 0);
patch(fv,'FaceColor',[0 0 .7],'EdgeColor',[0 0 1]);  title('Binary volume of a sphere');
view(45,45);
axis tight;
grid on;
xlabel('x-axis [pixels]'); ylabel('y-axis [pixels]'); zlabel('z-axis [pixels]')

Although, the solid bar is not diagonal as the figure attached below, I would expect the image to be something similar to this:

I do not know exactly what I am doing wrong here.

Answer 1

With regard to the problem in your code, it appears that you set points inside the sphere to 1, then set all the remaining points outside the sphere to 2, then a section through the y plane to 3. There is no value of 0 in the volume in this case, so trying to get an isosurface at the value of 0 isn't going to find anything.

However, if you'd rather create a "voxelated" Minecraft-like surface, like in your sample image showing the facets of your voxels, then I have another option for you...

First, I created a set of volume data as you did in your example, with the exception that I omitted the for loop that sets values to 2, and instead set the values of the solid bar to 2.

Next, I made use of a function build_voxels that I've used in a few 3D projects of mine:

function [X, Y, Z, C] = build_voxels(roiMask)
  maskSize = size(roiMask);

  % Create the ROI surface patches pointing toward -x:
  index = find(diff(padarray(roiMask, [1 0 0], 'pre'), 1, 1) > 0);
  [X1, Y1, Z1, C1] = make_patches([-1 -1 -1 -1], [1 1 -1 -1], [-1 1 1 -1]);

  % Create the ROI surface patches pointing toward +x:
  index = find(diff(padarray(roiMask, [1 0 0], 'post'), 1, 1) < 0);
  [X2, Y2, Z2, C2] = make_patches([1 1 1 1], [-1 -1 1 1], [-1 1 1 -1]);

  % Create the ROI surface patches pointing toward -y:
  index = find(diff(padarray(roiMask, [0 1 0], 'pre'), 1, 2) > 0);
  [X3, Y3, Z3, C3] = make_patches([-1 -1 1 1], [-1 -1 -1 -1], [-1 1 1 -1]);

  % Create the ROI surface patches pointing toward +y:
  index = find(diff(padarray(roiMask, [0 1 0], 'post'), 1, 2) < 0);
  [X4, Y4, Z4, C4] = make_patches([1 1 -1 -1], [1 1 1 1], [-1 1 1 -1]);

  % Create the ROI surface patches pointing toward -z:
  index = find(diff(padarray(roiMask, [0 0 1], 'pre'), 1, 3) > 0);
  [X5, Y5, Z5, C5] = make_patches([1 1 -1 -1], [-1 1 1 -1], [-1 -1 -1 -1]);

  % Create the ROI surface patches pointing toward +z:
  index = find(diff(padarray(roiMask, [0 0 1], 'post'), 1, 3) < 0);
  [X6, Y6, Z6, C6] = make_patches([-1 -1 1 1], [-1 1 1 -1], [1 1 1 1]);

  % Collect patch data:
  X = [X1 X2 X3 X4 X5 X6];
  Y = [Y1 Y2 Y3 Y4 Y5 Y6];
  Z = [Z1 Z2 Z3 Z4 Z5 Z6];
  C = [C1 C2 C3 C4 C5 C6];

  function [Xp, Yp, Zp, Cp] = make_patches(Xo, Yo, Zo)
    [Xp, Yp, Zp] = ind2sub(maskSize, index);
    Xp = bsxfun(@plus, Xp, Xo./2).';
    Yp = bsxfun(@plus, Yp, Yo./2).';
    Zp = bsxfun(@plus, Zp, Zo./2).';
    Cp = index(:).';
  end
end

This function accepts a 3D matrix, ideally a logical mask of the volume region(s) to create a surface for, and returns 4 4-by-N matrices: X/Y/Z matrices for the voxel face patches and an index matrix C that can be used to get values from the volume data matrix for use in coloring each surface.

Here's the code to render the surfaces:

[X, Y, Z, C] = build_voxels(Img > 0);
rgbData = reshape([1 0 0; 1 1 0], [2 1 3]);
hSurface = patch(X, Y, Z, rgbData(Img(C), :, :), ...
                 'AmbientStrength', 0.5, ...
                 'BackFaceLighting', 'unlit', ...
                 'EdgeColor', 'none', ...
                 'FaceLighting', 'flat');
axis equal;
axis tight;
view(45, 45);
grid on;
xlabel('x-axis (voxels)');
ylabel('y-axis (voxels)');
zlabel('z-axis (voxels)');
light('Position', get(gca, 'CameraPosition'), 'Style', 'local');

And here's the plot:

Note that the sphere and bar surfaces are colored differently since they are labeled with values 1 and 2, respectively, in the volume data Img. These values are extracted from Img using C and then used as an index into rgbData, which contains red (first row) and yellow (second row) RGB triplets. This will create an N-by-1-by-3 matrix of polygon face colors.