Fuse (blend) two 2D images to RGB image in SimpleITK

Dear all,

I am trying to write a python code that extracts specific 2D slices from a 3D image and its segmentation then generates RGB image that shows both image and its segmentation in the extracted slice.
So far everything works but I am wondering if I can enhance the for loops part in the fusion section e.g. by numpy indexing or something similar.

import numpy as np, SimpleITK as sitk
fnmA = "img.png" 
readerImg = sitk.ImageFileReader() ; readerImg.SetFileName(fnmA) ;
img =        readerImg.Execute();
a = sitk.GetArrayFromImage(img);

fnmB = "seg.png"    
readerSeg = sitk.ImageFileReader() ; readerSeg.SetFileName(fnmB) ; 
seg =  readerSeg.Execute();
b= sitk.GetArrayFromImage(seg);

maxJ=a.shape[0] ; print("max J: " +str(maxJ))
maxK=a.shape[0] ; print("max K: " +str(maxK))
img2d = sitk.Image([c.shape[1],c.shape[0]], sitk.sitkVectorUInt8, 3)
for j in range(0,maxJ):
     for k in range (0,maxK):
          v= [int(a[j][k]),int(a[j][k]),int(a[j][k])]
          vc=[150,80,100]
          img2d.SetPixel(k,j,v)           
          if b[j][k]== 0:  
             img2d.SetPixel(k,j,vc)
          #endif
      #endfor k
 #endfor j
 writerImg = sitk.ImageFileWriter()
 fnm ="rgb.png" 
 writerImg.SetFileName(fnm) 
 writerImg.Execute(img2d)

rgb

Hello @ibr,

This is not the way to go. There are a variety of filters in SimpleITK that will allow you to combine an image and segmentation for visualization purposes, the most straightforward is LabelOverlayImageFilter. Please take a look at this Jupyter notebook which illustrates various filters and approaches for visualizing results from segmentation and registration using SimpleITK.

Hopefully one of these is what you are looking for. If not, let us know what is missing.

3 Likes

Thanks a lot. This was helpful. Here is the modified code:

import numpy as np, SimpleITK as sitk
fnmA = "img.png" 
readerImg = sitk.ImageFileReader() ; readerImg.SetFileName(fnmA) ;
img =        readerImg.Execute();

fnmB = "seg.png"    
readerSeg = sitk.ImageFileReader() ; readerSeg.SetFileName(fnmB) ; 
seg =  readerSeg.Execute();

pink= [255,105,180] ;     green = [0,255,0] ;     gold = [255,215,0]
rgb = sitk.LabelOverlay(img, seg,opacity=0.5, backgroundValue= -1.0, colormap=pink+green+gold)   

writerImg = sitk.ImageFileWriter()
fnm ="rgb.png" 
writerImg.SetFileName(fnm) 
writerImg.Execute(rgb)

rgb

2 Likes

Another related question. How can we use the same concepts for two different images? e.g. fusion of the fixed and registered images to evaluate the registration visually.

It seems LabelOverlay supports only label images.

Please take a look at the Jupyter notebook I pointed to in the previous response. It includes a variety of techniques for visualizing results including combining two scalar images which is what you are asking for (includes alpha blending, checkerboard, combining channels to create color image).

3 Likes

Thanks, I will try it.

I just checked.What I want is to make the first image has a specific color e.g. Magenta and the second image has a different color e.g. green. The composed image should show the registration error (it will appear as magenta or green in our example). Assuming the two images bellow:

imgF
imgM

using this line from the notebook in the link:

rgb3 = sitk.Cast(sitk.Compose(   img1_255, 0.5 * img1_255 + 0.5 * img2_255, img2_255), sitk.sitkVectorUInt8)

I get this result:
imgR

What I want, is something similar to these results (images are generated by 3D Slicer). Note the part 1 in yellow where both images match and part 2 in yellow where there is a small registration error.


When aligned pixels are gray it means that the three channels have the same/similar values, so something along the lines of

#assuming image pixel type is sitkUInt8, otherwise Cast
rgb3 = sitk.Compose(img1_255, img1_255, img2_255)

Essentially one of the images is repeated in two channels. Which image and which channel is up to your visual preference.

The code bellow produced the green and magenta images but when I combine them, I don’t get the desired result.

import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import SimpleITK as sitk
reader1 = sitk.ImageFileReader() ;  reader1.SetFileName(imgF.png") ; img1_255=reader1.Execute()
reader2 = sitk.ImageFileReader() ;  reader2.SetFileName("imgM.png") ; 
img2_255=reader2.Execute()           
zeros = sitk.Image(img1_255.GetSize(), img1_255.GetPixelID()) ; zeros.CopyInformation(img1_255)
rgb1 = sitk.Cast(sitk.Compose(   img1_255 * 0,     img1_255, 0 * img1_255), sitk.sitkVectorUInt8)
rgb2 = sitk.Cast(sitk.Compose(   img2_255    , 0 * img2_255, img2_255), sitk.sitkVectorUInt8)
rgb  = rgb1+rgb2 
writer1 = sitk.ImageFileWriter() ;  writer1.SetFileName("imgFc.png") ;  writer1.Execute(rgb1)      
writer2 = sitk.ImageFileWriter() ;  writer2.SetFileName("imgMc.png") ;  writer2.Execute(rgb2)      
writer  = sitk.ImageFileWriter() ;  writer.SetFileName("imgR.png")   ;  writer.Execute(rgb)      
img1 = mpimg.imread('imgFc.png');    imgplot = plt.imshow(img1) ;plt.show()
img2 = mpimg.imread('imgMc.png');    imgplot = plt.imshow(img2) ;plt.show()
img3 = mpimg.imread('imgR.png');     imgplot = plt.imshow(img3) ;plt.show()

What I am missing?

imgFc
imgMc
imgR

The issue with your example images is that the intensity levels of your two images are very different.

The first image has a few bright spots that are washing out the rest of the image. So the gray levels in the main structures of that image are in the 80-90 range.

If you look at your second image, the gray levels the 220-230 range. So if you put img2 in the red/blue channels and img1 in the green channel, the red/blue channels dominate. So it looks primarily magenta.

Here’s a little SimpleITK script I wrote to try and balance out the levels. Basically I scale the intensities of img2 by a factor of .4, to try and match the levels.

import SimpleITK as sitk

img1 = sitk.ReadImage("imgF.png")
img2 = sitk.ReadImage("imgM.png")

img2a = sitk.Cast(sitk.Cast(img2, sitk.sitkFloat32)*.4, sitk.sitkUInt8)

c2img = sitk.Compose(img2a, img1, img2a)

It sorts of works but could probably be better. I also tried using SimpleITK’s HistogramMatching filter, but that didn’t quite do what I wanted.

Basically you have to go back and get the same gray levels in the original image acquisition. Trying to match them afterwards doesn’t work that well. Fixing bad data as a post process is never as good as having good data from the start.

1 Like