Is it possible to get a writable array view of a SimpleITK image?

Beginner Questions
, ,
1

Following up on this question: https://stackoverflow.com/questions/54557586/does-simpleitk-have-broadcasting/54577094

If GetArrayViewFromImage returned a writable view then you could use numpy broadcasting without issue. Is there a way to do this?

2

Just to emphasize, ignoring the spatial concept of an image in SimpleITK is not recommended. Having said that, you can still do what you want (would recommend that you add a comment in your code warning yourself that this is a dangerous thing to do):

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

img_view = sitk.GetArrayViewFromImage(img)
slc_view = sitk.GetArrayViewFromImage(slc)

new_vol = img_view*slc_view[None,:]

sitk.Show(sitk.GetImageFromArray(new_vol))
3

Unfortunately, that won’t work as I’m trying to avoid memory copies of large arrays.

4

Not sure that I understand the issue. Is the problem with the creation of new_vol instead of the overwriting of img_view? The only other copy is in sitk.GetImageFromArray(new_vol)).

5

Yeah the GetImageFromArray is the copy I’m worried about. The issues is some of the images a 50 GB+ so I’d like to make as few copies as possible.

6

Got it. Unfortunately, I don’t have a good solution off the top of my head. If I have an idea how to get around this will post it here.

7

You can look at the implementation of GetImageFromArray:

We have maintained a very conservative approach to managing array and references, which makes the interfaces quite stable. If you look at the implementation you can see it creates a new Image then fills it. I think we could very easily take an image as an argument and then “fill” it with the nparray. This would both be more efficient and very safe.

Would this help your situation?

8

When you say “fill” do you mean that both the ndarray object and the sitk.Image object would be looking at the same location in memory?

9

By fill I do mean copy. If you reuse the input image, you’ll still create the output numpy array but you can then copy back to the same location as the input.

10

Yes that would work for me.

11

Can you please create a feature request on GitHub to add an “outputImage” to the GetImageFromArray method.

I think the method could look something like this:

def GetImageFromArray( arr, isVector=None, outputImage=None):
“”“Get a SimpleITK Image from a numpy array. If isVector is True, then the Image will have a Vector pixel type, and the last dimension of the array will be considered the component index. By default when isVector is None, 4D images are automatically considered 3D vector images. If outpuImage is set then it will contain the output and the outputImage’s buffer may be reused. “””

if not HAVE_NUMPY:
    raise ImportError('Numpy not available.')

z = numpy.asarray( arr )

if isVector is None:
  if z.ndim == 4:
    isVector = True

if isVector:
  id = _get_sitk_vector_pixelid( z )
  if z.ndim > 2:
    number_of_components = z.shape[-1]
    shape = z.shape[-2::-1]
  else:
    number_of_components = 1
    shape = z.shape[::-1]
else:
  number_of_components = 1
  id = _get_sitk_pixelid( z )
  shape = z.shape[::-1]

if (outputImage and
    outputImage.GetSize() == shape and
    outputImage.GetNumberOfComponentsPerPixel() == number_of_components and
    outputImage.GetPixelIDValue() == id):
  img = outputImage
  print("Reusing outputImage's buffer!")
else:
  # SimpleITK throws an exception if the image dimension is not supported
  img = Image( shape, id, number_of_components )

sitk.SimpleITK._SetImageFromArray( z.tostring(), img )

return img


You might have to import some other things from SimpleITK to get it to work out of [this context](https://github.com/SimpleITK/SimpleITK/blob/4aabd77bddf508c1d55519fbf6002180a08f9208/Wrapping/Python/Python.i#L764-L794). But I think you can use it to get the job done.