<jupyter_start><jupyter_text>Medical Image Analysis workshop - IT-IST Image Filtering (edge preserving denoising, smoothing, etc.), Resampling and SegmentationLets load the image saved in previous notebook and view it.<jupyter_code>import SimpleITK as sitk
import itk
import itkwidgets as itkw
from ipywidgets import interactive
import ipywidgets as widgets
import numpy as np
image_itk = itk.imread('./cT1wNeuro.nrrd')

image_sitk = sitk.ReadImage('./cT1wNeuro.nrrd')
itkw.view(image_itk, cmap='Grayscale', mode='x')<jupyter_output><empty_output><jupyter_text>Filters Curvature Anisotropic Diffusion FilterThe Curvature Anisotropic Diffusion Filter ([Documentation SimpleITK](https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1CurvatureAnisotropicDiffusionImageFilter.html)) is widely used in medical images to denoise your image while preserving the edges. Note 1: ITK/SimpleITK often poses restriction into the pixel/voxel type to execute some operations. In this case we can use the Cast Image filter to convert our image of type short into float. Be careful when performing the cast operation from higher to lower pixel/voxel types (rescaling may be required)Note 2: There are several types of Anisotropic Diffusion Filters and another commonly used is the Gradient Anisotropic Diffusion filter (an example of that may be found on the Extras notebook).<jupyter_code>image_sitk_float = sitk.Cast(image_sitk, sitk.sitkFloat32)
smooth_cadf_sitk = sitk.CurvatureAnisotropicDiffusion(image_sitk_float)
itkw.view(smooth_cadf_sitk, cmap='Grayscale', mode='x')<jupyter_output><empty_output><jupyter_text>As previously mentioned, ITK is more verbose than SimpleITK but it more customizable and offers additional filters.The same result that was achieved by with three lines requires the following in ITK<jupyter_code>castImageFilter = itk.CastImageFilter[image_itk, itk.Image[itk.F,3]].New(image_itk)
castImageFilter.Update()
image_itk_float = castImageFilter.GetOutput()
curv_ani_dif_filter = itk.CurvatureAnisotropicDiffusionImageFilter[image_itk_float, itk.Image[itk.F,3]].New(image_itk_float)
curv_ani_dif_filter.Update()

smooth_cadf = curv_ani_dif_filter.GetOutput()
itkw.view(smooth_cadf, cmap='Grayscale', mode='x')<jupyter_output><empty_output><jupyter_text>Median Filter [Documentation SimpleITK](https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1MedianImageFilter.html)<jupyter_code># median_filter = sitk.MedianImageFilter()
# median_filter.SetRadius(2) # In pixels/voxels
# median_image = median_filter.Execute(image_sitk_float)

median_image = sitk.Median(image_sitk_float, [2, 2, 2])

itkw.view(median_image, cmap='Grayscale', mode='x')<jupyter_output><empty_output><jupyter_text>Sobel Filter [Documentation SimpleITK](https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1SobelEdgeDetectionImageFilter.html)<jupyter_code>sobel_edge_image = sitk.SobelEdgeDetection(image_sitk_float)
itkw.view(sobel_edge_image, cmap='Grayscale', mode='x')<jupyter_output><empty_output><jupyter_text>Laplacian Sharpening Image Filter[Documentation SimpleITK](https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1LaplacianSharpeningImageFilter.html)<jupyter_code>laplacian_sharped_image = sitk.LaplacianSharpening(image_sitk_float)
itkw.view(laplacian_sharped_image, cmap='Grayscale', mode='x')<jupyter_output><empty_output><jupyter_text>You can also compare two images using the following command<jupyter_code>itkw.checkerboard(median_image, laplacian_sharped_image, cmap='Grayscale', mode='x', pattern=5)<jupyter_output><empty_output><jupyter_text>Resampling an ImageThis image is isotropic with voxels of size 1mm x 1mm x 1mm. Often images come with different voxel sizes, and depending on the analysis we may have to normalize the voxel size accross the dataset. So lets learn how to change the voxel spacing/size from 1mm x 1mm x 1mm to 1.5mm x 2mm x 3mm.<jupyter_code>resample = sitk.ResampleImageFilter()
resample.SetInterpolator = sitk.sitkBSpline
resample.SetOutputDirection(image_sitk_float.GetDirection())
resample.SetOutputOrigin(image_sitk_float.GetOrigin())
new_spacing = [1.5, 2, 5]
resample.SetOutputSpacing(new_spacing)

