We will be using a library called SimpleITK to handle medical image files. It is a simplified interface around the Insight Segmentation and Registration Toolkit (ITK), one of the most popular C++ image processing libraries for medical imaging.
%% Cell type:code id: tags:
``` python
importSimpleITKassitk
```
%% Cell type:markdown id: tags:
We will use it primarily for reading and saving medical volumes encoded in NIfTI format, with the functions `sitk.ReadImage(<file_path>)` and `sitk.WriteImage(<image object>, <file_path>)`, respectively.
**Task:** Try loading the image `"ct-brain.nii.gz"` in our data directory `data_dir`:
%% Cell type:code id: tags:
``` python
img=sitk.ReadImage(data_dir+"ct-brain.nii.gz")
```
%% Cell type:markdown id: tags:
### Explore image information
Medical image files typically store raw volumetric data as a three-dimensional array. In addition, formats such as NIfTI and DICOM include a header containing a wealth of meta-information, such as:
- size: number of voxels in each dimension
- resolution/spacing: physical size of a voxel (e.g. 1mm x 1mm x 1mm)
- data type: e.g. `int32`, `uint8`, `float64`, vectors (+ number of components)
- scanner's origin and direction of coordinate axes
- voxel coordinate transformation matrices
- ... and [much more](https://brainder.org/2012/09/23/the-nifti-file-format/).
**Task:** Print the SimpleITK image object to see a summary of the loaded meta-information:
%% Cell type:code id: tags:
``` python
# Print image object
```
%% Cell type:code id: tags:
``` python
print(img)
```
%% Cell type:markdown id: tags:
SimpleITK also allows us to access each field directly.
**Task:** Let us have a look at the size and spacing of this image, with the methods `<image>.GetSize()` and `<image>.GetSpacing()`:
%% Cell type:code id: tags:
``` python
# Print image size and spacing
```
%% Cell type:code id: tags:
``` python
print(img.GetSize())
print(img.GetSpacing())
```
%% Cell type:markdown id: tags:
Although SimpleITK offers great functionality for manipulating the raw image data, we will often need to convert it to a NumPy array so it plays well with other Python packages, such as Matplotlib, which we will use for visualisation. This is done with the function `sitk.GetArrayFromImage(<image>)` (and vice-versa with `sitk.GetImageFromArray(<array>)`).
**Task:** Convert the SimpleITK image to a NumPy array
%% Cell type:code id: tags:
``` python
img_array=sitk.GetArrayFromImage(img)# Convert the SimpleITK image to a NumPy array
```
%% Cell type:markdown id: tags:
If you are curious about what else SimpleITK Image class offers, [this](http://simpleitk.github.io/SimpleITK-Notebooks/01_Image_Basics.html) is a good place to look. Additionally, you can run Python's `help` command on your image object to see all the available methods.
%% Cell type:markdown id: tags:
### Visualisation
The NumPy array you obtained from SimpleITK is three-dimensional, meaning we cannot visualise it directly. Fortunately, NumPy allows you to access entire 2D slices at a time, with `img_array[<z index>,:,:]`, `img_array[:,<y index>,:]` or `img_array[:,:,<x index>]`.
**Task:** Try printing a slice of your choice:
%% Cell type:code id: tags:
``` python
# Print an image slice
```
%% Cell type:code id: tags:
``` python
print(img_array[100])
```
%% Cell type:markdown id: tags:
As expected, we get a large two-dimensional array filled with numbers we cannot directly interpret.
To graphically display volume slices, we will use a library called Matplotlib. It is the most widely-used Python plotting library, and offers, beside more advanced functionality, a simple command-based interface, similar to Matlab's.
%% Cell type:code id: tags:
``` python
importmatplotlib.pyplotasplt
%matplotlibinline
```
%% Cell type:markdown id: tags:
To display images, Matplotlib offers the function `plt.imshow(<array>)`. It has [many options](https://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.imshow), including which colour map to use to visualise the given data (`cmap`). We will typically want to visualise medical images in greyscale (`cmap='gray'`), but feel free to play with [any of the ones available](http://matplotlib.org/examples/color/colormaps_reference.html).
By default, Matplotlib will place the origin of the y-axis at the top, increasing downwards. We can reverse it with the option `origin='lower'`.
**Task:** Investigate how to visualise axial (xy), coronal (xz) and sagittal (yz) slices according to the radiological convention.
%% Cell type:code id: tags:
``` python
# Display image slices
```
%% Cell type:code id: tags:
``` python
plt.imshow(img_array[100,:,:],cmap='gray');
```
%% Cell type:markdown id: tags:
If you tried coronal (xz) or sagittal (yz) slices, you will have noticed that the aspect ratio was wrong: the image seems 'squashed' in one dimension. This is because the resolution of this scan is *anisotropic*, i.e. the voxels do not have the same length in all directions.
`plt.imshow` has an option that lets you rescale the axes: `extent=(<left>, <right>, <bottom>, <top>)`.
**Task:** Using the `GetSize` and `GetSpacing` methods discussed earlier (note that SimpleITK's indexing convention is XYZ, whereas NumPy's is ZYX), try to set the display axes to the correct scale, in millimetres.
%% Cell type:code id: tags:
``` python
# Obtain the pysical dimensions of the image, in millimetres
Since the image data is numerical, it often makes sense to look at some of its statistics. Many basic statistics are readily available in NumPy, e.g. `np.min`, `np.max`, `np.mean`, `np.std` (standard deviation), `np.percentile` etc.
**Task:** Have a look at the minimum, maximum and mean values of your image array.
To get a better idea of the distribution of intensities, it is helpful to analyse the image *histogram*. It is a simple bar plot expressing the frequencies of pixel intensity values found in the data array, in regularly-spaced intervals ('bins').
**Task:** Use Matplotlib's `plt.hist` function to display the distribution of intensities in your image:
You will notice two main peaks in the histogram with the left corresponding to values between -1000 and -500. This is the "background" area of the image corresponding to non-tissue or air. This area often takes a significant part of the image domain.
**Task:** Try plotting the histogram again by excluding pixels from the background region, and once more plotting only values close to the second peak (in the range [-500,500]). Hint: use logical operators on the NumPy array.
In general, medical images have pixel intensity ranges beyond the typical 256 grayscales (8-bit) that we know from conventional digital images. The loaded CT scan, for example, has a range of about [-1000,2000].
By default, `plt.imshow` will display the entire intensity range, mapping the minimum and maximum values in the array to the extremes of the colour scale. However, we can manually specify the display range, setting `clim=(<low>, <high>)` or independently `vmin=<low>` and/or `vmax=<high>`.
**Task:** Using the concept of window/level, think about how to calculate parameters `clim=(<low>, <high>)` given a particular set of values for window and level. Also think about formulae for calculating the window and level that capture the full range of image intensities.
%% Cell type:code id: tags:
``` python
# Calculate parameters low and high from window and level
defwl_to_lh(window,level):
low=level-window/2
high=level+window/2
returnlow,high
print(wl_to_lh(160,70))
print(wl_to_lh(2000,300))
print(wl_to_lh(350,50))
```
%% Cell type:markdown id: tags:
**Task:** Explore displaying the loaded brain CT scan `"ct-brain.nii.gz"` with different window/level settings. For example, using a window = 120 and level = 40 should give a good setting for displaying brain tissue. Other useful settings for CT images are a window = 2000 and level = 300 for good bone contrast, or window = 350 and level = 50 for abdominal organs.
%% Cell type:code id: tags:
``` python
# Display image slices with different window/level settings
**Task:** Try also some of the other images provided in the data folder, such as `"ct-abdomen.nii.gz"` for an abdominal CT scan showing organs such as the liver or kidneys, or `"mri-brain.nii.gz"` for an example of a brain MRI scan. Try to find a good window/level setting for the MRI scan that shows good contrast for gray and white matter tissue (Hint: z-slice 130 shows a good cross-section through the ventricles).
%% Cell type:code id: tags:
``` python
importos
# Display other images with different window/level settings
**Task:** Now try writing a function to visualise an arbitrary medical volume, based on what you have done so far. It should take as input a SimpleITK image object (`img`), the slice indices (`x,y,z`) and window/level parameters (`window,level`), and display the specified axial (`z`), coronal (`y`) and sagittal (`x`) slices in grayscale and with correct axis dimensions and window/level contrast.
Note: If (`x,y,z`) are not specified (`=None`), the function should display the centre slices. If (`window,level`) are not specified, the function should calculate the window/level setting to cover the full intensity range.
The code below should give you an interactive way of displaying 3D medical images based on your `display_image` function.
%% Cell type:code id: tags:
``` python
fromipywidgetsimportinteract,fixed
fromIPython.displayimportdisplay
definteractive_view(img):
size=img.GetSize()
img_array=sitk.GetArrayFromImage(img)
interact(display_image,img=fixed(img),
x=(0,size[0]-1),
y=(0,size[1]-1),
z=(0,size[2]-1),
window=(0,np.max(img_array)-np.min(img_array)),
level=(np.min(img_array),np.max(img_array)));
interactive_view(img)
```
%% Cell type:markdown id: tags:
### Accessing voxel values
In Python, we can access a single element of a sequence with an integer index (`[index]`), or we can retrieve a contiguous sub-sequence with so-called *slice notation*: `[start:stop]`. Either or both arguments can be left blank: `[start:]` (`stop` defaults to the end), `[:stop]` (`start` defaults to `0`) or `[:]` (entire sequence).
A multi-dimensional array (such as our CT volume) can be indexed simultaneously in all dimensions, allowing us to access not only 2D slices, but also lines, rectangles, cuboids or individual voxels. A complete reference for NumPy array indexing can be found [here](https://docs.scipy.org/doc/numpy-1.11.0/reference/arrays.indexing.html).
**Task:** Explore the indexing options for extracting 1D, 2D and 3D sub-arrays, and check the resulting shapes, e.g. `img_array[50, :, 100:120].shape`, comparing with the shape of the original array:
%% Cell type:code id: tags:
``` python
print(img_array.shape)
```
%% Cell type:code id: tags:
``` python
# Print the shape of indexed sub-arrays
print(img_array[50,:,100:120].shape)
```
%% Cell type:markdown id: tags:
Indexing is not just useful for reading data in an array, we can also set values in entire regions without having to explicitly write iterative loops.
**Task:** Select any 2D slice, try setting a rectangular region to an arbitrary value and visualise the result with `plt.imshow`:
%% Cell type:code id: tags:
``` python
# Extract a 2D slice
# Set a rectangular region to a constant (e.g. 0)
# Visualise the result with plt.imshow
img=sitk.ReadImage(data_dir+'mri-brain.nii.gz')
slice=sitk.GetArrayFromImage(img)[120,:,:];
slice[100:150,50:150]=400
plt.imshow(slice,cmap='gray')
```
%% Cell type:markdown id: tags:
### Image arithmetic
Beside simply computing statistics, we can also perform arithmetic operations with image intensities, such as multiplication by scalars, contrast adjustment or any pointwise function (e.g. [gamma correction](https://en.wikipedia.org/wiki/Gamma_correction) also known as Power-Law Transform.
**Task:** Load the image `"mri-brain-t1-contrast.nii.gz"` and visualise slices (`x=105, y=170, z=95`) with window/level (`window=800, level=400`) using your `display_image` function.
**Task:** Implement a function for performing gamma correction, and apply it to the `"mri-brain-t1-contrast.nii.gz"` image (e.g., using (`c=10, gamma=3`)).
Note, for gamma correction you should first normalise the image intensity range to [0,1], apply gamma correction, and transform back to the original range. Some operations require the image data type to be floating point (`float`). NumPy arrays can be easily converted using `astype(<dtype>)`.
For the image above, an intravenous contrast agent has been used to highlight blood and tumourous tissue. In this case, we also have a corresponding brain scan taken before the contrast agent has been administered.
**Task:** Load the image `"mri-brain-t1.nii.gz"` and compare the two MRI scans with and without contrast using your `display_image` function. Hint: use the same window/level for both scans so you can better see the differences.
Note that these two images of the same patient are registered (i.e. spatially aligned).
**Task:** Now try displaying `(img1-img2)` (SimpleITK supports arithmetics with image objects) to see what was highlighted by the contrast agent. Hint: you might need to adjust the window/level for better visibility.
%% Cell type:code id: tags:
``` python
# Display the difference image
```
%% Cell type:code id: tags:
``` python
display_image(img1-img2,105,170,95,200,100)
```
%% Cell type:markdown id: tags:
**Task:** Try the same with gamma corrected version of the contrast image. Note, you might need to convert the data type of the non-contrast image to be able to subtract it from the gamma corrected image.
%% Cell type:code id: tags:
``` python
# Display the difference image of the gamma corrected images
When performing statistical analysis or running machine learning methods on images, it is often useful to first *standardise* the intensities, i.e. make them have zero mean and unit variance. This is achieved by subtracting the mean and dividing by the standard deviation of the intensities. In mathematical terms, the standardised image $\tilde{I}$ is computed as
We can also automatically enhance the contrast with [*histogram equalisation*](https://en.wikipedia.org/wiki/Histogram_equalization). It is a transformation applied to the image intensities which makes their distribution roughly uniform, i.e. all values occur approximately equally often in the image. Although it is not ideal for every application, the method is fast, simple to implement and useful for visualisation.
Below we define a function, `hist_eq`, which equalises the histogram of a given array, taking as input also the desired number of histogram bins (defaults to the maximum). The resulting array will have its values distributed almost uniformly between the original minimum and maximum values. The helper function `hist_eq_img` does exactly the same but can be called directly on SimpleITK images.
**Task:** Try to understand the individual steps of the implementation of `hist_eq`.
