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

Asked by Joy on 28 Jun 2012
Latest activity Commented on by Joy on 11 Jul 2012

Hi everyone.

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

Joy on 30 Jun 2012

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

Walter Roberson on 30 Jun 2012

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 on 2 Jul 2012

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

## Products

No products are associated with this question.

Answer by Image Analyst on 30 Jun 2012

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 on 2 Jul 2012

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 on 2 Jul 2012

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 on 11 Jul 2012

I got it! Thanks!! :D

Answer by Walter Roberson on 3 Jul 2012

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.

## 1 Comment

Joy on 6 Jul 2012

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