Pythonic Syntactic Sugar

The Image Basics Notebook was straight forward and closely follows ITK's C++ interface.

Sugar is great it gives your energy to get things done faster! SimpleITK has applied a generous about of syntactic sugar to help get things done faster too.

%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib as mpl

mpl.rc("image", aspect="equal")
import SimpleITK as sitk

# Download data to work on
%run update_path_to_download_script
from downloaddata import fetch_data as fdata

Let us begin by developing a convenient method for displaying images in our notebooks.

img = sitk.GaussianSource(size=[64] * 2)
plt.imshow(sitk.GetArrayViewFromImage(img))

<matplotlib.image.AxesImage at 0x1139279d0>

img = sitk.GaborSource(size=[64] * 2, frequency=0.03)
plt.imshow(sitk.GetArrayViewFromImage(img))

<matplotlib.image.AxesImage at 0x1139dd150>

def myshow(img):
    nda = sitk.GetArrayViewFromImage(img)
    plt.imshow(nda)


myshow(img)

Multi-dimension slice indexing

If you are familiar with numpy, sliced index then this should be cake for the SimpleITK image. The Python standard slice interface for 1-D object:

Operation	Result
d[i]	i-th item of d, starting index 0
d[i:j]	slice of d from i to j
d[i:j:k]	slice of d from i to j with step k

With this convenient syntax many basic tasks can be easily done.

img[24, 24]

0.048901304602622986

Cropping

myshow(img[16:48, :])

myshow(img[:, 16:-16])

myshow(img[:32, :32])

Flipping

img_corner = img[:32, :32]
myshow(img_corner)

myshow(img_corner[::-1, :])

myshow(
    sitk.Tile(
        img_corner,
        img_corner[::-1, ::],
        img_corner[::, ::-1],
        img_corner[::-1, ::-1],
        [2, 2],
    )
)

Slice Extraction

A 2D image can be extracted from a 3D one.

img = sitk.GaborSource(size=[64] * 3, frequency=0.05)

# Why does this produce an error?
myshow(img)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[12], line 4
      1 img = sitk.GaborSource(size=[64] * 3, frequency=0.05)
      3 # Why does this produce an error?
----> 4 myshow(img)

Cell In[4], line 3, in myshow(img)
      1 def myshow(img):
      2     nda = sitk.GetArrayViewFromImage(img)
----> 3     plt.imshow(nda)

File ~/toolkits/anaconda3/envs/sitkMaster3.11/lib/python3.11/site-packages/matplotlib/pyplot.py:3562, in imshow(X, cmap, norm, aspect, interpolation, alpha, vmin, vmax, origin, extent, interpolation_stage, filternorm, filterrad, resample, url, data, **kwargs)
   3541 @_copy_docstring_and_deprecators(Axes.imshow)
   3542 def imshow(
   3543     X: ArrayLike | PIL.Image.Image,
   (...)
   3560     **kwargs,
   3561 ) -> AxesImage:
-> 3562     __ret = gca().imshow(
   3563         X,
   3564         cmap=cmap,
   3565         norm=norm,
   3566         aspect=aspect,
   3567         interpolation=interpolation,
   3568         alpha=alpha,
   3569         vmin=vmin,
   3570         vmax=vmax,
   3571         origin=origin,
   3572         extent=extent,
   3573         interpolation_stage=interpolation_stage,
   3574         filternorm=filternorm,
   3575         filterrad=filterrad,
   3576         resample=resample,
   3577         url=url,
   3578         **({"data": data} if data is not None else {}),
   3579         **kwargs,
   3580     )
   3581     sci(__ret)
   3582     return __ret

File ~/toolkits/anaconda3/envs/sitkMaster3.11/lib/python3.11/site-packages/matplotlib/__init__.py:1473, in _preprocess_data.<locals>.inner(ax, data, *args, **kwargs)
   1470 @functools.wraps(func)
   1471 def inner(ax, *args, data=None, **kwargs):
   1472     if data is None:
-> 1473         return func(
   1474             ax,
   1475             *map(sanitize_sequence, args),
   1476             **{k: sanitize_sequence(v) for k, v in kwargs.items()})
   1478     bound = new_sig.bind(ax, *args, **kwargs)
   1479     auto_label = (bound.arguments.get(label_namer)
   1480                   or bound.kwargs.get(label_namer))

