<jupyter_start><jupyter_text>IntroHello! This rather quick and dirty kernel shows how to get started on segmenting nuclei using a neural network in Keras. The architecture used is the so-called [U-Net](https://arxiv.org/abs/1505.04597), which is very common for image segmentation problems such as this. I believe they also have a tendency to work quite well even on small datasets.Let's get started importing everything we need!<jupyter_code>import os
import sys
import random
import warnings

import numpy as np
import pandas as pd

import matplotlib.pyplot as plt

from tqdm import tqdm
from itertools import chain
from skimage.io import imread, imshow, imread_collection, concatenate_images
from skimage.transform import resize
from skimage.morphology import label

from keras.models import Model, load_model
from keras.layers import Input
from keras.layers.core import Lambda
from keras.layers.convolutional import Conv2D, Conv2DTranspose
from keras.layers.pooling import MaxPooling2D
from keras.layers.merge import concatenate
from keras.callbacks import EarlyStopping, ModelCheckpoint
from keras import backend as K

import tensorflow as tf

# Set some parameters
TRAIN_PATH = './stage1_train/'
TEST_PATH = './stage1_test/'
IMG_HEIGHT = 256
IMG_WIDTH = 256
IMG_CHANNELS = 3

warnings.filterwarnings('ignore', category=UserWarning, module='skimage')
# Get train and test IDs
train_ids = next(os.walk(TRAIN_PATH))[1]
test_ids = next(os.walk(TEST_PATH))[1]<jupyter_output><empty_output><jupyter_text>Get the dataLet's first import all the images and associated masks. I downsample both the training and test images to keep things light and manageable, but we need to keep a record of the original sizes of the test images to upsample our predicted masks and create correct run-length encodings later on. There are definitely better ways to handle this, but it works fine for now!  Data Generation<jupyter_code>from keras.utils import Sequence
import cv2
from sklearn.utils import shuffle
from skimage.transform import AffineTransform, warp
import copy

Y_train_label = []

class data_generator(Sequence):
    '''Generates data for Keras'''
    def __init__(self, ids, IMG_WIDTH, IMG_HEIGHT, IMG_CHANNELS, batch_size=10, training=True):
        '''Initialization'''
        self.ids = ids
        self.batch_size = batch_size
        self.training = training
        self.IMG_WIDTH = IMG_WIDTH
        self.IMG_HEIGHT = IMG_HEIGHT
        self.IMG_CHANNELS = IMG_CHANNELS
    def __len__(self):
        '''Denotes the number of batches per epoch'''
        return int(np.ceil(float(len(self.ids))/self.batch_size))
    def on_epoch_end(self):
        '''Updates indexes after each epoch'''
        if self.training: sel.fids= shuffle(self.ids)
    def __getitem__(self, i):
        l_bound =  i    * self.batch_size
        r_bound = (i+1) * self.batch_size
        if r_bound>len(self.data): # ensure every iteration has the same batch size
            r_bound = len(self.data)
            l_bound = r_bound - self.batch_size
        
        X = []
        Y = []
        for n, id_ in tqdm(enumerate(ids), total=len(ids)):
            path = TRAIN_PATH + id_
            img = imread(path + '/images/' + id_ + '.png')[:,:,:self.IMG_CHANNELS]
            img = resize(img, (self.IMG_HEIGHT, self.IMG_WIDTH), mode='constant', preserve_range=True) # skimage 0.13.1
#             img = resize(img, (IMG_HEIGHT, IMG_WIDTH), mode='constant', preserve_range=True, anti_aliasing=True) # wait for skimage 0.14.0
            img = color.rgb2hsv(img) # should I keep hue only? Maybe it can learn to invert purple images?
            mask_ = []
            for mask_file in next(os.walk(path + '/masks/'))[2]:
                mask_.append(imread(path + '/masks/' + mask_file))
                
            if self.training:
                # flip vertical
                if np.random.rand() < .5:
                    img = np.flip(img, 0)
                    for mask in mask_:
                        mask = np.flip(mask, 0)
                        
                # flip horizontal
                if np.random.rand() < .5:
                    img = np.flip(img, 1)
                    for mask in mask_:
                        mask = np.flip(mask, 1)
                        
                # rotation, shearing
                if np.random.rand() < 0.5:
                    y, x, _ = img.shape
                    h, w = y, x
                    img = np.pad(img, ((y//2, y//2), (x//2, x//2), (0, 0)), 'reflect')
                    for mask in mask_:
                        mask = np.pad(mask, ((y//2, y//2), (x//2, x//2), (0, 0)), 'reflect')
                    
                    rotT = np.random.uniform(-45,45)
                    shearT = np.random.uniform(-5,5)
                    translationX = np.random.uniform(-w,w) * .1
                    translationY = np.random.uniform(-h,h) * .1
                    atmtx = AffineTransform(rotation=rotT, shear=shearT, translation=(translationX, translationY))

                    img = warp(img, atmtx)
                    for mask in mask_:
                        mask = warp(mask, atmtx)

                    img = img[h//2:h//2+h,w//2:w//2+w,:]
                    for mask in mask_:
                        mask = mask[h//2:h//2+h,w//2:w//2+w,:]

            ### end of data augmentation ###
                        
            X.append(img)
            labeled_mask = np.zeros(img.shape[:2], dtype=np.uint8) # semantic segmentation ground truth
            dt = np.zeros(img.shape[:2], dtype=np.float32) # distance transform ground truth
            dir_grad = [np.zeros(img.shape[:2], dtype=np.float32)]
            dir_grad.append(np.zeros(img.shape[:2], dtype=np.float32))
            dir_grad = np.asarray(dir_grad) # direction ground truth (gradient of distance transform)
            for i, mask in enumerate(mask_):
                dt_   = -ndi.distance_transform_edt(mask).astype(np.float32) # basin
                dir_grad_y, dir_grad_x = np.array(mask).astype(int) * np.gradient(dt_)

        #         normalize direction
                dir_grad_ = np.sqrt(dir_grad_x * dir_grad_x + dir_grad_y * dir_grad_y)
                dir_grad_ = np.where(dir_grad_ < 1e-5, 1e-5, dir_grad_)
                dir_grad_x = dir_grad_x / dir_grad_
                dir_grad_y = dir_grad_y / dir_grad_
        #         print(dt_)
        #         sys.exit()

        #         superpositioning
                dir_grad[0] = np.where(dir_grad_x != 0, dir_grad_x, dir_grad[0])
                dir_grad[1] = np.where(dir_grad_y != 0, dir_grad_y, dir_grad[1])
                dt = np.where(dt_ != 0, dt_, dt)
                mask = np.where(mask > 0, i, 0)
                labeled_mask = labeled_mask | mask
        gt = np.maximum(np.where(dt == -1, 0, dt), -1) # ground truth
        Y_train_label.append([labeled_mask, dir_grad[0], dir_grad[1], dt, gt])
        
        
        return np.asarray(X), np.asarray(column(Y_train_label, 4))
    

# Get and resize test images
X_test = np.zeros((len(test_ids), IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS), dtype=np.uint8)
sizes_test = []
print('Getting and resizing test images ... ')
sys.stdout.flush()
for n, id_ in tqdm(enumerate(test_ids), total=len(test_ids)):
    path = TEST_PATH + id_
    img = imread(path + '/images/' + id_ + '.png')[:,:,:IMG_CHANNELS]
    sizes_test.append([img.shape[0], img.shape[1]])
    img = resize(img, (IMG_HEIGHT, IMG_WIDTH), mode='constant', preserve_range=True)
    X_test[n] = img<jupyter_output>Getting and resizing test images ...<jupyter_text>Let's see if things look all right by drawing some random images and their associated masks.<jupyter_code># # Check if training data looks all right
# ix = random.randint(0, len(train_ids))
# imshow(X_train[ix])
# plt.show()
# imshow(np.squeeze(Y_train[ix]))
# plt.show()<jupyter_output><empty_output><jupyter_text>Seems good! Create our Keras metricNow we try to define the *mean average precision at different intersection over union (IoU) thresholds* metric in Keras. TensorFlow has a mean IoU metric, but it doesn't have any native support for the mean over multiple thresholds, so I tried to implement this. **I'm by no means certain that this implementation is correct, though!** Any assistance in verifying this would be most welcome! *Update: This implementation is most definitely not correct due to the very large discrepancy between the results reported here and the LB results. It also seems to just increase over time no matter what when you train ... *<jupyter_code># Define IoU metric
def mean_iou_unet(y_true, y_pred):
    prec = []
    for t in np.arange(0.5, 1.0, 0.05):
        y_pred_ = tf.to_int32(y_pred > t)
        score, up_opt = tf.metrics.mean_iou(y_true, y_pred_, 2)
        K.get_session().run(tf.local_variables_initializer())
        with tf.control_dependencies([up_opt]):
            score = tf.identity(score)
        prec.append(score)
    return K.mean(K.stack(prec), axis=0)<jupyter_output><empty_output><jupyter_text>Build and train our neural networkNext we build our U-Net model, loosely based on [U-Net: Convolutional Networks for Biomedical Image Segmentation](https://arxiv.org/pdf/1505.04597.pdf) and very similar to [this repo](https://github.com/jocicmarko/ultrasound-nerve-segmentation) from the Kaggle Ultrasound Nerve Segmentation competition.![](https://lmb.informatik.uni-freiburg.de/people/ronneber/u-net/u-net-architecture.png)<jupyter_code># custom loss

def unet_custom(y_true, y_pred):
    y_true_ = y_true
    y_pred_ = y_pred
    return K.mean(K.binary_crossentropy(y_true_, y_pred_))    
# Build U-Net model
inputs = Input((IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS))
s = Lambda(lambda x: x / 255) (inputs)

c1 = Conv2D(8, (3, 3), activation='relu', padding='same') (s)
c1 = Conv2D(8, (3, 3), activation='relu', padding='same') (c1)
p1 = MaxPooling2D((2, 2)) (c1)

c2 = Conv2D(16, (3, 3), activation='relu', padding='same') (p1)
c2 = Conv2D(16, (3, 3), activation='relu', padding='same') (c2)
p2 = MaxPooling2D((2, 2)) (c2)

c3 = Conv2D(32, (3, 3), activation='relu', padding='same') (p2)
c3 = Conv2D(32, (3, 3), activation='relu', padding='same') (c3)
p3 = MaxPooling2D((2, 2)) (c3)

c4 = Conv2D(64, (3, 3), activation='relu', padding='same') (p3)
c4 = Conv2D(64, (3, 3), activation='relu', padding='same') (c4)
p4 = MaxPooling2D(pool_size=(2, 2)) (c4)

c5 = Conv2D(128, (3, 3), activation='relu', padding='same') (p4)
c5 = Conv2D(128, (3, 3), activation='relu', padding='same') (c5)

u6 = Conv2DTranspose(64, (2, 2), strides=(2, 2), padding='same') (c5)
u6 = concatenate([u6, c4])
c6 = Conv2D(64, (3, 3), activation='relu', padding='same') (u6)
c6 = Conv2D(64, (3, 3), activation='relu', padding='same') (c6)

u7 = Conv2DTranspose(32, (2, 2), strides=(2, 2), padding='same') (c6)
u7 = concatenate([u7, c3])
c7 = Conv2D(32, (3, 3), activation='relu', padding='same') (u7)
c7 = Conv2D(32, (3, 3), activation='relu', padding='same') (c7)

u8 = Conv2DTranspose(16, (2, 2), strides=(2, 2), padding='same') (c7)
u8 = concatenate([u8, c2])
c8 = Conv2D(16, (3, 3), activation='relu', padding='same') (u8)
c8 = Conv2D(16, (3, 3), activation='relu', padding='same') (c8)

u9 = Conv2DTranspose(8, (2, 2), strides=(2, 2), padding='same') (c8)
u9 = concatenate([u9, c1], axis=3)
c9 = Conv2D(8, (3, 3), activation='relu', padding='same') (u9)
c9 = Conv2D(8, (3, 3), activation='relu', padding='same') (c9)

outputs = Conv2D(1, (1, 1), activation='sigmoid') (c9)

model = Model(inputs=[inputs], outputs=[outputs])
model.compile(optimizer='adam', loss=unet_custom, metrics=[mean_iou_unet])
model.summary()<jupyter_output>__________________________________________________________________________________________________
Layer (type)                    Output Shape         Param #     Connected to                     
==================================================================================================
input_9 (InputLayer)            (None, 256, 256, 3)  0                                            
__________________________________________________________________________________________________
lambda_9 (Lambda)               (None, 256, 256, 3)  0           input_9[0][0]                    
__________________________________________________________________________________________________
conv2d_153 (Conv2D)             (None, 256, 256, 8)  224         lambda_9[0][0]                   
__________________________________________________________________________________________________
conv2d_154 (Conv2D)             (None, 256, 256, 8)  584         conv2d_153[0][0]                 
__________[...]<jupyter_text>Next we fit the model on the training data, using a validation split of 0.1. We use a small batch size because we have so little data. I recommend using checkpointing and early stopping when training your model. I won't do it here to make things a bit more reproducible (although it's very likely that your results will be different anyway). I'll just train for 10 epochs, which takes around 10 minutes in the Kaggle kernel with the current parameters. *Update: Added early stopping and checkpointing and increased to 30 epochs.*<jupyter_code># Fit model
earlystopper = EarlyStopping(patience=20, verbose=1)
checkpointer = ModelCheckpoint('model-dsbowl2018-unet.h5', verbose=1, save_best_only=True)

# partition
random.shuffle(train_ids)
train_ids_ = train_ids[:int(len(train_ids) * 0.9)]
validation_ids_ = train_ids[int(len(train_ids) * 0.9):]

# Parameters
params = {'IMG_WIDTH': IMG_WIDTH,
          'IMG_HEIGHT': IMG_HEIGHT,
          'IMG_CHANNELS': IMG_CHANNELS,
          'batch_size': 10,
          'training': True}

# Generators
training_generator = data_generator(train_ids_, **params)
params['training'] = False
validation_generator = data_generator(validation_ids_, **params)

results = model.fit_generator(generator=training_generator, use_multiprocessing=True, workers=8, 
                    callbacks=[earlystopper, checkpointer])<jupyter_output>Epoch 1/1<jupyter_text>All right, looks good! Loss seems to be a bit erratic, though. I'll leave it to you to improve the model architecture and parameters!  Make predictionsLet's make predictions both on the test set, the val set and the train set (as a sanity check). Remember to load the best saved model if you've used early stopping and checkpointing.<jupyter_code># Predict on train, val and test
model = load_model('model-dsbowl2018-unet.h5', custom_objects={'mean_iou': mean_iou_unet})
# preds_train = model.predict(X[:int(X.shape[0]*0.9)], verbose=1)
# preds_val = model.predict(X[int(X.shape[0]*0.9):], verbose=1)
preds_test = model.predict(X_test, verbose=1)

# Threshold predictions
# preds_train_t = (preds_train > 0.5).astype(np.uint8)
# preds_val_t = (preds_val > 0.5).astype(np.uint8)
preds_test_t = (preds_test > 0.5).astype(np.uint8)

# Create list of upsampled test masks
preds_test_upsampled = []
for i in range(len(preds_test)):
    preds_test_upsampled.append(resize(np.squeeze(preds_test[i]), 
                                       (sizes_test[i][0], sizes_test[i][1]), 
                                       mode='constant', preserve_range=True))
# # Perform a sanity check on some random training samples
# ix = random.randint(0, len(preds_train_t))
# imshow(X_train[ix])
# plt.show()
# imshow(np.squeeze(Y_train[ix]))
# plt.show()
# imshow(np.squeeze(preds_train_t[ix]))
# plt.show()<jupyter_output><empty_output><jupyter_text>The model is at least able to fit to the training data! Certainly a lot of room for improvement even here, but a decent start. How about the validation data?<jupyter_code># # Perform a sanity check on some random validation samples
# ix = random.randint(0, len(preds_val_t))
# imshow(X_train[int(X_train.shape[0]*0.9):][ix])
# plt.show()
# imshow(np.squeeze(Y_train[int(Y_train.shape[0]*0.9):][ix]))
# plt.show()
# imshow(np.squeeze(preds_val_t[ix]))
# plt.show()<jupyter_output><empty_output><jupyter_text>Not too shabby! Definitely needs some more training and tweaking. Encode and submit our resultsNow it's time to submit our results. I've stolen [this](https://www.kaggle.com/rakhlin/fast-run-length-encoding-python) excellent implementation of run-length encoding.<jupyter_code># Run-length encoding stolen from https://www.kaggle.com/rakhlin/fast-run-length-encoding-python
def rle_encoding(x):
    dots = np.where(x.T.flatten() == 1)[0]
    run_lengths = []
    prev = -2
    for b in dots:
        if (b>prev+1): run_lengths.extend((b + 1, 0))
        run_lengths[-1] += 1
        prev = b
    return run_lengths

def prob_to_rles(x, cutoff=0.5):
    lab_img = label(x > cutoff)
    for i in range(1, lab_img.max() + 1):
        yield rle_encoding(lab_img == i)<jupyter_output><empty_output><jupyter_text>Let's iterate over the test IDs and generate run-length encodings for each seperate mask identified by skimage ...<jupyter_code>new_test_ids = []
rles = []
for n, id_ in enumerate(test_ids):
    rle = list(prob_to_rles(preds_test_upsampled[n]))
    rles.extend(rle)
    new_test_ids.extend([id_] * len(rle))<jupyter_output><empty_output><jupyter_text>... and then finally create our submission!<jupyter_code># Create submission DataFrame
sub = pd.DataFrame()
sub['ImageId'] = new_test_ids
sub['EncodedPixels'] = pd.Series(rles).apply(lambda x: ' '.join(str(y) for y in x))
sub.to_csv('sub-dsbowl2018-1.csv', index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12856078408895513, 'fraction_numerical': 0.02659483005674328, 'mean_word_length': 3.2649718895135664, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 50, 'https://': 6, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '614830', 'license': 'Apache-2.0', 'n_tokens_mistral': 6018, 'n_tokens_neox': 5642, 'n_words': 1523, 'path': 'Nucli/dwt.ipynb', 'repo_name': 'Superdanby/Kaggle', 'size': '17447'}