Working with Images

In the previous tutorial you already saw that images are loaded as imfusion.SharedImageSet. As the name suggests, an imfusion.SharedImageSet is basically a list of images. You might think of it as a movie clip where each image in the list is a different frame.

Each image in an imfusion.SharedImageSet is an imfusion.SharedImage. The imfusion.SharedImage contains the pixel data of the image both in CPU memory but also in GPU memory (the same image is therefore shared between different locations). The framework takes care that the different locations stay synchronized. Additionally, a imfusion.SharedImage stores the position of the image in world coordinates. Also see Coordinate Systems for additional details.

A imfusion.SharedImageSet behaves similar to a list: you can iterate over it or index individual elements:

>>> (imageset,) = imfusion.load('us_sweep.dcm')
>>> imageset
imfusion.SharedImageSet(size: 20, [
        imfusion.SharedImage(UBYTE width: 164 height: 552 spacing: 0.228659x0.0724638x1 mm),
        imfusion.SharedImage(UBYTE width: 164 height: 552 spacing: 0.228659x0.0724638x1 mm),
        ...
        imfusion.SharedImage(UBYTE width: 164 height: 552 spacing: 0.228659x0.0724638x1 mm)
])
>>> image = imageset[0]
>>> image
imfusion.SharedImage(UBYTE width: 164 height: 552 spacing: 0.228659x0.0724638x1 mm)

Sometimes you might also come upon a imfusion.MemImage. This is a low-level container for the pixel data stored in CPU memory. Each imfusion.SharedImage contains a imfusion.MemImage.

When working with algorithms and other framework elements, you will most of the time work with imfusion.SharedImageSet since this is the most generic image container that can hold all kinds of medical image data: from single 3D CT volumes to a list of 2D ultrasound images.

Accessing image data with numpy

Note

Since numpy arrays can lead to direct image memory access, some of the conversions provided below require knowledge about how we internally store pixel values. In particular, shift and scale are an efficient mechanism that allow modification of the range of pixel values without actually updating all the pixels that are contained in the images. Please refer to Storage and Original values for further details about shift and scale.

The best way to work with the pixel data of an image is by converting it to a numpy.array. imfusion.MemImage, imfusion.SharedImage and imfusion.SharedImageSet allow for three slightly different ways to convert into a numpy array:

image_or_imageset.numpy(): creates a new copy of the data but shift and scale are automatically applied to the returned numpy array. The returned numpy array will have a dtype such that all the image values can be represented without information loss.
numpy.array(image_or_imageset): creates a new copy of the data but shift and scale are not applied to the returned numpy array. This means that, if needed, shift and scale need to be applied manually. The returned numpy array will have a dtype assigned by numpy’s own rules if the ‘dtype’ argument is not provided, or by the user if the ‘dtype’ argument is provided.
numpy.array(image_or_imageset, copy=False): does not create a new copy of the data but shift and scale are not applied to the returned numpy array. This means that any changes to the numpy array will be reflected to the memory pointed to by image_or_imageset and that, if needed, shift and scale need to be applied manually. The returned numpy array will have a dtype assigned by numpy’s own rules if the ‘dtype’ argument is not provided, or by the user if the ‘dtype’ argument is provided.

As an example, let’s start by looking at a imfusion.SharedImage first:

>>> import numpy as np
>>> (imageset,) = imfusion.load('ct_image.png')
>>> imageset
imfusion.SharedImageSet(size: 1, [imfusion.SharedImage(USHORT width: 512 height: 512 spacing: 0.661813x0.661813x1 mm)])
>>> image = imageset[0]
>>> arr = np.array(image)
>>> arr.shape
(512, 512, 1)

The number of dimensions can vary between 2 and 4. The shape should be read from right to left: the last dimension is always the number of channels. It is limited to values between 1 and 4. This is perceded by the width, height and slices. Dimensions of 1 on the left side are omitted. In this example the slices are omitted since the image is 2D.

