Setting certain pixels in a grayscale image to RGB (red) for MIP

Hey again. I have a grayscale image and would like to highlight certain pixels by setting them to red, e.g.

save_array = image;
save_array(Omega==1) = %RED

Is there an easy way to do this? I will be taking the image (3D array), calculating the MIP (max intensity projection) by

mip_final = max(save_array, [], 3);

and would then like the red pixels to appear on the MIP (2D). So if a pixel is red anywhere in the z-dimension, it will appear red in the 2D MIP

Answer 1

[Image Analyst]

Try this:

redChannel = grayImage; % Initialize
greenChannel = grayImage; % Initialize channel for green and blue channel.
redChannel(Omega) = 255;
greenChannel(Omega) = 0;
rgbImage = cat(3, redChannel, greenChannel, greenChannel);

[Alex G]

I don't think so. I'm looking at the .tiff file itself btw

    [recovery,~] = HaLRTC(im, obs_voxels, alpha, 1e-6, maxIter, 1e-5);
    recovery = uint16(recovery);
    redChannel = recovery; % Initialize
    greenChannel = recovery; % Initialize channel for green and blue channel.
    redChannel(Omega==1) = 255;
    greenChannel(Omega==1) = 0;
    rgbImage = cat(3, redChannel, greenChannel, greenChannel);
    saveastiff(rgbImage,fullfile(save_path, 'movie'));

Do you think it could be an issue with saving/ImageJ? If you have a way to do this with a MIP, that'd be great for troubleshooting (right now I'm going through all 600 z-slices looking for a red dot haha)

Also, Omega is a binary tensor (but it seems like you knew that)

[Image Analyst]

What I did will set all pixels in a 2-D gray scale image to red wherever binary image Omega is true. There is nothing to project.

Your mip_final will not do that, You said your save_array is "a grayscale image". If you did convert it to an RGB image using my code, then using your code for taking the max along the 3rd (color channel) dimension will simply give you a gray scale image that will be pure white wherever the RGB image is red, or wherever any pixel in the green or blue channel is white. Not sure why you want this since it's basically just the same as skipping making the RGB image with red altogether and simply setting all pixels in Omega to white. I.e., same as

save_array(Omega) = 255;

So, not sure what you're really after. A screenshot of starting and desired final image would be good.

[Alex G]

Of course, I'll try to get something to you soon. I'll try to explain it now though:

We have a 3D grayscale array in matlab. We have a binary array Omega of the same dimensions. I would like to make a MIP of that grayscale array, but whenever Omega is true in the xy-plane (irrespective of z), I would like the corresponding pixel to be red

...from writing this I guess I could take a MIP of the image, a MIP of Omega, and where Omega's MIP is true would be red in the final MIP?

[Image Analyst]

OK, so you have a volumetric image, like from CT or MRI. So just mask the image

masked3DImage = grayImage3d; % Initialize
masked3DImage(~Omega) = -inf;  % Mask
masked3DImage = max(masked3DImage , [], 3); % Do MIP
masked3DImage(~Omega) = 0;  % minus infinity set to zero or whatever value you want outside of the mask.

[Alex G]

This looks great. Also this might sound stupid, but where do I set the selected voxels to red (that's the issue I'm trying to resolve)? As a reminder, (Omega == 1) are the selected voxels

I believe I can just run masked3DImage through your first comment's code?

Thanks again.

P.S. I'm getting a "Attempt to grow array along ambiguous dimension" error. Probably because we reference Omega after making it a MIP. The MIP is 2D and Omega is 3D, so perhaps this is an issue?

Answer 2

[Thorsten]

From the above discussion, I came up with this solution:

 % fake some data
  X = rand(10, 10, 23);
  N = prod(size(X));
  idx = randi(N, 1, round(0.01*N)); % 1% of all pixels are 1
  X(idx) = 1;

Omega = X == 1; % find the values that are 1 in X
idx = sum(Omega, 3) > 0; % project to a 2D binary image

MIP = max(X, [], 3);

% create a color image from MIP
R = MIP; R(idx) = 1;
G = MIP; G(idx) = 0;
B = MIP; B(idx) = 0;
I = cat(3, R, G, B);
imshow(I)

[Alex G]

Awesome, that works perfectly. I'll try implementing it myself. The only concern is that my data is 512x1024x100 and I have 0.1% of the pixels being "1" (or chosen to be red). This results in ~10%, 0.1% * 100 z-stacks, of the MIP to be red, which is less than ideal

Is it trivial to get the RED onto the original 3D projection and imshow that? This would allow me to get to what I want (I'm trying to visualize an active learning project and showing this many red pixels doesn't really help my case)

Thanks!

[Image Analyst]

How are you showing the original 3-D projection? Do you mean like the MIP variable, which is a 2-D image created by projecting the 3-D image along the third dimension, or do you mean like some kind of perspective/volumetric view of the 3-D volume like it's a 3-D rectangular block with red pixels located at various locations inside that block?

[Alex G]

The latter. I've been showing the MIP until now, but for this representation, I'd like to do some 3D visualization on 2D (i.e. the rectangular block with red pixels located inside).

For context, we have some 3D image that is changing at every run of some algorithm and I would like to highlight the pixels that are causing this change (these pixels are also changing every run and are located where Omega == 1). In the end I would like to make a movie with each run as a frame (~10 frames)

Hopefully this all makes sense. Thanks for your patience thus far

[Image Analyst]

What if you used scatter3() to put a red marker at every location of the max in each column?

[Alex G]

I had considered this, but it's crashing on my 32GB machine for a simple simulation:

A = rand(512, 1024, 100);
idx = find(A);
[X, Y, Z] = ind2sub(size(A), idx);
pointsize = 2;
scatter3(X(:), Y(:), Z(:), pointsize, A(idx));

I'd imagine it'd do the same when I get to larger image (3000x1000x600). I can try working on one of my 256GB cluster nodes if this is the only way.

I tried it without the color map, e.g.:

A = rand(512, 1024, 100);
idx = find(A);
[X, Y, Z] = ind2sub(size(A), idx);
scatter3(X(:), Y(:), Z(:));

and it runs, but it just shows me a solid blue cube/box as the figure

[Image Analyst]

Yes, of course it does, because the other pixels are opaque. For true volume visualization, you can't use MATLAB, which is limited to isosurfaces and cutaway views last I checked. You'll need to use a true volume visualization program like Avizo.

[Alex G]

Thanks for the advice, I'll try out Avizo and see if my PI can get it for the lab. I do see your point - even if we increase transparency, it'll be difficult to see the voxels at the center of any 3D volume when viewed in 2D.

I realized I need to implement a density-weighted approach to my adaptive sampling, so hopefully the MIP visualization will work better after that :) I'm sure my PI doesn't want to pay for another license

[Thorsten]

You can visualize only the red blocks using plotcube from the FileExchange:

 X = rand(10, 10, 23);
  % X = rand(3000,1000,600);

N = prod(size(X));
idx = randi(N, 1, round(0.01*N)); % 1% of all pixels are 1
X(idx) = 1;

[x, y, z] = ind2sub(size(X), idx);
for i = 1:numel(x)
  plotcube([1, 1, 1], [x(i) y(i) z(i)], .8, 'r')
  if i == 1, hold on, end
end
axis equal
grid on