File ~/toolkits/anaconda3/envs/sitkMaster3.11/lib/python3.11/site-packages/matplotlib/axes/_axes.py:5895, in Axes.imshow(self, X, cmap, norm, aspect, interpolation, alpha, vmin, vmax, origin, extent, interpolation_stage, filternorm, filterrad, resample, url, **kwargs)
   5892 if aspect is not None:
   5893     self.set_aspect(aspect)
-> 5895 im.set_data(X)
   5896 im.set_alpha(alpha)
   5897 if im.get_clip_path() is None:
   5898     # image does not already have clipping set, clip to Axes patch

File ~/toolkits/anaconda3/envs/sitkMaster3.11/lib/python3.11/site-packages/matplotlib/image.py:729, in _ImageBase.set_data(self, A)
    727 if isinstance(A, PIL.Image.Image):
    728     A = pil_to_array(A)  # Needed e.g. to apply png palette.
--> 729 self._A = self._normalize_image_array(A)
    730 self._imcache = None
    731 self.stale = True

File ~/toolkits/anaconda3/envs/sitkMaster3.11/lib/python3.11/site-packages/matplotlib/image.py:697, in _ImageBase._normalize_image_array(A)
    695     A = A.squeeze(-1)  # If just (M, N, 1), assume scalar and apply colormap.
    696 if not (A.ndim == 2 or A.ndim == 3 and A.shape[-1] in [3, 4]):
--> 697     raise TypeError(f"Invalid shape {A.shape} for image data")
    698 if A.ndim == 3:
    699     # If the input data has values outside the valid range (after
    700     # normalisation), we issue a warning and then clip X to the bounds
    701     # - otherwise casting wraps extreme values, hiding outliers and
    702     # making reliable interpretation impossible.
    703     high = 255 if np.issubdtype(A.dtype, np.integer) else 1

TypeError: Invalid shape (64, 64, 64) for image data

myshow(img[:, :, 32])

myshow(img[16, :, :])

Subsampling

myshow(img[:, ::3, 32])

Mathematical Operators

Most python mathematical operators are overloaded to call the SimpleITK filter which does that same operation on a per-pixel basis. They can operate on a two images or an image and a scalar.

If two images are used then both must have the same pixel type. The output image type is usually the same.

As these operators basically call ITK filter, which just use raw C++ operators, care must be taken to prevent overflow, and divide by zero etc.

Operators

+

-

*

/

//

**

img = sitk.ReadImage(fdata("cthead1.png"))
img = sitk.Cast(img, sitk.sitkFloat32)
myshow(img)
img[150, 150]

Fetching cthead1.png

libpng warning: sCAL: invalid unit
libpng warning: sCAL: invalid unit
libpng warning: sCAL: invalid unit

140.0

timg = img**2
myshow(timg)
timg[150, 150]

19600.0

Division Operators

All three Python division operators are implemented __floordiv__, __truediv__, and __div__.

The true division's output is a double pixel type.

See PEP 238 to see why Python changed the division operator in Python 3.

Bitwise Logic Operators

Operators

&

|

^

~

img = sitk.ReadImage(fdata("cthead1.png"))
myshow(img)

Fetching cthead1.png

libpng warning: sCAL: invalid unit
libpng warning: sCAL: invalid unit
libpng warning: sCAL: invalid unit

Comparative Operators

Operators

>

>=

<

<=

==

These comparative operators follow the same convention as the reset of SimpleITK for binary images. They have the pixel type of sitkUInt8 with values of 0 and 1.

img = sitk.ReadImage(fdata("cthead1.png"))
myshow(img)

Fetching cthead1.png

libpng warning: sCAL: invalid unit
libpng warning: sCAL: invalid unit
libpng warning: sCAL: invalid unit

Amazingly make common trivial tasks really trivial

myshow(img > 90)

myshow(img > 150)

myshow((img > 90) + (img > 150))

Masking

Set the values in an image according to a given mask.

# Create a grid and use as mask
grid_image = sitk.GridSource(
    outputPixelType=sitk.sitkUInt16,
    size=img.GetSize(),
    sigma=(0.1, 0.1),
    gridSpacing=(20.0, 20.0),
)
# zero out the values in the original image that correspond to
# the grid lines in the grid_image
img[grid_image == 0] = 0
myshow(img)