<jupyter_start><jupyter_text>TP 2 : Computer Vision Part 1 : segmentationIn this part of the TP, we are going to implement two approaches to image segmentation :- algorithm of Otsu (histogram-based segmentation)- a simple region-growing based algorithmFirst, let us load some packages and define some helper functions (as in the previous TP)<jupyter_code>from matplotlib import pyplot as plt
import numpy as np 
import imageio
from skimage import color

is_colab = True

def read_image(file_name):
    img_color = imageio.imread(file_name)
    img_gray = color.rgb2gray(img_color)
    return img_gray,img_color
    
def write_image(img_in,file_name_out):
    imageio.imwrite(file_name_out, np.uint8(255.0*img_in))
    
def display_image(img_in):
    plt.figure(figsize=(10, 10))
    if (img_in.ndim == 2):
        plt.imshow(img_in,cmap='gray')
    elif (img_in.ndim == 3):
        # careful, in this case we supppose the pixel values are between 0 and 255
        plt.imshow(np.uint8(img_in))
    else:
        print('Error, unknown number of dimensions in image')
    return

file_dir = 'images/'
file_name = 'crossing'
file_ext = '.png'

if (is_colab == True):
  !wget "https://perso.telecom-paristech.fr/anewson/doc/images/crossing.png"
  img_gray,_ = read_image(file_name+file_ext)
else:
  img_gray,_ = read_image(file_dir+file_name+file_ext)

display_image(img_gray)
img_size = img_gray.shape
img_size<jupyter_output>--2021-01-31 17:17:17--  https://perso.telecom-paristech.fr/anewson/doc/images/crossing.png
Resolving perso.telecom-paristech.fr (perso.telecom-paristech.fr)... 172.16.58.3, fd00:c2b6:b24b:be67:2827:688d:e6a1:6a3b
Connecting to perso.telecom-paristech.fr (perso.telecom-paristech.fr)|172.16.58.3|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 638741 (624K) [image/png]
Saving to: ‘crossing.png’

crossing.png        100%[===================>] 623.77K  1.45MB/s    in 0.4s    

2021-01-31 17:17:18 (1.45 MB/s) - ‘crossing.png’ saved [638741/638741]<jupyter_text>Otsu's algorithmOtsu's method is based on finding a threshold to split the image into two separate segments, based simply on their grey-level values. Since the maximum number of grey-levels is finite (and small, 256), we can use an exhaustive search.The grey-level threshold is the one that minimises the following criterion :$argmin_{\tau} \omega_0(\tau)\sigma_0(\tau) + \omega_1(\tau)\sigma_1(t)$where- $\sigma_0(\tau)$ and $\sigma_1(\tau)$ are the variances of the pixels in the first and second regions- $\omega_0(\tau)$ and $\omega_1(\tau)$ are weights of the first and second regionsand - $\omega_0(\tau) = \sum_{p \in \Omega} \mathbb{1}_{\tau}(I(p)), \qquad$- $\omega_1(\tau) = \sum_{p \in \Omega} 1-\mathbb{1}_{\tau}(I(p))$and $\mathbb{1}_{\tau}(I(p)) =\begin{cases}        1 \; \text{if} \: I(p)<=\tau\\        0 \; \text{otherwise}\end{cases}$Implement this algorithm in a function called ``segmentation_otsu``. Note : only analyse thresholds which produce two distinct regions : you cannot have a region with no pixels. Furthemore, remember that the image is normalised between 0 and 1, with 256 discrete values. Choose the threshold values accordingly.<jupyter_code># This function returns the segmentation of the algorithm of otsu, as well as the threshold found
from skimage import exposure
def segmentation_otsu(image):

    # FILL IN CODE HERE
    counts, bins = np.histogram(image.ravel(), bins=256)
    bins = bins[:-1]

    weight1 = np.cumsum(counts)
    weight2 = np.cumsum(counts[::-1])[::-1]

    mean1 = np.cumsum(counts * bins) / weight1
    mean2 = (np.cumsum((counts * bins)[::-1]) / weight2[::-1])[::-1]

    # Understood via documentation that minimizing the within-class variance
    # is equivalent is equivalent to maximizing the between-class variance 
    variance12 = weight1[:-1] * weight2[1:] * (mean1[:-1] - mean2[1:]) ** 2

    idx = np.argmax(variance12)
    tau_out = bins[idx]

    seg_otstu_out = np.ones(shape=image.shape)

    for i in range(image.shape[0]):
      for j in range(image.shape[1]):
        if image[i,j] <= tau_out:
          seg_otstu_out[i,j] = 0

    return seg_otstu_out,tau_out<jupyter_output><empty_output><jupyter_text>Carry out the segmentation on the input image and display the result<jupyter_code>img_out_otsu,tau = segmentation_otsu(img_gray)
display_image(img_out_otsu)
write_image(img_out_otsu,file_name+'_segmentation_otsu.png')<jupyter_output><empty_output><jupyter_text>__Question__ What problem can you see with the segmentation result ? What morphological technique(s) (seen in the lesson) could you use to clean the result up ?  __Answer__ : It doesn't work well when there is noise and the distribution is not clear. We could simply denoise the image before. Carry out this technique to clean up the output image, and imwrite the result. You can find the necessary function(s) in the following package :skimage.morphologyChoose the neighbourhood parameter of the functions to obtain the best results.<jupyter_code>from skimage.morphology import opening, closing

