Main

import numpy as np
import matplotlib.pyplot as plt
from skimage import color 
import skimage.io
from skimage import morphology
from skimage import segmentation
from skimage import measure
from skimage.filters import sobel
from skimage.exposure import histogram
from skimage.color import label2rgb
from scipy import ndimage as ndi
from skimage import metrics
from skimage.metrics import (adapted_rand_error,
                              variation_of_information,
                             hausdorff_distance
                            )
from skimage.feature import canny
from skimage.util import img_as_float

%matplotlib inline

# Iput Image
image = skimage.io.imread('MRI.png')
img = color.rgb2gray(image)

# Edge enhancement
edge_sobel = sobel(img)
img1 = edge_sobel

# Histogram of MRI image
hist,hist_centers = histogram(img)
fig, axes = plt.subplots(1,2, figsize =(8, 3))
axes[0].imshow(img, cmap=plt.cm.gray)
axes[0].axis('off')
axes[1].plot(hist_centers, hist, lw=2)
axes[1].set_title('histogram of gray values')

# Initial Mask
mask_i = skimage.io.imread('InitialMask.png')

# GT
gt = skimage.io.imread('GT.png')
gt1 =sobel(gt)

# Histogram of the ground truth image
hist,hist_centers = histogram(gt)
fig, axes = plt.subplots(1,2, figsize =(8, 3))
axes[0].imshow(gt, cmap=plt.cm.gray)
axes[0].axis('off')
axes[1].plot(hist_centers, hist, lw=2)
axes[1].set_title('histogram of gray values')

# Redefine foreground and background
# Define the range of gray values
markers = np.zeros_like(img)
markers[img < 30] = 1
markers[img > 50] = 2

fig, ax = plt.subplots(figsize=(4, 3))
ax.imshow(markers, cmap=plt.cm.nipy_spectral)
ax.set_title('markers')
ax.axis('off')

# Watershed Segmentation of gray scale image
segmentation_img = segmentation.watershed(img1, markers)

fig, ax = plt.subplots(figsize=(4, 3))
ax.imshow(segmentation_img, cmap=plt.cm.gray)
ax.set_title('segmentation')
ax.axis('off')

# Draw the contour and color
segmentation_img = ndi.binary_fill_holes(segmentation_img - 1)
labeled_img, _ = ndi.label(segmentation_img)
image_label_overlay = label2rgb(labeled_img, image=img, bg_label=0)

fig, axes = plt.subplots(1, 2, figsize=(8, 3), sharey=True)
axes[0].imshow(img, cmap=plt.cm.gray)
axes[0].contour(segmentation_img, [0.5], linewidths=1.2, colors='y')
axes[1].imshow(image_label_overlay)

for a in axes:
    a.axis('off')

plt.tight_layout()

plt.show()

# Redefine foreground and background based on the ground truth
# Define the range of gray values
markers1 = np.zeros_like(img)
markers1[gt < 30] = 1
markers1[gt > 50] = 2

fig, ax = plt.subplots(figsize=(4, 3))
ax.imshow(markers1, cmap=plt.cm.nipy_spectral)
ax.set_title('markers')
ax.axis('off')

# Watershed based on the ground truth mask
segmentation_img1 = segmentation.watershed(img1, markers1)

fig, ax = plt.subplots(figsize=(4, 3))
ax.imshow(segmentation_img1, cmap=plt.cm.gray)
ax.set_title('segmentation')
ax.axis('off')

# Image based on ground truth 
im_true = segmentation_img1

# Watershed test
im_test1 = segmentation.watershed(img1, markers, compactness=0.0001)

fig, ax = plt.subplots(figsize=(4, 3))
ax.imshow(im_test1, cmap=plt.cm.gray)
ax.set_title('segmentation')
ax.axis('off')

# Hausdorff Distance between ground truth and image after watershed
hd1 = hausdorff_distance(gt1, im_true)
print(hd1)

# Hausdorff distance between watershed images after and before implemented ground truth
hd2 = hausdorff_distance(segmentation_img1, segmentation_img)
print(hd2)

image = img_as_float(gt1)
gradient = segmentation.inverse_gaussian_gradient(image)
init_ls = np.zeros(image.shape, dtype=np.int8)
init_ls[10:-10, 10:-10] = 1
im_test2 = segmentation.morphological_geodesic_active_contour(gradient, iterations=1000,
                                                 init_level_set=init_ls,
                                                 smoothing=4, balloon=-1,
                                                 threshold=0.8)
im_test2 = measure.label(im_test2)

fig, ax = plt.subplots(figsize=(4, 3))
ax.imshow(im_test2, cmap=plt.cm.gray)
ax.set_title('segmentation')
ax.axis('off')

method_names = ['Compact watershed',
                'Morphological Geodesic Active Contours']
short_method_names = ['Compact WS', 'GAC']

precision_list = []
recall_list = []
split_list = []
merge_list = []
for name, im_test in zip(method_names, [im_test1, im_test2]):
    error, precision, recall = adapted_rand_error(im_true, im_test)
    splits, merges = variation_of_information(im_true, im_test)
    split_list.append(splits)
    merge_list.append(merges)
    precision_list.append(precision)
    recall_list.append(recall)
    print(f"\n## Method: {name}")
    print(f"Adapted Rand error: {error}")
    print(f"Adapted Rand precision: {precision}")
    print(f"Adapted Rand recall: {recall}")
    print(f"False Splits: {splits}")
    print(f"False Merges: {merges}")

fig, axes = plt.subplots(2, 3, figsize=(9, 6), constrained_layout=True)
ax = axes.ravel()

ax[0].scatter(merge_list, split_list)
for i, txt in enumerate(short_method_names):
    ax[0].annotate(txt, (merge_list[i], split_list[i]),
                   verticalalignment='center')
ax[0].set_xlabel('False Merges (bits)')
ax[0].set_ylabel('False Splits (bits)')
ax[0].set_title('Split Variation of Information')

ax[1].scatter(precision_list, recall_list)
for i, txt in enumerate(short_method_names):
    ax[1].annotate(txt, (precision_list[i], recall_list[i]),
                   verticalalignment='center')
ax[1].set_xlabel('Precision')
ax[1].set_ylabel('Recall')
ax[1].set_title('Adapted Rand precision vs. recall')
ax[1].set_xlim(0, 1)
ax[1].set_ylim(0, 1)

ax[2].imshow(segmentation.mark_boundaries(img, im_true))
ax[2].set_title('True Segmentation')
ax[2].set_axis_off()

ax[3].imshow(segmentation.mark_boundaries(img, im_test1))
ax[3].set_title('Compact Watershed')
ax[3].set_axis_off()

ax[4].imshow(segmentation.mark_boundaries(img, im_test2))
ax[4].set_title('Morphological GAC')
ax[4].set_axis_off()

plt.show()

# Hausdorff distance between GAC images and watershed imaged after implemented ground truth
hd3 = hausdorff_distance(im_true, im_test2)
print(hd3)