import matplotlib.pyplot as plt
import numpy as np
from skimage.morphology import disk
from PIL import Image


# Read the original image

f = np.array(Image.open('hw62g.png').convert('L'))

def myzeropadding(f,m,n):
    
    [M,N] = f.shape
    
    a = int((m-1)/2)
    b = int((n-1)/2)
    
    padded_f = np.zeros((M+2*a,N+2*b))
    
    padded_f[a:M+a,b:N+b] = f
    
    return padded_f

Morphological image processing is a collection of nonlinear operations used to extract image components that are useful in the representation and description of region shape such as boundaries, skeletons, and the convex hull.

Erosion

# function for erosion

def myerosion(f,SE):
    
    [M,N] = f.shape
    [m,n] = SE.shape
    a = int((m-1)/2)
    b = int((n-1)/2)
    
    erroded_f = np.zeros((M+2*a,N+2*b))
    
    padded_f = myzeropadding(f,m,n)
    
    for i in range(a,M+a):
        for j in range(b,N+b):
            
            neighbourhood = padded_f[i-a:i+a+1,j-b:j+b+1]
            
            erroded_f[i,j] = int(np.amin(neighbourhood))
            
            
    erroded_f = erroded_f[a:M+a,b:N+b]
    
    return erroded_f

# structuring element

SE = disk(10)

erroded_f = myerosion(f,SE)

fig = plt.figure(figsize=(20,10))

ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

ax1.imshow(f, cmap='gray',vmin=0, vmax=255)
ax1.set_title('Original Image',fontsize=20)

ax2.imshow(erroded_f, cmap='gray',vmin=0, vmax=255)
ax2.set_title('Image after erosion',fontsize=20)

plt.show()

Dialation

# function for dilation

def mydialation(f,SE):
    
    [M,N] = f.shape
    [m,n] = SE.shape
    a = int((m-1)/2)
    b = int((n-1)/2)
    
    dialated_f = np.zeros((M+2*a,N+2*b))
    
    padded_f = myzeropadding(f,m,n)
    
    for i in range(a,M+a):
        for j in range(b,N+b):
            
            neighbourhood = padded_f[i-a:i+a+1,j-b:j+b+1]
            
            dialated_f[i,j] = np.amax(neighbourhood)
            
            
    dialated_f = dialated_f[a:M+a,b:N+b]
    
    return dialated_f

# structuring element

SE = disk(10)

dialted_f = mydialation(f,SE)

fig = plt.figure(figsize=(20,10))

ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

ax1.imshow(f, cmap='gray',vmin=0, vmax=255)
ax1.set_title('Original Image',fontsize=20)

ax2.imshow(dialted_f, cmap='gray',vmin=0, vmax=255)
ax2.set_title('Image after dialation',fontsize=20)

plt.show()

Opening

# function for opening

def myopening(f,SE):
    
    erroded_f = myerosion(f,SE)
    open_f = mydialation(erroded_f,SE)
    
    return open_f

# structuring element

SE = disk(10)


open_f = myopening(f,SE)

fig = plt.figure(figsize=(20,10))

ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

ax1.imshow(f, cmap='gray',vmin=0, vmax=255)
ax1.set_title('Original Image',fontsize=20)

ax2.imshow(open_f, cmap='gray',vmin=0, vmax=255)
ax2.set_title('Image after opening',fontsize=20)

plt.show()

Closing

# function for closing

def myclosing(f,SE):
    
    dialted_f = mydialation(f,SE)
    close_f = myerosion(dialted_f,SE)
    
    return close_f

# structuring element

SE = disk(10)


close_f = myclosing(f,SE)

fig = plt.figure(figsize=(20,10))

ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

ax1.imshow(f, cmap='gray',vmin=0, vmax=255)
ax1.set_title('Original Image',fontsize=20)

ax2.imshow(close_f, cmap='gray',vmin=0, vmax=255)
ax2.set_title('Image after closing',fontsize=20)

plt.show()

Boundary Extraction

# function for boundary extraction


def myboundaryextraction(f,SE):
    
    erroded_f = myerosion(f,SE)
    boundary_extracted_f = f - erroded_f
    
    return boundary_extracted_f

# structuring element

SE = disk(10)

boundary_extracted_f = myboundaryextraction(f,SE)

fig = plt.figure(figsize=(20,10))

ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

ax1.imshow(f, cmap='gray',vmin=0, vmax=255)
ax1.set_title('Original Image',fontsize=20)

ax2.imshow(boundary_extracted_f, cmap='gray',vmin=0, vmax=255)
ax2.set_title('Boundary Extraction',fontsize=20)

plt.show()