DEV Community: Sai Rajesh Vanimireddy

Brain Tumour Image Segmentation using U-Net Architecture.

Sai Rajesh Vanimireddy — Wed, 08 Jan 2020 22:05:21 +0000

Download the sample MRI dataset from my github . The input images files are in .jpg format and groundtruth is in .xml format. There are missing files to form a complete dataset for segmentation process.

Dataset:

The purpose of the dataset is to provide the research community with a resource to advance the state-of-the-art in image detection, segmentation, and classification as well as help evaluating shortcomings of existing methods. For a labeled region, we provide the location in terms of bounding boxes, classifications. We also provide a detailed description of our annotation pipeline. The results show that some methods achieve excellent detection precision and good transcription accuracy

Image Processing

Image processing is any form of signal processing for which the input is an image, such as a photograph or video frame, and the output of image processing may be either an image or a set of characteristics or parameters related to the image. Most image processing techniques involve treating the image as a two-dimensional signal and applying standard signal-processing techniques to it.

Image Segmentation

Image segmentation is the process of partitioning an image into multiple segments. Image segmentation is typically used to locate objects and boundaries in images. presents the segmenting result of a femur image. It shows the outer surface (red), the surface between compact bone and spongy bone (green) and the surface of the bone marrow (blue). The testing applied an example of image segmentation to demonstrate the PSO method to find the best clusters of image segmentation. The results showed that PSO runs 170% faster when it used GPU in a parallel mode other than that used CPU alone, for the number of particles 100. This speedup is growing as the number of particles gets higher.

Directory commands are !ls, !mkdir.

Getting Started:

Importing Packages that we need to execute


import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import tensorflow as tf
import re
import glob
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import ast
from PIL import Image, ImageDraw

Plotting Function:

the image is defined and arranged in rows and coloumns

plt.subplots()
def image_show(image, nrows=1, ncols=1, cmap='gray'):
    fig, ax = plt.subplots(nrows=nrows, ncols=ncols, figsize=(14, 14))
    ax.imshow(image, cmap='gray')
    ax.axis('off')
    return fig, ax

Function to read in Marks

def read_xml(filename):
    with open(filename, 'r') as f:
        lines = f.readlines()
    return lines

Get a list of files in Thyroid with Segmentatin Coordinates and create list of Image files, segments and Marks

filenames = glob.glob("thyroid/*")
images = [x for x in filenames if x.endswith(".jpg")]
segments = [x for x in filenames if x.endswith(".xml")]

Remove images in which the data does not make sense

pattern = r'"points"(.*?)"annotation"'
segments = segments[segments["mark"].apply(lambda x: len(re.findall(pattern, x))<=2)]

Prepare data for test plot

from ast import literal_eval

# test = segments.loc[segments['img_id']=="thyroid/88_1.jpg"]
test = segments.iloc[0,:]
l = test['mark']#.get_values()[0]
# l = literal_eval(test['mark'])[0]

i = 0
t = 0

test_1 = re.findall(r'\d+', l)

dims = len(test_1)//2

temp = np.empty([dims, 2])
temp

while i < dims:
    temp[i] = (test_1[t], test_1[t+1])
    i = i+1
    t = t+2

Create Plot

fig, ax = image_show(Image.open('thyroid/197_1.jpg'), cmap='gray')
ax.plot(temp[:, 0], temp[:, 1], '.r',lw=3)

Break up Marks for the segmented images

pattern = r'"points"(.*?)"annotation"'
segments["mark_1"] = segments["mark"].apply(lambda x: re.findall(pattern, x)[0])
segments["mark_2"] = segments["mark"].apply(lambda x: re.findall(pattern, x)[1] if len(re.findall(pattern, x)) > 1 else "")

Create temp dataframes and concat

temp_1 = segments[["img_id","mark_1"]].copy()
temp_2 = segments[["img_id","mark_2"]].copy()

temp_1=temp_1.rename(columns = {'mark_1':'mark'})
temp_2=temp_2.rename(columns = {'mark_2':'mark'})

frames = [temp_1, temp_2]

df_new = pd.concat(frames,ignore_index = True)
df_new.head()

Function to parse marks

def parse_mark(mark):
    i = 0
    t = 0
    test_1 = re.findall(r'\d+', mark)
    dims = len(test_1)//2
    temp = np.empty([dims, 2])

    while i < dims:
        temp[i] = (test_1[t], test_1[t+1])
        i = i+1
        t = t+2

    return temp

The final dataset is ready to contain images, masks images, masks inverted images, multi-segmented masks images, outlines images, overlay images, and single_segmented_masks images and it is composed of total 1965.

Unet Model

UNet was first designed especially for medical image segmentation. It showed such good results that it used in many other fields after. In this article, we'll talk about why and how UNet works.

The architecture looks like a ‘U’ which justifies its name. This architecture consists of three sections: The contraction, The bottleneck, and the expansion section. The contraction section is made of many contraction blocks. Each block takes an input that applies two 3X3 convolution layers followed by a 2X2 max pooling. The number of kernels or feature maps after each block doubles so that architecture can learn the complex structures effectively. The bottommost layer mediates between the contraction layer and the expansion layer. It uses two 3X3 CNN layers followed by a 2X2 up convolution layer. Similar to the contraction layer, it also consists of several expansion blocks. Each block passes the input to two 3X3 CNN layers followed by a 2X2 upsampling layer.

Unet Implementation

I implemented the UNet model using the Pytorch framework. You can check out the UNet module for my customized dataset.

