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

Hi everyone.

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

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.

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

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.

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

Answer 1

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

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?

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);

I got it! Thanks!! :D

Answer 2

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.

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