diff --git a/README.md b/README.md index c8d94094..b19e77fd 100644 --- a/README.md +++ b/README.md @@ -30,6 +30,7 @@ This project exists thanks to all mainteners. Sort list by last name (alphabetic ## Contributors This project exists thanks to all the people who contribute. Sort list by last name (alphabetically). * [Dave Cash](https://github.com/davecash75) +* [Clea Dronne](https://github.com/clead6) * [Milan Malfait](https://github.com/milanmlft) * [Jamie Mcclelland](https://github.com/jamie-mcclelland) * [Zakaria Senousy](https://github.com/zsenousy) diff --git a/config.yaml b/config.yaml index 3c3597c6..0a03c362 100644 --- a/config.yaml +++ b/config.yaml @@ -61,8 +61,10 @@ contact: 'team@carpentries.org' # Order of episodes in your lesson episodes: #- introduction.Rmd -- itksnap.Rmd -- excthreesteps.Rmd +- practical1.Rmd +- practical2.Rmd +- practical3.Rmd +- practical4.Rmd # Information for Learners learners: diff --git a/episodes/data/README.md b/episodes/data/README.md deleted file mode 100644 index 5baaa76d..00000000 --- a/episodes/data/README.md +++ /dev/null @@ -1,59 +0,0 @@ -# Medical Data Imaging - -## CT-vs-PET-Ventilation-Imaging -CT Ventilation as a Functional Imaging Modality for Lung Cancer Radiotherapy from [TCIA](https://www.cancerimagingarchive.net/collection/ct-vs-pet-ventilation-imaging/). -We recomend to focus on exhale/inhale breath hold CT (BHCT) Dicoms which belog to: -* The BHCT scans for CT-PET-VI-02 & CT-PET-VI-03 show very little motion between inhale and exhale -* The BHCT scans for CT-PET-VI-05 have a different number of slices between the inhale and exhale - -Data paths look like these: -``` -/CT-PET-VI-02$ tree -d -. -├── 1000.000000-PET SUV Factors-26496 [1 item, with size 42.1 kB] -├── 3.000000-CT Lung 3.0 B31f-61322 [175 items, totalling 92.4 MB] -├── 4.000000-Galligas Lung-03537 [159 items, totalling 51.5 MB] -├── 5.000000-Thorax Insp 2.0 B70f-29873 [167 items, totalling 88.1 MB] -├── 7.000000-Thorax Exp 2.0 B70f-73355 [167 items, totalling 88.1 MB] -└── 8.000000-RespLow 2.0 B30s 80 Ex - 100 In-44317 [810 items, totalling 426.9 MB] - -/CT-PET-VI-03$ tree -d -. -├── 06-25-1996-RTRTRespLowBreathRateRNS Adult-43080 -│ ├── 1000.000000-PET SUV Factors-06580 -│ ├── 3.000000-CT Lung 3.0 B31f-08354 -│ └── 4.000000-Galligas Lung-15379 -└── 06-25-1996-RTRTRespLowBreathRateRNS Adult-87093 - ├── 5.000000-Thorax Insp 2.0 B70f-42000 - ├── 7.000000-Thorax Exp 2.0 B70f-45310 - └── 8.000000-RespLow 2.0 B30s 80 Ex - 100 In-54790 - -/CT-PET-VI-05$ tree -d -└── 12-10-1996-RTRTRespLowBreathRateRNS Adult-37585 - ├── 1000.000000-PET SUV Factors-63910 - ├── 3.000000-CT Lung 3.0 B31f-90638 - ├── 4.000000-Galligas Lung-26361 - ├── 5.000000-Thorax Insp 2.0 B70f-45546 - ├── 7.000000-Thorax Exp 2.0 B70f-31579 - └── 8.000000-RespLow 2.0 B30s 80 Ex - 100 In-16454 -``` - - -### Citations & Data Usage Policy -Data Citation Required: Users must abide by [the TCIA Data Usage Policy and Restrictions](https://www.cancerimagingarchive.net/data-usage-policies-and-restrictions/). -Attribution must include the following citation, including the Digital Object Identifier: -> Eslick, E. M., Kipritidis, J., Gradinscak, D., Stevens, M. J., Bailey, D. L., Harris, B., Booth, J. T., & Keall, P. J. (2022). CT Ventilation as a functional imaging modality for lung cancer radiotherapy (CT-vs-PET-Ventilation-Imaging) (Version 1) [Data set]. The Cancer Imaging Archive. https://doi.org/10.7937/3ppx-7s22 - -### Data is converted from DICOM to nifti using `dicom2nifti`: -``` -from pathlib import Path -import dicom2nifti -dicoms = Path("~/3_000000-CT_Lung_30_B31f-61322") -dicom2nifti.convert_directory(dicoms, ".", compression=True, reorient=True) -``` - - -## Demons Image Registration Exercise Data - -Please check [`data`](https://github.com/HealthBioscienceIDEAS/Medical-Image-Registration-Short-Course/tree/main/episodes/data) folder inside `episodes` for accessing `cine_MR_img_1.png`, `cine_MR_img_2.png`, and `cine_MR_img_3.png` that are used for [`demonReg.py`](https://github.com/HealthBioscienceIDEAS/Medical-Image-Registration-Short-Course/blob/main/src/mirc/utils/demonsReg.py) and [`exampleSolution3.py`](https://github.com/HealthBioscienceIDEAS/Medical-Image-Registration-Short-Course/blob/main/src/mirc/examples/exampleSolution3.py). - diff --git a/episodes/data/cine_MR_1.png b/episodes/data/cine_MR_1.png deleted file mode 100644 index 36bc5897..00000000 Binary files a/episodes/data/cine_MR_1.png and /dev/null differ diff --git a/episodes/data/cine_MR_2.png b/episodes/data/cine_MR_2.png deleted file mode 100644 index 1209486c..00000000 Binary files a/episodes/data/cine_MR_2.png and /dev/null differ diff --git a/episodes/data/cine_MR_3.png b/episodes/data/cine_MR_3.png deleted file mode 100644 index d9f5848d..00000000 Binary files a/episodes/data/cine_MR_3.png and /dev/null differ diff --git a/episodes/excthreesteps.Rmd b/episodes/excthreesteps.Rmd deleted file mode 100644 index 9cdd9a47..00000000 --- a/episodes/excthreesteps.Rmd +++ /dev/null @@ -1,477 +0,0 @@ ---- -title: 'Demons Image Registration' -teaching: 10 -exercises: 2 ---- - -:::::::::::::::::::::::::::::::::::::: questions - -- How to understand and visualise three image in one pane for demons image registration algorithm? - -:::::::::::::::::::::::::::::::::::::::::::::::: - -::::::::::::::::::::::::::::::::::::: objectives - -- Understanding Demons Image Registration code for 3 images viewing pane - -:::::::::::::::::::::::::::::::::::::::::::::::: - -# Demons Image Registration Algorithm with Multi-Pane Display - -## Activate Conda environment -Open terminal, activate `mirVE` environment and run python - -```bash -conda activate mirVE -python -``` -:::::::::::::::::::::::::::::::::::::: spoiler - -## How to avoid `ModuleNotFoundError: No module named VARIABLE`? - -You need to export relative imports, using the following notation for either linux or windows (might be similar to mac) - -```bash -#WINDOWS -set PYTHONPATH=%PYTHONPATH%;C:\path\to\your\project\ -#LINUX -export PYTHONPATH="${PYTHONPATH}:/path/to/your/project/" -``` - -:::::::::::::::::::::::::::::::::::::: - -This document explains the implementation of a 2D image registration algorithm using the Demons algorithm. The implementation includes visualisation updates within a single window containing three panes for easier comparison of source, target, and warped images. - -## 1. Importing Libraries - -First, we import the necessary libraries for image processing, transformation, and visualisation. The `utils3` module is another Python file containing utility functions used for tasks such as displaying images and handling deformed and updated fields. Both `demonsReg.py` and `utils3.py` are available in `src/mirc/utils` - -```python -import matplotlib -matplotlib.use('TkAgg') # Set the backend for matplotlib -import matplotlib.pyplot as plt # Import matplotlib for plotting - -import numpy as np # Import numpy for numerical operations -from skimage.transform import rescale, resize # Import functions for image transformation from scikit-image -from scipy.ndimage import gaussian_filter # Import gaussian_filter function from scipy's ndimage module -from utils3 import dispImage, resampImageWithDefField, calcMSD, dispDefField # Import specific functions from a custom module -``` - -## 2. Function Definition - -The `demonsReg` function is designed for registering two 2D images using the Demons algorithm, transforming the source image to match the target. It offers flexibility through various optional parameters to customise the registration process: - -```python -def demonsReg(source, target, sigma_elastic=1, sigma_fluid=1, num_lev=3, use_composition=False, - use_target_grad=False, max_it=1000, check_MSD=True, disp_freq=3, disp_spacing=2, - scale_update_for_display=10, disp_method_df='grid', disp_method_up='arrows'): - """ - Perform a registration between the 2D source image and the 2D target - image using the demons algorithm. The source image is warped (resampled) - into the space of the target image. - - Parameters: - - source: 2D numpy array, the source image to be registered. - - target: 2D numpy array, the target image for registration. - - sigma_elastic: float, standard deviation for elastic regularisation. - - sigma_fluid: float, standard deviation for fluid regularisation. - - num_lev: int, number of levels in the multi-resolution scheme. - - use_composition: bool, whether to use composition in the update step. - - use_target_grad: bool, whether to use the target image gradient. - - max_it: int, maximum number of iterations for the registration. - - check_MSD: bool, whether to check Mean Squared Difference for improvements. - - disp_freq: int, frequency of display updates during registration. - - disp_spacing: int, spacing between grid lines or arrows in display. - - scale_update_for_display: int, scale factor for displaying the update field. - - disp_method_df: str, method for displaying the deformation field ('grid' or 'arrows'). - - disp_method_up: str, method for displaying the update field ('grid' or 'arrows'). - - Returns: - - warped_image: 2D numpy array, the source image warped into the target image space. - - def_field: 3D numpy array, the deformation field used to warp the source image. - """ -``` - - -## 3. Preparing Images and figure for demons image registration algorithm - -We start by making copies of the full-resolution images and initiating a figure for plotting during registration process. - -```python - # make copies of full resolution images - source_full = source - target_full = target - - #Preparing figure for live update during registration process - fig, axs = plt.subplots(1, 3, figsize=(12, 6)) - iteration_text = fig.text(0.5, 0.92, '', ha='center', va='top', fontsize=10, color='black') -``` - - - -## 4. Multi-Resolution Scheme and Initialisation - -In this step, we initialise variables and set up the multi-resolution scheme by looping over different resolution levels. - -```python -# Loop over resolution levels -for lev in range(1, num_lev + 1): - - # Resample images if not at the final level - if lev != num_lev: - resamp_factor = np.power(2, num_lev - lev) - target = rescale(target_full, 1.0 / resamp_factor, mode='edge', order=3, anti_aliasing=True) - source = rescale(source_full, 1.0 / resamp_factor, mode='edge', order=3, anti_aliasing=True) - else: - target = target_full - source = source_full -``` - - - -## 5. Deformation Field Initialisation -In this step, the code initialises the deformation field and related variables necessary for the registration process. The deformation field (def_field) is crucial as it defines the transformation that will be iteratively adjusted to align the source image with the target image. At the first resolution level (lev == 1), the initial deformation field is set up using the grid coordinates derived from the target image dimensions. Additionally, placeholders for the displacement field components (disp_field_x and disp_field_y) are created to track incremental changes in the deformation over iterations. - -```python - # If first level, initialise deformation and displacement fields - if lev == 1: - X, Y = np.mgrid[0:target.shape[0], 0:target.shape[1]] - def_field = np.zeros((X.shape[0], X.shape[1], 2)) - def_field[:, :, 0] = X - def_field[:, :, 1] = Y - disp_field_x = np.zeros(target.shape) - disp_field_y = np.zeros(target.shape) - else: - # Otherwise, upsample displacement field from previous level - disp_field_x = 2 * resize(disp_field_x, (target.shape[0], target.shape[1]), mode='edge', order=3) - disp_field_y = 2 * resize(disp_field_y, (target.shape[0], target.shape[1]), mode='edge', order=3) - # Recalculate deformation field for this level from displacement field - X, Y = np.mgrid[0:target.shape[0], 0:target.shape[1]] - def_field = np.zeros((X.shape[0], X.shape[1], 2)) # Clear def_field from previous level - def_field[:, :, 0] = X + disp_field_x - def_field[:, :, 1] = Y + disp_field_y -``` - - - -## 6. Update Initialisation -In this step, the code initialises the update fields and computes the initial warped image for the current resolution level. These updates are essential for applying iterative transformations to align the source image with the target image. - - -```python - # Initialise updates - update_x = np.zeros(target.shape) - update_y = np.zeros(target.shape) - update_def_field = np.zeros(def_field.shape) - - # Calculate the transformed image at the start of this level - warped_image = resampImageWithDefField(source, def_field) - - # Store the current def_field and MSD value to check for improvements at the end of iteration - def_field_prev = def_field.copy() - prev_MSD = calcMSD(target, warped_image) - - # If target image gradient is being used, this can be calculated now as it will not change during the registration - if use_target_grad: - img_grad_x, img_grad_y = np.gradient(target) -``` - - -## 7. Initial Transformation and Preparation - -In this step, the code performs the initial transformation of the source image using the current deformation field, calculates and stores the initial Mean Squared Difference (MSD), and optionally computes the gradient of the target image. This sets the stage for the iterative registration process. - - -```python -# calculate the transformed image at the start of this level -warped_image = resampImageWithDefField(source, def_field) - -# store the current def_field and MSD value to check for improvements at -# end of iteration -def_field_prev = def_field.copy() -prev_MSD = calcMSD(target, warped_image) - -# if target image gradient is being used this can be calculated now as it will -# not change during the registration -if use_target_grad: - [img_grad_x, img_grad_y] = np.gradient(target) -``` - - -## 8. Live update of Registration process - -Here, we introduce `live_update` function which works on updating warped_image, deformation field and update field subplots during the registration process. This function is used inside the main iterative loop of the registration process (See next step). - - ```python - # Function to update the display - def live_update(): - # Clear the axes - for ax in axs: - ax.clear() - - # Plotting the warped image in the first subplot - plt.sca(axs[0]) - dispImage(warped_image, title='Warped Image') - x_lims = plt.xlim() - y_lims = plt.ylim() - # Plotting the deformation field in the second subplot - plt.sca(axs[1]) - dispDefField(def_field, spacing=disp_spacing, plot_type=disp_method_df) - axs[1].set_xlim(x_lims) - axs[1].set_ylim(y_lims) - axs[1].set_title('Deformation Field') - - # Plotting the updated deformation field in the third subplot - plt.sca(axs[2]) - up_field_to_display = scale_update_for_display * np.dstack((update_x, update_y)) - up_field_to_display += np.dstack((X, Y)) - dispDefField(up_field_to_display, spacing=disp_spacing, plot_type=disp_method_up) - axs[2].set_xlim(x_lims) - axs[2].set_ylim(y_lims) - axs[2].set_title('Update Field') - - # Update the iteration text - iteration_text.set_text('Level {0:d}, Iteration {1:d}: MSD = {2:.6f}'.format(lev, it, prev_MSD)) - - # Adjust layout for better spacing - plt.tight_layout() - - # Redraw the figure - fig.canvas.draw() - fig.canvas.flush_events() # Ensure the figure updates immediately - - plt.pause(0.5) - ``` - -## 9. Main Iterative Loop for Image Registration -In this step, we perform the main iterative loop where the actual image registration occurs. The loop continues until the maximum number of iterations is reached or until convergence criteria are met. The loop updates the deformation field iteratively to minimise the Mean Squared Difference (MSD) between the target and warped images. - -```python -# main iterative loop - repeat until max number of iterations reached - for it in range(max_it): - - # calculate update from demons forces - # - # if the warped image gradient is used (instead of the target image gradient) - # this needs to be calculated - if not use_target_grad: - [img_grad_x, img_grad_y] = np.gradient(warped_image) - - # calculate difference image - diff = target - warped_image - # calculate denominator of demons forces - denom = np.power(img_grad_x, 2) + np.power(img_grad_y, 2) + np.power(diff, 2) - # calculate x and y components of numerator of demons forces - numer_x = diff * img_grad_x - numer_y = diff * img_grad_y - # calculate the x and y components of the update - #denom[denom < 0.01] = np.nan - update_x = numer_x / denom - update_y = numer_y / denom - - # set nan values to 0 - update_x[np.isnan(update_x)] = 0 - update_y[np.isnan(update_y)] = 0 - - # if fluid like regularisation used smooth the update - if sigma_fluid > 0: - update_x = gaussian_filter(update_x, sigma_fluid, mode='nearest') - update_y = gaussian_filter(update_y, sigma_fluid, mode='nearest') - - # update displacement field using addition (original demons) or - # composition (diffeomorphic demons) - if use_composition: - # compose update with current transformation - this is done by - # transforming (resampling) the current transformation using the - # update. we can use the same function as used for resampling - # images, and treat each component of the current deformation - # field as an image - # the update is a displacement field, but to resample an image - # we need a deformation field, so need to calculate deformation - # field corresponding to update. - update_def_field[:, :, 0] = update_x + X - update_def_field[:, :, 1] = update_y + Y - # use this to resample the current deformation field, storing - # the result in the same variable, i.e. we overwrite/update the - # current deformation field with the composed transformation - def_field = resampImageWithDefField(def_field, update_def_field) - # calculate the displacement field from the composed deformation field - disp_field_x = def_field[:, :, 0] - X - disp_field_y = def_field[:, :, 1] - Y - # replace nans in disp field with 0s - disp_field_x[np.isnan(disp_field_x)] = 0 - disp_field_y[np.isnan(disp_field_y)] = 0 - else: - # add the update to the current displacement field - disp_field_x = disp_field_x + update_x - disp_field_y = disp_field_y + update_y - - - # if elastic like regularisation used smooth the displacement field - if sigma_elastic > 0: - disp_field_x = gaussian_filter(disp_field_x, sigma_elastic, mode='nearest') - disp_field_y = gaussian_filter(disp_field_y, sigma_elastic, mode='nearest') - - # update deformation field from disp field - def_field[:, :, 0] = disp_field_x + X - def_field[:, :, 1] = disp_field_y + Y - - # transform the image using the updated deformation field - warped_image = resampImageWithDefField(source, def_field) - - # update images if required for this iteration - if disp_freq > 0 and it % disp_freq == 0: - # Create a single figure with 3 subplots in one row and three columns - live_update() - - # calculate MSD between target and warped image - MSD = calcMSD(target, warped_image) - - # display numerical results - print('Level {0:d}, Iteration {1:d}: MSD = {2:f}\n'.format(lev, it, MSD)) - - # check for improvement in MSD if required - if check_MSD and MSD >= prev_MSD: - # restore previous results and finish level - def_field = def_field_prev - warped_image = resampImageWithDefField(source, def_field) - print('No improvement in MSD') - break - - # update previous values of def_field and MSD - def_field_prev = def_field.copy() - prev_MSD = MSD.copy() - -``` - -### 9.1. Results after Level 1 iterations: -![Level 1 Results](fig/level1.png) - -### 9.2. Results after Level 2 iterations: -![Level 2 Results](fig/level2.png) - -### 9.3. Results after Level 3 iterations: -![Level 3 Results](fig/level3.png) - - - - -## 11. Displaying Final Results - -In this step, we visualise the final results of the image registration process. This code snippet sets up a visualisation interface using Matplotlib to display and interact with images and plots related to image processing tasks. Global variables are initialised to track the current indices for both images and display modes. Three images are defined: source, target, and warped_image, along with corresponding titles stored in image_titles. Two display modes, 'Deformation Field' and 'Jacobian', are defined in the modes list. The code defines an event handler function, on_key, which responds to key presses ('left', 'right', 'up', 'down') to navigate between images and switch display modes. The update_display function clears and updates three subplots within a single figure based on the current indices and modes, ensuring the visualisations are refreshed dynamically. Finally, a Matplotlib figure with subplots is created, initial images and plots are displayed, and instructions for navigation are added at the bottom. The figure is connected to the event handler to enable interactive navigation and mode switching. - -```python -if lev == num_lev: - # Initialize global variables for current index tracking - current_image_index = [0] - current_mode_index = [0] - - # Define the images and titles - images = [source, target, warped_image] - image_titles = ['Source Image', 'Target Image', 'Warped Image'] - modes = ['Deformation Field', 'Jacobian'] - - - def on_key(event): - if event.key == 'right': - current_image_index[0] = (current_image_index[0] + 1) % len(images) - elif event.key == 'left': - current_image_index[0] = (current_image_index[0] - 1) % len(images) - elif event.key == 'up' or event.key == 'down': - current_mode_index[0] = (current_mode_index[0] + 1) % len(modes) - update_display() - - - - def update_display(): - axs_combined[0].clear() - plt.sca(axs_combined[0]) - dispImage(images[current_image_index[0]], title=image_titles[current_image_index[0]]) - - axs_combined[1].clear() - plt.sca(axs_combined[1]) - if modes[current_mode_index[0]] == 'Deformation Field': - dispDefField(def_field, spacing=disp_spacing, plot_type=disp_method_df) - axs_combined[1].set_title('Deformation Field') - - else: - [jacobian, _] = calcJacobian(def_field) - dispImage(jacobian, title='Jacobian') - plt.set_cmap('jet') - #plt.colorbar() - - axs_combined[2].clear() - plt.sca(axs_combined[2]) - diff_image = images[current_image_index[0]] - target - dispImage(diff_image, title='Difference Image') - - fig_combined.canvas.draw() - - - # Create a single figure with 3 subplots - fig_combined, axs_combined = plt.subplots(1, 3, figsize=(12, 6)) - - # Display initial images - plt.sca(axs_combined[0]) - dispImage(images[current_image_index[0]], title=image_titles[current_image_index[0]]) - - plt.sca(axs_combined[1]) - dispDefField(def_field, spacing=disp_spacing, plot_type=disp_method_df) - axs_combined[1].set_title('Deformation Field') - - plt.sca(axs_combined[2]) - diff_image = images[current_image_index[0]] - target - dispImage(diff_image, title='Difference Image') - - # Add instructions for navigating images - fig_combined.text(0.5, 0.02, 'Press <- or -> to navigate between source, target and warped images, Press Up or Down to switch between deformation field and Jacobian', ha='center', va='top', fontsize=12, color='black') - - # Connect the key event handler to the figure - fig_combined.canvas.mpl_connect('key_press_event', on_key) - - plt.tight_layout() - plt.show() -``` - - -## 12. Preparing Images for calling the demonsReg function - -Follow `exampleSolution3.py` available in `src/mirc/example` for running the function. Below, we describe how we prepare images: `cine_MR_img_1.png`, `cine_MR_img_2.png`, and `cine_MR_img_3.png`. These images are available in `episodes/data`. - - -```python -import skimage.io -cine_MR_img_1 = skimage.io.imread('path/to/your/image/../cine_MR_1.png') -cine_MR_img_2 = skimage.io.imread('path/to/your/image/../cine_MR_2.png') -cine_MR_img_3 = skimage.io.imread('path/to/your/image/../cine_MR_3.png') -``` - -Converting all images into double data for standard orientation. - -```python -import numpy as np -# *************** -# ADD CODE HERE TO CONVERT ALL THE IMAGES TO DOUBLE DATA TYPE AND TO REORIENTATE THEM -# INTO 'STANDARD ORIENTATION' -cine_MR_img_1 = np.double(cine_MR_img_1) -cine_MR_img_2 = np.double(cine_MR_img_2) -cine_MR_img_3 = np.double(cine_MR_img_3) - -cine_MR_img_1 = np.flip(cine_MR_img_1.T, 1) -cine_MR_img_2 = np.flip(cine_MR_img_2.T, 1) -cine_MR_img_3 = np.flip(cine_MR_img_3.T, 1) -``` - - -### Calling the function: - -```python -demonsReg(cine_MR_img_1, cine_MR_img_2) -``` - -### Final Results Video (Navigation between images and modes) - - diff --git a/episodes/fig/itk-snap-additional-image-1.png b/episodes/fig/itk-snap-additional-image-1.png new file mode 100644 index 00000000..6b52533a Binary files /dev/null and b/episodes/fig/itk-snap-additional-image-1.png differ diff --git a/episodes/fig/itk-snap-additional-image-2.png b/episodes/fig/itk-snap-additional-image-2.png new file mode 100644 index 00000000..8586600a Binary files /dev/null and b/episodes/fig/itk-snap-additional-image-2.png differ diff --git a/episodes/fig/itk-snap-additional-image-3.png b/episodes/fig/itk-snap-additional-image-3.png new file mode 100644 index 00000000..287cac3f Binary files /dev/null and b/episodes/fig/itk-snap-additional-image-3.png differ diff --git a/episodes/fig/itk-snap-additional-image-4.png b/episodes/fig/itk-snap-additional-image-4.png new file mode 100644 index 00000000..2b8ab69e Binary files /dev/null and b/episodes/fig/itk-snap-additional-image-4.png differ diff --git a/episodes/fig/itk-snap-additional-image-5.png b/episodes/fig/itk-snap-additional-image-5.png new file mode 100644 index 00000000..ad3e5b1f Binary files /dev/null and b/episodes/fig/itk-snap-additional-image-5.png differ diff --git a/episodes/fig/itk-snap-additional-image-6.png b/episodes/fig/itk-snap-additional-image-6.png new file mode 100644 index 00000000..ed144ecc Binary files /dev/null and b/episodes/fig/itk-snap-additional-image-6.png differ diff --git a/episodes/fig/itk-snap-additional-image.png b/episodes/fig/itk-snap-additional-image.png new file mode 100644 index 00000000..1ee08f11 Binary files /dev/null and b/episodes/fig/itk-snap-additional-image.png differ diff --git a/episodes/fig/itk-snap-additional-image.svg b/episodes/fig/itk-snap-additional-image.svg deleted file mode 100644 index f2f7f7da..00000000 --- a/episodes/fig/itk-snap-additional-image.svg +++ /dev/null @@ -1,3979 +0,0 @@ - -image/svg+xml diff --git a/episodes/fig/itk-snap-aligned-images-1.png b/episodes/fig/itk-snap-aligned-images-1.png new file mode 100644 index 00000000..6f01227a Binary files /dev/null and b/episodes/fig/itk-snap-aligned-images-1.png differ diff --git a/episodes/fig/itk-snap-aligned-images-2.png b/episodes/fig/itk-snap-aligned-images-2.png new file mode 100644 index 00000000..92881ff6 Binary files /dev/null and b/episodes/fig/itk-snap-aligned-images-2.png differ diff --git a/episodes/fig/itk-snap-cropped-ct-aligned-1.png b/episodes/fig/itk-snap-cropped-ct-aligned-1.png new file mode 100644 index 00000000..5daa0d88 Binary files /dev/null and b/episodes/fig/itk-snap-cropped-ct-aligned-1.png differ diff --git a/episodes/fig/itk-snap-cropped-ct-aligned-2.png b/episodes/fig/itk-snap-cropped-ct-aligned-2.png new file mode 100644 index 00000000..e1f70f01 Binary files /dev/null and b/episodes/fig/itk-snap-cropped-ct-aligned-2.png differ diff --git a/episodes/fig/itk-snap-cropped-ct.png b/episodes/fig/itk-snap-cropped-ct.png new file mode 100644 index 00000000..67c478c2 Binary files /dev/null and b/episodes/fig/itk-snap-cropped-ct.png differ diff --git a/episodes/fig/itk-snap-ct-pet-overlay.png b/episodes/fig/itk-snap-ct-pet-overlay.png new file mode 100644 index 00000000..d95711f0 Binary files /dev/null and b/episodes/fig/itk-snap-ct-pet-overlay.png differ diff --git a/episodes/fig/itk-snap-image-layer-inspector-color-map.png b/episodes/fig/itk-snap-image-layer-inspector-color-map.png new file mode 100644 index 00000000..207178d5 Binary files /dev/null and b/episodes/fig/itk-snap-image-layer-inspector-color-map.png differ diff --git a/episodes/fig/itk-snap-image-layer-inspector-contrast.png b/episodes/fig/itk-snap-image-layer-inspector-contrast.png new file mode 100644 index 00000000..3ef5d2f0 Binary files /dev/null and b/episodes/fig/itk-snap-image-layer-inspector-contrast.png differ diff --git a/episodes/fig/itk-snap-image-layer-inspector-general.png b/episodes/fig/itk-snap-image-layer-inspector-general.png new file mode 100644 index 00000000..6cd8649a Binary files /dev/null and b/episodes/fig/itk-snap-image-layer-inspector-general.png differ diff --git a/episodes/fig/itk-snap-image-layer-inspector-info.png b/episodes/fig/itk-snap-image-layer-inspector-info.png new file mode 100644 index 00000000..44763bae Binary files /dev/null and b/episodes/fig/itk-snap-image-layer-inspector-info.png differ diff --git a/episodes/fig/itk-snap-image-layer-inspector-metadata.png b/episodes/fig/itk-snap-image-layer-inspector-metadata.png new file mode 100644 index 00000000..8e079d03 Binary files /dev/null and b/episodes/fig/itk-snap-image-layer-inspector-metadata.png differ diff --git a/episodes/fig/itk-snap-info-16.png b/episodes/fig/itk-snap-info-16.png new file mode 100644 index 00000000..d94c1b98 Binary files /dev/null and b/episodes/fig/itk-snap-info-16.png differ diff --git a/episodes/fig/itk-snap-info-32.png b/episodes/fig/itk-snap-info-32.png new file mode 100644 index 00000000..416e3bcd Binary files /dev/null and b/episodes/fig/itk-snap-info-32.png differ diff --git a/episodes/fig/itk-snap-inhale-and-exhale-scans.png b/episodes/fig/itk-snap-inhale-and-exhale-scans.png new file mode 100644 index 00000000..2cdb9999 Binary files /dev/null and b/episodes/fig/itk-snap-inhale-and-exhale-scans.png differ diff --git a/episodes/fig/itk-snap-inhale-scan.png b/episodes/fig/itk-snap-inhale-scan.png new file mode 100644 index 00000000..b23cd646 Binary files /dev/null and b/episodes/fig/itk-snap-inhale-scan.png differ diff --git a/episodes/fig/itk-snap-install-1.png b/episodes/fig/itk-snap-install-1.png new file mode 100644 index 00000000..80cc8bd2 Binary files /dev/null and b/episodes/fig/itk-snap-install-1.png differ diff --git a/episodes/fig/itk-snap-install-2.png b/episodes/fig/itk-snap-install-2.png new file mode 100644 index 00000000..e1c83d96 Binary files /dev/null and b/episodes/fig/itk-snap-install-2.png differ diff --git a/episodes/fig/itk-snap-layout-pref-reminder.png b/episodes/fig/itk-snap-layout-pref-reminder.png new file mode 100644 index 00000000..da6e6378 Binary files /dev/null and b/episodes/fig/itk-snap-layout-pref-reminder.png differ diff --git a/episodes/fig/itk-snap-metadata-16.png b/episodes/fig/itk-snap-metadata-16.png new file mode 100644 index 00000000..8f749209 Binary files /dev/null and b/episodes/fig/itk-snap-metadata-16.png differ diff --git a/episodes/fig/itk-snap-metadata-32.png b/episodes/fig/itk-snap-metadata-32.png new file mode 100644 index 00000000..c5912a62 Binary files /dev/null and b/episodes/fig/itk-snap-metadata-32.png differ diff --git a/episodes/fig/itk-snap-misaligned-images.png b/episodes/fig/itk-snap-misaligned-images.png new file mode 100644 index 00000000..6cc0737f Binary files /dev/null and b/episodes/fig/itk-snap-misaligned-images.png differ diff --git a/episodes/fig/itk-snap-nifti-header.png b/episodes/fig/itk-snap-nifti-header.png new file mode 100644 index 00000000..451de5a7 Binary files /dev/null and b/episodes/fig/itk-snap-nifti-header.png differ diff --git a/episodes/fig/itk-snap-open-image-1.png b/episodes/fig/itk-snap-open-image-1.png new file mode 100644 index 00000000..39f8ace3 Binary files /dev/null and b/episodes/fig/itk-snap-open-image-1.png differ diff --git a/episodes/fig/itk-snap-open-image-2.png b/episodes/fig/itk-snap-open-image-2.png new file mode 100644 index 00000000..7fda4e42 Binary files /dev/null and b/episodes/fig/itk-snap-open-image-2.png differ diff --git a/episodes/fig/itk-snap-open-image-3.png b/episodes/fig/itk-snap-open-image-3.png new file mode 100644 index 00000000..a053277c Binary files /dev/null and b/episodes/fig/itk-snap-open-image-3.png differ diff --git a/episodes/fig/itk-snap-recent-images.png b/episodes/fig/itk-snap-recent-images.png new file mode 100644 index 00000000..18e32e93 Binary files /dev/null and b/episodes/fig/itk-snap-recent-images.png differ diff --git a/episodes/fig/itk-snap-save-ct-for-pet.png b/episodes/fig/itk-snap-save-ct-for-pet.png new file mode 100644 index 00000000..415616e3 Binary files /dev/null and b/episodes/fig/itk-snap-save-ct-for-pet.png differ diff --git a/episodes/fig/itk-snap-save-image-1.png b/episodes/fig/itk-snap-save-image-1.png new file mode 100644 index 00000000..72cc412a Binary files /dev/null and b/episodes/fig/itk-snap-save-image-1.png differ diff --git a/episodes/fig/itk-snap-tools-reorient-image.png b/episodes/fig/itk-snap-tools-reorient-image.png index 5a028137..c05a3934 100644 Binary files a/episodes/fig/itk-snap-tools-reorient-image.png and b/episodes/fig/itk-snap-tools-reorient-image.png differ diff --git a/episodes/fig/p3-entropies-values.png b/episodes/fig/p3-entropies-values.png new file mode 100644 index 00000000..0a98da7f Binary files /dev/null and b/episodes/fig/p3-entropies-values.png differ diff --git a/episodes/fig/p3-images.png b/episodes/fig/p3-images.png new file mode 100644 index 00000000..fdec15b0 Binary files /dev/null and b/episodes/fig/p3-images.png differ diff --git a/episodes/fig/p3-int8-intensities.png b/episodes/fig/p3-int8-intensities.png new file mode 100644 index 00000000..08cb49a4 Binary files /dev/null and b/episodes/fig/p3-int8-intensities.png differ diff --git a/episodes/fig/p3-mi-values.png b/episodes/fig/p3-mi-values.png new file mode 100644 index 00000000..37798838 Binary files /dev/null and b/episodes/fig/p3-mi-values.png differ diff --git a/episodes/fig/p3-msd-values.png b/episodes/fig/p3-msd-values.png new file mode 100644 index 00000000..eb99fa9c Binary files /dev/null and b/episodes/fig/p3-msd-values.png differ diff --git a/episodes/fig/p3-ncc-values.png b/episodes/fig/p3-ncc-values.png new file mode 100644 index 00000000..882e4dda Binary files /dev/null and b/episodes/fig/p3-ncc-values.png differ diff --git a/episodes/fig/p3-nmi-values.png b/episodes/fig/p3-nmi-values.png new file mode 100644 index 00000000..494eeddc Binary files /dev/null and b/episodes/fig/p3-nmi-values.png differ diff --git a/episodes/fig/p3-rotated-image.png b/episodes/fig/p3-rotated-image.png new file mode 100644 index 00000000..6b021b5a Binary files /dev/null and b/episodes/fig/p3-rotated-image.png differ diff --git a/episodes/fig/p3-ssd-values.png b/episodes/fig/p3-ssd-values.png new file mode 100644 index 00000000..6171a407 Binary files /dev/null and b/episodes/fig/p3-ssd-values.png differ diff --git a/episodes/fig/p4-comp-jac.png b/episodes/fig/p4-comp-jac.png new file mode 100644 index 00000000..664f9bb6 Binary files /dev/null and b/episodes/fig/p4-comp-jac.png differ diff --git a/episodes/fig/p4-def-field-comp.png b/episodes/fig/p4-def-field-comp.png new file mode 100644 index 00000000..412b2678 Binary files /dev/null and b/episodes/fig/p4-def-field-comp.png differ diff --git a/episodes/fig/p4-def-field-zoom.png b/episodes/fig/p4-def-field-zoom.png new file mode 100644 index 00000000..8a78024f Binary files /dev/null and b/episodes/fig/p4-def-field-zoom.png differ diff --git a/episodes/fig/p4-diff-images.png b/episodes/fig/p4-diff-images.png new file mode 100644 index 00000000..40aac712 Binary files /dev/null and b/episodes/fig/p4-diff-images.png differ diff --git a/episodes/fig/p4-displayed-images.png b/episodes/fig/p4-displayed-images.png new file mode 100644 index 00000000..5cbaf6db Binary files /dev/null and b/episodes/fig/p4-displayed-images.png differ diff --git a/episodes/fig/p4-jac-binarymask.png b/episodes/fig/p4-jac-binarymask.png new file mode 100644 index 00000000..03b92f73 Binary files /dev/null and b/episodes/fig/p4-jac-binarymask.png differ diff --git a/episodes/fig/p4-reg-first-step.png b/episodes/fig/p4-reg-first-step.png new file mode 100644 index 00000000..2d90caa2 Binary files /dev/null and b/episodes/fig/p4-reg-first-step.png differ diff --git a/episodes/fig/p4-reg-output-numlev1.png b/episodes/fig/p4-reg-output-numlev1.png new file mode 100644 index 00000000..6ff54299 Binary files /dev/null and b/episodes/fig/p4-reg-output-numlev1.png differ diff --git a/episodes/fig/p4-reg-output-sigelastic0.5.png b/episodes/fig/p4-reg-output-sigelastic0.5.png new file mode 100644 index 00000000..33dc26d6 Binary files /dev/null and b/episodes/fig/p4-reg-output-sigelastic0.5.png differ diff --git a/episodes/fig/p4-reg-output-sigelastic0.png b/episodes/fig/p4-reg-output-sigelastic0.png new file mode 100644 index 00000000..66e11eda Binary files /dev/null and b/episodes/fig/p4-reg-output-sigelastic0.png differ diff --git a/episodes/fig/p4-reg-output-sigfluid0.png b/episodes/fig/p4-reg-output-sigfluid0.png new file mode 100644 index 00000000..300b6f17 Binary files /dev/null and b/episodes/fig/p4-reg-output-sigfluid0.png differ diff --git a/episodes/fig/p4-reg-output.png b/episodes/fig/p4-reg-output.png new file mode 100644 index 00000000..81e52cf5 Binary files /dev/null and b/episodes/fig/p4-reg-output.png differ diff --git a/episodes/fig/p4-reg-output2-numlev1.png b/episodes/fig/p4-reg-output2-numlev1.png new file mode 100644 index 00000000..e539fe18 Binary files /dev/null and b/episodes/fig/p4-reg-output2-numlev1.png differ diff --git a/episodes/fig/p4-reg-output2-sigelastic0.5.png b/episodes/fig/p4-reg-output2-sigelastic0.5.png new file mode 100644 index 00000000..7c8df108 Binary files /dev/null and b/episodes/fig/p4-reg-output2-sigelastic0.5.png differ diff --git a/episodes/fig/p4-reg-output2-sigelastic0.png b/episodes/fig/p4-reg-output2-sigelastic0.png new file mode 100644 index 00000000..d28011cd Binary files /dev/null and b/episodes/fig/p4-reg-output2-sigelastic0.png differ diff --git a/episodes/fig/p4-reg-output2-sigfluid0.png b/episodes/fig/p4-reg-output2-sigfluid0.png new file mode 100644 index 00000000..13d155c3 Binary files /dev/null and b/episodes/fig/p4-reg-output2-sigfluid0.png differ diff --git a/episodes/fig/p4-reg-output2.png b/episodes/fig/p4-reg-output2.png new file mode 100644 index 00000000..8b7bf789 Binary files /dev/null and b/episodes/fig/p4-reg-output2.png differ diff --git a/episodes/fig/registration_animation.gif b/episodes/fig/registration_animation.gif new file mode 100644 index 00000000..5481ed5e Binary files /dev/null and b/episodes/fig/registration_animation.gif differ diff --git a/episodes/fig/registration_animation_comp.gif b/episodes/fig/registration_animation_comp.gif new file mode 100644 index 00000000..95708e63 Binary files /dev/null and b/episodes/fig/registration_animation_comp.gif differ diff --git a/episodes/fig/registration_animation_elastic0.5fluid1.gif b/episodes/fig/registration_animation_elastic0.5fluid1.gif new file mode 100644 index 00000000..3ffce01e Binary files /dev/null and b/episodes/fig/registration_animation_elastic0.5fluid1.gif differ diff --git a/episodes/fig/registration_animation_elastic0.gif b/episodes/fig/registration_animation_elastic0.gif new file mode 100644 index 00000000..d5c7ae7f Binary files /dev/null and b/episodes/fig/registration_animation_elastic0.gif differ diff --git a/episodes/fig/registration_animation_fluid0.gif b/episodes/fig/registration_animation_fluid0.gif new file mode 100644 index 00000000..222a7fb7 Binary files /dev/null and b/episodes/fig/registration_animation_fluid0.gif differ diff --git a/episodes/fig/registration_animation_level1.gif b/episodes/fig/registration_animation_level1.gif new file mode 100644 index 00000000..1957520c Binary files /dev/null and b/episodes/fig/registration_animation_level1.gif differ diff --git a/episodes/fig/rotated_images.gif b/episodes/fig/rotated_images.gif new file mode 100644 index 00000000..19ae5322 Binary files /dev/null and b/episodes/fig/rotated_images.gif differ diff --git a/episodes/fig/trans_diff_image.png b/episodes/fig/trans_diff_image.png new file mode 100644 index 00000000..bc0ebabf Binary files /dev/null and b/episodes/fig/trans_diff_image.png differ diff --git a/episodes/fig/trans_diff_image_10.5.png b/episodes/fig/trans_diff_image_10.5.png new file mode 100644 index 00000000..1673f5dd Binary files /dev/null and b/episodes/fig/trans_diff_image_10.5.png differ diff --git a/episodes/fig/trans_diff_image_intlims.png b/episodes/fig/trans_diff_image_intlims.png new file mode 100644 index 00000000..8b9807fc Binary files /dev/null and b/episodes/fig/trans_diff_image_intlims.png differ diff --git a/episodes/fig/trans_disp_image.png b/episodes/fig/trans_disp_image.png new file mode 100644 index 00000000..fa5ca34a Binary files /dev/null and b/episodes/fig/trans_disp_image.png differ diff --git a/episodes/fig/trans_disp_image_translated.png b/episodes/fig/trans_disp_image_translated.png new file mode 100644 index 00000000..4e38222d Binary files /dev/null and b/episodes/fig/trans_disp_image_translated.png differ diff --git a/episodes/fig/trans_rot_final_image.gif b/episodes/fig/trans_rot_final_image.gif new file mode 100644 index 00000000..10b1617b Binary files /dev/null and b/episodes/fig/trans_rot_final_image.gif differ diff --git a/episodes/fig/trans_rot_final_image.png b/episodes/fig/trans_rot_final_image.png new file mode 100644 index 00000000..0af7436b Binary files /dev/null and b/episodes/fig/trans_rot_final_image.png differ diff --git a/episodes/fig/trans_rot_final_image_compose_linear.gif b/episodes/fig/trans_rot_final_image_compose_linear.gif new file mode 100644 index 00000000..32168819 Binary files /dev/null and b/episodes/fig/trans_rot_final_image_compose_linear.gif differ diff --git a/episodes/fig/trans_rot_final_image_compose_nearest.gif b/episodes/fig/trans_rot_final_image_compose_nearest.gif new file mode 100644 index 00000000..89fe377a Binary files /dev/null and b/episodes/fig/trans_rot_final_image_compose_nearest.gif differ diff --git a/episodes/fig/trans_rot_final_image_compose_splinef2d.gif b/episodes/fig/trans_rot_final_image_compose_splinef2d.gif new file mode 100644 index 00000000..ff5633d7 Binary files /dev/null and b/episodes/fig/trans_rot_final_image_compose_splinef2d.gif differ diff --git a/episodes/fig/trans_rot_final_image_pad0.gif b/episodes/fig/trans_rot_final_image_pad0.gif new file mode 100644 index 00000000..5c12dc36 Binary files /dev/null and b/episodes/fig/trans_rot_final_image_pad0.gif differ diff --git a/episodes/fig/trans_rot_final_image_pad0.png b/episodes/fig/trans_rot_final_image_pad0.png new file mode 100644 index 00000000..2954b69e Binary files /dev/null and b/episodes/fig/trans_rot_final_image_pad0.png differ diff --git a/episodes/fig/trans_rot_final_image_pad0_nn.gif b/episodes/fig/trans_rot_final_image_pad0_nn.gif new file mode 100644 index 00000000..bc9d6713 Binary files /dev/null and b/episodes/fig/trans_rot_final_image_pad0_nn.gif differ diff --git a/episodes/fig/trans_rot_final_image_pad0_spline.gif b/episodes/fig/trans_rot_final_image_pad0_spline.gif new file mode 100644 index 00000000..41db5e92 Binary files /dev/null and b/episodes/fig/trans_rot_final_image_pad0_spline.gif differ diff --git a/episodes/fig/trans_rot_final_image_push.gif b/episodes/fig/trans_rot_final_image_push.gif new file mode 100644 index 00000000..8a3cbbce Binary files /dev/null and b/episodes/fig/trans_rot_final_image_push.gif differ diff --git a/episodes/fig/trans_rot_image.png b/episodes/fig/trans_rot_image.png new file mode 100644 index 00000000..672bd083 Binary files /dev/null and b/episodes/fig/trans_rot_image.png differ diff --git a/episodes/fig/trans_rot_image_nn_spline.png b/episodes/fig/trans_rot_image_nn_spline.png new file mode 100644 index 00000000..141d9f92 Binary files /dev/null and b/episodes/fig/trans_rot_image_nn_spline.png differ diff --git a/episodes/fig/vscode-ideas-reg.png b/episodes/fig/vscode-ideas-reg.png new file mode 100644 index 00000000..d6e1106a Binary files /dev/null and b/episodes/fig/vscode-ideas-reg.png differ diff --git a/episodes/fig/vscode-python-env.png b/episodes/fig/vscode-python-env.png new file mode 100644 index 00000000..f448826a Binary files /dev/null and b/episodes/fig/vscode-python-env.png differ diff --git a/episodes/fig/vscode-select-kernel.png b/episodes/fig/vscode-select-kernel.png new file mode 100644 index 00000000..0bda1512 Binary files /dev/null and b/episodes/fig/vscode-select-kernel.png differ diff --git a/episodes/fig/vscode.png b/episodes/fig/vscode.png new file mode 100644 index 00000000..623b9e83 Binary files /dev/null and b/episodes/fig/vscode.png differ diff --git a/episodes/itksnap.Rmd b/episodes/itksnap.Rmd deleted file mode 100644 index b91686e2..00000000 --- a/episodes/itksnap.Rmd +++ /dev/null @@ -1,938 +0,0 @@ ---- -title: 'A Practical Introduction to Working with Medical Image Data' -teaching: 10 -exercises: 2 ---- - -:::::::::::::::::::::::::::::::::::::: questions - -- How can we visualize volumetric data? -- How do locations in an image get mapped to real-world coordinates? -- How can we register different types of images together using ITK-SNAP? - -:::::::::::::::::::::::::::::::::::::::::::::::: - -::::::::::::::::::::::::::::::::::::: objectives - -- Describe the structure of medical imaging data and popular formats (DICOM and NifTi) -- Discover and display medical imaging data in ITK-SNAP -- Demonstrate how to convert between different medical image formats -- Inspecting the NifTi header with `nibabel` -- Manipulate image data (cropping images) with `nibabel` and learn how it is displayed on screen with ITK-SNAP - -:::::::::::::::::::::::::::::::::::::::::::::::: - -## Medical Imaging Data - -### The Cancer Imaging Archive (TCIA) -In these exercises we will be working with real-world medical imaging data from the Cancer Imaging Archive [(TCIA)](https://www.cancerimagingarchive.net/collection/ct-vs-pet-ventilation-imaging/). -TCIA is a resource of public datasets hosts a large archive of medical images of cancer accessible for public download. -We are using CT Ventilation as a Functional Imaging Modality for Lung Cancer Radiotherapy (also known as `CT-vs-PET-Ventilation-Imaging dataset`) from TCIA. -The CT-vs-PET-Ventilation-Imaging collection is distributed under the Creative Commons Attribution 4.0 International License (https://creativecommons.org/licenses/by/4.0/). -The `CT-vs-PET-Ventilation-Imaging` dataset contains 20 lung cancer patients underwent exhale/inhale breath hold CT (BHCT), free-breathing four-dimensional CT (4DCT) and Galligas PET ventilation scans in a single session on a combined 4DPET/CT scanner. - -We recommend using the BHCT scans for participants `CT-PET-VI-02` and `CT-PET-VI-03` showing very little motion between inhale and exhale (for PET scan and accompanying CT scan), and BHCT scans for participants `CT-PET-VI-05` and `CT-PET-VI-07` (the inhale and exhale CT scans), where 05 has a different number of slices between the inhale and exhale. - -Data paths contain: (A) `inhale_BH_CT` and `exhale_BH_CT` contain CT scans acquired during an inhalation breath hold and exhalation breath hold respectively, -(B) `PET` contains a PET scan that measures the local lung function, and (C) `CT_for_PET` contains a CT scan acquired at the same time as the PET scan for attenuation correction of the PET scan and to provide anatomical reference for the PET data. - - -* CT-PET-VI-02 -```bash -├── CT-PET-VI-02 -│   ├── CT_for_PET -│   ├── exhale_BH_CT -│   ├── inhale_BH_CT -│   └── PET -``` - -* CT-PET-VI-03 -```bash -/CT-PET-VI-03$ tree -d -├── CT-PET-VI-03 -│   ├── CT_for_PET -│   └── PET -``` - -* CT-PET-VI-05 - -```bash -CT-PET-VI-05$ tree -d -. -├── CT-PET-VI-05 -│   ├── CT_for_PET -│   ├── exhale_BH_CT -│   ├── inhale_BH_CT -│   └── PET -``` - -* CT-PET-VI-07 -```bash -/CT-PET-VI-07$ tree -d -. -└── CT-PET-VI-07 - ├── exhale_BH_CT - ├── inhale_BH_CT - └── PET -``` - -For example, datasets in ITK-SNAP are illustrated in the following figures. -![Figure. CT-PET-VI-02 dataset](fig/CT-PET-VI-02.svg) - -Refer to [`Summary and Setup`](https://healthbioscienceideas.github.io/Medical-Image-Registration-Short-Course/#data-sets) section for further details on CT-vs-PET-Ventilation-Imaging dataset. - - - -:::::::::::::::::::::::::::::::::::::: spoiler - -### Data sharing ethics and policies -1. Data may require extensive cleaning and pre-processing - before it is suitable to be used. -2. When using open datasets that contain images of humans, - such as those from TCIA, for your own research, - you will still need to get ethical approval - as you do for any non-open datasets you use). - This can usually be obtained from the local ethics - committee (both Medical Physics and Computer Science have - their own local ethics committee) or from - UCL central ethics if your department does not have one. -3. Although there are many formats that medical imaging - data can come in, we are going to focus on DICOM and NifTi - as they are two of the most common formats. - Most data that comes from a hospital will be in - DICOM format, whereas NifTi is a very popular format - in the medical image analysis community. - -:::::::::::::::::::::::::::::::::::::: - -### DICOM format -Digital Imaging and Communications in Medicine (DICOM) is a technical standard for the digital storage and transmission of medical images and related information. -While DICOM images typically have a separate file for every slice, more modern DICOM images can come with all slices in a single file. - -If you look at the data `CT-PET-VI-02/CT_for_PET` you will see there are 175 individual files corresponding to 175 slices in the volume. -In addition to the image data for each slice, each file contains a header which can contain an extensive amount of extra information relating to the scan and subject. - -In a clinical setting this will include patient identifiable information such as their name, address, and other relevant details. -Such information should be removed before the data is transferred from the clinical network for use in research. -If you ever discover patient identifiable information in the header of data you are using you should immediately alert your supervisor, manager or collaborator - -The majority of the information in the DICOM header is not directly useful for typical image processing and analysis tasks. -Furthermore, there are complicated ‘links’ (provided by unique identifiers, UIDs) between the DICOM headers of different files belonging to the same scan or subject. Together with DICOM routinely storing each slice as a separate files, it makes processing an entire imaging volume stored as DICOM format rather cumbersome, and the extra housekeeping required could lead to a greater chance of an error being made. -Therefore, a common first step of any image processing pipeline is to convert the DICOM image to a more suitable format such as `NifTi`. -Generally, most conversions go from DICOM to NIfTI. There are scenarios when you might want to convert from NIfTI **back** to DICOM, for example, if you need to import them into a clinical system that only works with DICOM. - - - -:::::::::::::::::::::::::::::::::::::: spoiler - -### Warning on converting images back to DICOM -Converting images back to DICOM such that -they are correctly interpreted by a clinical system -can be very tricky and requires a good understanding -of the DICOM standard. -More information on the DICOM standard -can be found here: https://www.dicomstandard.org - -:::::::::::::::::::::::::::::::::::::: - -### NifTi format -The Neuroimaging Informatics Technology Initiative (NIfTI) image format are usually stored as a single file containing the imaging data and header information. -While the NifTI file format was originally developed by the neuroimaging community it is not specific to neuroimaging and is now widely used for many different medical imaging applications outside of the brain. -During these exercises you will learn about some of the key information stored in the NifTi header. - -NifTI files can also store the header and image data in separate files but this is not very common. -For more information, please see Anderson Winkler's blog on the [NIfTI-1]( https://brainder.org/2012/09/23/the-nifti-file-format/) and[NIfTI-2]( https://brainder.org/2015/04/03/the-nifti-2-file-format/) formats. - -## Visualisation of volumetric data with ITK-SNAP - -### Getting Started with ITK-SNAP -#### Introduction -ITK-SNAP started in 1999 with SNAP (SNake Automatic Partitioning) and developed by Paul Yushkevich with the guidance of Guido Gerig. -ITK-SNAP is open-source software distributed under the GNU General Public License. -ITK-SNAP is written in C++ and it leverages the Insight Segmentation and Registration Toolkit (ITK) library. -See `Summary and Setup` section for requirements and installation of ITK-SNAP - -### ITK-SNAP application -ITK-SNAP application shows three orthogonal slices and a fourth window for three-dimensional view segmentation. - -![Figure. Description of ITK-SNAP tools in the main window using [MRI-crop example datasets](http://www.nitrc.org/frs/download.php/750/MRI-crop.zip) for plane and 3D views.](fig/itk-snap-main-window.svg) - - -## Viewing and understanding DICOM data -We are using CT Lung DICOM dataset from `CT-PET-VI-02/CT_for_PET`, containing 175 dicom items and totalling 92.4 MB. - -### Opening and Viewing DICOM images -* **Open ITK-SNAP application** -* **Load inhale scan using file browser** - - Open the itk-snap application. If you have used it before it will display recent images when you open it, and you can select one of these to open it. - - If the image you want to open is not listed under recent images, you can open it either by clicking the 'Open Image...' button at the bottom right of the window, or File->Open Main Image. - When you do this a small window opens where you can enter the path of the file you want to open, or click Browse in your file explorer (exFinder in Mac and nautilus in Ubuntu) to open a file browser (for instance, "CT-PET-VI-02/CT_for_PET"). - Do this and select the first slice from the Inhale_BH_CT folder, which is called 1-001.dcm. The set the File Format to DICOM Image Series and then click `Next`. - You will then be asked to Select DICOM series to open, but as this folder only includes one series you can just click `Next`. - Click `Finish` - -![Figure. ITK-SNAP opening image file](fig/itk-snap-opening-dicoms.svg) - -![Figure. ITK-SNAP main window](fig/itk-snap-viewing-dicoms.svg) - - -### Navigating DICOM images -* Navigate around image, see intensity value, and adjust intensity window -* To navigate around the image, i.e. change the displayed slices, you can left clicking in the displayed slices, using the mouse wheel, or the slider next to the displayed slices. -* To zoom, you can use the right mouse button (or control + button 1), and pan using the centre mouse button (or alt + button 1). - -![Figure. ITK-SNAP cursor inspector](fig/itk-snap-cursor-inspector.png) - -### Handling multiple DICOM images -We provide instructions for handling multiple DICOM images, applying overlays, viewing image information, and using color maps in ITK-SNAP. -For this exmaple, we are using `CT-PET-VI-02/inhale_BH_CT` and `CT-PET-VI-02/exhale_BH_CT` dicoms. - -#### Load additional image using `exhale_BH_CT` scan by drag and drop -- Images can also be opened by simply dragging them on to the ITK-SNAP window. -In a file explorer (Finder in Mac) go to the `exhale_BH_CT` folder and drag any of the dicom files on to the itksnap window. -A small window will pop up asking What should itksnap do with this image? -If you select Load as Main Image it will replace the current image. -Instead select Load as Additional Image - this will enable easy comparison of the two images. - -![Figure. ITK-SNAP additional image](fig/itk-snap-additional-image.svg) - -#### Using thumbnails to see image differences -* Once the second image is loaded you will see two thumbnails in the top corner of each slice display -* Click on these to swap the displayed slices from one image to the other - this will enable you to easily see the similarities and differences between the images. -* You will also see that both images now appear in the Cursor Inspector, where you can see the intensity value at the cursor for both images. -* You can also change the displayed image by clicking on it in the Cursor Inspector - -![Figure. ITK-SNAP thumbnails, handlding two DICOM images (`inhale_BH_CT` and `exhale_BH_CT` for `CT-PET-VI-02`)](fig/itk-snap-thumbnails.png) - -#### Colour overlay -In some cases you may want to display an image as a colour overlay on top of another image, e.g. displaying a PET image overlaid on the corresponding CT. - -##### Colour overlay as a separate image (shown besides other images) -* **Load Primary Image**: - - Open your primary DICOM series of the `CT-PET-VI-02/CT_for_PET` image either using the drag and drop method or `File->Open Image` as described above. - - Make sure File Format is DICOM Image Series. - - This time select Load as Main Image. This will load in the new image replacing the two that were previously loaded. -* **Load Overlay Image**: - - Load the `CT-PET-VI-02/PET` image either using the drag and drop method or `File->Open Image` - -![Figure. ITK-SNAP thumbnails](fig/itk-snap-colour-overlay-as-separate-image.svg) - -##### Colour overlay as a semi-transparent overlay (shown on the top of other images) -* **Load Primary Image**: - - Open your primary DICOM series of the `CT-PET-VI-02/CT_for_PET` image either using the drag and drop method or `File->Open Image` as described above. - - Make sure File Format is DICOM Image Series. - - Select Load as Main Image. -* **Load Overlay Image**: - - Load the `CT-PET-VI-02/PET` image either using the drag and drop method or `File->Open Image` - - When the image has loaded select As a semi-transparent overlay. You can also set the overlay color map here - select `Hot` from the drop-down menu and click Next. - - The image summary will then be displayed, along with a warning about possible loss of precision (which can be ignored). Click Finish - -![Figure. ITK-SNAP semi-transparent overlays](fig/itk-snap-colour-overlay-as-semi-transparent-image.svg) - - -### Layer inspector tools -* Load `CT-PET-VI-02/CT_for_PET` and `CT-PET-VI-02/PET` scan datasets, using hot semi-transparency (as shown above). -* Go to `Tools->Layer Inspector`. This opens the Image Layer Inspector window. You will see the two images on the left, and five tabs along the top of the window. -* PET image appears in Cursor Inspector (called Galligas Lung) in italic and CT image as CT Lung in bold. - -* **The General tab** displays the filename and the nickname for the selected image. - * The Nickname can be edited. For the PET image the general tab also displays a slider to set the opacity, and the option to display it as a separate image or semi-transparent overlay. - -![Figure. ITK-SNAP: Layer inspector with general tab](fig/itk-snap-layer-inspector-general-tab.png) - -* **The Contrast tab** enables you to adjust how the image intensities are displayed for the different images, e.g. Select the lung image and set the maximum value to 1000 (push enter after typing the number) and this will make the structures in the lung easier to see. - -![Figure. ITK-SNAP: Layer inspector with contrast tab](fig/itk-snap-layer-inspector-contrast-tab.png) - -* **The color map tab** can be used to select and manipulate the colour map - * One of the key elements of an imaging viewer is to provide different means to map those values into grey scales or different colors. - * We will show you how to apply different color maps or lookup tables to your data and how this affects how the images are presented to the user. - - 1. **Load the Image**: - - Open the desired DICOM series. - - 2. **Open Color Map Settings**: - - Go to `Tools` > `Color Maps Editor`. - - 3. **Apply and Adjust Color Map**: - - Choose a predefined color map from the list. You might also create a customised map. - - Adjust the intensity and transparency settings as needed to enhance the visualisation of the images. - -![Figure. ITK-SNAP: Layer inspector with color map tab](fig/itk-snap-layer-inspector-color-map-tab.png) - -* **The info tab** displays the Image Header information for the images and also gives the Cursor Coordinates in both Voxel and World units. - * A dialog box will appear displaying metadata and other relevant information about the loaded DICOM images (e.g. values for image header and cursor coordinates). - * If you look at the info tab for both the CT image and PET image you will see that the header information for the images is different, and the cursor coordinates in voxel units for the PET image are not whole numbers. - * This is because the displayed slices are from the main (CT) image, and the overlaid (PET) image slices are interpolated at the corresponding locations. - -![Figure. ITK-SNAP: Layer inspector with info tab](fig/itk-snap-layer-inspector-info-tab.png) - -* **The metadata tab** contains (some of) the information from the DICOM header of the images - -![Figure. ITK-SNAP: Layer inspector with metadata tab](fig/itk-snap-layer-inspector-metadata-tab.png) - - - - - - - -### Converting DICOM images to NifTi -As mentioned earlier in the exercise, the NIfTI image format tends to be much easier to work with when processing and analysing medical image data. -We will now work on converting DICOM images to a NIfTI image volume. - -#### Using ITK-SNAP -In the ITK-SNAP application, save the file by clicking 'File' --> 'Save image' --> rename image and choose format as 'NiFTI'. - - - -:::::::::::::::::::::::::::::::::::::: spoiler - -## Warning for Linux users. -You might get this error `Error: exception occurred during image IO. Exception: std::expedition` to which we suggest using dicom2nifti library. - -:::::::::::::::::::::::::::::::::::::: - -#### Using [dicom2nifti](https://github.com/icometrix/dicom2nifti) -1. Open terminal, activate `mirVE` environment and run python -```bash -conda activate mirVE -python -``` -2. Run the following commands to convert your DICOM scans to nifti format. -```python -from pathlib import Path -import dicom2nifti -dicom_path = Path("CT_for_PET") -dicom2nifti.convert_directory(dicom_path, ".", compression=True, reorient=True) -pet_dicom_path = Path("PET") -dicom2nifti.convert_directory(pet_dicom_path, ".", compression=True, reorient=True) -``` -Creating the following files: -``` -`3_ct_lung__30__b31f.nii.gz` [47M] and renamed as `ct_for_pet.nii.gz` -`4_galligas_lung.nii.gz` [12M] and renamed as `pet.nii.gz` -``` - -#### Using others packages -You might also be interested to look other packages: -* https://github.com/rordenlab/dcm2niix -* https://nipy.org/nibabel/dicom/dicom.html#dicom -* https://neuroimaging-cookbook.github.io/recipes/dcm2nii_recipe/ - -## Viewing and understanding the NifTi header with NiBabel -We are going to use the python package [NiBabel](https://nipy.org/nibabel/) to upload `NifTI` images and learn some basic properties of the image using Nibabel. -[`nibabel` api](https://nipy.org/nibabel/api.html#api) provides various image utilities, conversions methods, helpers. -We recommend to checking [documentation](https://nipy.org/nibabel/#documentation) for further deatils on using nibabel library. -Please see [instructions](https://github.com/HealthBioscienceIDEAS/Medical-Image-Registration-Short-Course/tree/main/_dependencies) to install `NiBabel` package and other dependecies. - -### Loading images -1. Open terminal, activate `mirVE` environment and run python -```bash -conda activate mirVE -python -``` - -2. Importing python packages under `*/episodes` path - -A NifTi image with extension `*.nii.gz` can be read using `nibabel`’s `load` function. -```python -import numpy as np -import nibabel as nib -import matplotlib.pyplot as plt -nii3ct = nib.load("data/ct_for_pet.nii.gz") -``` - -```callout -`NiBabel`'s `load` function does not actually read the image data itself from disk, but -does read and interpret the header, and provides functions for accessing the image data. -Calling these functions reads the data from disk (unless it has already been read, in which -case it may be cached in memory already. -``` - -2. `NifTI` object type, shape/size -The Nifti1Image object allows us to see the size/shape of the image and its data type. -```python -print(type(nii3ct)) -# -``` - -```python -print(nii3ct.get_data_dtype()) -#int16 -``` - -```python -nii3ct.shape -# (512, 512, 175) -``` - -3. Affine transform, mapping from voxel space to world space -There are two common ways of specifying the affine transformation mapping from voxel coordinates -to world coordinates in the nifti header, called the `sform` and the `qform'. -Another transformation is that neither the sform or qform is used the mapping is performed based just on the voxel dimensions. -The **sform** directly stores the affine transformation in the header -whereas the **qform** stores the affine transformation using [quaternions](https://en.wikipedia.org/wiki/Quaternion). - -:::::::::::::::::::::::::::::::::::::: spoiler - -### Recommending `sform` over `qform` -We strongly advise always using the `sform` over the `qform`, as: -* It is easier to understand (at least for instructors). -* It can contain shears which cannot be represented using the `qform`. - -:::::::::::::::::::::::::::::::::::::: - -```python -nii3ct.affine -#array([[ -0.9765625 , 0. , 0. , 249.51171875], -# [ -0. , 0.9765625 , 0. , -33.01171875], -# [ 0. , -0. , 2. , -553.5 ], -# [ 0. , 0. , 0. , 1. ]]) -``` - -### Header features -1. Printing `nii3ct.header` outputs -```python -print(nii3ct.header) - -# object, endian='<' -#sizeof_hdr : 348 -#data_type : b'' -#db_name : b'' -#extents : 0 -#session_error : 0 -#regular : b'' -#dim_info : 0 -#dim : [ 3 512 512 175 1 1 1 1] -#intent_p1 : 0.0 -#intent_p2 : 0.0 -#intent_p3 : 0.0 -#intent_code : none -#datatype : int16 -#bitpix : 16 -#slice_start : 0 -#pixdim : [-1. 0.9765625 0.9765625 2. 1. 1. -# 1. 1. ] -#vox_offset : 0.0 -#scl_slope : nan -#scl_inter : nan -#slice_end : 0 -#slice_code : unknown -#xyzt_units : 2 -#cal_max : 0.0 -#cal_min : 0.0 -#slice_duration : 0.0 -#toffset : 0.0 -#glmax : 0 -#glmin : 0 -#descrip : b'' -#aux_file : b'' -#qform_code : unknown -#sform_code : aligned -#quatern_b : 0.0 -#quatern_c : 1.0 -#quatern_d : 0.0 -#qoffset_x : 249.51172 -#qoffset_y : -33.01172 -#qoffset_z : -553.5 -#srow_x : [ -0.9765625 0. 0. 249.51172 ] -#srow_y : [ -0. 0.9765625 0. -33.01172 ] -#srow_z : [ 0. -0. 2. -553.5] -#intent_name : b'' -#magic : b'n+1' -``` -* `sform` transformation matrix -You can see that the `sform` is stored in the `srow_x`, `srow_y`, and `srow_z` fields of the header. -These specify the top three rows of a 4x4 matrix representing an affine transformation using homogeneous coordinates. -The fourth row is not stored as it is always `[0 0 0 1]`. -You will notice that the sform matrix matches the affine transform in the `Nifti1Image` object. - -* `sform_code` and `qform_code` -It can also be seen from the `sform_code` and `qform_code` that both the sform and qform are set for this image. -There are 5 different sform/qform codes defined: - -| Code | Label | Meaning | -| --- | --- | --- | -| 0 | unknown | sform not defined | -| 1 | scanner | RAS+ is scanner coordinates | -| 2 | aligned | RAS+ aligned to some other scan | -| 3 | talairach | RAS+ in Talairach atlas space | -| 4 | mni | RAS+ in MNI atlas space | - - -:::::::::::::::::::::::::::::::::::::: spoiler - -### Which code to use? -But for most purposes all that matters is whether the code is unknown (numerical value 0) or one of the other valid values. -If the code is unknown then the sform/qform is not specified (and any values provided in the sform/qform fields of the header will be ignored). -For any other valid values the sform/qform is specified and the affine transform will be determined from the corresponding header values. - -:::::::::::::::::::::::::::::::::::::: - -### `qform` transform -You can also see `qform` transform using the `get_qform` function, which returns the affine matrix represented by the qform: -```python -print(nii3ct.get_qform()) -#[[ -0.9765625 0. 0. 249.51171875] -# [ 0. 0.9765625 0. -33.01171875] -# [ 0. 0. 2. -553.5 ] -# [ 0. 0. 0. 1. ]] -``` - - - -:::::::::::::::::::::::::::::::::::::: spoiler - -### Advincing on using `sform` -However, the nifti format does not require that the sform and qform specify the same matrix. -It is also not well defined what should be done when they are both provided and are different from each other. -Often the `qform` is simply ignored and the `sform` is used, so the general advice is to only use the `sform` and set the `qform` to unknown to avoid any confusion. - -:::::::::::::::::::::::::::::::::::::: - - -### Coordinate system -#### LPS -The original DICOM images assume the world coordinate system is LPS (i.e. values increase when moving to the left, posterior, and superior). -See this helpful information on [DICOM orientations](https://nipy.org/nibabel/dicom/dicom_orientation.html), and determine if you obtained a similar affine transform from the DICOM headers. -The values in the first two rows corresponding to the x and y dimensions would be the negative of those in the NifTi header. - -#### RAS -You will notice that the diagonal elements of the affine matrix match the voxel dimensions -(as the affine matrix does not contain any rotations), but the values for the x and y -dimensions are negative. This is because the voxel indices increase as you move to the left -and posterior of the image/patient, but the nifti format assumes a RAS world coordinate -system (i.e. the values increase as you move to the right and anterior of the patient). -Therefore, the world coordinates decrease as the voxel indices increase. - -#### ITK-SNAP visualisation coordinate system -ITK-SNAP provides the option, `Tools>Reorient Image`, to set image orientation using three-letter code (**RAI** code), describing the image rows, columns, and slice map to the anatomical coordinates. -For example, the code **ASL** means that image rows (X) run from **A**nterior to **P**osterior, image columns (Y) run from **S**uperior to **I**nferior, and images slices (Z) run from **L**eft to **R**ight. - -![Figure. ITK-SNAP Reorient image](fig/itk-snap-tools-reorient-image.png) - -## Modifying `NifTi` images and headers with Nibabel - -### Obtaining image pixel data -The image data for a Nifti1Image object can be accessed using the `get_fdata` function. -This will load the data from disk and cast it to float64 type before returning the data as a `numpy` array. -```python -image_data=nii3ct.get_fdata() -print(type(image_data)) -# - -print(image_data.dtype, image_data.shape) -#float64 (512, 512, 175) -``` - -You can also plot a particular slide, (e.g. 140) for such data. -```python -plt.imshow(image_data[:,:,140], cmap="gray"); plt.show() -#launch plot window -``` - -### Using float64 data type -Casting image data to float64 prevents any integer-related errors and problems in downstream processing when using the data read by NiBabel. -However, this does not change the type of the image stored on disk, which as we have seen is int16 for this image, and this could still lead to downstream errors or unexpected behaviour for any processing that works directly on the images stored on disk. -We are going to convert the image saved on disk to a floating point image, using the `set_data_dtype` function to set the data type to float32. -Essentially, we are going to use 32 bit rather than 64 bit floating point, as this is usually sufficiently accurate. - -```python -print(nii3ct.get_data_dtype()) -#int16 -``` - -```python -nii3ct.set_data_dtype(np.float32) -print(nii3ct.get_data_dtype()) -#float32 -``` - -We are also going to set the `qform` to `unknown`, updating the corresponding fields in the header. -We use the `set_qform` function to set the `qform` code. -To just set the code and not modify the other `qform` values in the header, the first input should be `None`. -It is not necessary to modify the other `qform` values as they should be ignored when the code is set to unknown. - -```python -nii3ct.set_qform(None, code = 'unknown') -print(nii3ct.header) -... -qform_code : unknown -... - -``` - -The floating point image can now be written to disk using the save function. - - -:::::::::::::::::::::::::::::::::::::: spoiler - -### We do not need to manually modify the type of the image data? -Setting the data type for the Nifti1Image object tells it which type to use when writing the image data to disk. - -:::::::::::::::::::::::::::::::::::::: - -```python -nib.save(nii3ct, 'data/ct_for_pet_float32.nii.gz') -``` - -### 16-bit integers vs 32 bit floats - -The file size of `ct_for_pet_float32.nii.gz` is larger than `ct_for_pet.nii.gz` (61 MB, 47 MB, respectively) because of the conversion from 16 bit integers to 32 bit floats. -Both images are compressed, but the compression is more efficient for the floating point image: while a 32-bit floating point variable is twice the size of a 16-bit integer, the floating image is less than double the size of the original integer image. - -* int16 -```python -nii3ct_int16 = nib.load("data/ct_for_pet.nii.gz") -``` - -```python -print(nii3ct_int16.get_data_dtype()) -#int16 -``` - -```python -print(nii3ct_int16.header) -# object, endian='<' -#sizeof_hdr : 348 -#data_type : b'' -#db_name : b'' -#extents : 0 -#session_error : 0 -#regular : b'' -#dim_info : 0 -#dim : [ 3 512 512 175 1 1 1 1] -#intent_p1 : 0.0 -#intent_p2 : 0.0 -#intent_p3 : 0.0 -#intent_code : none -#datatype : int16 -#bitpix : 16 -#slice_start : 0 -#pixdim : [-1. 0.9765625 0.9765625 2. 1. 1. -# 1. 1. ] -#vox_offset : 0.0 -#scl_slope : nan -#scl_inter : nan -#slice_end : 0 -#slice_code : unknown -#xyzt_units : 2 -#cal_max : 0.0 -#cal_min : 0.0 -#slice_duration : 0.0 -#toffset : 0.0 -#glmax : 0 -#glmin : 0 -#descrip : b'' -#aux_file : b'' -#qform_code : unknown -#sform_code : aligned -#quatern_b : 0.0 -#quatern_c : 1.0 -#quatern_d : 0.0 -#qoffset_x : 249.51172 -#qoffset_y : -33.01172 -#qoffset_z : -553.5 -#srow_x : [ -0.9765625 0. 0. 249.51172 ] -#srow_y : [ -0. 0.9765625 0. -33.01172 ] -#srow_z : [ 0. -0. 2. -553.5] -#intent_name : b'' -#magic : b'n+1' -``` - -* float32 -```python -nii3ct_float32 = nib.load("data/ct_for_pet_float32.nii.gz") -``` - -```python -print(nii3ct_float32.get_data_dtype()) -#float32 -``` - -```python -print(nii3ct_float32.header) -# object, endian='<' -#sizeof_hdr : 348 -#data_type : b'' -#db_name : b'' -#extents : 0 -#session_error : 0 -#regular : b'' -#dim_info : 0 -#dim : [ 3 512 512 175 1 1 1 1] -#intent_p1 : 0.0 -#intent_p2 : 0.0 -#intent_p3 : 0.0 -#intent_code : none -#datatype : float32 -#bitpix : 32 -#slice_start : 0 -#pixdim : [-1. 0.9765625 0.9765625 2. 1. 1. -# 1. 1. ] -#vox_offset : 0.0 -#scl_slope : nan -#scl_inter : nan -#slice_end : 0 -#slice_code : unknown -#xyzt_units : 2 -#cal_max : 0.0 -#cal_min : 0.0 -#slice_duration : 0.0 -#toffset : 0.0 -#glmax : 0 -#glmin : 0 -#descrip : b'' -#aux_file : b'' -#qform_code : unknown -#sform_code : aligned -#quatern_b : 0.0 -#quatern_c : 1.0 -#quatern_d : 0.0 -#qoffset_x : 249.51172 -#qoffset_y : -33.01172 -#qoffset_z : -553.5 -#srow_x : [ -0.9765625 0. 0. 249.51172 ] -#srow_y : [ -0. 0.9765625 0. -33.01172 ] -#srow_z : [ 0. -0. 2. -553.5] -#intent_name : b'' -#magic : b'n+1' -``` - - -#### ITK-SNAP visualistaion of 16-bit integers vs 32 bit floats - -![Figure. ITK-SNAP: Metadata for 16-bit integers vs 32-bit floats](fig/itk-snap-16-32-bits.svg) - - -### Cropping data -Cropping data might often contain background voxels which can be useful in some instances but might take up valuable RAM memory, especially for large images. -We are going to crop the image from slice 103 to slice 381 in the x dimension (Sagittal slices) and from slice 160 to 340 in the y dimension (Coronal slices). -Use ITK-SNAP to confirm that cropping the image to these slices will only remove slices that do not contain the patient. - -:::::::::::::::::::::::::::::::::::::: spoiler - -### Notes on `Nifti1Image` object -It is not possible to change the image data in a Nifti1Image object, e.g. to use a smaller array corresponding to a cropped image. -Therefore, if you want to modify the image data you need to create a new array with the modified data, then create a new Nifti1Image object using this new array which is used to save the modified image to disk. - -:::::::::::::::::::::::::::::::::::::: - -* Read the image data from disk using `get_fdata` and then create a new array containing a copy of the desired slices: -```python -nii3ct_float32 = nib.load("data/ct_for_pet_float32.nii.gz") -nii3ct_float32.set_data_dtype(np.float32) -print(nii3ct_float32.get_data_dtype()) -#float32 -``` - -```python -image_data=nii3ct_float32.get_fdata() -print(image_data.dtype, image_data.shape) -#float64 (512, 512, 175) -``` - -```python -plt.imshow(image_data[:,:,100], cmap="gray"); plt.show() -#launch plot window -``` - -```python -image_data_cropped = image_data[103:381,160:340,:].copy() -#x[sagittal], y[coronal], z[axial] #voxel units -``` - -:::::::::::::::::::::::::::::::::::::: spoiler - -### Notes on copying data -It is important to create a copy of the data (using the copy function, as show above). -Otherwise the new image will reference the data in the original image, and if you make any changes to the values in the new image they will also be changed in the original image. - -:::::::::::::::::::::::::::::::::::::: - - - -When creating a new Nifti1Image object you can provide an affine transform and/or a Nifti1Header object (e.g. nii_1.header). -When you just provide an affine, a default header is created with the sform set according to the provided affine transform and the qform set to unknown. -When you provide a header, the values in it will be copied into the header for the new image, except that the dim values will be updated based on the image data provided. -If you provide an affine as well as a header, and the affine is different to the sform or qform in the header they will be set to the provided affine and unknown respectively. -If the affine is not provided, it is not set from the header but left empty. -So, if any of your downstream processing will make use of the affine, make sure you manually set it from the header either when creating the Nifti1Image object or by calling the set_sform or set_qform functions for the object (which update the affine as well as setting the sform/qform). -Here we will provide a Nifti1Header object as we want the header for the new image to based on the uncropped image and will use the sform from the header as the affine for the new Nifti1Image. - -```python -nii3ct_float32_cropped = nib.nifti1.Nifti1Image(image_data_cropped, nii3ct_float32.get_sform(), nii3ct_float32.header) -nii3ct_float32_cropped.shape -#(278, 180, 175) -``` - -```python -print(nii3ct_float32_cropped.header) -# object, endian='<' -#sizeof_hdr : 348 -#data_type : b'' -#db_name : b'' -#extents : 0 -#session_error : 0 -#regular : b'' -#dim_info : 0 -#dim : [ 3 278 180 175 1 1 1 1] -#intent_p1 : 0.0 -#intent_p2 : 0.0 -#intent_p3 : 0.0 -#intent_code : none -#datatype : float32 -#bitpix : 32 -#slice_start : 0 -#pixdim : [-1. 0.9765625 0.9765625 2. 1. 1. -# 1. 1. ] -#vox_offset : 0.0 -#scl_slope : nan -#scl_inter : nan -#slice_end : 0 -#slice_code : unknown -#xyzt_units : 2 -#cal_max : 0.0 -#cal_min : 0.0 -#slice_duration : 0.0 -#toffset : 0.0 -#glmax : 0 -#glmin : 0 -#descrip : b'' -#aux_file : b'' -#qform_code : unknown -#sform_code : aligned -#quatern_b : 0.0 -#quatern_c : 1.0 -#quatern_d : 0.0 -#qoffset_x : 249.51172 -#qoffset_y : -33.01172 -#qoffset_z : -553.5 -#srow_x : [ -0.9765625 0. 0. 249.51172 ] -#srow_y : [ -0. 0.9765625 0. -33.01172 ] -#srow_z : [ 0. -0. 2. -553.5] -#intent_name : b'' -#magic : b'n+1' -``` - -The cropped image can now be saved using Nibabel’s save function. - -```python -nib.save(nii3ct_float32_cropped, "data/ct_for_pet_cropped.nii.gz") -``` - -Now load the cropped image into ITK-SNAP. -You will see that the image has been cropped in the x and y dimensions as expected, but it is shifted relative to the uncropped image: - -![Figure. ITK-SNAP: Cropped image in ITK-SNAP](fig/itk-snap_ct_for_pet_cropped_nii_gz.png) - -This is because we did not update the origin/offset in the nifti header to account for the slices we removed. -The origin gives the world coordinates of the voxel (0, 0, 0). -So, if we want the images to remain aligned, we need to set the origin for the cropped image to be the same as the world coordinates of voxel (103, 160, 0) in the uncropped image. -This can be calculated using the affine from the uncropped image, and used to create a new affine matrix with this as the origin (the origin is the top 3 values in the 4th column of the matrix). -The new affine matrix can then be used to set the sform (which also sets the affine) for the cropped image, which can then be saved. - -```python -cropped_origin = nii3ct_float32.affine@np.array([103,160,0,1]) -cropped_origin -#array([ 148.92578125, 123.23828125, -553.5 , 1. ]) -``` - - -```python -aff_mat_cropped = nii3ct_float32_cropped.get_sform() -print(aff_mat_cropped) -#[[ -0.9765625 0. 0. 249.51171875] -# [ -0. 0.9765625 0. -33.01171875] -# [ 0. -0. 2. -553.5 ] -# [ 0. 0. 0. 1. ]] -``` - -```python -aff_mat_cropped[:, 3] = cropped_origin -print(aff_mat_cropped) -#[[ -0.9765625 0. 0. 148.92578125] -# [ -0. 0.9765625 0. 123.23828125] -# [ 0. -0. 2. -553.5 ] -# [ 0. 0. 0. 1. ]] -``` - -Note, NiBabel provides handy functionality for slicing nifti images and updating the affine transforms and header accordingly using the slicer attribute: - -```python -nii3ct_float32_cropped.set_sform(aff_mat_cropped) -nii3ct_float32_cropped.affine -#array([[ -0.9765625 , 0. , 0. , 148.92578125], -# [ -0. , 0.9765625 , 0. , 123.23828125], -# [ 0. , -0. , 2. , -553.5 ], -# [ 0. , 0. , 0. , 1. ]]) -``` - -```python -nib.save(nii3ct_float32_cropped, "data/ct_for_pet_float32_cropped_affine.nii.gz") -``` - -* updating the affine transforms with `nibabel.slicer` -```python -nii3ct_float32_cropped_with_nibabel_slicer = nii3ct_float32.slicer[103:381,160:340,:] -nii3ct_float32_cropped_with_nibabel_slicer.affine -#array([[ -0.9765625 , 0. , 0. , 148.92578125], -# [ 0. , 0.9765625 , 0. , 123.23828125], -# [ 0. , 0. , 2. , -553.5 ], -# [ 0. , 0. , 0. , 1. ]]) -``` - -But you wouldn’t have learnt so much if we’d simply done that to start with :) - -If you load the cropped and aligned image into ITK-SNAP using overlay semi-transparent feature, you will see that it is now aligned with the original image but has fewer slices in the x and y dimensions. - -![Figure. ITK-SNAP: Cropped image in ITK-SNAP](fig/itk-snap-overlay-semi-transparent-cropped-affine.png) - - -If we try to perform a registration between these images as they are, it will likely have difficulties as the images are so far out of alignment to start with. -One way to roughly align the images prior to performing a registration is to align the centres of the images by modifying the origin for the 2nd subject such that the image centres for both images have the same world coordinates. - -This can be done by following these steps: -* Load `ct_for_pet_cropped.nii.gz` -* Calculate the centre of this image in voxel coordinates -* Calculate the centre of this image in world coordinates -* Calculate the centre of the cropped image from subject 1 in voxel coordinates -* Calculate the centre of the cropped image from subject 1 in world coordinates -* Calculate the translation (in world coordinates) required to align the images: Centre image 1 = centre image 2 + translation - -Add this translation to the origin for image 2 -Save the image for image 2 using the header with the updated origin -Try and write code to implement these steps yourself, based on what you have learnt so far. -If you get stuck ask me or one of the PGTAs for help. - -If you have implemented this correctly when you load the aligned image from image 2 into ITK-SNAP is should appear roughly aligned with the images from image 1. -You can use `Color Map` feature tab and the scroll bar to see the differences between two images. - -![Figure. ITK-SNAP main window](fig/itk-snap-calculating-registrations.png) - -### Cropped and aligned image -Automated image registrations can be prone to be failure if there are very large initial differences between the images. A manual alignment is subjective, but can provide a good enough starting guess that keeps the objective, automated registration more reliable. -* Manually aligning images -For manual image registration in ITK-SNAP go to Tools > image registration - -![Figure. ITK-SNAP main window](fig/itk-snap-image-manual-registration.png) - -* Automatic registration (affine) -You can use [niftyreg](https://nipype.readthedocs.io/en/latest/api/generated/nipype.interfaces.niftyreg.reg.html), [Reg_aladin](http://cmictig.cs.ucl.ac.uk/wiki/index.php/Reg_aladin) or [DL methods](https://github.com/Project-MONAI/tutorials/tree/main/3d_registration). - -## References -* NiBabel: "The NiBabel documentation also contains some useful information and tutorials for working with NifTi images and understanding world coordinate systems and radiological vs neurological view." - * https://nipy.org/nibabel/nifti_images.html - * https://nipy.org/nibabel/coordinate_systems.html - * https://nipy.org/nibabel/neuro_radio_conventions.html -* ITK-SNAP: - * [Tutorial: Getting Started with ITK-SnAP](http://www.itksnap.org/docs/viewtutorial.php) - * [ITK-SNAP 3.x Training Class Final Program](http://itksnap.org/files/handout_201409.pdf) - -* AI-based medical image registration: - * https://github.com/Project-MONAI/tutorials/tree/main/2d_registration - * https://github.com/Project-MONAI/tutorials/tree/main/3d_registration - * https://github.com/Project-MONAI/tutorials/tree/main/3d_regression - diff --git a/episodes/practical1.Rmd b/episodes/practical1.Rmd new file mode 100644 index 00000000..d97770e0 --- /dev/null +++ b/episodes/practical1.Rmd @@ -0,0 +1,617 @@ +--- +title: "An introduction to medical imaging data" +teaching: 10 +exercises: 2 +output: pdf_document +--- + +:::::::::::::::::::::::::::::::::::::: questions + +- How is medical imaging data typically stored? +- What tools and methods can be used to visualise medical imaging data? +- How can medical imaging data be manipulated and modified consistently? + +:::::::::::::::::::::::::::::::::::::::::::::::: + +::::::::::::::::::::::::::::::::::::: objectives + +- Explore medical imaging data formats, focusing on DICOM and NIfTI. +- Learn to use ITK-SNAP software to display medical imaging data. +- Understand the structure and content of the NIfTI header. +- Modify and manipulate NIfTI images and headers using Python. + +:::::::::::::::::::::::::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::: prereq + +You should have followed the instructions from [the Setup page](../learners/setup.md) to download the course material, setup Python, and install ITK-SNAP. + +:::::::::::::::::::::::::::::::::::::: + +The first part of this tutorial makes use of the ITK-SNAP software. The second part is based in Python, and uses the Jupyter notebook `practical1-exercises.ipynb` downloaded with the course material. + +### Data + +In these exercises, we will be working with real-world medical imaging data from The Cancer Imaging Archive [(TCIA)](https://www.cancerimagingarchive.net/collection/ct-vs-pet-ventilation-imaging/). +TCIA is a resource of public datasets that hosts a large archive of medical images of cancer accessible for public download. +We are going to use data from the CT Ventilation as a Functional Imaging Modality for Lung Cancer Radiotherapy dataset (also known as CT-vs-PET-Ventilation-Imaging dataset). +It contains 20 lung cancer patients who underwent exhale/inhale breath-hold CT (BHCT), free-breathing four-dimensional CT (4DCT) and Galligas PET ventilation scans in a single session on a combined 4DPET/CT scanner. We won't make use of the 4DCT. + +We used the scans from the patient CT-PET-VI-02. In the unzipped `data` folder, you can find them in the `practical1` folder. You should see four folders: + +* `inhale_BH_CT` and `exhale_BH_CT`: CT scans acquired during an inhalation breath hold and exhalation breath hold, respectively. +* `PET`: PET scan measuring local lung function. +* `CT_for_PET`: CT scan acquired at the same time as the PET scan for attenuation correction of the PET scan and to provide anatomical reference for the PET data. + + +:::::::::::::::::::::::::::::::::::::: spoiler + +#### Ethical approval for using open datasets +When using open datasets that contain images of humans, such as those from TCIA, for your own research, you will still need to get ethical approval as you do for any non-open datasets you use. You should contact your local Research Ethics Committee at your institution for details on how to obtain ethical approval for your research. + +Ethical approval has been obtained from UCL for using these datasets as part of this teaching course. + +:::::::::::::::::::::::::::::::::::::: + +### Medical imaging formats + +Although there are many different file formats used in medical imaging, we will focus on two of the most common formats: + +* [DICOM](https://www.dicomstandard.org/) (Digital Imaging and Communications in Medicine) +* [NIfTI](https://nifti.nimh.nih.gov/) (Neuroimaging Informatics Technology Initiative) + +Most data from a hospital will be in DICOM format, whereas NIfTI is a very popular format in the medical image analysis community. The files from the TCIA are in DICOM format. + +:::::::::::::::::::::::::::::::::::::: callout + +#### DICOM format +DICOM is a widely used technical standard for digitally storing and transmitting medical images and related information. +DICOM images typically have a separate file for every slice. +If you look in the `CT_for_PET` folder, you will see that there are 175 individual files corresponding to 175 slices in the volume. More modern DICOM images can come with all slices in a single file. + +In addition to the image data for each slice, each file contains a header which can include extensive information relating to the scan and subject. +In a clinical setting, this will include identifiable information such as the patient's name, address, and other relevant details. +Such information should be removed before the data is transferred from the clinical network for use in research. +If you discover patient-identifiable information in the header of the data you are using, you should immediately alert your supervisor, manager or collaborator. + +Most of the information in the DICOM header is not directly useful for typical image processing and analysis tasks. +Furthermore, there are complicated ‘links’ (provided by unique identifiers, UIDs) between the DICOM headers of different files belonging to the same scan or subject. With DICOM routinely storing each slice as a separate file, it makes processing an entire imaging volume stored in DICOM format rather cumbersome, and the extra housekeeping required could lead to a greater chance of an error being made. +Therefore, a common first step of any image processing pipeline is to convert the DICOM image to a more suitable format, such as NIfTI. +Generally, most conversions go from DICOM to NIfTI. There are scenarios when you might want to convert from NIfTI back to DICOM, for example, if you need to import them into a clinical system that only works with DICOM. + + +:::::::::::::::::::::::::::::::::::::::::: spoiler + +#### Warning on converting images back to DICOM +Converting images back to DICOM such that they are correctly interpreted by a clinical system can be very tricky and requires a good understanding of the DICOM standard. +More information on the DICOM standard can be found here: https://www.dicomstandard.org. + +:::::::::::::::::::::::::::::::::::::::::: + +#### NIfTI format +The NIfTI image format usually stores all of the imaging data and header information within a single file. +While the neuroimaging community originally developed the NIfTI file format, it is not specific to neuroimaging. It is now widely used for many medical imaging applications outside the brain. +The NIfTI header contains much less information than DICOM, but it includes all the key information required to interpret, manipulate, and process the image. +During these exercises, you will learn about some of the key information stored in the NIfTI header. +For more information, please see Anderson Winkler's blog on the [NIfTI file format]( https://brainder.org/2012/09/23/the-nifti-file-format/). + +:::::::::::::::::::::::::::::::::::::: + + +## 1. Visualising images with ITK-SNAP + +### 1.1. Getting started with ITK-SNAP + +:::::::::::::::::::::::::::::::::::::: callout + +#### ITK-SNAP + +ITK-SNAP is a free, open-source, multi-platform software application used to segment structures in 3D and 4D biomedical images. Although ITK-SNAP has been developed to segment structures, it can also be used to view 3D and 4D medical image data in a wide variety of formats. In this practical we will just be using ITK-SNAP to view 3D images and will not be exploring any of it's segmentation functionality. For more information on ITK-SNAP, including documentation and tutorial videos, see the [ITK-SNAP website](http://www.itksnap.org/pmwiki/pmwiki.php?n=Main.HomePage). + +:::::::::::::::::::::::::::::::::::::: + +* Open the ITK-SNAP application. + * When it opens it may display a 'Layout Preference Reminder' popout window: + + ![](fig/itk-snap-layout-pref-reminder.png) + + * If this popout appears and you have not changed the settings in ITK-SNAP it should say the same as above, i.e. that radiological convention is used. In this case you can check 'Don't remind me again' and click 'No'. + +* When ITK-SNAP opens you can select between the 'Getting Started' pane, the 'Recent Images' pane, or the 'Recent Workspaces' pane using the tabs near the top of the window. You should select the 'Recent Images' pane if it is not already showing. + * If you have used ITK-SNAP before it will display recent images which you can select to open. Note, for DICOM images the name displayed is the name of the file containing the first DICOM slice. This can be confusing when the folder names rather than the file names specify which scan is which, as is the case in this practical. + +![Figure. ITK-SNAP Recent Images](fig/itk-snap-recent-images.png) + +### 1.2. Opening and viewing DICOM images + +* Load the inhale scan. + * If the image you want to open is not listed under recent images (or you are not sure which one it is), you can open it by clicking the 'Open Image...' button at the bottom right of the window or by clicking on 'File -> Open Main Image...' in the top ribbon. When you do this, a small window opens where you can enter the path of the file you want to open or click 'Browse...' to search your file explorer. + * Select the first slice from the `inhale_BH_CT` folder, which is called `1-001.dcm`. Make sure the 'File Format:' is set to 'DICOM Image Series' and click 'Next >'. + * :warning: **Warning!** some versions of ITK-SNAP incorrectly set the 'File Format:' for the images used in this practical to '4D CTA DICOM Series' + * If this happens use the drop-down menu to set the 'File Format:' to 'DICOM Image Series' instead before clicking 'Next >' + + ![](fig/itk-snap-open-image-1.png) + + * You will then be asked to select the DICOM series to open, but as this folder only includes one series, you can click 'Next >'. + + ![](fig/itk-snap-open-image-2.png) + + * Once the image has loaded a small window will display some summary information about the image. This can be reviewed and then the window can be closed by clicking 'Finish'. + + ![](fig/itk-snap-open-image-3.png) + +* The inhale scan will now be displayed in ITK-SNAP, as shown below. On the left you will see the 'Main Toolbar', 'Cursor Inspector', 'Segmentation Labels', and '3D Toolbar'. In the rest of the window you will see four panes, three of which show orthogonal slices through the 3D volume: + * Top-left shows an axial slice. + * Top-right shows a sagittal slice. + * Bottom-right shows a coronal slice. + * Bottom-left is currently empty (but would show a 3D reconstruction of the structures that have been segmented in the image when using ITK-SNAP for segmentation) + +![Figure. ITK-SNAP after loading the inhale scan](fig/itk-snap-inhale-scan.png) + +### 1.3. Navigating around 3D images +* To navigate around the image, you can left-click in any of the three orthogonal slices to reposition the cursor to where you clicked. This will also change the slices displayed in the other two panes so that they show the slices that pass through the cursor (as indicated by the two dotted blue lines). + * If you left-click and drag the cursor the other panes will scroll through the slices as you move the mouse. + * If you use the mouse wheel (or two fingers on the touchpad) you will scroll through the slices in the pane where the mouse pointer currently is. You can also use the slider next to the displayed slice. + * As you move the cursor around or scroll through the slices you should notice that the information displayed in the 'Cursor Inspector' is updated accordingly. Note, the cursor position displayed here is in voxel coordinates. + + ![](fig/itk-snap-cursor-inspector.png) + +* To zoom in , you can use the right click and drag the cursor upwards (to zoom out, drag the cursor downwards). +* You can pan, or move the entire image, by clicking and dragging the mouse centre button (or alt + left click on Windows, option key + left click on Mac OS). + +### 1.4. Loading and displaying multiple images +If you open a new main image in ITK-SNAP it will replace the current image, making it difficult to directly compare the images. However, it is possible to load multiple images into ITK-SNAP so that they can be easily compared. The additional images can either be loaded as separate images, enabling you to quickly switch the displayed image from one to another, or as semi-transparent overlay images, enabling you to see one image overlaid on the other. + +* Load an additional image as a separate image + * Click 'File -> Add Another Image...' + * This will open a similar window to the one for opening the main image, where you can enter the name of the file you want to open or use the 'Browse...' button to locate it. + * This time select the first slice from the `exhale_BH_CT` folder, which is also called `1-001.dcm`. As before, make sure the 'File Format:' is set to 'DICOM Image Series' and click 'Next >'. + + ![](fig/itk-snap-additional-image-1.png) + + * As before you will be asked to select the DICOM series to open, but there is only one series in this folder, so just click 'Next >'. + * A window will then be displayed asking 'How should the image be displayed?' + * Select 'As a separate image (shown beside other images)' and then click 'Next >'. + + ![](fig/itk-snap-additional-image-2.png) + + * Once the second image is loaded, you will see two thumbnails in the top corner of each slice display. + * Click on these to swap the displayed slices from one image to the other, allowing you to easily see the similarities and differences between the images. + * You can also swap quickly by using the square bracket shortcut keys '[' and ']'. + * You will also see that both images now appear in the 'Cursor Inspector', where you can see the intensity value at the cursor for both images. + * You can also change the displayed image by clicking on it in the 'Cursor Inspector'. + * If you compare the inhale and exhale image you will see they are very similar, although the patient's chest is further out in the inhale scan, and the rest of the anatomy has shifted accordingly. + + ![Figure. ITK-SNAP after loading the inhale and exhale scans](fig/itk-snap-inhale-and-exhale-scans.png) + +* Load an additional image as a colour overlay + * To demonstrate the use of colour overlays we are going to load the PET scan as a colour overlay on top of the corresponding CT scan. You should first load the CT scan as the main image (which will close any other images that are currently open, e.g. the inhale and exhale scans). + * Click 'File -> Open Main Image...'. + * Select the first slice from the `CT_for_PET` folder (also called `1-001.dcm`) and again make sure the 'File Format:' is set to 'DICOM Image Series' before clicking 'Next >'. + + ![](fig/itk-snap-additional-image-3.png) + + * As before you will be asked to select the DICOM series to open, but there is only one series in this folder, so just click 'Next >'. + * And then the summary information about the image will be displayed, which can be reviewed before closing the window by clicking 'Finish'. + * Now click 'File -> Add Another Image...' + * Select the first slice from the `PET` folder (also called `1-001.dcm`) and check the 'File Format:' is set to 'DICOM Image Series' before clicking 'Next >'. + + ![](fig/itk-snap-additional-image-4.png) + + * Once again there is only one DICOM series in the folder so just click 'Next >'. + * The window will now be displayed asking 'How should the image be displayed?' + * This time select 'As a semi-transparent overlay (shown op top of other images)'. + * Use the dropdown menu to set the 'Overlay color map:' to 'Hot'. + * Click 'Next >'. + + ![](fig/itk-snap-additional-image-5.png) + + * The 'Image Summary' will now be displayed. This time it will include a warning about loss of precision, but this can be ignored. Click 'Finish' to close the window. + + ![](fig/itk-snap-additional-image-6.png) + + * The PET scan will now be displayed in ITK-SNAP using the hot colourmap (so the lungs mostly appear in red). Depending on the version of ITK-SNAP used, the opacity of the PET scan may initially be set to 0% (so all you will see is the CT scan), 100% (so all you will see is the PET scan), or some other value (so you will see both the CT scan and the PET scan overlaid on it). + * The opacity of the PET scan can be adjusted by right-clicking on the PET scan in the 'Cursor Inspector' (where it is called *Galligas Lung*) and dragging the 'Opacity:' slider. + + ![Figure. ITK-SNAP after loading the PET scan as a colour overlay on top of the corresponding CT scan](fig/itk-snap-ct-pet-overlay.png) + +### 1.5. Image Layer Inspector tool +The 'Image Layer Inspector' tool displays further information about the images and allows you to modify and control how they are displayed. We will now briefly look at some of the functionality provided by the 'Image Layer Inspector'. + +* Click 'Tools -> Layer Inspector...' in the top ribbon. This will open the 'Image Layer Inspector' in a separate window. + +You will see the open images listed on the left, and five tabs across the top: 'General', 'Contrast', 'Color Map', 'Info', and 'Metadata'. You can click on the images listed on the left to select a different image. + +#### General tab +The 'General' tab displays the filename and the nickname for the selected image. For any 'Additional Images' the 'General' tab also enables you to select if the image is displayed as a 'Separate image' or a 'Semi-transparent overlay', and to change the 'Overlay opacity' if the image is being displayed as an overlay. + +![](fig/itk-snap-image-layer-inspector-general.png) + +#### Contrast tab +The 'Contrast' tab enables you to adjust how the image intensities are displayed for the different images. You can use the text boxes in the 'Linear Image Contrast Adjustment:' panel to modify the range of intensity values displayed, e.g.: + +* Select the CT image (called 'CT Lung 3.0 B31f' and listed as the 'Main Image') from the list on the left. +* Use the text box to set the 'Maximum:' value to '1000' (push enter after typing the number). The CT image should now appear brighter than it previously did. + +The 'Reset' button can be used to return the intensity range to the original values. The 'Auto' button can be used to automatically set the intensity range based on the distribution of the intensity values in the image. +The 'Curve-Based Image Contrast Adjustment:' panel can also be used to provide more fine grained control of how the image intensities are mapped to displayed intensities. The details of how this works are beyond the scope of this practical, but feel free to have a play with it in your own time (and remember the 'Reset' button resets the curve as well). + +![](fig/itk-snap-image-layer-inspector-contrast.png) + +#### Color Map tab +The 'Color Map' tab can be used to select and manipulate the colour map used for each image. The colour map can be changed for any of the images, not just those displayed as semi-transparent overlays. You can select from predefined colour maps using the dropdown menu in the 'Presets:' panel. +The 'Color Map Editor:' panel can be used to modify the colour maps, which can then be saved as new presets (using the '+' button in the 'Presets:' panel). The details of how this works are beyond the scope of this practical, but feel free to have a play with it in your own time. You can always restore the original colour map by reselecting it from the dropdown menu in the 'Presets:' panel. + +![](fig/itk-snap-image-layer-inspector-color-map.png) + +#### Info tab +The 'Info' tab displays some of the important information about the selected image under 'Image Header'. It also displays the 'Cursor Coordinates' in both voxel and world units (ITK-SNAP uses the same world coordinate system as NIfTI images - see section 4.3 below for more information on NIfTI and DICOM world coordinate systems) as well as the intensity at the cursor for the selected image. + +* First select the CT image. Then select the PET image. Notice how the values change for both the 'Image Header' and 'Cursor Coordinates'. + * The resolution ('Spacing') of the images is different, with the PET image having coarser voxels (larger spacing) than the CT. + * The voxel coordinates are always integers for the main image (here the CT image) but may be non-integer for any additional images. + * The world coordinates are the same for all images. + +In ITK-SNAP the cursor can be located at the centre of any voxel in the main image. However, the voxel locations in the additional images may not be aligned with the voxels in the main image, e.g. if they have a different resolution, as for the CT and PET images here. In this case the cursor will not be located at the centre of a voxel in the additional images, hence the non-integer voxel coorindates. + +![](fig/itk-snap-image-layer-inspector-info.png) + +#### Metadata tab +The 'Metadata' tab contains information stored in the headers of the image files, so in this case (some of) the DICOM header information for these images. As can be seen, a lot of the information stored in the DICOM headers is not required or useful for processing the images. + +![](fig/itk-snap-image-layer-inspector-metadata.png) + +### 1.6 Converting DICOM images to NIfTI +As mentioned earlier, converting DICOM images to another format, such as NIfTI, is often one of the first steps when processing and analysing medical image data. ITK-SNAP can save images in NIfTI format (and a few others), so can be used to convert the DICOM images to NIfTI. To save the images in NIfTI format: + +* Open the 'General' tab in the 'Image Layer Inspector'. +* Select the image you want to save from the list on the left, in this case select the CT image. Then click on the disk icon just above the list of images on the left. This will open the 'Save Image' window. + * Click 'Browse...' and navigate up one folder to the `practical1` folder and enter the 'File name:' as `CT_for_PET.nii.gz`. + * Note, using the `.nii.gz` extension saves the image as a compressed NIfTI file. Many softwares and libraries, including ITK-SNAP, NiBabel, and NiftyReg, can work directly with the compressed `.nii.gz` files, so these are often used to save storage space, and will be used in these exercises. + * Check the 'File Format:' is set to 'NiFTI' and click 'Finish'. + + ![](fig/itk-snap-save-image-1.png) + +* If you look in the `practical1` folder you should now see a file called `CT_for_PET.nii.gz`. +* If you load this file as an additional image in ITK-SNAP you should see that the new NIfTI image looks exactly the same as the original DICOM image, and the 'Info' tab in the 'Image Layer Inspector' has exactly the same information for both images, as we would expect. +* However, the information in the 'Metadata' tab for the two images is different, with the information from the NIfTI file header being displayed for the NIfTI image and the DICOM header information being displayed for the DICOM image. + +![](fig/itk-snap-nifti-header.png) + +:::::::::::::::::::::::::::::::::::::: spoiler + +### Converting large numbers of DICOM images to NIfTI +While ITK-SNAP provides a convenient way to convert DICOM images to NIfTI format, it is not very practical to use if you need to convert more than a few images. There are many free tools available that can be used to convert large numbers of images from DICOM to NIfTI, such as: + +* [dicom2nifti](https://dicom2nifti.readthedocs.io/en/latest/) +* [dcm2niix](https://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage) +* [Scikit-rt](https://scikit-rt.github.io/scikit-rt/index.html) + +:::::::::::::::::::::::::::::::::::::: + +## 2. Viewing and understanding the NIfTI header with NiBabel + +:::::::::::::::::::::::::::::::::::::: callout + +#### NiBabel + +We will use the Python package [NiBabel](https://nipy.org/nibabel/) for viewing and manipulating NIfTI images and their headers. +We will demonstrate some of the functionality of NiBabel, but recommend checking the [documentation](https://nipy.org/nibabel/#documentation) for further details on using the NiBabel library. This also contains some useful information and tutorials on [working with NIfTI images](https://nipy.org/nibabel/nifti_images.html), [understanding world coordinate systems](https://nipy.org/nibabel/coordinate_systems.html), and [neurological vs radiological views](https://nipy.org/nibabel/neuro_radio_conventions.html). + +:::::::::::::::::::::::::::::::::::::: + +### 2.1. Reading and displaying the NIfTI header +The NIfTI header can be read using NiBabel’s `load` function. + +* Run the first two cells from `practical1_exercises.ipynb` to import the required libraries and `load` the NIfTI image saved in the previous section, `CT_for_PET.nii.gz`. + +This creates a Nifti1Image object. This has a number of useful attributes and functions for accessing and manipulating the NIfTI header and image data. E.g.: + +* Run the next cell to display size/shape of the image, the data types used to store the image data, and the affine transform that maps from voxel space to world space. +```output +(512, 512, 175) +int16 +[[ -0.9765625 0. 0. 249.51171875] + [ 0. -0.9765625 0. 466.01171875] + [ 0. 0. 2. -553.5 ] + [ 0. 0. 0. 1. ]] + ``` + +:::::::::::::::::::::::::::::::::::::: callout + +#### Note + +The `load` function does not read the image data itself from disk, it just reads and interprets the header, but the Nifti1Image object it returns provides functions for accessing the image data, e.g. `get_fdata` (as we'll see in the next section). + +:::::::::::::::::::::::::::::::::::::: + + +* You can also display the full header information, as seen by running the next cell. +```output + object, endian='<' +sizeof_hdr : 348 +data_type : np.bytes_(b'') +db_name : np.bytes_(b'') +extents : 0 +session_error : 0 +regular : np.bytes_(b'r') +dim_info : 0 +dim : [ 3 512 512 175 1 1 1 1] +intent_p1 : 0.0 +intent_p2 : 0.0 +intent_p3 : 0.0 +intent_code : none +datatype : int16 +bitpix : 16 +slice_start : 0 +pixdim : [1. 0.9765625 0.9765625 2. 0. 0. 0. + 0. ] +vox_offset : 0.0 +scl_slope : nan +scl_inter : nan +slice_end : 0 +slice_code : unknown +xyzt_units : 10 +cal_max : 0.0 +cal_min : 0.0 +slice_duration : 0.0 +toffset : 0.0 +glmax : 0 +glmin : 0 +descrip : np.bytes_(b'') +aux_file : np.bytes_(b'') +qform_code : scanner +sform_code : scanner +quatern_b : 0.0 +quatern_c : 0.0 +quatern_d : 1.0 +qoffset_x : 249.51172 +qoffset_y : 466.01172 +qoffset_z : -553.5 +srow_x : [ -0.9765625 0. 0. 249.51172 ] +srow_y : [ 0. -0.9765625 0. 466.01172 ] +srow_z : [ 0. 0. 2. -553.5] +intent_name : np.bytes_(b'') +magic : np.bytes_(b'n+1') +``` + +### 2.2 Specifying the affine transform + +There are two common ways of specifying the affine transformation mapping from voxel coordinates to world coordinates in the NIfTI header, called the *sform* and the *qform*. +The sform directly stores the affine matrix in the header, whereas the qform stores the affine transformation using [quaternions](https://en.wikipedia.org/wiki/Quaternion). + +You can see that the sform is stored in the `srow_x`, `srow_y`, and `srow_z` fields in the header. These specify the top three rows of a 4x4 matrix representing the affine transformation using [homogeneous coordinates](https://en.wikipedia.org/wiki/Homogeneous_coordinates). The fourth row is not stored as it is always `[0 0 0 1]`. You can see that the sform matrix matches the affine transformation from the Nifti1Image object. + +The qform is stored in the `quartern_b`, `quartern_c`, and `quartern_d` fields, which represent 3D rotations, and in the `qoffset_x`, `qoffset_y`, and `qoffset_z` fields, which give the offset between the first voxel in the image, i.e. voxel coordinates (0,0,0), and the origin in world coordinate space, i.e. world coordinates (0,0,0). + +The `sform_code` and `qform_code` fields indicate whether the sform and/or qform have been used to specify the affine transformation in the NIfTI header. There are five different sform/qform codes defined: + +Table: Sform/qform codes. + +| Code | Label | Meaning | +| ---- | --------- | ------------------------------- | +| 0 | unknown | not defined | +| 1 | scanner | RAS+ in scanner coordinates | +| 2 | aligned | RAS+ aligned to some other scan | +| 3 | talairach | RAS+ in Talairach atlas space | +| 4 | mni | RAS+ in MNI atlas space | + +The different codes are meant to indicate which space the affine transform aligns the image to, but they are not always used/set correctly, and for most purposes all that matters is whether the code is 0 (unknown), in which case the sform/qform should not be used to specify the affine transformation, or the code is set to any of the other values, in which case the sform/qform should be used to specify the affine transformation. + +You will see that when ITK-SNAP saved the NIfTI image it set both the sform and qform code to 1 (scanner). It also set the values of the sform and qform fields in the header so that they represent the same affine transformation. You can see this using the `get_qform` function, which returns the affine matrix specified by the qform. + +* Run the next cell of the Juypter notebook to display the affine matrix represented by the qform. +```output +[[ -0.9765625 0. 0. 249.51171875] + [ 0. -0.9765625 0. 466.01171875] + [ 0. 0. 2. -553.5 ] + [ 0. 0. 0. 1. ]] + ``` + +:warning: **Warning**, the NIfTI format does not require that the sform and qform specify the same transformation, and in general it is not well defined what should be done when they both have a non-zero code and specify different transformations. Therefore, **we recommend only using one of the sform or qform to specify the affine transformation**, and setting the code for the other to 0 (as we will do later). And **we recommend using the sform rather than the qform** for various reasons: + +* It is easier to interpret as it directly provides the affine matrix. +* The qform cannot represent shears but the sform can represent any arbitrary affine transformation. +* NiBabel (and many other libraries and software tools) gives precedence to the sform and if both the sform and qform are provided it will ignore the qform. + +:::::::::::::::::::::::::::::::::::::: callout + +#### Using the sform to align images + +As the sform can store any arbitrary affine transformation it can be used to store the result of an affine registration between two images, effectively aligning the images without needing to resample either of them. + +:::::::::::::::::::::::::::::::::::::: + +### 2.3. World coordinate systems assumed by the DICOM and NIfTI formats +#### NIfTI - RAS +In the above affine matrix, you will notice that the diagonal elements of the matrix match the voxel dimensions (stored in the `pix_dim` field in the NIfTI header), but the values for the X and Y dimensions are negative. This is because the voxel indices increase as you move to the left and posterior of the image/patient, but the NIfTI format assumes a RAS world coordinate system (i.e. the values increase as you move to the **Right**, **Anterior**, and **Superior** of the patient). Therefore, the world coordinates decrease as the voxel indices increase. + +#### DICOM - LPS +On the other hand, the original DICOM images assume the world coordinate system is LPS (i.e. values increase when moving to the **Left**, **Posterior**, and **Superior**). So if you calculated the affine transform from the DICOM headers (see [here](https://nipy.org/nibabel/dicom/dicom_orientation.html) if you want to know how to do this) you would obtain the following affine matrix: +```output +[[ 0.9765625 0. 0. -249.51171875] + [ 0. 0.9765625 0. -466.01171875] + [ 0. 0. 3. -205.5 ] + [ 0. 0. 0. 1. ]] +``` +As you can see, the values in the first two rows corresponding to the X and Y dimensions are the negative of those in the NIfTI header. + +It should be noted that neither format assumes the voxels are stored in a particular order/orientation on disk. It is quite common that images originally saved as DICOMs are written to disk with the voxels also stored in LPS orientation, as is the case here, but there are many times when this is not the case. + +:::::::::::::::::::::::::::::::::::::: callout + +#### What about ITK-SNAP? + +ITK-SNAP uses a RAS coordinate system, like NIfTI images. However, it knows that DICOM images use a LPS coordinate system, so when a DICOM image is loaded into ITK-SNAP it calculates the affine transformation from the DICOM header, and then takes the negative of the X and Y rows, so that the image is displayed correctly. + +* Look at the 'Info' tab in the 'Image Layer Inspector' for the DICOM CT image while you drag the cursor around the image. You will see that the coordinates in 'Voxel units' increase as you move the cursor to the left and posterior of the image, but the coordinates in 'World units (NIfTI)' increase as you move to the right and anterior of the image. Both the voxel and world coordinates increase as you move superior in the images - at least DICOM and NIfTI can agree on up and down! :smile: + +However, it gets even more confusing! :weary: + +While ITK-SNAP uses a RAS coordinate system, it calls it a LPI coordinate system as it names coordinate systems according to where they start, not the direction in which they increase. Naming coordinate systems by the direction they increase is more common, but it is not universal, as we have seen with ITK-SNAP. Therefore, sometimes it will be written RAS+ (as in the table above), to be clear that this means the axis increase in the specified direction. + +You will see that the 'Info' tab lists the 'Orientation:' of the images used in these exercises as RAI. This corresponds to the orientation of the voxels (as determined from the affine transformation), and means that the voxel coordinates start on the right, anterior, and inferior of the patient, and increase as you move to the left, posterior, and superior of the patient (i.e. LPS+). + +:::::::::::::::::::::::::::::::::::::: + +## 3. Modifying NIfTI images and headers with NiBabel + +### 3.1. Changing the data type and qform code +The image data for a Nifti1Image object can be accessed using the `get_fdata` function. +This will load the data from the disk and cast it to float64 type before returning it as a numpy array. + +* Run the next cell in the Juypter notebook to read the image data and check it is returned as a numpy array. + +Casting it to float64 prevents any integer-related errors and problems in downstream processing when using the data read by NiBabel. +However, this does not change the type of the image stored on disk, which we have previously seen is int16 for this image. We are going to save a new image to disk where the image is stored as 32-bit floating point numbers. We are also going to set the qform code to 0, as suggested above. + +* Run the next cell to change the data type and qform code of the Nifti1Image object using its `set_data_dtype` and `set_qform` functions. This also updates the corresponding values in the header. The image is then written to disk using the `save` function. + +We do not need to manually modify the type of the image data. Setting the data type for the Nifti1Image object tells it which type to use when writing the image data to disk. + +The new file `CT_for_PET_float32.nii.gz` is larger than `CT_for_PET.nii.gz` (62 MB and 46 MB, respectively) as the image is stored as 32-bit floats instead of 16-bit integers. However, it is not twice the size, as may be expected (a 32-bit floating point number is twice the size of a 16-bit integer). This is because both images are compressed, but the compression is more efficient for the floating point image. + +#### Check data types with NiBabel +* Run the next cell in the Juypter notebook to load `CT_for_PET_float32.nii.gz` and check that the data type is `float32`, and to reload `CT_for_PET.nii.gz` and check that the data type is still `int16`. + +#### Check the data types with ITK-SNAP +* Load `CT_for_PET_float32.nii.gz` as an additional image in ITK-SNAP. You should see that the new NIfTI image looks exactly the same as the previous NIfTI image (and the original DICOM image) and has the same intensity values. +* However, if you look at the 'Info' tab in the 'Image Layer Inspector' you will see that the 'Pixel Format:' for the new image is float32, but is int16 for the previous image. + +![](fig/itk-snap-info-32.png) + +![](fig/itk-snap-info-16.png) + +* And if you look at he 'Metadata' tab you will see that the 'bitpix' value is 32 for the new image and 16 for the previous image, and the 'datatype' is 16 for the new image and 4 for the previous image. + +![](fig/itk-snap-metadata-32.png) +![](fig/itk-snap-metadata-16.png) + + +### 3.2. Cropping images +We are going to crop the CT image to remove slices that do not contain the region of interest, in this case is the lungs. This can make downstream processing of the images more efficient. However, there are a few considerations to pay careful attention to when cropping to ensure it is performed correctly and to preserve the mapping of the image to world coordinates. + +You usually want to leave 5-10 slices on each side of the region of interest when cropping. Therefore, we are going to crop the image so that we keep slices 91 to 390 in the x dimension (sagittal slices), slices 131 to 375 in the y dimension (coronal slices), and slices 21 to 155 in the z dimension (axial slices). + +* Use ITK-SNAP to confirm that cropping the image to these slices will only remove slices that do not contain the lungs. + +Changing the image data of a Nifti1Image object is impossible, e.g., to set it to a smaller array corresponding to a cropped image. +Therefore, if you want to modify the image data, you need to create a new array containing the modified data and then create a new Nifti1Image object using this new array, which is used to save the modified image to disk. + +* Run the next cell of the Juypter notebook to create a new array containing a copy of the desired slices. + +It is important to create a copy of the data, otherwise, the new image will reference the data in the original image, and if you make any changes to the values in the new image, they will also be changed in the original image. + +When creating a new Nifti1Image object, you can provide an affine transform and/or a Nifti1Header object (e.g. `ct_for_pet_nii.header`) as well as the image data. The affine transform and the header for the new Nifti1Image object will be based on those provided, but may be updated so that they are consistent with the image data and each other, e.g. the `dim` field in the new header will be updated based on the image data provided. + +* Run the next cell to create a new Nifti1Image object using the header and affine from the uncropped image, and check the shape and header are as expected. It then saves the cropped image using NiBabel’s save function. +* Load the uncropped CT image (`CT_for_PET.nii.gz`) into ITK-SNAP as the main image. Now load the cropped image as an additional image. +You will see that the cropped image has been cropped in all three dimensions so that it just contains the lungs, but it is shifted relative to the uncropped image. + +![Figure. ITK-SNAP after loading the cropped image (and repositioning the cursor to see the image more clearly)](fig/itk-snap-cropped-ct.png) + +:::::::::::::::::::::::::::::::::::::: callout + +#### Why are voxels outside the image grey? + +ITK-SNAP sets voxels outside the image to have a value of 0. For CT images in Hounsfield Units (HU) this corresponds to soft tissue, so the voxels outside the image have a similar intensity to soft tissue. Unfortunately, it does not appear to be possible to change the value that ITK-SNAP assigns to the voxels outside the image. + +:::::::::::::::::::::::::::::::::::::: + +The cropped image is shifted relative to the original image because we did not update the image origin in the NIfTI header to account for the slices we removed. +The image origin gives the world coordinates of the voxel (0, 0, 0). +As we did not update the image origin, voxel (0,0,0) in the cropped image is aligned with voxel (0,0,0) in the uncropped image. However, it should be aligned with voxel (91, 131, 21) in the uncropped image. So, we need to set the image origin for the cropped image to be the same as the world coordinates of the voxel (91, 131, 21) in the uncropped image. + +* Run the next cell in the Juypter notebook. This does the following: + * Uses the affine from the uncropped image to calculate the cropped image origin, i.e. the world coordinates of voxel (91, 131, 21). + * Uses the cropped image origin to update the corresponding values affine matrix for the cropped image. + * Saves the cropped image. + * This updates the sform in the header of the cropped image to be the same as the affine before saving the image. + +:::::::::::::::::::::::::::::::::::::: callout + +#### NiBabel's `slicer` attribute + +NiBabel provides handy functionality for cropping NIfTI images and updating the affine transforma and header accordingly using the `slicer` attribute. Using the `slicer` attribute you can do the same thing we did above with a single line of code: + +```python +ct_for_pet_slicer_nii = ct_for_pet_nii.slicer[x_first:x_last+1, y_first:y_last+1, z_first:z_last+1] +``` + +But you wouldn’t have learnt so much if we’d simply done that to start with! :smile: + +:::::::::::::::::::::::::::::::::::::: + +If you load the new image into ITK-SNAP you will see that it is now aligned with the original image. + +![Figure. ITK-SNAP after loading the cropped and aligned image)](fig/itk-snap-cropped-ct-aligned-1.png) + +![Figure. ITK-SNAP showing the same slices through the uncropped image)](fig/itk-snap-cropped-ct-aligned-2.png) + +### 3.3. Aligning images +If you load `CT_for_PET.nii.gz` and the DICOM image in the inhale_BH_CT` folder in to ITK-SNAP, you can see that there is a large misalignment between the images. + +![Figure. ITK-SNAP after loading `CT_for_PET.nii.gz` and the DICOM image in the inhale_BH_CT` folder](fig/itk-snap-misaligned-images.png) + +If we perform a registration between these images as they are, it may have difficulties as the images are so far out of alignment to start with. One way to roughly align the images before performing a registration is to align the centres of the images. + +:::::::::::::::::::::::::::::::::::::: challenge + +#### Align the centre of the `inhale_BH_CT` image with the centre of `CT_for_PET.nii.gz` + +This can be achieved by following these steps: + +* Save the `inhale_BH_CT` DICOM image as a NIfTI image, `inhale_BH_CT.nii.gz`. +* Calculate the centre of `CT_for_PET.nii.gz` in voxel coordinates. +* Calculate the centre of `CT_for_PET.nii.gz` in world coordinates. +* Calculate the centre of `inhale_BH_CT.nii.gz` in voxel coordinates. +* Calculate the centre of `inhale_BH_CT.nii.gz` in world coordinates. +* Calculate the translation (in world coordinates) required to align the images. +* Use this translation to modify the affine transform for `inhale_BH_CT.nii.gz`. +* Save the aligned image as `inhale_BH_CT_aligned.nii.gz` + +Try writing code to implement these steps yourself, based on what you have learned from this practical. + +:::::::::::::::::::::::::::::::::: hint + +The centre of the image in voxel coordinates can be calculated as $\frac{NumVox - 1}{2}$ + +You can use the affine transform to transform the voxel coordinates to world coordinates, as seen earlier in the practical. Don't forget that you must use homogeneous coordinates, i.e. the coordinates are stored in a 4x1 vector, with the last element set to 1. + +The translation should not replace the image origin values in the affine header - it should be added to (or subtracted from) from them. + +:::::::::::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::::: + +Use ITK-SNAP to verify that you have correctly aligned the centres of the images. + +![Figure. ITK-SNAP after loading `inhale_BH_CT_aligned.nii.gz`](fig/itk-snap-aligned-images-1.png) +![Figure. ITK-SNAP showing the cprresponding slices of `CT_for_PET.nii.gz`](fig/itk-snap-aligned-images-2.png) + +## 4. References +#### Data +* [The Cancer Imaging Archive (TCIA)](https://www.cancerimagingarchive.net/collection/ct-vs-pet-ventilation-imaging/) +* [DICOM (Digital Imaging and Communications in Medicine)](https://www.dicomstandard.org/) +* [NIfTI (Neuroimaging Informatics Technology Initiative)](https://nifti.nimh.nih.gov/) +* [Anderson Winkler's blog on the NIfTI file format]( https://brainder.org/2012/09/23/the-nifti-file-format/) + +#### ITK-SNAP +* [ITK-SNAP website](http://www.itksnap.org/pmwiki/pmwiki.php?n=Main.HomePage). +* [Tutorial: Getting Started with ITK-SNAP](http://www.itksnap.org/docs/viewtutorial.php) +* [ITK-SNAP 3.x Training Class Final Program](http://itksnap.org/files/handout_201409.pdf) + +#### NiBabel +* [NiBabel website](https://nipy.org/nibabel/) +* [Documentation](https://nipy.org/nibabel/#documentation) +* [Working with NIfTI images](https://nipy.org/nibabel/nifti_images.html) +* [Understanding world coordinate systems](https://nipy.org/nibabel/coordinate_systems.html) +* [Neurological vs radiological views](https://nipy.org/nibabel/neuro_radio_conventions.html) +* [How to calculate the affine for DICOM images](https://nipy.org/nibabel/dicom/dicom_orientation.html) + +#### Misc. +* [Quaternions](https://en.wikipedia.org/wiki/Quaternion) +* [Homogeneous coordinates](https://en.wikipedia.org/wiki/Homogeneous_coordinates) + + + diff --git a/episodes/practical2.Rmd b/episodes/practical2.Rmd new file mode 100644 index 00000000..e4cec1c6 --- /dev/null +++ b/episodes/practical2.Rmd @@ -0,0 +1,232 @@ +--- +title: 'Image transformations in medical imaging' +teaching: 10 +exercises: 2 +--- + +:::::::::::::::::::::::::::::::::::::: questions + +- What are the steps for applying transformations to medical images? +- How can multiple transformations be combined and applied to the same image? +- What are the advantages of pull interpolation over push interpolation? + +:::::::::::::::::::::::::::::::::::::::::::::::: + +::::::::::::::::::::::::::::::::::::: objectives + +- Learn how to apply various transformations, such as translations and rotations, to medical images. +- Understand the process of composing multiple transformations. +- Differentiate between push and pull interpolations. + +:::::::::::::::::::::::::::::::::::::::::::::::: + +In these exercises, we will use a 2D slice from a lung MRI image, *lung_MRI_slice.png*, which you can find in the zipped data folder under `practical2`. + +This tutorial uses Python. There is a template Jupyter notebook for you to work on `practical2-exercises.ipynb`. You have also been provided with some utility functions in `utils.py`. When using a function from this file, you should read it and make sure you understand what it does. + +* Run the first cell of the Juypter notebook. This imports the required libraries, and sets the *matplotlib* library to display figures in separate windows. This is required for the animations in these exercises to run correctly. + +## 1. Loading and displaying images +You can load the 2D lung MRI image using the `imread` function from the *scikit-image* python library. If you check the data type of the image you will see it is stored as 8-bit integers. As you know, you should convert the image to double precision so that errors do not occur when processing the image due to the limited precision. + +* Run the next cell of the notebook to load the image and convert it to double precision. + +The *scikit-image* and *matplotlib* libraries used in these exercises index images using matrix coordinates. The first coordinate is the row number (i.e. the y coordinate of the image) and the second coordinate is the column number (i.e. the x coordinate of the image). The rows are also numbered from top to bottom. + +It is possible to do image processing using matrix coordinates, but it can get confusing, especially when working with deformation fields or in 3D. Therefore, the approach used in these exercises is to store the image memory in 'standard orientation', so that the first coordinate indexes the x (horizontal) dimension and the second coordinate indexes the y (vertical) dimension, and the first pixel (0,0) is at the bottom-left of the image. Then, when displaying the image a wrapper function such as the `dispImage` function provided in *utils.py* can be used to reorientate the image into matrix coordinates as required for the `imshow` function from *matplotlib*. + +Therefore, before proceeding you should reorientate the image into ‘standard orientation’. This can be done by first taking the transpose of the matrix (switching x and y dimensions) and then flipping along the second dimension (moving the first pixel from the top to the bottom of the image). + +* Write your own code in the next cell of the Juypter notebook to reorientate the image. It will then display the image using the `dispImage` function from *utils.py*. + * If you have done this correctly, the image should look like this: +![Figure. Reorientated image](fig/trans_disp_image.png) + +## 2. Translating and resampling images +* Edit the code in the next cell of the notebook to create an affine matrix representing a translation by 10 pixels in the x direction and 20 pixels in the y direction. The cell already contains the code to create a 3 x 3 identity matrix, which can be edited as required. + +:::::::::::: discussion + +#### Question + +* Can you remember Why we use a 3 x 3 matrix for a 2D affine transformation? + +If you are not sure of the answer to this question (or any of the others in the practicals) look back at the lecture material or ask one of the assistants to explain. + +:::::::::::: + +:::::::::::: callout + +#### Note + +*numpy* has a [matrix class](https://numpy.org/doc/stable/reference/generated/numpy.matrix.html), but the documentation recommends not using it. We are using *numpy* arrays for this exercise. +Matrix multiplication can be performed between two arrays using the @ operator or the `numpy.matmul` function + +:::::::::::: + +* Edit the next cell to create a deformation field from the affine matrix, and resample the image with the deformation field. + * You should use the provided `defFieldFromAffineMatrix` and `resampImageWithDefField function` for this. The function definitions in *utils.py* include comments explaining what the inputs and outputs of the functions should be. + +* Now run the next cell to display the transformed image. + * If you have done this correctly it should appear like this: +![Figure. Translated image](fig/trans_disp_image_translated.png) + +:::::::::::: discussion + +#### Question + +* Does the image appear as expected? +* Has the image been translated in the direction you expected? + +The `resampImageWithDefField` function uses pull-interpolation, so the image will appear to have been transformed by the inverse of the transformation in the affine matrix (i.e. it has been translated by -10 and -20 pixels in the x and y directions respectively). We will try using push-interpolation instead of pull-interpolation later in the exercises to see the difference. + +:::::::::::: + +* Run the next cell to check what value has been assigned to pixels that were originally outside the image by printing the value of the top right pixel (255,255). + +This is known as the ‘padding value’ or ‘extrapolation value’. A value of NaN (not a number) is often used to indicate that the true value for these pixels is unknown, and therefore they should be ignored when calculating similarity measures during image registration. + +The `resampImageWithDefField` function uses linear interpolation by default. It is implemented using *scipy*’s `interpn` function, so can also use the other interpolation methods available for this function. These include *nearest neighbour* and *splinef2d*, which is an efficient implementation of cubic interpolation using splines (note – it also has an interpolation method called *cubic*, which gives very similar results, but is much slower). + +* Write code to resample the image using *nearest neighbour* and *splinef2d* interpolations, and display the results in separate figures. + +The results should look the same as the result from using *linear* interpolation. However, the effects of different interpolation methods can be subtle and difficult to spot. To make any differences between the results clear we can display difference images, i.e. one image minus another, between the original result with *linear* interpolation and the results using *nearest neighbour* and *splinef2d* interpolations. + +* Write the code to display the difference images. + * They should appear like this: +![Figure. Difference images](fig/trans_diff_image.png) +(the difference image for *nearest neighbour* is on the left, and *splinef2d* on the right) + +:::::::::::: discussion + +#### Question + +The values in the difference image for *splinef2d* are very small, in the order of 10-14, and can be attributed to numerical errors. So, the different interpolation all give (very almost) the same results. + +* Do you understand why this is the case? + +:::::::::::: + +* Now write code that repeats the steps above using a translation of 10.5 pixels in the x direction and 20.5 pixels in the y direction. + * The difference images should look like this: +![Figure. Difference images](fig/trans_diff_image_10.5.png) +(the difference image for nearest neighbour is on the left, and splinef2d on the right) + +By default, the `dispImage` function displays images by scaling the values so that they use the full intensity range, i.e. the lowest value in the image is set to black and the highest to white. However, this can be misleading when comparing images, as the grey values seen in the images will not correspond to the same intensity values. Therefore, it is often a good idea to ensure exactly the same intensity range is used when displaying and comparing different images by using the `int_lims` input to the `dispImage` function. + +* Write code to redisplay the difference images, in both cases using an intensity limits of [-20, 20]. + * The difference images should now appear like this: +![Figure. Difference images](fig/trans_diff_image_intlims.png) +(the difference image for nearest neighbour is on the left, and splinef2d on the right) + +::::::::::::::::::: discussion + +#### Question + +* Why do the different interpolation methods now give different results? + +* And why is the difference larger for the nearest neighbour interpolation? + +::::::::::::::::::: + + +## 3. Rotating images +* Add code to the next cell of the notebook to implement a function to calculate the affine matrix corresponding to a rotation about a point, P. + * The inputs to the function should be the angle of rotation (in degrees) and the coordinates of the point and the output should be the affine matrix. + +::::::::::::::::::: spoiler + +#### Hints + +You will need to convert the angle from degrees to radians, as *numpy* `sin` and `cos` functions expect the input in radians. + +You will need to create 3 affine matrices - one representing the rotation, and two representing the translations from the origin to the point p and from p to the origin. + +You will then need to compose the transformations (using matrix multiplication) in the correct order and return the result. + +::::::::::::::::::: + +* Now write code that uses the function to calculate the affine matrix representing an anticlockwise rotation of 5 degrees about the centre of the image. Then transforms the original image using the rotation you just created and display the result. + * Linear interpolation should be used when resampling the image, and the intensity limits from the original image used when displaying the results. + * The result should look like this: +![Figure. Rotated image](fig/trans_rot_image.png) +* Apply the same transformation again to the resampled image and display the result. This should be repeated 71 times, so that the image appears to rotate a full 360 degrees. + * It is necessary to add a short pause (i.e. 0.05 seconds) using `plt.pause(0.05)` so that the figure updates each time the new image is displayed. + * The final image should look this this: +![Figure. Rotated image](fig/trans_rot_final_image.png) + * Here is the animation: +![Figure. Rotated image](fig/trans_rot_final_image.gif) + +You will notice that the image gets smaller and smaller as it rotates. This is because of the NaN padding values – when a pixel value is interpolated from one or more NaN values it also gets set to NaN, so the pixels at the edge of the image keep getting set to NaN, and the image gets smaller after each rotation. + +* To prevent this, repeat your code so that it uses a padding value of 0 rather than NaN and run your new code. + * Before the for loop You will need to first resample the original image (also using a padding value of 0), display it, and add a short pause so that the displayed image is updated. + * This time the final image should look like this: +![Figure. Rotated image](fig/trans_rot_final_image_pad0.png) + * Here is the animation: +![Figure. Rotated image](fig/trans_rot_final_image_pad0.gif) + +::::::::::::::::::: discussion + +#### Question + +You will notice that the corners of the image still get ‘rounded off’ as it rotates so that it has become a circle after rotating 90 degrees. + +* Do you understand why this happens? + +::::::::::::::::::: + +* Now add code to the template script to repeat the above, but first using *nearest neighbour* interpolation, and then using *splinef2d*. + * The final images should look like this: +![Figure. Rotated image](fig/trans_rot_image_nn_spline.png) + * Here is the animation for *nearest neighbour* interpolation: +![Figure. Rotated image](fig/trans_rot_final_image_pad0_nn.gif) + * And for *splinef2d*: +![Figure. Rotated image](fig/trans_rot_final_image_pad0_spline.gif) + +The blurring artefacts and the ‘rounding off’ of the images seen above are caused by multiple resamplings of the image. Resampling essentially makes a copy of the image, and the interpolation slightly degrades that image. Multiple resamplings are like copies of copies of copies - you start to notice the cumulative degradation of the images. This can be prevented by composing the rotations into a single transformation and then applying the resulting composed transformation to the original image instead of resampling the result image each time. + +## 4. Composing transformations +* Write code that makes an animation of the rotating image as above (using *linear* interpolation), but composes the transformations to avoid multiple resamplings of the image. + +::::::::::::::::::: spoiler + +#### Hints + +You will need to store the current rotation, `R_current` as well as the original rotation. + +`R_current` will initially be equal to `R` and should be used to resample the original image and display the result prior to starting the loop. + +At each iteration of the loop `R_current` will be updated by composing it with `R`, and then used to resample the original image and display the result. + +::::::::::::::::::: + +* Repeat the animation using *nearest neighbour* and *splinef2d* interpolation. + * The resulting animation should look like this for *linear*: +![Figure. Rotated image](fig/trans_rot_final_image_compose_linear.gif) + * And this for *nearest neighbour*: +![Figure. Rotated image](fig/trans_rot_final_image_compose_nearest.gif) + * And this for *splinef2d*: +![Figure. Rotated image](fig/trans_rot_final_image_compose_splinef2d.gif) + +You should notice that the corners of the images do not get ‘rounded off’ during these animations. Furthermore, you should notice that the intermediate images are different for the different interpolation methods (although this is only really noticeable for nearest neighbour interpolation), but the blurring artefacts do not get worse as the animation progresses, and the final image should be the same as the original image. + +## 5. Push interpolation +As discussed in the lectures, it is possible to resample an image using push-interpolation, but it is far less computationally efficient than using pull-interpolation. + +* Copy your code from above (that composes the transformations and uses *linear* interpolation) and modify it to perform push interpolation instead of pull interpolation by using the `resampImageWithDefFieldPushInterp` function instead of the `resampImageWithDefField` function. + * Note, you cannot provide a `pad_value` for `resampImageWithDefFieldPushInterp` + * The resulting animation should look like this: +![Figure. Rotated image](fig/trans_rot_final_image_push.gif) + +::::::::::::::::::: discussion + +#### Question + +* Other than the animation being much slower, what other difference do you notice? + +* If you have time experiment with using different angles (smaller/larger, positive/negative) and rotating about a different point. Try this for both pull and push interpolation, and for different interpolation methods. Make sure you understand all the results you get. + +::::::::::::::::::: + + + diff --git a/episodes/practical3.Rmd b/episodes/practical3.Rmd new file mode 100644 index 00000000..a8b0ba4f --- /dev/null +++ b/episodes/practical3.Rmd @@ -0,0 +1,116 @@ +--- +title: 'Exploring similarity measures for cost functions' +teaching: 10 +exercises: 2 +--- + +:::::::::::::::::::::::::::::::::::::: questions + +- What measures are used to calculate similarity between two images? +- What are some commonly used similarity measures and how do they differ in implementation and functionality? + +:::::::::::::::::::::::::::::::::::::::::::::::: + +::::::::::::::::::::::::::::::::::::: objectives + +- Understand how common similarity measures for image registration are calculated and used. +- Learn about the advantages and disadvantages of these metrics. + +:::::::::::::::::::::::::::::::::::::::::::::::: + +In these exercises, we use a three 2D images: + +* `ct_slice_int8.png`: an axial CT slice of a head, pixels are stored as unsigned 8-bit integers. +* `ct_slice_int16.png`: the same axial slice CT slice, pixels are stored as unsigned 16-bit integers. +* `mr_slice_int16.png`: the corresponding axial MR slice, pixels are stored as unsigned 16-bit integers. + +You can find in the zipped data folder under `practical3`. + +This tutorial uses Python. There is a template Jupyter notebook for you to work on *practical3-exercises.ipynb*. You have also been provided with some utility functions in `utils.py`. When using a function from this file, you should read it and make sure you understand what it does. + +## 1. Loading and displaying the images +Load the three 2D images. Display their data type and check these match the expected data types. + +Convert these images to double, reorient them into 'standard orientation' and display each image in a separate figure using the `dispImage` function from *utils2.py*. + +The images should appear like this: + +![Figure. Displayed images](fig/p3-images.png) + +The left and middle images are the CT images, while the image on the right is an MRI. +You should notice that the two CT images appear exactly the same, but if you examine the actual values in the images you will see that the images use different intensity ranges. If you move your cursor over the image the intensity values can be seen next to the figure: + +![Figure. Intensity values](fig/p3-int8-intensities.png) +In the 16 bit version, the intensity values correspond to Hounsfield Units (or rather Hounsfield Units + 1000, since air has a value of -1000, but the image cannot contain negative values as it is stored as unsigned integers), but in the 8 bit image the values have been scaled as the image can only contain values from 0 to 255. + +## 2. Rotations + +The template script contains some code to rotate each of the images between -90 degrees and 90 degrees, in steps of 1 degree. Edit the code so that on each iteration of the loop it uses the `affineMatrixFromRotationAboutPoint` function from *utils2.py* to create an affine matrix representing an anti-clockwise rotation by theta degrees about the point 10,10, and uses the `defFieldFromAffineMatrix` function to create the corresponding deformation field. Then resample each of the 3 images using the `resampImageWithDefField` function and display the transformed 8-bit CT image using the `dispImage` function. + +If this has been implemented correctly the image should appear to rotate clockwise, starting with a rotation of -90 degrees (so mostly being ‘off to the left’ of the original image) and finishing with a rotation of +90 degrees (so mostly being ‘off the bottom’ of the original image), as shown in the intermediate images below: + +![Figure. Rotated image](fig/p3-rotated-image.png) + +![Figure. Rotated image](fig/rotated_images.gif) + +Make sure you understand why the image appears to rotate clockwise when the function produces an affine matrix representing an anti-clockwise rotation. + +Note, the default padding value of NaN should be used when resampling the image so that pixels from outside the original images are ignored when calculating the similarity measures +below. + +### 2.1. Sum of Squared DIfferences (SSD) + +In this section, we are going to observe how Sum of Squared Differences (SSD) varies between our original image and our rotated image. Edit the code so that on each iteration of the loop it calculates and stores the SSD between: + +* 1) The original 16-bit CT image and the transformed 16-bit CT image +* 2) The original 16-bit CT image and the transformed 8-bit CT image +* 3) The original 16-bit CT image and the transformed MR image +* 4) The original 8-bit CT image and the transformed 8-bit CT image + +Now rerun the cell with the for loop (you may want to comment out the lines that display the image and pause so that the code runs faster). + +Plot the SSD values on the y-axis against the theta on the x-axis for each of the four cases above. + +The plots should look like this: + +![Figure. SSD values](fig/p3-ssd-values.png) + +Note that: + +* The SSD reaches a minimum (of 0) for cases 1 and 4 when the images are in alignment. +* For cases 2 and 3 SSD does not have a minimum when the images are aligned. +* For all cases the SSD decreases as the overlap between the images decreases. +* Although the shape of the SSD curve for cases 1 and 4 is the same, the values of the SSD are different by 2 orders of magnitude. + +Make sure you understand why you get these results. + +### 2.2. Mean of Squared Differences (MSD) +Now edit the code so that it rotates the images as above but calculates the MSD instead of the SSD at each iteration, and then plots the MSD values. + +The plots should look like this: + +![Figure. MSD values](fig/p3-msd-values.png) + +Note that: + +* The MSD for cases 1 and 4 does not decrease as the amount of overlap decreases. +* The shape of the MSD curves for cases 1 and 4 are the same but the values are larger for case 1. +* The MSD values for cases 2 and 3 are lower for negative values of theta and higher for positive values (this one is a bit tricky!). + +Make sure you understand why you get these results. + +### 2.3. Normalized Cross Correlation (NCC) and Normalized Mutual Information (NMI) +Implement the `calcNCC` function so that it calculates the Normalised Cross Correlation (NCC) between two images. + +Edit the code so that it uses your `calcNCC` function and the `calcEntropies` function from `utils.py` to calculate the (NCC) and the joint and marginal entropies ($H_{AB}$, $H_A$, and $H_B$) instead of the MSD/SSD at each iteration. Also calculate the Mutual Information (MI) and Normalised Mutual Information (NMI) from the entropy values, and plot the results for NCC, $H_{AB}$, MI, and NMI. + +The plots should look like this: + +![Figure. NCC values](fig/p3-ncc-values.png) +![Figure. entropy values](fig/p3-entropies-values.png) +![Figure. MI values](fig/p3-mi-values.png) +![Figure. NMI values](fig/p3-nmi-values.png) + +Do the different measures perform as expected for the different cases? Based on these results, which measures are suitable for registering the different pairs of images? Does this agree with what you were taught in the lecture? + +Make sure you understand all the results you get. diff --git a/episodes/practical4.Rmd b/episodes/practical4.Rmd new file mode 100644 index 00000000..c243935e --- /dev/null +++ b/episodes/practical4.Rmd @@ -0,0 +1,196 @@ +--- +title: 'Understanding the Demons registration algorithm' +teaching: 10 +exercises: 2 +--- + +:::::::::::::::::::::::::::::::::::::: questions + +- How does the Demons registration algorithm achieve image registration? +- What parameters should be adjusted to optimise the performance of the Demons registration algorithm? +- What are the key strengths and potential limitations of using the Demons registration algorithm? + +:::::::::::::::::::::::::::::::::::::::::::::::: + +::::::::::::::::::::::::::::::::::::: objectives + +- Explore the underlying principles and methods of the Demons registration algorithm. +- Understand the key parameters that influence the performance of the algorithm and how to tune them effectively. +- Visualise and interpret the results by generating plots of deformation fields and Jacobian determinant maps. + +:::::::::::::::::::::::::::::::::::::::::::::::: + +In these exercises, we use three 2D sagittal MR slices `cine_MR_1.png`, `cine_MR_2.png`, and `cine_MR_3.png` showing the lung, liver, and surrounding anatomy, at different points in time during free-breathing. Image 1 is at end-exhalation, image 2 at the end-inhalation for a normal breath, and image 3 at end-inhalation for a deep breath. You can find them in the zipped data folder under `practical4`. + +This tutorial uses Python. There is a template Jupyter notebook for you to work on `practical4-exercises.ipynb`. You have also been provided with some utility functions in `utils.py`. When using a function from this file, you should read it and make sure you understand what it does. +There is also a function that implements the Demons registration algorithm in `demonsReg.py`. You should look at the code in `demonsReg.py` and check you can follow what it is doing. You should be able to understand what each section of code is doing, even if you don’t fully understand every line. + +## 1. Loading and displaying images +Load and display the three images. Convert the images to doubles and reorientate them as in the previous exercises. + +The images should look like this: + +![Figure. Displayed images](fig/p4-displayed-images.png) + +Compare the images to each other and observe the motion that occurs between the images. Difference images (one image minus another, e.g. as shown below) are a good way seeing where there are differences between the images, but it can be difficult to tell how things have moved between the images from the difference image. + +![Figure. Difference images](fig/p4-diff-images.png) + +The image on the left is the difference image between the first and second images, and the one on the right is the difference image between the first and third images. + +It is often easier to appreciate how things have moved/changed between images if you 'flip' between the images. To achieve this, we can use the `dispImageFlip` function. This function takes in the two images you want to look at as arguments, and you can flip between the two using the arrow keys on your keyboard. + +See how the images differ from each other. Pay particular attention to motion that you think could be challenging for the demons algorithm to recover, such as very large motion and deformation, or sliding motion that can occur between the lung/liver and surrounding anatomy. + +## 2. Demons algorithm +### 2.1. `demonsReg` function +Have a look at the comments at the start of `demonsReg` function in `demonsReg.py` (which are also displayed as the help info for the function). There are two required input parameters: the **source** and **target** images for the registration. There are also a number of optional input parameters, which have default values assigned to them that will be used if an alternative value is not specified. We will explore the effects of these parameters as we go through these exercises. The outputs of the function are the final warped image and the corresponding deformation field. + +Start by running a registration with `cine_MR_1` as the source image and `cine_MR_2` as the target image and using default values for all the other parameters. +When you start the registration, a figure appears with 3 subplots. The first subplot on the left shows the warped image (transformed source image), the second subplot in the middle shows the current deformation field and the third subplot on the right shows the update to the deformation field for the current iteration. The deformation field (middle plot) is displayed as a deformed grid whereas the update is displayed as a vector field. All three images update during the registration. + +![Figure. Initialisation of the demons algorithm](fig/p4-reg-first-step.png) + +While the registration is running, you will notice that the figures are updated and that text is outputted in the notebook stating the current resolution level, iteration number, and value of the Mean Squared Difference (MSD) between the target image and the current warped image. The MSD is the similarity metric used by this demons registration. At the end of the registration, the figure should look like this: + +![Figure. Output of the demons algorithm](fig/p4-reg-output.png) + +Here is a summary of the registration: +![Figure. Output of the demons algorithm](fig/registration_animation.gif) + + +At the end of the registration, a second figure appears. The figure allows us to flip between the source, target and warped image and shows the deformation field and difference image. We can also display the Jacobian (more details later). + +![Figure. Output of the demons algorithm](fig/p4-reg-output2.png) + +If you look at the text output you will see that the final MSD value is 80.20 (rounded to 2 decimal places – more are shown in the console), and that 74 iterations were performed at the final (3rd) resolution level. If you scroll up in the console you will see that 24 iterations were also performed at the 2nd resolution level, and 17 iterations at the 1st level. +:::::::::::::::challenge +Within each resolution level, the MSD value decreases with each iteration. However, the MSD increases between resolution levels. Do you know why it does this? +::::::::::::::::hint +It’s not just because there are more pixels in the higher resolution images, as using the mean rather than sum of squared differences should account for different numbers of pixels in the images. +:::::::::::::::::::: +:::::::::::::::::::: + +Compare the final warped image with the target and source images. How well do you think the registration did at aligning the images? Which parts of the images are well aligned, and which parts are less well aligned? + +If you look at the deformation field you will see that it corresponds to the deformation that would be applied to the warped image to get back to the source image (rather than the deformation that would be applied to the source image to get to the warped image). Make sure you understand why this is the case (hint – pull-interpolation). + +### 2.2. Regularisation + +You can specify the values of the optional input parameters when calling the function, e.g. `demonsReg(cine_MR_1, cine_MR_2, sigma_elastic=0)` will run the registration with sigma_elastic set to 0, i.e. with no elastic like regularisation applied. Try running this registration. At the end of the registration, the figures should look like this: + +![Figure. Output of the demons algorithm](fig/p4-reg-output-sigelastic0.png) +![Figure. Output of the demons algorithm](fig/p4-reg-output2-sigelastic0.png) + +Here is a summary of the registration: +![Figure. Output of the demons algorithm](fig/registration_animation_elastic0.gif) + +You will see that parts of the final warped image looks very similar to the target image, and the final MSD is 79.01, which is a little lower than for the default parameters. However, you can see that the final deformation field is very ‘crumpled’ and does not look very physically plausible. Therefore, for most applications this would not be considered a good registration result. The reason that the result is like this is that without any elastic-like regularisation applied to the deformation field it is free to ‘evolve’ in any way that better aligns the intensities, even if this results in implausible transformations, and highlights the need for using suitable regularisation. + +Now rerun the registration with sigma_elastic at its default value again but this time setting sigma_fluid to 0. At the end of the registration, the figures should look like this: + +![Figure. Output of the demons algorithm](fig/p4-reg-output-sigfluid0.png) +![Figure. Output of the demons algorithm](fig/p4-reg-output2-sigfluid0.png) + +Here is a summary of the registration: +![Figure. Output of the demons algorithm](fig/registration_animation_fluid0.gif) + +You will now see that the updates are much larger and noisier as no regularisation is being applied to them, but the deformation field is still relatively smooth as it does have regularisation applied. If you compare the result from this registration to the result obtained with the default parameters (remember, you can store the outputs of the function so that you can directly compare the results from different registrations) you will see that the results are very similar. However, the final MSD value is lower (72.73), but some regions of the deformation field are slightly less smooth. Which of the results should be considered ‘best’ will depend on exactly what the application is. But it should be clear that the elastic-like regularisation usually has a larger effect on the result than the fluid-like regularisation does. Now try running several more registrations with different levels of elastic-like and fluid-like regularisation to get a good feel for how the combination of these parameters affects the results, and which ranges of values produce good/bad results. + +:::::::::::::::::::::::callout +You can control how often the plots are updated by setting the value of the `disp_freq` input parameter. The registrations themselves run relatively quickly, but updating the images takes some time. This makes the registrations run considerably slower but allows you to see how the deformation field and warped image change during the course of the registration, which can be very useful for understanding why the registration produces the result it does. If you set `disp_freq` to 0 then the images will only update at the start of each level and to show the final result. If you do this you will notice that the registration runs much quicker, but you miss all the details of how the registration produces the result. For most of these exercises, keep the `disp_freq` value relatively low, e.g. at the default value of 5 iterations or less, so that you can see how the deformation field and warped image change as the registration is progressing. +::::::::::::::::::::::::::: + +### 2.3. Resolution levels +You should have noticed that the registration is using a multi-resolution approach, with 3 resolution levels by default. You can specify the number of resolution levels to use with the `num_lev` input parameter. Try running a registration with `num_lev` set to 1, i.e. so that only a single resolution level is used, and the multi-resolution approach is not applied. At the end of the registration, the figures should look like this: + +![Figure. Output of the demons algorithm](fig/p4-reg-output-numlev1.png) +![Figure. Output of the demons algorithm](fig/p4-reg-output2-numlev1.png) + +Here is a summary of the registration: +![Figure. Output of the demons algorithm](fig/registration_animation_level1.gif) + +::::::::::::::::::::discussion +Do you think using a multi-resolution approach helps? Why? +::::::::::::::::::::: + +Now try running the registration with 6 resolution levels (if you use 7 levels then the first level only contains 4 pixels, so there is not very much to change and improve. Anything greater than 7 levels will result in an error, since no further down-sampling can be done). The final MSD with six levels is worse (MSD = 86.24 ) then with the default 3 resolution levels (MSD = 80.20). However, visually the results from six resolution levels appears to align the structures within the lung better than the result from using three levels. It is the structures outside the lung, e.g. the ribs, that are not aligned as well. Again, assessing the "best" registration will be dependent on the exact application. + +Now experiment with varying the values for the regularisation parameters as well as the number of resolution levels, and see how the results are affected. + +### 2.4. Jacobian determinant + +When investigating the effects of the regularisation parameters you should have noticed that some parameters will result in deformation fields do not look very plausible. In particular, sometimes the deformation field will appear to fold over itself, indicating that folding has occurred, which is often undesirable when performing registrations. + +For example, if you perform a registration with the following parameters: `sigma_elastic = 0.5, sigma_fluid = 1, num_lev = 3`, the result is the deformation field shown below, and as can be seen in the zoomed-in image on the right, the grid has folded over itself in some regions. + +![Figure. Deformation field](fig/p4-def-field-zoom.png) +You have been provided with a function, `calcJacobian`, to calculate the Jacobian determinant (and optionally the full Jacobian matrix) at each pixel in the deformation field, and return a map (image) of the Jacobian determinants. Have a look at this function and make sure you understand how it works. + +You can display the Jacobian by using the up and down arrows on the second output figure at the end of the registration. + +![Figure. Output of the demons algorithm](fig/p4-reg-output2-sigelastic0.5.png) + +You can also plot the Jacobian yourself by calling the `calcJacobian` and `dispImage` functions. It is useful to add a colourbar and to use the `jet` colourmap. You can also plot a binary mask of the Jacobian by setting all values <= 0 to 1s. + +![Figure. Jacobian](fig/p4-jac-binarymask.png) + +You can notice that the image is very ‘patchy’, and does not show smooth expansions/contractions as would be expected during breathing. The colourbar shows us that the Jacobian contains some very large values (>4) as well as negative values, indicating that folding has occurred. + +If you compare the location of the pixels with folding in the binary image to the deformation field you will notice that they do not align with where you can see folding occurring in the grid. This is because the grid shows the deformation in the source image space whereas the binary image corresponds to the pixels in the warped image space. + +### 2.5. Compositions +With the `demonsReg` function it is possible to compose the updates rather than add them, so that the deformation field will always represent a diffeomorphic transformation (as the updates are themselves always diffeomorphic), and hence will not contain any folding. If you look at lines 150-157 of the `demonsReg` function and you will see that the composition of two deformation fields can be implemented by treating the current deformation field as an image, and resampling it with (a deformation field derived from) the update to the transformation. To see why, we can write the equation for resampling an image as: + +$I_{resamp}(\vec{x})=I_{orig}(T(\vec{x}))$ + +i.e. the intensity in $I_{resamp}$ at $\vec{x}$ is the intensity in $I_{orig}$ at $T(\vec{x})$ ($\vec{x}$ transformed by $T$). As $T(\vec{x})$) will typically not be the centre of a pixel in $I_{orig}$, the intensities in $I_{orig}$ must be interpolated to calculate $I_{resamp}$. + +The equation for composing the update, $T_{up}$, with the current transformation, $T_t$, is: + +$T_{t+1}(\vec{x}) = T_t(T_{up}(\vec{x}))$ + +i.e. to calculate the composed transformation, $T_{t+1}$, we first apply the update and then the current transformation (make sure you understand why they are applied in this order – remember order matters for composition!). + +You can see that the equation for composing the transformation is the same as the equation for resampling the image, with the current transformation, $T_t$, replacing the image, $I_{orig}$, and $T_{up}$ used as the transformation to resample the image/transformation. + +Now perform a registration using the same parameters as above but use composition rather than addition to update the transformation, i.e. use the following parameters: `sigma_elastic = 0.5, sigma_fluid = 1, num_lev = 3, use_composition = True`. + +Here is a summary of the registration: +![Figure. Output of the demons algorithm](fig/registration_animation_comp.gif) + +You will notice that the deformation field no longer contains any folding: + +![Figure. Deformation field](fig/p4-def-field-comp.png) + +but the Jacobian determinant map is still very patchy + +![Figure. Jacobian](fig/p4-comp-jac.png) + +and the max value is larger (note – a different intensity scale is used compared to the Jacobian determinant image on the previous page to show there is a different range of values). + +This indicates that the deformation is still not very physically plausible. And if you look at the warped image you will see that some of the features in the lungs are better aligned when addition was used for updating the transformation, although the MSD value is lower when composition is used (63.10) than when addition is used (64.07). But there is no folding when using composition and there is folding when using addition. + +So which registration do you think is best? + +### 2.6. Other parameters +As you will have seen, there are a number of other optional parameters that can be set. The parameters `disp_spacing`, `scale_update_for_display`, `disp_method_df`, and `disp_method_up` all effect how the deformation field and update are displayed but have no effect on the actual registration itself. The comments/help should explain how these parameters work. + +The `max_it` parameter sets the maximum number of iterations to be performed. The default value is 1000, but you should have noticed that usually the registration stops well before reaching this number of iterations. This is because the registration checks for an improvement in the MSD at each iteration, and if there has not been an improvement it will move on to the next resolution level or finish if it is already on the last level. This behaviour can be disabled by setting the `check_MSD` parameter to false. In this case the registration will run for the full max_it number of iterations at each level. + +Try running a registration with the default parameters and `check_MSD = False`. It is also recommended to set `disp_freq = 50` so that the registration does not take too long to run. By looking at the MSD values output to the console you will notice that some iterations now cause (small) increases to the MSD. You will also notice that the registration reaches a stage at each level where the result starts to ‘oscillate’, with a few iterations making minor increases to the MSD and then a few making minor decreases to the MSD, and barely any noticeable changes to the warped image or deformation field. The final MSD value is lower than when `check_MSD` was set to true (74.53 compared to 80.20), but at the expense of running many more iterations. We get a similar result (MSD = 75.21) by setting `max_it` to 500 and only running half as many iterations, indicating that for at least the last 500 iterations the result was mostly oscillating and not improving, but the problem is it is difficult to know how many iterations will be required before the results start to ‘oscillate’ as this will depend on both the images being registered and the other parameters. + +The last parameter that can be set is `use_target_grad`. This is false by default, and setting it to true will cause the registration to use the spatial gradient of the target image when calculating the demons forces rather than the spatial gradient of the source image, i.e. it will use the original demons forces rather than the common demons forces. If you run a registration with the default parameters and `use_target_grad` set to true and compare it to the results when `use_target_grad` is false, you will see there are differences in the results. However, each result aligns some regions/structures better than the other, and neither of the results are clearly the best overall (the final MSC is lower when `use_target_grad` is false, but this is not always the case when different parameters are used). + +### 2.7. Deeper breath image + +Now try registering `cine_MR_1.png` to `cine_MR_3.png`, which is a more challenging registration due to the deep breath taken when `cine_MR_3.png` was acquired. + +Experiment with different parameters – can you find some that give a good result? + +Try swapping the source and target image over. Do you think you get better results with `cine_MR_3.png` as the source or target image (or are they equally bad for both!)? + + + + + diff --git a/index.md b/index.md index af662764..89e087c6 100644 --- a/index.md +++ b/index.md @@ -2,8 +2,9 @@ site: sandpaper::sandpaper_site --- -This is a new lesson built with [The Carpentries Workbench][workbench]. + +Add course summary, lecture slides and timetable. [workbench]: https://carpentries.github.io/sandpaper-docs diff --git a/learners/setup.md b/learners/setup.md index 2445ea3a..731b47c2 100644 --- a/learners/setup.md +++ b/learners/setup.md @@ -1,8 +1,12 @@ --- title: Setup +editor_options: + markdown: + wrap: sentence --- - +````{=html} + +```` -## Dataset +## Course material - -Download the [data zip file](https://example.com/FIXME) and unzip it to your Desktop - -### CT-vs-PET-Ventilation-Imaging -CT Ventilation as a Functional Imaging Modality for Lung Cancer Radiotherapy from [TCIA](https://www.cancerimagingarchive.net/collection/ct-vs-pet-ventilation-imaging/). -We recommend to focus on exhale/inhale breath hold CT (BHCT) Dicoms. -Bare in mind that: -* The BHCT scans for CT-PET-VI-02 & CT-PET-VI-03 show very little motion between inhale and exhale -* The BHCT scans for CT-PET-VI-05 have a different number of slices between the inhale and exhale - -#### Extracting and renaming data paths -Original data directory paths for this work were renamed for easy interpretability as follows: -* `5.000000-Thorax Insp 2.0 B*` renamed as `inhale_BH_CT` -* `7.000000-Thorax Exp 2.0 B*` renamed as `exhale_BH_CT` -* `3.000000-CT Lung 3.0 B*` renamed as `CT_for_PET` -* `4.000000-Galligas Lung-0*` renamed as `PET` - -See how the original data paths and renamed paths look like -```bash -#Original directory paths -/CT-PET-VI-02$ tree -d -. -├── 1000.000000-PET SUV Factors-26496 [1 item, with size 42.1 kB] -├── 3.000000-CT Lung 3.0 B31f-61322 [175 items, totalling 92.4 MB] -├── 4.000000-Galligas Lung-03537 [159 items, totalling 51.5 MB] -├── 5.000000-Thorax Insp 2.0 B70f-29873 [167 items, totalling 88.1 MB] -├── 7.000000-Thorax Exp 2.0 B70f-73355 [167 items, totalling 88.1 MB] -└── 8.000000-RespLow 2.0 B30s 80 Ex - 100 In-44317 [810 items, totalling 426.9 MB] -#Renamed directory paths -. -├── CT-PET-VI-02 -│   ├── CT_for_PET -│   ├── exhale_BH_CT -│   ├── inhale_BH_CT -│   └── PET - - - - -#Original directory paths -/CT-PET-VI-03$ tree -d -. -├── 06-25-1996-RTRTRespLowBreathRateRNS Adult-43080 -│ ├── 1000.000000-PET SUV Factors-06580 -│ ├── 3.000000-CT Lung 3.0 B31f-08354 -│ └── 4.000000-Galligas Lung-15379 -└── 06-25-1996-RTRTRespLowBreathRateRNS Adult-87093 - ├── 5.000000-Thorax Insp 2.0 B70f-42000 - ├── 7.000000-Thorax Exp 2.0 B70f-45310 - └── 8.000000-RespLow 2.0 B30s 80 Ex - 100 In-54790 -#Renamed directory paths for `06-25-1996-RTRTRespLowBreathRateRNS Adult-43080` -├── CT-PET-VI-03 -│   ├── CT_for_PET -│   └── PET - - - - - -#Original directory paths -/CT-PET-VI-05$ tree -d -└── 12-10-1996-RTRTRespLowBreathRateRNS Adult-37585 - ├── 1000.000000-PET SUV Factors-63910 - ├── 3.000000-CT Lung 3.0 B31f-90638 - ├── 4.000000-Galligas Lung-26361 - ├── 5.000000-Thorax Insp 2.0 B70f-45546 - ├── 7.000000-Thorax Exp 2.0 B70f-31579 - └── 8.000000-RespLow 2.0 B30s 80 Ex - 100 In-16454 -#Renamed directory paths -├── CT-PET-VI-05 -│   ├── CT_for_PET -│   ├── exhale_BH_CT -│   ├── inhale_BH_CT -│   └── PET - - - - - -#Original directory paths -/CT-PET-VI-07$ tree -d -. -└── 02-14-1997-RTRTRespLowBreathRateRNS Adult-80512 - ├── 1000.000000-PET SUV Factors-99730 - ├── 4.000000-Galligas Lung-82849 - ├── 5.000000-Thorax Insp 2.0 B70f-20440 - ├── 7.000000-Thorax Exp 2.0 B70f-03567 - └── 8.000000-RespLow 2.0 B30s 80 Ex - 100 In-01782 -#Renamed directory paths -└── CT-PET-VI-07 - ├── exhale_BH_CT - ├── inhale_BH_CT - └── PET -``` +Start by downloading the course material [here](https://liveuclac-my.sharepoint.com/:u:/g/personal/rmapcdr_ucl_ac_uk/ETtoAQ1qXedMiU69mRdFXOIBcnKQ-nxMQ1c74gsPm7a6Sg?download=1). +Unzip it and make a note of where the folder is saved. +In this folder, you will find: +* A data directory which contains the datasets for each practical. +* 6 Jupyter notebooks with extensions `.ipynb` for each practical. +* 2 Python scripts with extensions `.py`. +* The `python-environment.yml` file which we will use in the steps below. -## Software Setup +### Dataset +We provide a summary of the different datasets for each practical here. Datasets are described at the beginning of each practical. -### Installing Dependencies in Python Virtual environments -Installing python virtual environment with dependencies (e.g., nibabel, etc). -See further instructions [here](https://github.com/HealthBioscienceIDEAS/Medical-Image-Registration-Short-Course/tree/main/_dependencies). +| Practical | Dataset | +|-------------------------------|-----------------------------------------| +| Practical 1 | [TCIA CT-vs-PET-Ventilation-Imaging](https://www.cancerimagingarchive.net/collection/ct-vs-pet-ventilation-imaging/) | +| Practical 2 | Lung MRI slice | +| Practical 3 | Head CT and MRI slices | +| Practical 4 | Lung MRI slices | +| Practical 5 | Dataset 5 | +| Practical 6 | Dataset 6 | -#### Using conda -```bash -conda create -n "mirVE" python=3.11 pip -c conda-forge -conda config --set pip_interop_enabled True -conda activate mirVE -cd ~/_dependencies -pip install -r requirements.txt -``` +## Python Setup -#### Using python virtual environment -```bash -# Installing dependencies in Ubuntu 22.04 -sudo apt-get install -y python3-pip -sudo apt-get install -y python3-venv -# Create path for venv -cd $HOME -mkdir mirVE -cd mirVE -# Create virtual environment -python3 -m venv mirVE -source mirVE/bin/activate -cd ~/_dependencies -python3 -m pip install -r requirements.txt -``` +### Installing Python -### ITK-Snap +For Python beginners, we recommend installing miniconda [here](https://www.anaconda.com/download/success). +This provides a basic installation of Python and conda, which allows us to create virtual environments. +Once this has been installed, you should be able to open a terminal. -#### Hardware Requirements -* Memory: "Overall 2Gb of memory should be enough for most users; 4Gb should suffice for users with large images." [:link:](http://www.itksnap.org/pmwiki/pmwiki.php?n=Documentation.HardwareRequirements#:~:text=Memory%20usage%20in%20SNAP%20is,for%20a%20512%20cubed%20image.) -* Multi - Operating System: - * ITK-SNAP binaries are distributed for Linux (32 and 64 bit)), - * Windows (64 bit), - * and MacOS Leopard (64 bit). +::::::::::::::::::::::::::::::::::::::::::: spoiler +### Opening a Terminal -::::::::::::::::::::::::::::::::::::::: discussion +* Windows: Click Start > Search for Anaconda Prompt > Click to Open +* macOS: Launchpad > Other Application > Terminal +* Linux: Open a terminal window -### Installation +::::::::::::::::::::::::::::::::::::::::::: -Installing ITK-SNAP in differen operating systems +### Creating an environment -::::::::::::::::::::::::::::::::::::::::::::::::::: +Once you have opened the terminal, move to the location of the course files +(use the `cd` command to move to the correct location - have a quick look [here](https://tutorials.codebar.io/command-line/introduction/tutorial.html) if you've never used the command line before). -:::::::::::::::: solution +:::::::::::::::::::::::::: discussion +### What does this look like? +Say you have saved the extracted zip file on your Desktop. +:::::::::::::::::::::::::: +:::::::::::::::::::::::::::::::::::: solution ### Windows +``` +cd Desktop\ideas-reg\ +``` +Note: the above tutorial for the command line is for Mac/Linux users. +In your Anaconda prompt in Windows, you can use `dir` instead of `ls` and `cd` without any arguments instead of `pwd`. +:::::::::::::::::::::::::::::::::::: + +:::::::::::::::::::::::::::::::::::: solution +### Mac/Linux +```bash +cd Desktop/ideas-reg/ +``` +:::::::::::::::::::::::::::::::::::: -- Open the following link: https://sourceforge.net/projects/itk-snap/files/itk-snap/4.2.0/ -- Download `itksnap-4.2.0-20240422-win64-AMD64.exe` -- Install it by following default settings. +Create the environment from the existing `python-environment.yml` file. +``` bash +conda env create --name ideas-reg -f python-environment.yml +``` +This environment has been tested on all the tutorials and should help avoiding issues when running them. -::::::::::::::::::::::::: +### VSCode +If you do not already have a Python IDE installed, we recommend using VSCode. +It is very convenient as it allows us to open both Python scripts and Jupyter notebooks in the same interface. +You can download it [here](https://code.visualstudio.com/download). -:::::::::::::::: solution +We recommend looking at the following instructions to get started: -### Linux (tested on Ubuntu 22.04 x64 GNU/Linux) +* [Python tutorial](https://code.visualstudio.com/docs/python/python-tutorial) +* [Instructions for Jupyter notebooks](https://code.visualstudio.com/docs/datascience/jupyter-notebooks) -- Open your Terminal to install it -```bash -sudo apt-get install tree # to visualise files and paths -FILENAME=itksnap-nightly-rel_4.0-Linux-gcc64.tar.gz #Length: 200943059 (192M) [application/x-gzip] -cd ~/Downloads/ -mkdir -p itksnap && cd itksnap -wget https://sourceforge.net/projects/itk-snap/files/itk-snap/Nightly/$FILENAME/download -O $FILENAME -tar -xvzf $FILENAME -``` +#### Open the `ideas-reg` folder +Start by opening the folder with the course content. Go on `File --> Open Folder` and navigate to the `ideas-reg` unzipped folder. +![Figure. VSCode with ideas-reg opened](fig/vscode.png) -::::::::::::::::::::::::: +On the left, you can see the files available in the directory. Open one of the practicals (`.ipynb`). This is a Jupyter notebook. +If this is the first time you open a Jupyter notebook in VSCode, you should be prompted to install some extensions (Python, Pylance, Jupyter and a few more). +Install them. +#### Select the environment +You now need to select the conda environment we created above. At the top right of the notebook, there should be a `Select Kernel` button. +![Figure. VSCode: selecting an environment](fig/vscode-select-kernel.png) -:::::::::::::::: solution +Click on Python Environment. You should see a dropdown with your existing environments. Select `ideas-reg`. +![Figure. VSCode: selecting an environment](fig/vscode-ideas-reg.png) -### MacOS +## Software Setup -- Open the following link: https://sourceforge.net/projects/itk-snap/files/itk-snap/4.2.0/ and download `itksnap-4.2.0-20240422-Darwin-x86_64.dmg` -- Double click the icon for the file you downloaded (e.g., itksnap-3.2.0-rc2-20140919-MacOS-x86_64.dmg). -- Drag the icon ITK-SNAP.app to the Applications folder. -- Drag the icon itksnap on top of the icon `usr_local_bin` in that folder. -- To add the ITK-SNAP launcher to your Dock, open the Applications folder and drag the ITK-SNAP.app icon onto the dock. +### ITK-SNAP +ITK-SNAP started in 1999 with SNAP (SNake Automatic Partitioning) and was first developed by Paul Yushkevich with the guidance of Guido Gerig. +It is open-source software distributed under the GNU General Public License. +It is written in C++ and leverages the [Insight Segmentation and Registration Toolkit (ITK)](https://itk.org/). +ITK-SNAP will be used in Practical 1. + +* The ITK-SNAP installer should be downloaded from [here](http://www.itksnap.org/pmwiki/pmwiki.php?n=Downloads.SNAP4). + * Version 4.2.2 was used when preparing the practical exercises. Other versions can also be used for the practical, although they may behave slightly differently to what is described in the practical. +* Once you have downloaded the installer you should run it to install ITK-SNAP. + * When you run the installer on Windows a window may pop up as shown below, with 'Don't run' as the only option given. + + ![](fig/itk-snap-install-1.png) + + * However, If you click on 'More info' a new 'Run anyway' button will appear which you can click on to run the installer. + + ![](fig/itk-snap-install-2.png) + + * Then all of the default options can be selected during installation. -:::::::::::::::::::::::::