Unet_Model

import os
import numpy as np

from skimage.io import imread, imshow, concatenate_images
from skimage.transform import resize
from skimage.morphology import label

import tensorflow as tf

from keras.models import Model, load_model
from keras.layers import Input, BatchNormalization, Activation, Dense, Dropout
from keras.layers.core import Lambda, RepeatVector, Reshape
from keras.layers.convolutional import Conv2D, Conv2DTranspose
from keras.layers.pooling import MaxPooling2D, GlobalMaxPool2D
from keras.layers.merge import concatenate, add
from keras.callbacks import EarlyStopping, ModelCheckpoint, ReduceLROnPlateau
from keras.optimizers import Adam
from keras.preprocessing.image import ImageDataGenerator, array_to_img, img_to_array, load_img
import os
import random
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
plt.style.use("ggplot")
%matplotlib inline

from tqdm import tqdm_notebook, tnrange
from itertools import chain
from sklearn.model_selection import train_test_split



def conv2d_block(input_tensor, n_filters, kernel_size = 3, batchnorm = True):
    """Function to add 2 convolutional layers with the parameters passed to it"""
    # first layer
    x = Conv2D(filters = n_filters, kernel_size = (kernel_size, kernel_size),\
              kernel_initializer = 'he_normal', padding = 'same')(input_tensor)
    if batchnorm:
        x = BatchNormalization()(x)
    x = Activation('relu')(x)

    # second layer
    x = Conv2D(filters = n_filters, kernel_size = (kernel_size, kernel_size),\
              kernel_initializer = 'he_normal', padding = 'same')(input_tensor)
    if batchnorm:
        x = BatchNormalization()(x)
    x = Activation('relu')(x)

    return x

def get_unet(input_img, n_filters = 8, dropout = 0.2, batchnorm = True):
    """Function to define the UNET Model"""
    # Contracting Path
    c1 = conv2d_block(input_img, n_filters * 1, kernel_size = 3, batchnorm = batchnorm)
    p1 = MaxPooling2D((2, 2))(c1)
    p1 = Dropout(dropout)(p1)

    c2 = conv2d_block(p1, n_filters * 2, kernel_size = 3, batchnorm = batchnorm)
    p2 = MaxPooling2D((2, 2))(c2)
    p2 = Dropout(dropout)(p2)

    c3 = conv2d_block(p2, n_filters * 4, kernel_size = 3, batchnorm = batchnorm)
    p3 = MaxPooling2D((2, 2))(c3)
    p3 = Dropout(dropout)(p3)

    c4 = conv2d_block(p3, n_filters * 8, kernel_size = 3, batchnorm = batchnorm)
    p4 = MaxPooling2D((2, 2))(c4)
    p4 = Dropout(dropout)(p4)

    c5 = conv2d_block(p4, n_filters = n_filters * 16, kernel_size = 3, batchnorm = batchnorm)

    # Expansive Path
    u6 = Conv2DTranspose(n_filters * 8, (3, 3), strides = (2, 2), padding = 'same')(c5)
    u6 = concatenate([u6, c4])
    u6 = Dropout(dropout)(u6)
    c6 = conv2d_block(u6, n_filters * 8, kernel_size = 3, batchnorm = batchnorm)

    u7 = Conv2DTranspose(n_filters * 4, (3, 3), strides = (2, 2), padding = 'same')(c6)
    u7 = concatenate([u7, c3])
    u7 = Dropout(dropout)(u7)
    c7 = conv2d_block(u7, n_filters * 4, kernel_size = 3, batchnorm = batchnorm)

    u8 = Conv2DTranspose(n_filters * 2, (3, 3), strides = (2, 2), padding = 'same')(c7)
    u8 = concatenate([u8, c2])
    u8 = Dropout(dropout)(u8)
    c8 = conv2d_block(u8, n_filters * 2, kernel_size = 3, batchnorm = batchnorm)

    u9 = Conv2DTranspose(n_filters * 1, (3, 3), strides = (2, 2), padding = 'same')(c8)
    u9 = concatenate([u9, c1])
    u9 = Dropout(dropout)(u9)
    c9 = conv2d_block(u9, n_filters * 1, kernel_size = 3, batchnorm = batchnorm)

    outputs = Conv2D(1, (1, 1), activation='sigmoid')(c9)
    model = Model(inputs=[input_img], outputs=[outputs])
    return model

Set Parameters(resize all images in height and width)

im_width = 128
im_height = 128
border = 5

Convert images & masks into arrays

for n, id_ in tqdm_notebook(enumerate(ids), total=len(ids)):
    # Load images
    img = load_img("./data/images/"+id_, grayscale=True)
    x_img = img_to_array(img)
    x_img = resize(x_img, (128, 128, 1), mode = 'constant', preserve_range = True)
    # Load masks
    mask_orig = img_to_array(load_img("./data/masks_inverted/"+id_, grayscale=True))
    mask = resize(mask_orig, (128, 128, 1), mode = 'constant', preserve_range = True)
    # Save images
    X[n] = x_img/255.0
    y[n] = mask/255.0

Split train and valid the dataset

X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.1, random_state=42)
#Calculate test size ratio
test_size = (X_valid.shape[0]/X_train.shape[0])
# Split train and test
X_train, X_test, y_train, y_test = train_test_split(X_train, y_train, test_size=test_size, random_state=42)
y_train.shape
y_train_plt = y_train.reshape(1571, 128, 128)
import matplotlib.pyplot as plt
plt.imshow(y_train_plt[179, :,:], cmap='gray')