orig_size = np.array(image_sitk_float.GetSize(), dtype=np.int)
orig_spacing = np.array(image_sitk_float.GetSpacing(), dtype=np.float)
new_size = orig_size * (orig_spacing / new_spacing)
new_size = np.ceil(new_size).astype(np.int)  # Image dimensions are in integers
new_size = [int(s) for s in new_size]
resample.SetSize(new_size)

resampled_image = resample.Execute(image_sitk_float)

print('Resample image spacing:', list(resampled_image.GetSpacing()))
print('Resample image size:', list(resampled_image.GetSize()))
print('Original image spacing:', list(image_sitk_float.GetSpacing()))
print('Original image size:', list(image_sitk_float.GetSize()))<jupyter_output>Resample image spacing: [1.5, 2.0, 5.0]
Resample image size: [171, 128, 36]
Original image spacing: [1.0, 1.0, 1.0]
Original image size: [256, 256, 176]<jupyter_text>Segmentation FiltersAnother common task on the medical image domain is the segmentation.In this example we will work with a different image.<jupyter_code>image_brainT1 = sitk.ReadImage('./brain_T1.nii.gz')
image_brainT1_mask_float = sitk.ReadImage('./brain_T1_mask.nii.gz')
# it is possible to force reading an image with a specific pixel/voxel type
image_brainT1_mask_uc = sitk.ReadImage('./brain_T1_mask.nii.gz', sitk.sitkUInt8)

itkw.view(image_brainT1, cmap='Grayscale', mode='z')<jupyter_output><empty_output><jupyter_text>It is possible to observe that for this case we also have a binary mask corresponding to the brain.Using the multiply filter we can obtain an image with just the brain.<jupyter_code>brain_image = sitk.Multiply(image_brainT1, image_brainT1_mask_float)# image_atlas * image_atlas_mask
itkw.view(brain_image, cmap='Grayscale', mode='z')<jupyter_output><empty_output><jupyter_text>NOTE: show multiplication by using *  Thresholding Filters  Binary Threshold with lower and upper threshold[Documentation SimpleITK](https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1BinaryThresholdImageFilter.html)<jupyter_code>binary_mask = sitk.BinaryThreshold(brain_image, lowerThreshold=300, upperThreshold=600, insideValue=1, outsideValue=0)
itkw.view(binary_mask, cmap='Grayscale', mode='z')<jupyter_output><empty_output><jupyter_text>Interactively discretize your image using Otsu Multiple Threshold method[Documentation SimpleITK](https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1OtsuMultipleThresholdsImageFilter.html)In this example you can see how to use SimpleITK and ipywidgets interactively change parameters of a filter (Otsu Multiple Thresholds) and check the results immediately!!<jupyter_code>from ipywidgets import interactive
import ipywidgets as widgets
viewer_int = itkw.view(brain_image, cmap='Grayscale', mode='z', annotations=False)

# Create an itk smoother filter object. By re-using the object, output memory-reallocation is avoided
otsu_filter = sitk.OtsuMultipleThresholdsImageFilter() # [brain_image, itk.Image[itk.F,3]].New(brain_image)
otsu_filter.SetNumberOfHistogramBins(64)

def otsu_and_view(thresholds=2):
    otsu_filter.SetNumberOfThresholds(thresholds)
    # Execute and Update the image used in the viewer
    viewer_int.image = otsu_filter.Execute(brain_image)
slider = interactive(otsu_and_view, thresholds=(1, 5, 1))

widgets.VBox([viewer_int, slider])<jupyter_output><empty_output><jupyter_text>Region Growing Filters[Documentation SimpleITK](https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1ConfidenceConnectedImageFilter.html)These are a type filters that require a seed point to perform the segmentation. We will try to segment the white matter using the Confidence Connected algorithm but other implementations like the Connected Threshold and the Neighborhood Connected.<jupyter_code>confidence_connected_filter = sitk.ConfidenceConnectedImageFilter()
confidence_connected_filter.SetMultiplier(1.80)
confidence_connected_filter.SetSeed([30, 61, 77])
confidence_connected_filter.SetNumberOfIterations( 10 );
confidence_connected_filter.SetReplaceValue( 255 );
confidence_connected_filter.SetInitialNeighborhoodRadius( 3 );
white_matter_image = confidence_connected_filter.Execute(brain_image)

