CODEWITHSUNDEEP
pythonnumpyimage-processingsimpleitk
David Hoffman
David Hoffman

Reputation: 2343

Does SimpleITK have broadcasting?

If I want to multiply a 3D volume plane by plane with a 2D slice in numpy I can use broadcasting:

import numpy as np
vol = np.random.rand(10, 20, 30)
slc = np.random.rand(10, 30)
new_vol = vol * slc[:, None]

If I try something similar in SimpleITK I get an error

import SimpleITK as sitk
vol_img = sitk.GetImageFromArray(vol)
slc_img = sitk.GetImageFromArray(slc[:, None])
new_vol_img = vol_img * slc_img

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-1-7d2c0160b591> in <module>
      9 vol_img = sitk.GetImageFromArray(vol)
     10 slc_img = sitk.GetImageFromArray(slc[:, None])
---> 11 new_vol_img = vol_img * slc_img

~\AppData\Local\Continuum\anaconda3\lib\site-packages\SimpleITK\SimpleITK.py in __mul__(self, other)
   4273     def __mul__( self, other ):
   4274         if isinstance( other, Image ):
-> 4275            return Multiply( self, other )
   4276         try:
   4277            return Multiply( self, float(other) )

~\AppData\Local\Continuum\anaconda3\lib\site-packages\SimpleITK\SimpleITK.py in Multiply(*args)
  50874
  50875     """
> 50876     return _SimpleITK.Multiply(*args)
  50877 class N4BiasFieldCorrectionImageFilter(ImageFilter_0):
  50878     """

RuntimeError: Exception thrown in SimpleITK Multiply: C:\b\3.6-64\ITK\Modules\Core\Common\src\itkDataObject.cxx:393:
Requested region is (at least partially) outside the largest possible region.

Upvotes: 0

Views: 372

Answers (1)

zivy
zivy

Reputation: 541

This cannot be done directly in SimpleITK because the concept of an image is not equivalent to an array of intensities, it has a physical spatial extent (see this read-the-docs explanation). The spacing, origin and direction of two multiplied images must be the same.

To do what you want you will need to iterate over the slices and then recompose back into a volume.

Here is sample code that does just that:

import SimpleITK as sitk

img = sitk.ReadImage('training_001_ct.mha')
slc = sitk.GridSource(outputPixelType=img.GetPixelID(), size=img.GetSize()[0:2], 
                             sigma=(0.1,0.1), gridSpacing=(20.0,20.0))
slc.SetSpacing(img.GetSpacing()[0:2])

modified_slices = []
for i in range(img.GetDepth()):
    current_img_slc = img[:,:,i]
    slc.SetOrigin(current_img_slc.GetOrigin())
    slc.SetDirection(current_img_slc.GetDirection())
    modified_slices.append(current_img_slc*slc)

sitk.Show(sitk.JoinSeries(modified_slices))

Please post future questions on the ITK discourse forum, and use the simpleitk tag.

Upvotes: 2

Related Questions