The output has shown below

Plot Loss vs Epoch

plt.figure(figsize=(8, 8))
plt.title("Learning curve")
plt.plot(results.history["loss"], label="loss")
plt.plot(results.history["val_loss"], label="val_loss")
plt.plot( np.argmin(results.history["val_loss"]), np.min(results.history["val_loss"]), marker="x", color="r", label="best model")
plt.xlabel("Epochs")
plt.ylabel("log_loss")
plt.legend();

The output of learning curve is plotted loss vs epoch as shown below

Predictions on test dataset

ix = random.randint(0, len(preds_val))
print(ix)
plot_sample(X_test, y_test, preds_test, preds_test_t,ix=ix)
threshold =.4082
binarize = .1
intersection = np.logical_and(y_test[ix].squeeze() > binarize, preds_test[ix].squeeze() > threshold)
union = np.logical_or(y_test[ix].squeeze() > binarize, preds_test[ix].squeeze() > threshold)
iou=np.sum(intersection) / np.sum(union)
print('IOU:',iou)

Results:

The accuracy of the result is 89.4% from the test dataset.

results = model.fit(X_train, y_train, batch_size=32, epochs=50, callbacks=callbacks,\validation_data=(X_valid, y_valid))

References

https://www.sciencedirect.com/science/article/pii/B978012809362750011X

https://towardsdatascience.com/u-net-b229b32b4a71

https://www.sciencedirect.com/science/article/pii/B9780128155530000021

https://www.elsevier.com/books/soft-computing-based-medical-image-analysis/dey/978-0-12-813087-2

https://arxiv.org/abs/1505.04597

Brain Tumour Image Segmentation using U-Net Architecture.

Sai Rajesh Vanimireddy — Wed, 08 Jan 2020 21:32:11 +0000

Dataset:

Image Processing

Image Segmentation

Directory commands are !ls, !mkdir.

Getting Started:

Importing Packages that we need to execute


import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import tensorflow as tf
import re
import glob
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import ast
from PIL import Image, ImageDraw

Plotting Function:

the image is defined and arranged in rows and coloumns

plt.subplots()
def image_show(image, nrows=1, ncols=1, cmap='gray'):
    fig, ax = plt.subplots(nrows=nrows, ncols=ncols, figsize=(14, 14))
    ax.imshow(image, cmap='gray')
    ax.axis('off')
    return fig, ax

Function to read in Marks

def read_xml(filename):
    with open(filename, 'r') as f:
        lines = f.readlines()
    return lines

Get a list of files in Thyroid with Segmentatin Coordinates and create list of Image files, segments and Marks

filenames = glob.glob("thyroid/*")
images = [x for x in filenames if x.endswith(".jpg")]
segments = [x for x in filenames if x.endswith(".xml")]

Remove images in which the data does not make sense

pattern = r'"points"(.*?)"annotation"'
segments = segments[segments["mark"].apply(lambda x: len(re.findall(pattern, x))&lt;=2)]

Prepare data for test plot

from ast import literal_eval

# test = segments.loc[segments['img_id']=="thyroid/88_1.jpg"]
test = segments.iloc[0,:]
l = test['mark']#.get_values()[0]
# l = literal_eval(test['mark'])[0]

i = 0
t = 0

test_1 = re.findall(r'\d+', l)

dims = len(test_1)//2

temp = np.empty([dims, 2])
temp

while i &lt; dims:
    temp[i] = (test_1[t], test_1[t+1])
    i = i+1
    t = t+2

Create Plot

fig, ax = image_show(Image.open('thyroid/197_1.jpg'), cmap='gray')
ax.plot(temp[:, 0], temp[:, 1], '.r',lw=3)

Break up Marks for the segmented images

pattern = r'"points"(.*?)"annotation"'
segments["mark_1"] = segments["mark"].apply(lambda x: re.findall(pattern, x)[0])
segments["mark_2"] = segments["mark"].apply(lambda x: re.findall(pattern, x)[1] if len(re.findall(pattern, x)) &gt; 1 else "")

Create temp dataframes and concat

temp_1 = segments[["img_id","mark_1"]].copy()
temp_2 = segments[["img_id","mark_2"]].copy()

temp_1=temp_1.rename(columns = {'mark_1':'mark'})
temp_2=temp_2.rename(columns = {'mark_2':'mark'})

frames = [temp_1, temp_2]

df_new = pd.concat(frames,ignore_index = True)
df_new.head()

Function to parse marks

def parse_mark(mark):
    i = 0
    t = 0
    test_1 = re.findall(r'\d+', mark)
    dims = len(test_1)//2
    temp = np.empty([dims, 2])

    while i &lt; dims:
        temp[i] = (test_1[t], test_1[t+1])
        i = i+1
        t = t+2

    return temp

Unet Model

Unet Implementation

I implemented the UNet model using the Pytorch framework. You can check out the UNet module for my customized dataset.

Unet_Model

import os
import numpy as np

from skimage.io import imread, imshow, concatenate_images
from skimage.transform import resize
from skimage.morphology import label

import tensorflow as tf

from keras.models import Model, load_model
from keras.layers import Input, BatchNormalization, Activation, Dense, Dropout
from keras.layers.core import Lambda, RepeatVector, Reshape
from keras.layers.convolutional import Conv2D, Conv2DTranspose
from keras.layers.pooling import MaxPooling2D, GlobalMaxPool2D
from keras.layers.merge import concatenate, add
from keras.callbacks import EarlyStopping, ModelCheckpoint, ReduceLROnPlateau
from keras.optimizers import Adam
from keras.preprocessing.image import ImageDataGenerator, array_to_img, img_to_array, load_img
import os
import random
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
plt.style.use("ggplot")
%matplotlib inline

