0

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

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

1 回答 1

2

这不能直接在 SimpleITK 中完成，因为图像的概念不等同于强度数组，它具有物理空间范围（请参阅此阅读文档说明）。两个相乘图像的间距、原点和方向必须相同。

要执行您想要的操作，您需要迭代切片，然后重新组合回一个卷。

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

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

请在ITK 讨论论坛上发布未来的问题，并使用 simpleitk 标签。

于 2019-02-07T15:46:45.890 回答