How to calculate Maximum Intensity Projection(MIP) of a 3D image ????

Hi everyone.

Does anyone know how to calculate MIP of a 3D image????

[Ryan]

I know this is not in the spirit of the forum, but I believe the free NIH software, ImageJ, has a feature to visualize 3D datasets using MIP.

[Joy]

But I need to use Matlab to do it. Do you know how? Thanks for your suggestion anyway.

[Walter Roberson]

Would this be similar to Projection Persuit ? If so then it becomes a global minimization over a very large search space, which takes a very long time to calculate.

[Joy]

Oh no. This is too complicated. I just need to know how to calculate a simple 3D image's mip that's all. I found my answer. Thanks for the help though :D

[supermanu]

I still do not know how to create a maximum intesity projection using matlab. Could you please exaplin?

[Image Analyst]

Did you see my answer? It has code. If that's not good for you, then start a new question and attach your images.

[Cosmas Mwikirize]

For an image axbxc, MIP=max(image,[],dim) where dim is the dimension along which you want to calculate the MIP (either 1,2 or 3). This result is similar to what you would get with "viewer3d".

[Walter Roberson]

That does not give Maximum Intensity. See my Answer, where I give code that calculates intensity based upon RGB.

[Mohamed Amine Henchir]

@Walter Roberson, I couldn't find the code in the link you provided for projection pursuit.

Answer 1

[Image Analyst]

This may sound obvious, but did you look at the max() function? It can take a 3D array and give you the max along any dimension you specify.

% Generate sample 3D array
A = rand(2,2,2)
% Get maximum intensity projection.
mip = max(A, [], 3)

In the command window:

A(:,:,1) =
    0.6352    0.9889
    0.4384    0.9632
A(:,:,2) =
    0.2460    0.6582
    0.4045    0.9417
mip =
    0.6352    0.9889
    0.4384    0.9632

[Joy]

Hi

Thanks for the answer. It's very useful. However, may i ask you, how come the output dimension is only 2D and not 3D after the mip function has been performed?

I have an image which is interpolated, and it's an color image. I've interpolated the image into 3 seperated channels (red, green, blue). And the output dimension is 128x170x91. So if I perform mip to each channel, how do i combine them back into one image and if the output of the mip resulted in 2D, what should I do so that I will be able to visualize the full 3D image?

As for the combining back into 1 color channel, can i use the cat() function to do it?

[Image Analyst]

A projection will reduce the dimensionality by 1. How could it not? If you have a 3D color image, that's a 4D image. You need to do it on each color channel at a time, then combine again.

RGBImage = cat(3, redChannel, greenChannel, blueChannel);

[Joy]

I got it! Thanks!! :D

[Walter Roberson]

A Maximum Intensity Projection over each of the channels individually and recombined is not the same as a Maximum Intensity Projection over the entire array. Intensity is a property of the three channels combined. "Intensity" for color images is the same as "brightness", and the various bit planes do not contribute equally to brightness.

When using (0 to 255) as the intensities of the channels, the maximum that the blue channel can contribute to the intensity, when the blue value is 255, is 29.1840. The red channel only needs to be 49.71720614 to have an equal contribution to the brightness.

The code I provided in my answer finds, for each pixel location, the slice that has the maximum brightness, and extracts the R, G, and B that went into making that brightest pixel in order to create the intensity projection.

Answer 2

[Walter Roberson]

I will assume for this discussion that your RGB image stack is arranged as (X, Y, RGB, Z) and that it is named Stack

grayStack = squeeze( 0.2989 * Stack(:,:,1,:) + 0.5870 * Stack(:,:,2,:) + 0.1140 * Stack(:,:,3,:) ); %X, Y, Z after squeeze

[maxgray, maxidx] = max(grayStack, [], 3);  %along the 3rd axis (Z)

nX = size(Stack,1);
nY = size(Stack,2);
[X, Y] = ndgrid( 1 : nX, 1 : nY );

MIP_R = Stack( sub2idx( size(Stack), X(:), Y(:), 1, maxidx(:) ) );
MIP_G = Stack( sub2idx( size(Stack), X(:), Y(:), 2, maxidx(:) ) );
MIP_B = Stack( sub2idx( size(Stack), X(:), Y(:), 3, maxidx(:) ) );

MIP = cat(3, reshape(MIP_R, nX, nY), reshape(MIP_G, nX, nY), reshape(MIP_B, nX, nY) );

image(MIP);

There are performance tweaks that can be done.

The code would be slightly different if the RGB is the 4th dimension instead of the third.

The first line of the code is effectively doing an rgb2gray() but for all of the image layer simultaneously. This conversion calculates the intensity (brightness) of each pixel in each slice. In the discussion with IA, the conversion from RGB to brightness (intensity) is not done, so the code there is not doing a Maximum Intensity Mapping. You should not be calculating the maximum R along the stack and the maximum G along the stack independently: you need the R, G, and B information combined at each pixel in order to calculate intensity there.

With intensity in hand, the code finds the maximum intensity along the Z, and the Z slice number that corresponds to the brightest point.

Once the slice number is done, there is some code that, in a vectorized way, pulls out the R, G, and B pixel values of the appropriate Z layer. And once it has those, it combines the three channels into a single RGB image.

[Joy]

I've figured it out myself. But anyway thanks :D

[JoaquinB]

Hello Walter, How would it be for a grayscale image. Im doing it with the max function, but it does give a lot of white pixels and the image is really saturared compared to ImageJ Z proyection. Am i skipping any type of math procedure? Thanks!

Im doing Z proyection every 50 stacks (my image dimensions are 280,1000,4000).

PD: I do have an error on s=41 because of crop dimensions.

STACK_pass=50;
INPH_filt_pass=0:STACK_pass:4000; 
INPH_filt_pass(1)=1;
NUM_stacks= length(INPH_filt_pass);
BOT_filt=zeros(280,1000,NUM_stacks);


    for s=1:1:(NUM_stacks-1) 
    
    BOT_stack=imcrop3(BOT_crop,[1,1,INPH_filt_pass(s),999,279,(INPH_filt_pass(s+1)-1)]);
    BOT_max=max(BOT_stack,[],3);
    BOT_filt(:,:,s)=BOT_max;
    
    end 

Answer 3

[navid alemi]

Hi Joy

would you please let me know how you figure it out I have the same problem.

thank you in advance

[Image Analyst]

Is there something about the max() function you don't understand?

[Walter Roberson]

I gave complete code in my answer. The method Joy followed of doing each color plane independently and then combining the results does not give the maximum intensity because the different color planes do not contribute equally to intensity.