from tqdm import tqdm_notebook, tnrange
from itertools import chain
from sklearn.model_selection import train_test_split



def conv2d_block(input_tensor, n_filters, kernel_size = 3, batchnorm = True):
    """Function to add 2 convolutional layers with the parameters passed to it"""
    # first layer
    x = Conv2D(filters = n_filters, kernel_size = (kernel_size, kernel_size),\
              kernel_initializer = 'he_normal', padding = 'same')(input_tensor)
    if batchnorm:
        x = BatchNormalization()(x)
    x = Activation('relu')(x)

    # second layer
    x = Conv2D(filters = n_filters, kernel_size = (kernel_size, kernel_size),\
              kernel_initializer = 'he_normal', padding = 'same')(input_tensor)
    if batchnorm:
        x = BatchNormalization()(x)
    x = Activation('relu')(x)

    return x

def get_unet(input_img, n_filters = 8, dropout = 0.2, batchnorm = True):
    """Function to define the UNET Model"""
    # Contracting Path
    c1 = conv2d_block(input_img, n_filters * 1, kernel_size = 3, batchnorm = batchnorm)
    p1 = MaxPooling2D((2, 2))(c1)
    p1 = Dropout(dropout)(p1)

    c2 = conv2d_block(p1, n_filters * 2, kernel_size = 3, batchnorm = batchnorm)
    p2 = MaxPooling2D((2, 2))(c2)
    p2 = Dropout(dropout)(p2)

    c3 = conv2d_block(p2, n_filters * 4, kernel_size = 3, batchnorm = batchnorm)
    p3 = MaxPooling2D((2, 2))(c3)
    p3 = Dropout(dropout)(p3)

    c4 = conv2d_block(p3, n_filters * 8, kernel_size = 3, batchnorm = batchnorm)
    p4 = MaxPooling2D((2, 2))(c4)
    p4 = Dropout(dropout)(p4)

    c5 = conv2d_block(p4, n_filters = n_filters * 16, kernel_size = 3, batchnorm = batchnorm)

    # Expansive Path
    u6 = Conv2DTranspose(n_filters * 8, (3, 3), strides = (2, 2), padding = 'same')(c5)
    u6 = concatenate([u6, c4])
    u6 = Dropout(dropout)(u6)
    c6 = conv2d_block(u6, n_filters * 8, kernel_size = 3, batchnorm = batchnorm)

    u7 = Conv2DTranspose(n_filters * 4, (3, 3), strides = (2, 2), padding = 'same')(c6)
    u7 = concatenate([u7, c3])
    u7 = Dropout(dropout)(u7)
    c7 = conv2d_block(u7, n_filters * 4, kernel_size = 3, batchnorm = batchnorm)

    u8 = Conv2DTranspose(n_filters * 2, (3, 3), strides = (2, 2), padding = 'same')(c7)
    u8 = concatenate([u8, c2])
    u8 = Dropout(dropout)(u8)
    c8 = conv2d_block(u8, n_filters * 2, kernel_size = 3, batchnorm = batchnorm)

    u9 = Conv2DTranspose(n_filters * 1, (3, 3), strides = (2, 2), padding = 'same')(c8)
    u9 = concatenate([u9, c1])
    u9 = Dropout(dropout)(u9)
    c9 = conv2d_block(u9, n_filters * 1, kernel_size = 3, batchnorm = batchnorm)

    outputs = Conv2D(1, (1, 1), activation='sigmoid')(c9)
    model = Model(inputs=[input_img], outputs=[outputs])
    return model

Set Parameters(resize all images in height and width)

im_width = 128
im_height = 128
border = 5

Convert images & masks into arrays

for n, id_ in tqdm_notebook(enumerate(ids), total=len(ids)):
    # Load images
    img = load_img("./data/images/"+id_, grayscale=True)
    x_img = img_to_array(img)
    x_img = resize(x_img, (128, 128, 1), mode = 'constant', preserve_range = True)
    # Load masks
    mask_orig = img_to_array(load_img("./data/masks_inverted/"+id_, grayscale=True))
    mask = resize(mask_orig, (128, 128, 1), mode = 'constant', preserve_range = True)
    # Save images
    X[n] = x_img/255.0
    y[n] = mask/255.0

Split train and valid the dataset

X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.1, random_state=42)
#Calculate test size ratio
test_size = (X_valid.shape[0]/X_train.shape[0])
# Split train and test
X_train, X_test, y_train, y_test = train_test_split(X_train, y_train, test_size=test_size, random_state=42)
y_train.shape
y_train_plt = y_train.reshape(1571, 128, 128)
import matplotlib.pyplot as plt
plt.imshow(y_train_plt[179, :,:], cmap='gray')

The output has shown below

Plot Loss vs Epoch

plt.figure(figsize=(8, 8))
plt.title("Learning curve")
plt.plot(results.history["loss"], label="loss")
plt.plot(results.history["val_loss"], label="val_loss")
plt.plot( np.argmin(results.history["val_loss"]), np.min(results.history["val_loss"]), marker="x", color="r", label="best model")
plt.xlabel("Epochs")
plt.ylabel("log_loss")
plt.legend();

The output of learning curve is plotted loss vs epoch as shown below

