Working with Images
===================

In the previous tutorial, you already saw that images are loaded as :class:`imfusion.SharedImageSet`.
As the name suggests, an :class:`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 :class:`imfusion.SharedImageSet` is an :class:`imfusion.SharedImage`. The :class:`imfusion.SharedImage` stores the image's pixel data in both the CPU and the GPU memory, enabling *shared* access of the same image across these locations. The framework takes care that they always stay synchronized. In addition, an :class:`imfusion.SharedImage` stores the position of the image in world coordinates. For more details, see :doc:`coordinatesystems`.

An :class:`imfusion.SharedImageSet` behaves similarly to a Python ``list``: you can iterate over it and index individual elements:

>>> imageset, *_ = imfusion.io.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)

Occasionally, you might come across an :class:`imfusion.MemImage`, which is a low-level container for the pixel data stored in CPU memory. Therefore, each :class:`imfusion.SharedImage` contains an :class:`imfusion.MemImage`.

When working with algorithms and other framework components, you will typically deal with an :class:`imfusion.SharedImageSet`, as it is the most generic image container. It can hold various types of medical image data, ranging from single 3D CT volumes to lists 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 efficient mechanisms 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``.
:class:`imfusion.MemImage`, :class:`imfusion.SharedImage`, and :class:`imfusion.SharedImageSet` offer three slightly different ways to convert them into a numpy array:

#. `image_or_imageset.numpy()`: Creates a new copy of the data, whereby shift and scale are automatically applied to the returned numpy array. The `numpy` array's `dtype` is set such that all image values can be represented without information loss.
#. `numpy.array(image_or_imageset)`: Creates a new copy of the data, whereby shift and scale are not automatically applied to the returned numpy array. Thus, if needed, shift and scale must be applied manually. The numpy array's `dtype` can be specified by passing it as an additional argument. Otherwise, it will be assigned by numpy's internal rules.
#. `numpy.array(image_or_imageset, copy=False)`: Does not create a new copy of the data, and shift and scale are not applied to the returned numpy array. This means that any changes to the numpy array will be reflected in the memory pointed to by `image_or_imageset`. If needed, shift and scale must be applied manually. The numpy array's `dtype` can be specified by passing it as an additional argument. Otherwise, it will be assigned by numpy's internal rules.

As a first example, let's take a closer look at an :class:`imfusion.SharedImage`:

>>> import numpy as np
>>> (imageset,) = imfusion.io.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, limited to values between 1 and 4. This is preceded by the width, height, and the number of slices. Dimensions of 1 on the left side are omitted. In this example, the number of slices is omitted as the given 2D image only has one slice.

Unlike the `Python Imaging Library (PIL) <https://en.wikipedia.org/wiki/Python_Imaging_Library>`_, the color dimension is also 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...>

.. image:: images/ct_image_matplotlib.png
   :align: center

Accordingly, you have to add the color dimension if you want to load images with PIL into an :class:`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_col = np.expand_dims(arr, axis=-1)  # append color dimension
>>> arr_col.shape
(512, 512, 1)
>>> imfusion.SharedImage(arr_col)
imfusion.SharedImage(USHORT width: 512 height: 512)

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

>>> imfusion.SharedImage(arr, greyscale=True)
imfusion.SharedImage(USHORT width: 512 height: 512)

Modifying image data
--------------------

Using ``numpy``, we can now access the pixel data, but how do we modify it? You will notice that changes made to the ``numpy`` array are not reflected in the corresponding :class:`imfusion.SharedImage`. This is because ``numpy`` creates a copy by default. Therefore, we have to copy the altered data back into the :class:`imfusion.SharedImage` using :meth:`imfusion.SharedImage.assign_array` (see :ref:`direct-memory-access` for how to avoid copies in the first place). Let's assume we want to apply a simple threshold:

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

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

SharedImageSet
--------------

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

>>> imageset, *_ = imfusion.io.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 an :class:`imfusion.SharedImageSet`, will always be the number of frames. It is important to remember that the number of slices and number of frames will be omitted if they are 1. Consequently, 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 how each dimension in a ``numpy.shape`` relates to the corresponding dimension in an :class:`imfusion.SharedImageSet`:

.. image:: images/dimensions.png
   :align: center

Aside from the additional dimension, an :class:`imfusion.SharedImageSet` can be used similarly to an :class:`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:

Direct memory access
--------------------

:func:`numpy.array` always creates a copy of the input array as long as the optional *copy* argument is not set (default: copy=True). 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 :class:`imfusion.SharedImage` and *not* on :class:`imfusion.SharedImageSet`.

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

>>> (imageset,) = imfusion.io.load(path_to / '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 :class:`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 :meth:`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 pixel values, whereby their relation is defined by the :attr:`imfusion.SharedImage.shift` and :attr:`imfusion.SharedImage.scale` properties.

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

The following conversion rules apply:

.. math::

   original = (storage / scale) - shift

   storage = (original + shift) * scale

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

    >>> imageset, *_ = imfusion.io.load(ct_image_png)
    >>> imageset[0].shift = 1024.0
    >>> arr = np.array(imageset)
    >>> arr.dtype
    dtype('uint16')
    >>> int(np.min(arr))
    0
    >>> arr = imageset.apply_shift_and_scale(arr)
    >>> arr.dtype
    dtype('float64')
    >>> float(np.min(arr))
    -1024.0

Type conversions
---------------------------
:class:`imfusion.MemImage`, :class:`imfusion.SharedImage` and :class:`imfusion.SharedImageSet` allow to convert the image or images into a new format. It works similarly to conversions in ``numpy``:

>>> from imfusion import PixelType
>>> imageset, *_ = imfusion.io.load(ct_image_png)
>>> image = imageset[0]
>>> image_as_float = image.astype(PixelType.FLOAT)
>>> image.descriptor.pixel_type, image_as_float.descriptor.pixel_type
(<PixelType.USHORT: 5123>, <PixelType.FLOAT: 5126>)

Besides the above example, 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
----------

.. autoclass:: imfusion.Selection
	:noindex: