how to fuse image

Dear all,

this is my code to view CT image by slice

P = zeros(256, 256, 72);
for K = 1 : 72
    petname = sprintf('I4%03d.dcm', K);
    P(:,:,K) = dicomread(petname);
end
imshow3D(P)

then, this is my code for view SPECT image by slice,

Noted: all my 512 slice SPECT image stored in one file.

[spect map]=dicomread('128x128');
info = dicominfo('128x128');
gp=info.SliceThickness;
spect=(squeeze(spect));%smooth3

aa=size(spect);aa=aa(3);

imshow3D(spect);

Anybody can help me to fuse both SPECT and CT images?

Answer 1

[Prathamesh]

Hi @mohd akmal masud

I understand that you have a code to view CT image by slice and then there is another code to view SPECT image by slice. So you are trying to fuse the CT and SPECT images.

So you can follow some steps which will help you fuse the CT and SPECT images.

Below are the steps:

1. Normalize both CT and SPECT images to the range [0, 1] using ‘mat2gray’.

2. If necessary, resize the SPECT image to match the dimensions of the CT image using ‘imresize3’.

3. Use different color maps for CT (e.g., grayscale) and SPECT (e.g., hot) for distinction.

4. Overlay the images using ‘imfuse’ with the 'blend' method.

Please refer below resources for more help:

• mat2gray documentation : https://www.mathworks.com/help/images/ref/mat2gray.html
• imresize3 documentation : https://www.mathworks.com/help/images/ref/imresize3.html

[DGM]

How exactly do you propose to apply different colormaps to two (ostensibly) volumetric images? There are various ways this could be done, some worse than others.

Since it's not clear how the images are intended to be colormapped, it's not clear whether normalizing the images is necessary.

What is the particular benefit of using imfuse() here, other than the fact that "fuse" is in the name?

In what way is #389248 relevant to this question? In that answer, I describe how to explicitly avoid using imfuse() for a composition involving binarized images.

I admit that "how do I fuse images" is terribly vague, but this still doesn't really answer the question.

Answer 2

[DGM]

I've posted other answers regarding "fusion". While it's a vague question, a fixed 50% weighted mean (the only "blend" imfuse() can do) doesn't ever seem like an appropriate answer.

For this example, I'm ignoring the fact that these are volumetric images in a particular format. I only assume that the data is arbitrarily-scaled intensity data. So long as the slices of the volume are oriented and scaled such that they are valid 2D images of compatible geometry, the volumes can simply be processed in a loop.

If the goal is to apply a colormap to such arbitrarily-scaled data, that's relatively easy. There are other ways of doing it, but I use MIMT gray2pcolor(), because it's convenient, and it can replicate the quantization used by the display tools. I've posted standalone copies elsewhere on the forum, but I'm not including one here, since I'm also going to use other MIMT tools.

As I said, a 50% opacity blend between two mostly-black images will just create an image which is unnecessarily dark. It's an average; of course it'll be dark if the image backgrounds are black. It's probably better to use some sort of unidirectional lightening blend mode. Again, I'm going to use MIMT imblend(), because that's what it's for. See the link below for a handful of examples of doing simple blends the manual way.

% these are both arbitrarily-scaled intensity images in uint16
A = imread('mrarb.png');
B = imread('tcarb.png');
imshow([A B]) % not normalized yet

% composing two colormapped images only makes an unreadable mess
% if we're keeping one as grayscale, no map needs to be applied
% but it still needs to be rescaled
A = mat2gray(A); % normalize to extrema

% the other one needs to be rescaled _and_ have an applied map
CT = magma(256); % some colormap (preferably one that starts at black)
B = gray2pcolor(B,CT); % this will do both for default settings

% combine them somehow. that's literally as specific as "fuse images"
outpict = imblend(B,A,1,'screen'); % complementary multiplication

imshow(outpict)

See also:

https://www.mathworks.com/matlabcentral/answers/1838343-fusion-of-colour-pet-mri-image-of-brain

EDIT: I forgot I used magma(), which isn't a built-in map either. It doesn't really matter, but there aren't many good built-in maps that start at black.

https://www.mathworks.com/matlabcentral/fileexchange/62729-matplotlib-perceptually-uniform-colormaps