Predictions on test dataset

ix = random.randint(0, len(preds_val))
print(ix)
plot_sample(X_test, y_test, preds_test, preds_test_t,ix=ix)
threshold =.4082
binarize = .1
intersection = np.logical_and(y_test[ix].squeeze() &gt; binarize, preds_test[ix].squeeze() &gt; threshold)
union = np.logical_or(y_test[ix].squeeze() &gt; binarize, preds_test[ix].squeeze() &gt; threshold)
iou=np.sum(intersection) / np.sum(union)
print('IOU:',iou)

Plot Accuracy vs Epoch

plt.figure(figsize=(8, 8))
acc = results.history['acc']
val_acc = results.history['val_acc']
epochs = range(len(acc))
plt.plot(epochs, acc, 'r', label='Training acc')
plt.plot(epochs, val_acc, 'b', label='Validation acc')
plt.title('Training and validation accuracy')
plt.legend()
plt.figure()

Results

results = model.fit(X_train, y_train, batch_size=32, epochs=50, callbacks=callbacks,\validation_data=(X_valid, y_valid))

Conclusion

Image segmentation is an important problem and every day some new research papers are published. UNet contributed significantly to such research. Many new architectures are inspired by UNet. But still, there is so much to explore. There are so many variants of this architecture in the industry and hence it is necessary to understand the first one to understand them better. So if you have any doubts please comment below.

References

https://towardsdatascience.com/u-net-b229b32b4a71

https://www.sciencedirect.com/book/9780128155530/intelligent-data-analysis-for-biomedical-applications

https://www.sciencedirect.com/science/article/pii/B978012809362750011X

https://arxiv.org/abs/1601.07140

https://www.elsevier.com/books/soft-computing-based-medical-image-analysis/dey/978-0-12-813087-2

https://www.sciencedirect.com/science/article/pii/S2213158214000394

https://github.com/sairajes

Brain Tumour Image Segmentation using U-Net Architecture.

Sai Rajesh Vanimireddy — Wed, 08 Jan 2020 20:18:29 +0000

Dataset:

Image Processing

Image Segmentation

Getting Started:

Directory commands are !ls./folder name/, !mkdir.

Importing Packages that we need to execute

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import tensorflow as tf
import re
import glob
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import ast
from PIL import Image, ImageDraw

Plotting Function:

the image is defined and arranged in rows and coloumns

plt.subplots()
def image_show(image, nrows=1, ncols=1, cmap='gray'):
    fig, ax = plt.subplots(nrows=nrows, ncols=ncols, figsize=(14, 14))
    ax.imshow(image, cmap='gray')
    ax.axis('off')
    return fig, ax

Function to read in Marks

def read_xml(filename):
    with open(filename, 'r') as f:
        lines = f.readlines()
    return lines

Get a list of files in Thyroid with Segmentatin Coordinates and create list of Image files, segments and Marks

filenames = glob.glob("thyroid/*")
images = [x for x in filenames if x.endswith(".jpg")]
segments = [x for x in filenames if x.endswith(".xml")]
marks = [read_xml(file) for file in segments]

Remove images in which the data does not make sense

pattern = r'"points"(.*?)"annotation"'
segments = segments[segments["mark"].apply(lambda x: len(re.findall(pattern, x))<=2)]

Prepare data for test plot

from ast import literal_eval

# test = segments.loc[segments['img_id']=="thyroid/88_1.jpg"]
test = segments.iloc[0,:]
l = test['mark']#.get_values()[0]
# l = literal_eval(test['mark'])[0]

i = 0
t = 0

test_1 = re.findall(r'\d+', l)

dims = len(test_1)//2

temp = np.empty([dims, 2])
temp

while i < dims:
    temp[i] = (test_1[t], test_1[t+1])
    i = i+1
    t = t+2

Create Plot

fig, ax = image_show(Image.open('thyroid/197_1.jpg'), cmap='gray')
ax.plot(temp[:, 0], temp[:, 1], '.r',lw=3)

Break up Marks for the segmented images

pattern = r'"points"(.*?)"annotation"'
segments["mark_1"] = segments["mark"].apply(lambda x: re.findall(pattern, x)[0])
segments["mark_2"] = segments["mark"].apply(lambda x: re.findall(pattern, x)[1] if len(re.findall(pattern, x)) > 1 else "")

Create temp dataframes and concat

temp_1 = segments[["img_id","mark_1"]].copy()
temp_2 = segments[["img_id","mark_2"]].copy()

temp_1=temp_1.rename(columns = {'mark_1':'mark'})
temp_2=temp_2.rename(columns = {'mark_2':'mark'})

frames = [temp_1, temp_2]

df_new = pd.concat(frames,ignore_index = True)
df_new.head()

Function to parse marks

def parse_mark(mark):
    i = 0
    t = 0
    test_1 = re.findall(r'\d+', mark)
    dims = len(test_1)//2
    temp = np.empty([dims, 2])

    while i < dims:
        temp[i] = (test_1[t], test_1[t+1])
        i = i+1
        t = t+2

    return temp

Unet Model

Unet Implementation

I have worked on the the UNet model using the Pytorch framework. You can check out the UNet module for my customized dataset.

Unet_Model

import os
import numpy as np

from skimage.io import imread, imshow, concatenate_images
from skimage.transform import resize
from skimage.morphology import label

import tensorflow as tf