struct_size = 7
struct_el = np.ones((struct_size,struct_size))
img_out_otsu_cleaned = closing(opening(img_out_otsu,struct_el),struct_el) # FILL IN CODE HERE

display_image(img_out_otsu_cleaned)
write_image(img_out_otsu_cleaned,file_name+'_segmentation_otsu_cleaned_1.png')<jupyter_output><empty_output><jupyter_text>__Question__ In which order did you propose to carry out the morphological operations ? Why ? __Answer__ : The opening and closing was performed on a binary image. It is quicker, and the objective being to connect black or white cracks, applying the filter on a binary image makes more sense in this case. Let's change the image, and carry out the segmentation again.<jupyter_code>file_dir = 'images/'
file_name = 'town'
file_ext = '.png'

if (is_colab == True):
  !wget "https://perso.telecom-paristech.fr/anewson/doc/images/town.png"
  img_gray,_ = read_image(file_name+file_ext)
else:
  img_gray,_ = read_image(file_dir+file_name+file_ext)

display_image(img_gray)
img_size = img_gray.shape
img_size
img_out_otsu,tau = segmentation_otsu(img_gray)
display_image(img_out_otsu)
write_image(img_out_otsu,file_name+'_segmentation_otsu.png')<jupyter_output>--2021-01-31 17:33:31--  https://perso.telecom-paristech.fr/anewson/doc/images/town.png
Resolving perso.telecom-paristech.fr (perso.telecom-paristech.fr)... 172.16.58.3, fd00:c2b6:b24b:be67:2827:688d:e6a1:6a3b
Connecting to perso.telecom-paristech.fr (perso.telecom-paristech.fr)|172.16.58.3|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 806256 (787K) [image/png]
Saving to: ‘town.png’

town.png            100%[===================>] 787.36K  1.46MB/s    in 0.5s    

2021-01-31 17:33:32 (1.46 MB/s) - ‘town.png’ saved [806256/806256]<jupyter_text>As we can see, there is no notion of connectivity in the algorithm of Otsu (two regions far away can be in the same segment). For such a notion, we turn to region merging.  Region mergingThe region merging algorithm is quite simple : we start with a small region, and add new pixels progressively to the region. More precisely, at each iteration of the algorithm, we attempt to add the neighbours of the pixels in the current segmentation. We add new pixels according to the criterion :- $I(p)$ is added to the region if the absolute difference between the pixel value and the region's average value is less than a threshold value (which you must set).First, create a function called ``initialize_segmentation`` which initialises the segmentation with a small circular region centred on a given seed point (the seed point is an input). For this, the following function :- scipy.ndimage.morphology.distance_transform_edtcan be useful.<jupyter_code>from scipy.ndimage.morphology import distance_transform_edt

def initialize_segmentation(seed_pt,img_shape):
    # FILL IN CODE HERE 
    n,m = seed_pt

    seg_init=np.ones(shape=img_shape)
    seg_init[n,m] = 0
    seg_init = distance_transform_edt(seg_init)
    seg_init[seg_init <= 10] = 1
    seg_init[seg_init > 10] = 0

    return seg_init<jupyter_output><empty_output><jupyter_text>Now, create a function called ``segmentation_region_growing`` which carries out the region growing algorithm, with the criterion above, and which uses the ``initialize_segmentation`` function. Again, the functions from ``scipy.ndimage.morphology`` can be useful here.__IMPORTANT NOTE__ Once a pixel has been added to the region, it is in the region permanently : it is not checked again in further iterations (otherwise we could remove pixels which were added).<jupyter_code>from scipy.ndimage import binary_dilation


def segmentation_region_growing(img,seed_pt,tau):
  seg_init = initialize_segmentation(seed_pt,img.shape)
  seg = seg_init
  hold = True
  while hold:
    hold = False
    mean_reg = np.mean(img[seg == 1])
    new_points_idx = np.where(binary_dilation(seg) - seg == 1)
    for idx in range(len(new_points_idx[0])):
      i = new_points_idx[0][idx]
      j = new_points_idx[1][idx]
      if (np.abs(img[i,j] - mean_reg) < tau):
        seg[i,j] = 1
        hold = True
  return seg,seg_init<jupyter_output><empty_output><jupyter_text>Carry out the segmentation on the input image with the seed point and threshold specified in the following code :<jupyter_code>seed_pt = np.asarray([130,350])
tau = 20.0 * (1.0/256.0)
img_out_region_growing,seg_init = segmentation_region_growing(img_gray,seed_pt,tau)
display_image(img_out_region_growing)
write_image(seg_init,file_name+'_segmentation_region_growing_init.png')
write_image(img_out_region_growing,file_name+'_segmentation_region_growing.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1069435833849969, 'fraction_numerical': 0.03265137425087828, 'mean_word_length': 4.809723889555823, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 34, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 4, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '419189', 'license': 'MIT', 'n_tokens_mistral': 3450, 'n_tokens_neox': 3173, 'n_words': 935, 'path': 'segmentation_for_students.ipynb', 'repo_name': 'cadmoska/Computer_Visions_TP2', 'size': '9646'}