How can I color encode depth information of a greyscale 3-D image?

% Get useful information about the stack
[height, width, depth] = size(stack);

% Useful for Color encoding the stack in 3D
colorMap = jet(depth); % This is a 121x3 matrix of RGB values

% Convert each frame in the stack into a unit8 image with each frame corresponding to the colorMap values
stack_color = stack; 
for i = 1:depth
    maxVal = max(max(stack_color(:,:,i))); % Find the maximum value in the frame
    minVal = min(min(stack_color(:,:,i))); % Find the minimum value in the frame

    % Normalize the frame to the range [0,1]
    stack_color(:,:,i) = (stack_color(:,:,i) - minVal)/(maxVal - minVal);

    % Color encode the stack to a unit8 image
    stackRGB(:,:,1,i) = uint8(colorMap(i,1)*255*stack_color(:,:,i)); % Red channel
    stackRGB(:,:,2,i) = uint8(colorMap(i,2)*255*stack_color(:,:,i)); % Green channel
    stackRGB(:,:,3,i) = uint8(colorMap(i,3)*255*stack_color(:,:,i)); % Blue channel
end

% Create maximum intensity projections
stackMIP = max(stackRGB,[],4); % This is the MIP of the color encoded stack
MIP = max(stack,[],3); % This is the MIP of the grayscale stack

This is the code I've writen but I have an issue with it.

The output is this.

Dark intensity values seem to just go away entiirely. How can change this code, so that the brightness of each slice becomes a brightness for that color? Essentially I want to create something like the depth encoding feature in ImageJ.

[DGM]

What is the class of the variable stack?

How are you combining/displaying the two images?

[Aaron Woods]

I am using imshow()

Answer 1

[Walter Roberson]

Suppose you were to take the furthest slice and select a Hue for it and use the grayscale information as the Saturation, and use a constant Value, and image() that with an AlphaData property that is less than 1. Then likewise for the second furthest with a different hue, image() it on top of the previous image, and so on onto you have overlaid all of the images. Or, instead of constructing so many image() objects, you could mathematically calculate the result of overlaying the images and then image() only the final result. As a modification of the process you might do some kind of background detection on each slice and set the alpha information to 0 for the background locations.

Or instead of using image(), you could use a surface() with texture mapping; warp makes that convenient. If you do that with different Z values for the different slices, the result you get would be something you could rotate in 3 space. (You cannot rotate image() objects in 3-space, they will disappear.)

[Aaron Woods]

I'm excited to try the three methods you wrote, but can you provide a bit more with how I'd go about this?

Answer 2

[Image Analyst]

If the "stack" image is a volumetric image with zeros except where there is material, you can use min() to find the slice number of the first non-zero slice by setting zeros to inf and using

volumetricImage3d(volumetricImage3d==0) = inf;
[values, topSlices] = min(volumetricImage3d, [], 3);

Now it's simply an indexed image and you can apply a colormap as usual.

[Aaron Woods]

Can you explain this a bit more?

[Image Analyst]

Let's say you had a 2 voxel by 2 voxel by 3 slice image, with the top voxel in slice 1 in the upper left, the top voxel in slice 2 in the upper right, and the top voxel in slice 3 in the lower left,

slice1 = [130, 0; 
    0, 0];
slice2 = [0, 34; 
    0, 0];
slice3 = [0, 0; 
    190, 0];
volumetricImage3d = cat(3, slice1, slice2, slice3)

volumetricImage3d = 
volumetricImage3d(:,:,1) =

   130     0
     0     0


volumetricImage3d(:,:,2) =

     0    34
     0     0


volumetricImage3d(:,:,3) =

     0     0
   190     0

mask = volumetricImage3d == 0

mask = 2×2×3 logical array
mask(:,:,1) =

   0   1
   1   1


mask(:,:,2) =

   1   0
   1   1


mask(:,:,3) =

   1   1
   0   1

volumetricImage3d(mask) = inf

volumetricImage3d = 
volumetricImage3d(:,:,1) =

   130   Inf
   Inf   Inf


volumetricImage3d(:,:,2) =

   Inf    34
   Inf   Inf


volumetricImage3d(:,:,3) =

   Inf   Inf
   190   Inf

[~, topSlices] = min(volumetricImage3d, [], 3)

topSlices = 2×2

1 2 3 1

topSlices tells you which slice contained the first non-zero voxel value. Now topSlices has 1 in the upper left, 2 in the upper right, and 3 in the lower left. The lower right column was all zeros, so the "top" slice in that case is slice 1, though there actually in no voxel in that column with a non-zero value.