from keras.models import Model, load_model
from keras.layers import Input, BatchNormalization, Activation, Dense, Dropout
from keras.layers.core import Lambda, RepeatVector, Reshape
from keras.layers.convolutional import Conv2D, Conv2DTranspose
from keras.layers.pooling import MaxPooling2D, GlobalMaxPool2D
from keras.layers.merge import concatenate, add
from keras.callbacks import EarlyStopping, ModelCheckpoint, ReduceLROnPlateau
from keras.optimizers import Adam
from keras.preprocessing.image import ImageDataGenerator, array_to_img, img_to_array, load_img
import os
import random
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
plt.style.use("ggplot")
%matplotlib inline

from tqdm import tqdm_notebook, tnrange
from itertools import chain
from sklearn.model_selection import train_test_split



def conv2d_block(input_tensor, n_filters, kernel_size = 3, batchnorm = True):
    """Function to add 2 convolutional layers with the parameters passed to it"""
    # first layer
    x = Conv2D(filters = n_filters, kernel_size = (kernel_size, kernel_size),\
              kernel_initializer = 'he_normal', padding = 'same')(input_tensor)
    if batchnorm:
        x = BatchNormalization()(x)
    x = Activation('relu')(x)

    # second layer
    x = Conv2D(filters = n_filters, kernel_size = (kernel_size, kernel_size),\
              kernel_initializer = 'he_normal', padding = 'same')(input_tensor)
    if batchnorm:
        x = BatchNormalization()(x)
    x = Activation('relu')(x)

    return x

def get_unet(input_img, n_filters = 8, dropout = 0.2, batchnorm = True):
    """Function to define the UNET Model"""
    # Contracting Path
    c1 = conv2d_block(input_img, n_filters * 1, kernel_size = 3, batchnorm = batchnorm)
    p1 = MaxPooling2D((2, 2))(c1)
    p1 = Dropout(dropout)(p1)

    c2 = conv2d_block(p1, n_filters * 2, kernel_size = 3, batchnorm = batchnorm)
    p2 = MaxPooling2D((2, 2))(c2)
    p2 = Dropout(dropout)(p2)

    c3 = conv2d_block(p2, n_filters * 4, kernel_size = 3, batchnorm = batchnorm)
    p3 = MaxPooling2D((2, 2))(c3)
    p3 = Dropout(dropout)(p3)

    c4 = conv2d_block(p3, n_filters * 8, kernel_size = 3, batchnorm = batchnorm)
    p4 = MaxPooling2D((2, 2))(c4)
    p4 = Dropout(dropout)(p4)

    c5 = conv2d_block(p4, n_filters = n_filters * 16, kernel_size = 3, batchnorm = batchnorm)

    # Expansive Path
    u6 = Conv2DTranspose(n_filters * 8, (3, 3), strides = (2, 2), padding = 'same')(c5)
    u6 = concatenate([u6, c4])
    u6 = Dropout(dropout)(u6)
    c6 = conv2d_block(u6, n_filters * 8, kernel_size = 3, batchnorm = batchnorm)

    u7 = Conv2DTranspose(n_filters * 4, (3, 3), strides = (2, 2), padding = 'same')(c6)
    u7 = concatenate([u7, c3])
    u7 = Dropout(dropout)(u7)
    c7 = conv2d_block(u7, n_filters * 4, kernel_size = 3, batchnorm = batchnorm)

    u8 = Conv2DTranspose(n_filters * 2, (3, 3), strides = (2, 2), padding = 'same')(c7)
    u8 = concatenate([u8, c2])
    u8 = Dropout(dropout)(u8)
    c8 = conv2d_block(u8, n_filters * 2, kernel_size = 3, batchnorm = batchnorm)

    u9 = Conv2DTranspose(n_filters * 1, (3, 3), strides = (2, 2), padding = 'same')(c8)
    u9 = concatenate([u9, c1])
    u9 = Dropout(dropout)(u9)
    c9 = conv2d_block(u9, n_filters * 1, kernel_size = 3, batchnorm = batchnorm)

    outputs = Conv2D(1, (1, 1), activation='sigmoid')(c9)
    model = Model(inputs=[input_img], outputs=[outputs])
    return model

Set Parameters(resize all images in height and width)

im_width = 128
im_height = 128
border = 5

Convert images & masks into arrays

for n, id_ in tqdm_notebook(enumerate(ids), total=len(ids)):
    # Load images
    img = load_img("./data/images/"+id_, grayscale=True)
    x_img = img_to_array(img)
    x_img = resize(x_img, (128, 128, 1), mode = 'constant', preserve_range = True)
    # Load masks
    mask_orig = img_to_array(load_img("./data/masks_inverted/"+id_, grayscale=True))
    mask = resize(mask_orig, (128, 128, 1), mode = 'constant', preserve_range = True)
    # Save images
    X[n] = x_img/255.0
    y[n] = mask/255.0

Split train and valid the dataset

X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.1, random_state=42)
#Calculate test size ratio
test_size = (X_valid.shape[0]/X_train.shape[0])
# Split train and test
X_train, X_test, y_train, y_test = train_test_split(X_train, y_train, test_size=test_size, random_state=42)
y_train.shape
y_train_plt = y_train.reshape(1571, 128, 128)
import matplotlib.pyplot as plt
plt.imshow(y_train_plt[179, :,:], cmap='gray')

The output has shown below

Plot Loss vs Epoch

plt.figure(figsize=(8, 8))
plt.title("Learning curve")
plt.plot(results.history["loss"], label="loss")
plt.plot(results.history["val_loss"], label="val_loss")
plt.plot( np.argmin(results.history["val_loss"]), np.min(results.history["val_loss"]), marker="x", color="r", label="best model")
plt.xlabel("Epochs")
plt.ylabel("log_loss")
plt.legend();