Unlike the Python Imaging Library (PIL), the color dimension will even be present on grayscale images (like in the example above). This way you can always differentiate between a 2D color image (ndim of 3) and a grayscale 3D volume (ndim of 4). To interact with PIL or matplotlib you might have to remove the color dimension if it is 1:

>>> import matplotlib.pyplot as plt
>>> plt.imshow(arr[:,:,0], cmap='gray')
<matplotlib.image.AxesImage object at 0x...>

Likewise you have to add the color dimension if you want to load images with PIL into a imfusion.SharedImage:

>>> from PIL import Image
>>> arr = np.array(Image.open('ct_image.png'))
>>> arr.shape
(512, 512)
>>> imfusion.SharedImage(arr)
imfusion.SharedImage(USHORT width: 512 channels: 512)
>>> arr = np.expand_dims(arr, axis=-1)  # append color dimension
>>> arr.shape
(512, 512, 1)
>>> imfusion.SharedImage(arr)
imfusion.SharedImage(USHORT width: 512 height: 512)

Alternatively you can also use the optional greyscale argument to tell the imfusion.SharedImage to use an implicit color dimension of 1:

>>> arr = np.array(Image.open('ct_image.png'))
>>> imfusion.SharedImage(arr, greyscale=True)
imfusion.SharedImage(USHORT width: 512 height: 512)

Modifying image data

Through numpy we can now access the pixel data, but how do we modify it? You will notice that when you change values in the numpy array, the corresponding imfusion.SharedImage won’t change. This is because numpy creates a copy by default. Therefore we have to copy the changed data back into the imfusion.SharedImage using imfusion.SharedImage.assign_array() (see Direct memory access for how to avoid copies). Let’s assume we want to apply a simple thresholding:

>>> arr = np.array(image)
>>> threshold = 2000
>>> arr[arr < threshold] = 0
>>> arr[arr >= threshold] = 1
>>> image.assign_array(arr)

Note that you cannot resize or reshape the image through this way. Neither can you change the type of the image.

SharedImageSet

Since a imfusion.SharedImageSet is basically a list of imfusion.SharedImage, it has an additional dimension. Let’s load an ultrasound sweep, i.e. a time series of 2D ultrasound frames:

>>> (imageset,) = imfusion.load('us_sweep.dcm')
>>> imageset
imfusion.SharedImageSet(size: 20, [
        imfusion.SharedImage(UBYTE width: 164 height: 552 spacing: 0.228659x0.0724638x1 mm),
        imfusion.SharedImage(UBYTE width: 164 height: 552 spacing: 0.228659x0.0724638x1 mm),
        ...
        imfusion.SharedImage(UBYTE width: 164 height: 552 spacing: 0.228659x0.0724638x1 mm)
])
>>> arr = np.array(imageset)
>>> arr.shape
(20, 552, 164, 1)

The first dimension of an array which was created from a imfusion.SharedImageSet will always be the number of frames. The number of slices and number of frames will be omitted if they are 1. For example a set of 3D volumes will have 5 dimensions, a set of 2D images 4 and a set of 1D images 3. The following image shows which dimension in a numpy shape corresponds to which dimension in a imfusion.SharedImageSet:

Aside from the additional dimension, a imfusion.SharedImageSet can be used similar to a imfusion.SharedImage:

>>> imageset = imfusion.SharedImageSet(np.zeros([10, 600, 800, 1], dtype='uint8'))
>>> imageset
imfusion.SharedImageSet(size: 10, [
        imfusion.SharedImage(UBYTE width: 800 height: 600),
        imfusion.SharedImage(UBYTE width: 800 height: 600),
        ...
        imfusion.SharedImage(UBYTE width: 800 height: 600)
])
>>> imageset.assign_array(np.ones([10, 600, 800, 1], dtype='uint8'))

Direct memory access

numpy.array() always creates a copy of the input array as long as the optional copy argument is not set. While it is possible and generally supported to avoid the copy, you have to be very careful with variables referring to such arrays. Since numpy and the C++ backend are accessing the same memory block, your program will crash if the memory gets deleted while numpy is still accessing it. Direct memory access is only available on imfusion.SharedImage and not on imfusion.SharedImageSet.

