<jupyter_start><jupyter_text>Brain Image Segmentation by Clustering  Data and imports<jupyter_code>import os
import numpy as np
import SimpleITK as sitk
import matplotlib.pyplot as plt
%matplotlib inline

data_dir = "path to data"
print(os.listdir(data_dir))<jupyter_output>['mri-brain-tissues.nii.gz', 'mri-brain-mask.nii.gz', 'ct-brain.nii.gz', 'spine-data.txt', 'ct-brain-mask.nii.gz', 'mri-brain.nii.gz', 'ct-brain-lesion.nii.gz']<jupyter_text>Set up the image viewer<jupyter_code>from ipywidgets import interact, fixed
from IPython.display import display
from utils.image_viewer import display_image<jupyter_output><empty_output><jupyter_text>Image segmentation via clustering<jupyter_code>img = sitk.ReadImage(data_dir + 'mri-brain.nii.gz')
msk = sitk.ReadImage(data_dir + 'mri-brain-mask.nii.gz')
print('MR image')
display_image(img, window=400, level=200)
print('Brain mask')
display_image(msk)<jupyter_output>MR image<jupyter_text>use the brain mask to mask out non-brain regions by setting them to zero<jupyter_code>img_array = sitk.GetArrayFromImage(img)
msk_array = sitk.GetArrayFromImage(msk)

masked_array = img_array
masked_array[msk_array==0] = 0

img_masked = sitk.GetImageFromArray(masked_array)
img_masked.CopyInformation(img)

print('Masked image')
display_image(img_masked, window=400, level=200)<jupyter_output>Masked image<jupyter_text>Data preparation for clustering<jupyter_code># Take all non-zero voxels and flatten the data into a 1D numpy array
X = img_array[msk_array > 0].flatten().reshape(-1, 1)

# Get the number of points
num_pts = len(X.flatten())

# Extract the minimum and maximum intensity values and calculate the number of bins for the histogram
lim_low = np.min(X).astype(np.int)
lim_high = np.max(X).astype(np.int)
num_bins = (lim_high - lim_low + 1)

plt.figure(figsize=(10, 4), dpi=100)
plt.hist(X, bins=num_bins, density=True, range=(lim_low, lim_high), color='lightgray');
plt.xlim([0,350]);
plt.show()

print('Number of points ' + str(num_pts))<jupyter_output><empty_output><jupyter_text>Random subsamplinguse a random subset of the original data for efficiency.<jupyter_code>sampling = 0.05
X_subset = np.random.choice(X.flatten(),int(num_pts*sampling)).reshape(-1, 1)

plt.figure(figsize=(10, 4), dpi=100)
plt.hist(X_subset, bins=num_bins, density=True, range=(lim_low, lim_high), color='lightgray');
plt.xlim([0,350]);
plt.show()

print('Number of points ' + str(len(X_subset)))<jupyter_output><empty_output><jupyter_text>Number of clustersAssume the brain consists of mostly three tissue types, grey matter (GM), white matter (WM), and cerebrospinal fluid (CSF).<jupyter_code>num_clusters = 3<jupyter_output><empty_output><jupyter_text>K-means<jupyter_code>import sklearn.cluster as cluster

# Create a k-means instance
kmeans = cluster.KMeans(n_clusters=num_clusters)

# Running k-means via the fit function
kmeans.fit(X_subset)

# Produce the clustering result for all image points
y = kmeans.predict(img_array.flatten().reshape(-1, 1))

# K-means will produce labels between 0 and (k-1), we want 0 to be background, so we shift the labels by one
y = y + 1 # shift labels
y[(msk_array == 0).flatten()] = 0 # zero background

# Construct a 3D label map
lab_array = y.reshape(img_array.shape).astype('uint8')
seg_kmeans = sitk.GetImageFromArray(lab_array)
seg_kmeans.CopyInformation(img)

# Display the results using SimpleITK mapping of label maps to colours
display_image(sitk.LabelToRGB(seg_kmeans))<jupyter_output><empty_output><jupyter_text>Gaussian Mixture Models<jupyter_code>import sklearn.mixture as mixture

# Create a GaussianMixture instance
gmm = mixture.GaussianMixture(n_components=num_clusters)
gmm.fit(X_subset)
y = gmm.predict(img_array.flatten().reshape(-1, 1))
y = y + 1 # shift labels
y[(msk_array == 0).flatten()] = 0 # zero background