The output of learning curve is plotted loss vs epoch as shown below

Predictions on test dataset

ix = random.randint(0, len(preds_val))
print(ix)
plot_sample(X_test, y_test, preds_test, preds_test_t,ix=ix)
threshold =.4082
binarize = .1
intersection = np.logical_and(y_test[ix].squeeze() > binarize, preds_test[ix].squeeze() > threshold)
union = np.logical_or(y_test[ix].squeeze() > binarize, preds_test[ix].squeeze() > threshold)
iou=np.sum(intersection) / np.sum(union)
print('IOU:',iou)

Now run the model and fit your dataset

results = model.fit(X_train, y_train, batch_size=32, epochs=50, callbacks=callbacks,\validation_data=(X_valid, y_valid))

inference

I have achieved the results 89.4% accuracy after 50 epochs output and it is valid to use for image segmentation process.

References

https://towardsdatascience.com/u-net-b229b32b4a71

https://www.sciencedirect.com/book/9780128093627/gpu-based-parallel-implementation-of-swarm-intelligence-algorithms

https://www.sciencedirect.com/book/9780128155530/intelligent-data-analysis-for-biomedical-applications

https://github.com/sairajes/unet_segmentation

Data Scientist Roles and Responsibilities(My Opinion)

Sai Rajesh Vanimireddy — Fri, 03 Jan 2020 01:30:04 +0000

There is a difference in Data Engineer, Data Analysist and Data Scientist. If you want to become a Data Scientist the main moto is to think, learn and analyze globally with your Neuron.It is possible when you work on as many various datasets as you can. Focus Of Role

Create Hypothesis:

Test

Analyse

Translate Using Clean Data

In order to get a job in an industry as a data scientist position, the following skills are needed.

Dealing with various data operations.

Involved with several data underlying procedures.

Can able to handle Structure and Unstructured Data

possesses knowledge of mathematics, statistics, and machine learning algorithms

Proficient in Python, R, SAS, Matlab, Hadoop and Spark

Data Preprocessing for MRI Datasets

Sai Rajesh Vanimireddy — Thu, 02 Jan 2020 05:14:42 +0000

A lot of effort in solving any machine learning problem goes into preparing the data. In this tutorial, we will see how to load and preprocess/augment data from a non-trivial dataset. The various image format classes give full or selective access to the header (meta) information and access to the image data is made available via NumPy arrays.

To run this tutorial, please make sure the following packages are installed:

scikit-image: For image io and transforms.

nibabel: Read / write access to some common neuroimaging file formats

matplotlib: which allows mathematical representations of all kinds to be created and visualized.

Parameter:

Path: A path-like object representing a file system path.

Return Type: This method returns a string value that represents the directory name from the specified path.

Download the sample nibabel files from my github which is in nii.giz form

Import Libraries:

import os
import nibabel as nib
import numpy as np
import matplotlib.pyplot as plt
}

Directory for MRI data
Importing Datasets: dirName = r"file_path"

function to read data from dirName

def FileRead(file_path):
    nii = nib.load(file_path)
    data = nii.get_data()
    #if (np.shape(np.shape(data))[0] == 3): #If channel data not present
     #   data = np.expand_dims(data, 3)
    return data

function returns a 2 dimensional numpy array from a 3 dimensional data

def Nifti3Dto2D(Nifti3D):
    Nifti3DWOChannel = Nifti3D#[:,:,:,0] #Considering there is only one channel info
    Nifti2D = Nifti3DWOChannel.reshape(np.shape(Nifti3DWOChannel)[0], np.shape(Nifti3DWOChannel)[1] * np.shape(Nifti3DWOChannel)[2])
    return Nifti2D

function returns a 1 dimensional numpy array from a 2 dimensional array

def Nifti2Dto1D(Nifti2D):
    Nifti1D = Nifti2D.reshape(np.shape(Nifti2D)[0] * np.shape(Nifti2D)[1])
    return Nifti1D

reshapes a 1 dimensional array to 2 dimensional data

def Nifti1Dto2D(Nifti1D, height):
    Nifti2D = Nifti1D.reshape(height,int(np.shape(Nifti1D)[0]/height))
    return Nifti2D

reshapes 2 dimensional array to 3 dimensional data