Without the implicit copy, our thresholding example from before now looks like this:

>>> (imageset,) = imfusion.load('ct_image.png')
>>> image = imageset[0]
>>> arr = np.array(image, copy=False)
>>> threshold = 2000
>>> arr[arr < threshold] = 0
>>> arr[arr >= threshold] = 1
>>> image.set_dirty_mem()

Note the last line: this tells the imfusion.SharedImage that the data has changed and other shared memory locations (like the one on the GPU) must be updated. A missing call to imfusion.SharedImage.setDirtyMem() is a common source of bugs, e.g. GPU algorithms working on outdated data or overwriting changes.

A general recommendation is to write your algorithm by copying the data and only avoid the copy if you notice performance issues. This tends to avoid bugs and leads to less headaches when developing.

Storage and Original values

After converting an image to a numpy array, you might have noticed that the values are not in the range you expect them to be. Especially if you loaded some data with negative values, all the values in the array might be positive. The ImFusionLib distinguishes between Storage- and Original values, which are indicated by the imfusion.SharedImage.shift and imfusion.SharedImage.scale properties.

Original: Pixel values are the same as in their original source (e.g. when loaded from a file). Same as the storage pixel values if the image’s scale is 1 and the shift is 0
Storage: Pixel values as they are stored in a imfusion.MemImage. The user may decide to apply such a rescaling in order to better use the available limits of the underlying type.

The following conversion rules apply:

\[ \begin{align}\begin{aligned}original = (storage / scale) - shift\\storage = (original + shift) * scale\end{aligned}\end{align} \]

Since this is a common operation, there are imfusion.SharedImage.apply_shift_and_scale() and imfusion.SharedImageSet.apply_shift_and_scale() that convert a numpy array from storage to original values:

>>> (imageset,) = imfusion.load('ct_image.png')
>>> imageset[0].shift = 1024.0
>>> arr = np.array(imageset)
>>> arr.dtype
dtype('uint16')
>>> np.min(arr)
0
>>> arr = imageset.apply_shift_and_scale(arr)
>>> arr.dtype
dtype('float64')
>>> np.min(arr)
-1024.0

Type conversions

imfusion.MemImage, imfusion.SharedImage and imfusion.SharedImageSet allow to convert the image or images into a new format. It works similarly as with numpy:

>>> from imfusion import PixelType
>>> (imageset,) = imfusion.load('ct_image.png')
>>> image = imageset[0]
>>> image_as_float = imageset[0].astype(PixelType.FLOAT)
>>> image.descriptor.pixel_type, image_as_float.descriptor.pixel_type
(<PixelType.USHORT: 5123>, <PixelType.FLOAT: 5126>)

Besides, the astype method supports as input:

a PixelType (e.g. imfusion.PixelType.UINT);
most of the numpy’s dtypes (e.g. np.uint);
python’s float or int types.

Selections

class imfusion.Selection

Utility class for describing a selection of elements out of a set. Conceptually, a Selection pairs a list of bools describing selected items with the index of a “focus” item and provides syntactic sugar on top. For instance, the set of selected items could define which ones to show in general while the focus item is additionally highlighted. The class is fully separate from the item set of which it describes the selection. This means for instance that it cannot know the actual number of items in the set and the user/parent class must manually make sure that they match. Also, a Selection only manages indices and offers no way of accessing the underlying elements. In order to iterate over all selected indices, you can do for instance the following:

for index in range(selection.start, selection.stop):
        if selection[index]:
                ...

The same effect can also be achieved in a much more terse fashion:

for selected_index in selection.selected_indices():
...

For convenience, the selection can also be converted to a slice object (if the selection has a regular spacing, see below):

selected_subset = container[selection.slice()]

Sometimes it can be more convenient to “thin out” a selection by only selecting every N-th element. To this end, the Selection constructor takes the arguments start, stop and step. Setting step to N will only select every N-th element, mimicking the signature of range, slice, etc.