lab_array = y.reshape(img_array.shape).astype('uint8')
seg_gmm = sitk.GetImageFromArray(lab_array)
seg_gmm.CopyInformation(img)

display_image(sitk.LabelToRGB(seg_gmm))<jupyter_output><empty_output><jupyter_text>Plotting the GMM result<jupyter_code>import matplotlib.mlab as mlab
from scipy.stats import norm

def plot_gmm(x, gmm):
    omega = gmm.weights_
    mu = gmm.means_
    sigma = np.sqrt(gmm.covariances_)
    for ind in range(0,omega.shape[0]): 
        plt.plot(x,omega[ind]*norm.pdf(x, mu[ind], sigma[ind]), linewidth=2, label='GMM Component '+str(ind))

plt.figure(figsize=(10, 4), dpi=100)
plt.hist(X, bins=num_bins, density=True, range=(lim_low, lim_high), label='Intensity histogram', color='lightgray');
x = np.linspace(lim_low,lim_high,num_bins).reshape(-1,1)
plot_gmm(x,gmm)
plt.plot(x,np.exp(gmm.score_samples(x)), linewidth=2, color='k', label='Gaussian Mixture Model')
plt.xlim([0,350])
plt.legend(loc=0, shadow=True, fontsize=12)<jupyter_output><empty_output><jupyter_text>Lesion segmentation using clustering<jupyter_code>img = sitk.ReadImage(data_dir + 'ct-brain.nii.gz')
msk = sitk.ReadImage(data_dir + 'ct-brain-mask.nii.gz')

print('CT image')
display_image(img, x=70, y=100, z=90, window=120, level=40)

print('Brain mask')
display_image(msk, x=70, y=100, z=90)
img_array = sitk.GetArrayFromImage(img)
msk_array = sitk.GetArrayFromImage(msk)

masked_array = img_array
masked_array[msk_array==0] = 0

img_masked = sitk.GetImageFromArray(masked_array)
img_masked.CopyInformation(img)

print('Masked image')
display_image(img_masked, x=70, y=100, z=90, window=120, level=40)
# Take all non-zero voxels and flatten the data into a 1D numpy array
X = img_array[msk_array > 0].flatten().reshape(-1, 1)

# Get the number of points
num_pts = len(X.flatten())

# Extract the minimum and maximum intensity values and calculate the number of bins for the histogram
lim_low = -20 # manually set intensity range of interest
lim_high = 100 # manually set intensity range of interest
num_bins = (lim_high - lim_low + 1)

plt.figure(figsize=(10, 4), dpi=100)
plt.hist(X, bins=num_bins, density=True, range=(lim_low,lim_high), color='lightgray');
plt.xlim([0,80]) # we limit the x-axis to the range of interest
plt.show()

print('Number of points ' + str(num_pts))
sampling = 0.05
X_subset = np.random.choice(X.flatten(),int(num_pts*sampling)).reshape(-1, 1)

plt.figure(figsize=(10, 4), dpi=100)
plt.hist(X_subset, bins=num_bins, density=True, range=(lim_low, lim_high), color='lightgray');
plt.xlim([0,80]);
plt.show()

print('Number of points ' + str(len(X_subset)))
num_clusters = 5
import sklearn.cluster as cluster

# Create a k-means instance
kmeans = cluster.KMeans(n_clusters=num_clusters)

# Running k-means via the fit function
kmeans.fit(X_subset)

# Produce the clustering result for all image points
y = kmeans.predict(img_array.flatten().reshape(-1, 1))

# K-means will produce labels between 0 and (k-1), we want 0 to be background, so we shift the labels by one
y = y + 1 # shift labels
y[(msk_array == 0).flatten()] = 0 # zero background

# Construct a 3D label map
lab_array = y.reshape(img_array.shape).astype('uint8')
seg_kmeans = sitk.GetImageFromArray(lab_array)
seg_kmeans.CopyInformation(img)

# Display the results using SimpleITK mapping of label maps to colours
display_image(sitk.LabelToRGB(seg_kmeans))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1290367483296214, 'fraction_numerical': 0.02185412026726058, 'mean_word_length': 6.587117212249208, 'pattern_counts': {'":': 0, '<': 42, '<?xml version=': 0, '>': 44, 'https://': 0, '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': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '89996', 'license': 'MIT', 'n_tokens_mistral': 2829, 'n_tokens_neox': 2621, 'n_words': 604, 'path': 'brain_image_segmentation_by_clustering.ipynb', 'repo_name': 'iyezhiyu/Medical_Image_Processing', 'size': '7184'}