def Nifti2Dto3D(Nifti2D):
    Nifti3DWOChannel = Nifti2D.reshape(np.shape(Nifti2D)[0],np.shape(Nifti2D)[0],np.shape(Nifti2D)[1]//np.shape(Nifti2D)[0])
    #Nifti3D = np.expand_dims(Nifti3DWOChannel, 3)
    return Nifti3DWOChannel

normalize data between o-1

def normalize(x):
    min_val = np.min(x)
    max_val = np.max(x)
    x = (x-min_val) / (max_val-min_val)
    return x

To copy data to local storage

def FileSave(data, file_path):
    nii = nib.Nifti1Image(data, np.eye(4))
    nib.save(nii, file_path)

returns list of files from the given dirName

def getListOfFiles(dirName):
    # creating a list of file and sub directories 
    listOfFile = os.listdir(dirName)
    allFiles = []
# Iterate over entries
    for entry in listOfFile:
        # Creating full path
        fullPath = os.path.join(dirName, entry)
        if os.path.isdir(fullPath):
            allFiles = allFiles + getListOfFiles(fullPath)
        else:
            allFiles.append(fullPath)
    return allFiles

Takes list of files and creates 1 dimension and 2 dimension data from MR data

def main():
    listOfFiles = getListOfFiles(dirName)
    twod_data = []# to save two dimensional data
    oned_data = []
    listOfFiles = []
    for (dirpath, dirnames, filenames) in os.walk(dirName):
        listOfFiles += [os.path.join(dirpath, file) for file in filenames]
       # listOfFiles_under.append(listOfFiles_under)
        #print(listOfFiles_under)
    for element in listOfFiles:
        print(element)
        red_files = FileRead(element)
        twod_array = Nifti3Dto2D(red_files)
        twod_data.append(twod_array)
        #print(twod_unsam)
        oned_array= Nifti2Dto1D(twod_array)
        oned_data.append(oned_array)
    oned_data = np.asarray(oned_data, dtype=np.float64, order='C')
    twod_data = np.asarray(twod_data, dtype=np.float64, order='C')
    return oned_data, twod_data

oned_data, twod_data = main()

for visualization

data_vis = Nifti2Dto3D(Nifti1Dto2D(oned_data[0,:], 256))
#plt.imshow(data_vis[20,:,:], cmap='gray')
plt.imshow(data_vis[:,:,20], cmap='gray')

Output for one slice image

Datatypes in array object

Sai Rajesh Vanimireddy — Tue, 31 Dec 2019 20:57:46 +0000

1D Image:

One-dimensional pictures are those containing only one dimension. This is only possible when you're dealing with a line, as the only dimension you have is length, defined by a single array.When data objects are stored in an array, individual objects are selected by an index that is usually a non-negative scalar integer. Indexes are also called subscripts. An index maps the array value to a stored object. There are three ways in which the elements of an array can be indexed:

0 (zero-based indexing): The first element of the array is indexed by subscript of 0.

1 (one-based indexing): The first element of the array is indexed by subscript of 1.

n (n-based indexing): The base index of an array can be freely chosen. Usually programming languages allowing n-based indexing also allow negative index values and other scalar data types like enumerations, or characters may be used as an array index.

import numpy as np
a = np.array([1, 2, 3,4,5])
a
a.ndim
a.shape
len(a)
Out[69]: 5
}

2D Image:

The two dimensions depicted are length and width and the objects on the picture are flat. one of the examples of 2D images are gray scale images, a kind of black-and-white or gray monochrome.

import numpy as np
b = np.array([[0, 1, 2], [3, 4, 5]])    
b
b.ndim
b.shape
len(b)
Out[67]: 2

3D Image:

Three-dimensional pictures contain yet another dimension can be a color image. Information in 3-D coordinates can also be objects that have less than three spatial dimensions and in which the third coordinate describes another property, such as time, color or a color channel.

import numpy as np
c = np.array([[[1], [2]], [[3], [4]]])
c
c.shape
Out[70]: (2, 2, 1)

Brain Tumour, MRI Images and CT Images

Sai Rajesh Vanimireddy — Thu, 26 Dec 2019 07:59:47 +0000

Brain Tumor:

A brain tumor occurs when infrequent cells shape appears inside the cerebrum, which has two main types, namely, benign tumors and malignant tumors. Malignant tumors can be categorized into basic tumors, and secondary tumors that spread elsewhere. Automated brain tumor detection through MRI can offer a valued outlook and earlier accurate detection of the brain tumor. The tumor detection is performed in several stages, namely, enhancement, segmentation, classification. Depending on the tumorous tissue origin and its behavior, it is classified from less aggressive (benign: grades I and II) to aggressive (malignant: grades III and IV). The clinicians diagnose brain tumors using medical imaging techniques like magnetic resonance imaging(MRI) and computerized tomography (CT).

Magnetic Resonance Imaging:

Magnetic resonance imaging or MRI scanning uses magnetism, radio waves, and a computer to produce images of body structures. MRI imaging gives structural details of the brain and gives excellent results even during the early onset of disease. images. In general, the artifacts present in the MRI are partial volume, RF noise, intensity overlaps, motion, and gradient. The different MRI modalities are shown below. There are numerous challenges affecting automated or semiautomated brain tumor localization in MRI. The challenges are as follows.

The tumorous tissues may originate anywhere in the brain; therefore, their localization is very challenging.
All normal brain tissues have probable and prefixed locations. The tumor growth compresses normal tissues.
In MRI, image sequences like benigs(T1-weighted, T2-weighted) malignant(T1c, fluid-attenuated inversion recovery(FLAIR) )are generated. Each sequence gives different biological information about brain tissues. So using only some of these sequences may not guarantee accurate tumor detection.
The accuracy of segmentation also affects growth rate prediction, which is a very important issue for post-treatment planning.

CT Images:

CT imaging is preferred over MRI due to wider availability, lower cost, and sensitiveness to early stoke. Computed tomography (CT) and magnetic resonance imaging (MRI) are the two modalities that are regularly used for brain imaging. While MRI is the preferred clinical modality for imaging the anatomy of the brain, due to its lack of ionizing radiation and a higher contrast resolution1, CT plays an essential role in diagnostic imaging due to a much faster acquisition speed, which reduces movement artifacts, and a higher spatial resolution. In addition, CT can be used in the presence of metallic implants and tattoos with metallic ink, for which MRI is contraindicated due to the dangers imposed by the high magnetic fields required. CT is also better suited for the evaluation of penetrating brain injuries and other trauma and acute neurological emergencies.