SimpleITK 有广播吗？

Question

如果我想在 numpy 中将一个 3D 体积平面与一个 2D 切片相乘，我可以使用广播：

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

如果我在 SimpleITK 中尝试类似的操作，我会收到错误消息

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.6-64\ITK\Modules\Core\Common\src\itkDataObject.cxx:393:
Requested region is (at least partially) outside the largest possible region.

Answer 1

这不能在 SimpleITK 中直接完成，因为图像的概念不等同于强度数组，它具有物理空间范围（参见 this read-the-docs explanation）。两张相乘图片的间距、原点、方向必须相同

要执行您想要的操作，您需要遍历切片，然后重新组合成一个卷。

下面是执行此操作的示例代码：

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

请 post 以后的问题 ITK discourse forum，并使用 simpleitk 标签。

SimpleITK 有广播吗？

Does SimpleITK have broadcasting?

python

numpy

image-processing

simpleitk