itkw.view(white_matter_image, cmap='Grayscale', mode='z')<jupyter_output><empty_output><jupyter_text>This segmentation is far from perfect. If we go back to the image, it is possible to observe that this MRI image was affected by bias field inhomogeneity. This artifact affects the performance of segmentation filters.  Filtering MRI bias field inhomogeneity [Image Source](http://digitool.library.mcgill.ca/webclient/StreamGate?folder_id=0&dvs=1579874409466~728) N4 Bias Field Correction Image Filter[Documentation SimpleITK](https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1N4BiasFieldCorrectionImageFilter.html)<jupyter_code>bias_field_filter = sitk.N4BiasFieldCorrectionImageFilter()
bias_field_filter.SetNumberOfControlPoints([4,4,4])
bias_corrected_image = bias_field_filter.Execute(sitk.Cast(brain_image,sitk.sitkFloat32), image_brainT1_mask_uc)

itkw.checkerboard(bias_corrected_image, brain_image, cmap='Grayscale', mode='z', pattern=4)<jupyter_output><empty_output><jupyter_text>Let's now try to segment the white matter again using the same filter and parameters<jupyter_code>confidence_connected_filter = sitk.ConfidenceConnectedImageFilter()
# confidence_connected_filter.SetInput(bias_corrected_image)
confidence_connected_filter.SetSeed([30, 61, 77])
confidence_connected_filter.SetMultiplier( 1.80 );
confidence_connected_filter.SetNumberOfIterations( 10 );
confidence_connected_filter.SetReplaceValue( 255 );
confidence_connected_filter.SetInitialNeighborhoodRadius( 3 );
white_matter_mask = confidence_connected_filter.Execute(bias_corrected_image)

itkw.view(white_matter_mask, cmap='Grayscale', mode='z')<jupyter_output><empty_output><jupyter_text>Checkboard of Brain Image  White Matter Segmentation<jupyter_code>white_matter_mask_float = sitk.Cast(white_matter_mask, sitk.sitkFloat32) * 5 ## * 5 -> bring the mask to an intensity ~ of the image
itkw.checkerboard(white_matter_mask_float, image_brainT1, cmap='Grayscale', mode='z', pattern=6)<jupyter_output><empty_output><jupyter_text>More Advanced Segmentation Filters Bayesian Classifier Image Filter SegmentationThis filterThis filter does not exist in SimpleITK, only in ITK. Using the numpy (python scientific computing package providing powerful N-dimensional array object) interface of both libraries we will recreate the bias field corrected simple itk images as an itk image.The conversion of SimpleITK image to and ITK image will be your first exercise.  Exercise 1Convert the bias_corrected_image, which is a SimpleITK image object to ITK image object called bias_corrected_itk_image. Tips (functions required):- sitk.GetArrayFromImage() - converts the SimpleITK image into a numpy array (physical properties of image like spacing, origin and direction are lost)- itk.GetImageFromArray() - converts the numpy array into a SimpleITK image (physical properties of image like spacing, origin and direction are set to default values)- bias_corrected_image.GetOrigin() - get SimpleITK image world coordinates origin vector- bias_corrected_image.GetSpacing() - get SimpleITK image spacing vector- bias_corrected_image.GetDirection() - get SimpleITK image direction vector- bias_corrected_itk_image.SetOrigin() - set the correct world coordinates origin using the values prodived in the vector- bias_corrected_itk_image.SetSpacing() - set the correct spacing using the values prodived in the vector- Set direction is a little bit trickier so here is the example of the code that you can use:```python The interface for the direction is a little trickiernp_dir_vnl = itk.vnl_matrix_from_array(np.array(original_direction).reshape(3,3))DirectionType = type(bias_corrected_itk_image.GetDirection())direction = DirectionType(np_dir_vnl)bias_corrected_itk_image.SetDirection(direction)```<jupyter_code>## use this code box to write your solution




















bias_corrected_np_image = sitk.GetArrayFromImage(bias_corrected_image)
original_spacing = bias_corrected_image.GetSpacing()
original_origin = bias_corrected_image.GetOrigin()
original_direction = bias_corrected_image.GetDirection()
bias_corrected_itk_image = itk.GetImageFromArray(bias_corrected_np_image)
bias_corrected_itk_image.SetSpacing(original_spacing)
bias_corrected_itk_image.SetOrigin(original_origin)

# The interface for the direction is a little trickier
np_dir_vnl = itk.vnl_matrix_from_array(np.array(original_direction).reshape(3,3))
DirectionType = type(bias_corrected_itk_image.GetDirection())
direction = DirectionType(np_dir_vnl)
bias_corrected_itk_image.SetDirection(direction)

itkw.view(bias_corrected_itk_image, cmap='Grayscale', mode='z')
# The goal of this filter is to generate a membership image that indicates the membership of each pixel to each class. 
# These membership images are fed as input to the Bayesian classifier filter.
# BayesianClassifierInitializationImageFilter runs K-means on the image to determine k Gaussian density functions centered 
# around 'n' pixel intensity values. This is equivalent to generate Gaussian mixture model for the input image.
instance = itk.BayesianClassifierInitializationImageFilter[bias_corrected_itk_image, itk.F].New(bias_corrected_itk_image)
instance.SetNumberOfClasses(5)
instance.Update()
# The output to this filter is an itk::VectorImage that represents pixel memberships to 'n' classes. 
image_bayesian_class = instance.GetOutput() 

# Performs Bayesian Classification on an image using the membership Vector Image obtained by the itk.BayesianClassifierInitializationImageFilter.
bc = itk.BayesianClassifierImageFilter[image_bayesian_class, itk.SS, itk.F, itk.F].New(image_bayesian_class)
bc.Update()
labelled_bayesian = bc.GetOutput()
itkw.view(labelled_bayesian, cmap='Grayscale', mode='z')<jupyter_output><empty_output><jupyter_text>Gaussian Misture Models (GMMs) SegmentationTo execute this we will make use of the interface ITK  numpy and use scikit-learn use Gaussian Mixture Models (GMMs) to perform the segmentation.We will start by getting some statistics of the image using the segmentation labels.These will be used in the GMM.<jupyter_code>import sklearn
from sklearn.mixture import GaussianMixture
## Copy of itk.Image to numpy.ndarray
np_copy = itk.array_from_image(bias_corrected_itk_image)
np_mask_copy = itk.array_from_image(labelled_bayesian)
middle_slice = int(np.floor(np_copy.shape[0]/2))
gmm = GaussianMixture(n_components = 4)

## Fit the GMM on a single slice of the MRI data
data = np_copy[middle_slice]
gmm.fit(np.reshape(data, (data.size, 1)))

## Classify the all images according to the GMM
for j in range(2):
    im = np_copy[j]
    cls = gmm.predict(im.reshape((im.size, 1)))
    seg = np.zeros(np_copy[0].shape)
    seg = cls.reshape(np_copy[0].shape)
    np_mask_copy[j] = seg

## Copy of numpy.ndarray to itk.Image
mask_itk = itk.image_from_array(np_mask_copy)

## The conversion  itk -> numpy -> itk will change axis orientation. Correct it with the following filter
flipFilter = itk.FlipImageFilter[itk.Image.SS3].New()
flipFilter.SetInput(mask_itk)
flipAxes = (False, True, False)
flipFilter.SetFlipAxes(flipAxes)
flipFilter.Update()

corrected_mask = flipFilter.GetOutput()
itkw.view(corrected_mask, cmap='Grayscale', mode='z')<jupyter_output><empty_output><jupyter_text>Save mask for next notebook on Quantification<jupyter_code>write_filter = itk.ImageFileWriter[labelled_bayesian].New(labelled_bayesian)
write_filter.SetFileName('bayesian_mask.nii.gz')
write_filter.Update()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08353964324255697, 'fraction_numerical': 0.012124674665143147, 'mean_word_length': 7.86550365785031, 'pattern_counts': {'":': 0, '<': 76, '<?xml version=': 0, '>': 79, 'https://': 8, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '29881', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5340, 'n_tokens_neox': 4910, 'n_words': 1244, 'path': '02_Filtering_And_Segmentation_MIAW.ipynb', 'repo_name': 'JoaoSantinha/Medical_Image_Analysis_Workshop', 'size': '15753'}