{"cells":[{"metadata":{"collapsed":true,"trusted":true,"_uuid":"1eab3d459f095d4e32cd12d762f0e137a4572324"},"cell_type":"code","source":"# fork from https://github.com/jfpuget/DSB_2018 kaggler:CPMP\nfname = 'keras_67'","execution_count":1,"outputs":[]},{"metadata":{"_cell_guid":"c332549b-8d23-4bb5-8497-e7a8eb8b21d2","_uuid":"5c38504af3a84bee68c66d3cde74443c58df422f","trusted":false},"cell_type":"code","source":"import os\nos.environ[\"CUDA_VISIBLE_DEVICES\"]=\"1\"\n\nimport sys\nimport random\nimport warnings\n\nimport numpy as np\nimport pandas as pd\n\nimport matplotlib.pyplot as plt\n\nfrom tqdm import tqdm\nfrom itertools import chain\nfrom skimage.io import imread, imshow, imread_collection, concatenate_images\nfrom skimage.transform import resize\nfrom skimage.morphology import label, binary_erosion, binary_dilation, disk\nfrom skimage.morphology import square, watershed, closing, binary_closing\nfrom skimage.morphology import remove_small_holes, remove_small_objects\nfrom skimage.filters.rank import gradient\nfrom skimage.exposure import rescale_intensity\nfrom skimage.segmentation import random_walker\n\nfrom sklearn.model_selection import KFold\n\nfrom scipy.ndimage.morphology import binary_fill_holes\n\nfrom keras.models import Model, load_model\nfrom keras.layers import Input\nfrom keras.layers.core import Dropout, Lambda\nfrom keras.layers.convolutional import Conv2D, Conv2DTranspose\nfrom keras.layers.pooling import MaxPooling2D\nfrom keras.layers.merge import concatenate\nfrom keras.callbacks import EarlyStopping, ModelCheckpoint\nfrom keras import backend as K\nfrom keras.optimizers import Adam\n\nimport tensorflow as tf\n\nimport pickle as pkl\nimport gc\n\n# Set some parameters\nIMG_WIDTH = 256\nIMG_HEIGHT = 256\nIMG_CHANNELS = 3\nTRAIN_PATH = '../input/stage1_train/'\nTEST_PATH = '../input/stage1_test/'\n\nwarnings.filterwarnings('ignore', category=UserWarning, module='skimage')","execution_count":2,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"b5490387a85e6d1019c6fed05f83e327e66b729f"},"cell_type":"code","source":"import cv2","execution_count":3,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"eae7519ecd4efbf758b3dba9ee4c80d973f597c9"},"cell_type":"code","source":"import random as rn\ndef init_seeds(seed):\n    os.environ['PYTHONHASHSEED'] = '0'\n\n    # The below is necessary for starting Numpy generated random numbers\n    # in a well-defined initial state.\n\n    np.random.seed(seed)\n\n    # The below is necessary for starting core Python generated random numbers\n    # in a well-defined state.\n\n    rn.seed(seed)\n\n    # Force TensorFlow to use single thread.\n    # Multiple threads are a potential source of\n    # non-reproducible results.\n    # For further details, see: https://stackoverflow.com/questions/42022950/which-seeds-have-to-be-set-where-to-realize-100-reproducibility-of-training-res\n\n    session_conf = tf.ConfigProto(intra_op_parallelism_threads=1, inter_op_parallelism_threads=1)\n\n    from keras import backend as K\n\n    # The below tf.set_random_seed() will make random number generation\n    # in the TensorFlow backend have a well-defined initial state.\n    # For further details, see: https://www.tensorflow.org/api_docs/python/tf/set_random_seed\n\n    tf.set_random_seed(seed)\n\n    sess = tf.Session(graph=tf.get_default_graph(), config=session_conf)\n    K.set_session(sess)\n    return sess","execution_count":4,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"57405ebd7bb4045a183eeaedf2ad26f766f59a97"},"cell_type":"code","source":"# Get train and test IDs\ntrain_ids = next(os.walk(TRAIN_PATH))[1]\ntest_ids = next(os.walk(TEST_PATH))[1]\n\nnp.random.seed(0)\nperm = np.random.permutation(len(train_ids))\ntrain_ids = [train_ids[i] for i in perm]","execution_count":5,"outputs":[]},{"metadata":{"_cell_guid":"59c4a25d-645f-4b74-9c53-145ac78cc481","_uuid":"875af74f980236825de3a650825b46e25632422c"},"cell_type":"markdown","source":"# Get the data\nLet'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!"},{"metadata":{"trusted":false,"_uuid":"724228fb8c746453756086e5b99323ecd6093d4c"},"cell_type":"code","source":"train = []\ntrain_mask = []\ntest = []\n\nfor n, id_ in tqdm(enumerate(train_ids), total=len(train_ids)):\n    file = \"../input/stage1_train/{}/images/{}.png\".format(id_,id_)\n    mfile = \"../input/stage1_train/{}/masks/*.png\".format(id_)\n    image = cv2.imread(file)\n    image = rescale_intensity(image, out_range=np.uint8)\n    masks = imread_collection(mfile).concatenate()\n    train.append(image)\n    train_mask.append(masks)\n    \nfor n, id_ in tqdm(enumerate(test_ids), total=len(test_ids)):\n    file = \"../input/stage1_test/{}/images/{}.png\".format(id_,id_)\n    image = cv2.imread(file)\n    image = rescale_intensity(image, out_range=np.uint8)\n    test.append((image))\n    ","execution_count":6,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"4e3a6d32ae299ed44058adb2fe0175a37b65a60c"},"cell_type":"code","source":"def to_flip(img_rgb):\n    # do not flip colored images\n    if (img_rgb[:,:,0] != img_rgb[:,:,1]).any():\n        return img_rgb\n    #green channel happends to produce slightly better results\n    #than the grayscale image and other channels\n    img_gray=img_rgb[:,:,1]#cv2.cvtColor(img_rgb, cv2.COLOR_BGR2GRAY)\n    #morphological opening (size tuned on training data)\n    circle7=cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(7,7))\n    img_open=cv2.morphologyEx(img_gray, cv2.MORPH_OPEN, circle7)\n    #Otsu thresholding\n    img_th=cv2.threshold(img_open,0,255,cv2.THRESH_OTSU)[1]\n    #Invert the image in case the objects of interest are in the dark side\n    if (np.sum(img_th==255)>np.sum(img_th==0)):\n        return img_rgb\n    else:\n        return 255 - img_rgb","execution_count":7,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"65b5af1be82dd954e836e1d08120ad60326c0a32"},"cell_type":"code","source":"train = [to_flip(img_rgb) for img_rgb in train]\ntest = [to_flip(img_rgb) for img_rgb in test]","execution_count":8,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"f589e82c27cdfa6b9e0e77d2e8f4eafe44d590e0"},"cell_type":"code","source":"def split_aux(img):\n    height = img.shape[0]\n    width = img.shape[1]\n    if height > 2*width:\n        half = int(height//2)\n        return [img[:half, :, :], img[half:, :, :]]\n    elif height > width:\n        return [img[:width, :, :], img[height-width:, :, :]]\n    elif width > 2*height:\n        half = int(width//2)\n        return [img[:, :half, :], img[:, half:, :]]\n    else:\n        return [img[:, :height, :], img[:, width-height:, :]]\n\ndef split(img):\n    s = split_aux(img)\n    return s\n\ndef split_mask_aux(img):\n    height = img.shape[1]\n    width = img.shape[2]\n    if height > 2*width:\n        half = int(height//2)\n        return [img[:, :half, :], img[:, half:, :]]\n    elif height > width:\n        return [img[:, :width, :], img[:, height-width:, :]]\n    elif width > 2*height:\n        half = int(width//2)\n        return [img[:, :, :half], img[:, :, half:]]\n    else:\n        return [img[:, :, :height], img[:, :, width-height:]]\n\ndef split_mask(img):\n    s = split_mask_aux(img)\n    return s","execution_count":9,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"2b9e33b56a4de52a05a46320e1bb554c8fe894a6"},"cell_type":"code","source":"train_split = [split(img) for img in train]\ntrain_split = [t_split[i] for t_split in train_split for i in [0, 1] ]\n\ntrain_mask_split = [split_mask(img) for img in train_mask]\ntrain_mask_split = [t_split[i] for t_split in train_mask_split for i in [0, 1] ]\n\ntest_split = [split(img) for img in test]\ntest_split = [t_split[i] for t_split in test_split for i in [0, 1] ]","execution_count":10,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"0fa48836d355243b8f9e087a8115e6cc2a6a31f6"},"cell_type":"code","source":"# Get and resize train images and masks\nX_train = np.zeros((len(train_split) * 4, IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS), np.uint8)\nY_train = np.zeros((len(train_split) * 4, IMG_HEIGHT, IMG_WIDTH, 1), np.uint8)\nZ_train = np.zeros((len(train_split) * 4, IMG_HEIGHT, IMG_WIDTH, 3), np.uint8)\nprint('Getting and resizing train images and masks ... ')\nsys.stdout.flush()\nfor n, (img, masks) in enumerate(zip(tqdm(train_split), train_mask_split)):\n    img = img[:, :, :IMG_CHANNELS]\n    img = resize(img, (IMG_HEIGHT, IMG_WIDTH), mode='constant', preserve_range=True)\n    img_mean = np.mean(img, axis=2).astype(np.uint8)\n    for c in range(IMG_CHANNELS):\n        img[:,:,c] = img_mean\n    X_train[n * 4 + 0] = img\n    X_train[n * 4 + 1] = np.fliplr(img)\n    X_train[n * 4 + 2] = np.flipud(img)\n    X_train[n * 4 + 3] = np.flipud(np.fliplr(img))\n\n    mask = np.zeros((IMG_HEIGHT, IMG_WIDTH), np.uint8)\n    mask_lr = np.zeros((IMG_HEIGHT, IMG_WIDTH), np.uint8)\n    mask_ud = np.zeros((IMG_HEIGHT, IMG_WIDTH), np.uint8)\n    mask_lr_ud = np.zeros((IMG_HEIGHT, IMG_WIDTH), np.uint8)\n    for mask_id in range(masks.shape[0]):\n        mask_ = masks[mask_id, :, :]\n        mask_ = mask_ // 255\n        mask_ = resize(mask_, (IMG_HEIGHT, IMG_WIDTH), mode='constant', \n                       preserve_range=True).astype(np.uint8)\n        mask = np.maximum(mask, mask_)\n        mask_lr = np.maximum(mask_lr, np.fliplr(mask_))\n        mask_ud = np.maximum(mask_ud, np.flipud(mask_))\n        mask_lr_ud = np.maximum(mask_lr_ud, np.flipud(np.fliplr(mask_)))\n        \n    Y_train[4*n + 0, :, :, 0] = mask\n    Y_train[4*n + 1, :, :, 0] = mask_lr\n    Y_train[4*n + 2, :, :, 0] = mask_ud\n    Y_train[4*n + 3, :, :, 0] = mask_lr_ud\n  \n    mask = np.zeros((IMG_HEIGHT, IMG_WIDTH), np.uint8)\n    mask_lr = np.zeros((IMG_HEIGHT, IMG_WIDTH), np.uint8)\n    mask_ud = np.zeros((IMG_HEIGHT, IMG_WIDTH), np.uint8)\n    mask_lr_ud = np.zeros((IMG_HEIGHT, IMG_WIDTH), np.uint8)\n    for mask_id in range(masks.shape[0]):\n        mask_ = masks[mask_id, :, :]\n        mask_ = mask_ // 255\n        mask_ = resize(mask_, (IMG_HEIGHT, IMG_WIDTH), mode='constant', \n                       preserve_range=True).astype(np.uint8)\n        mask_ = binary_dilation(mask_, selem=square(3))\n        mask += mask_\n        mask_lr += np.fliplr(mask_)\n        mask_ud += np.flipud(mask_)\n        mask_lr_ud += np.flipud(np.fliplr(mask_))\n        \n    Z_train[4*n + 0, :, :, 0] = (mask > 1)\n    Z_train[4*n + 1, :, :, 0] = (mask_lr > 1)\n    Z_train[4*n + 2, :, :, 0] = (mask_ud > 1)\n    Z_train[4*n + 3, :, :, 0] = (mask_lr_ud > 1)\n","execution_count":12,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"e3bfb19b6023ce066afe86b678a761e55afe789c"},"cell_type":"code","source":"for i in tqdm(range(len(train_split) * 4)):\n    Z_train[i, :, :, 1] = Y_train[i, :, :, 0]\n    Z_train[i, :, :, 2] = np.where(Z_train[i, :, :, 0] == 1, 0, 1 - Y_train[i, :, :, 0])","execution_count":13,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"11c961de98bfc4e10a01c4450e89a40092e61e2c"},"cell_type":"code","source":"# Get and resize test images and masks\nX_test = np.zeros((len(test_split) * 4, IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS), np.uint8)\nprint('Getting and resizing test images ... ')\nsys.stdout.flush()\nfor n, img in enumerate(tqdm(test_split)):\n    img = img[:, :, :IMG_CHANNELS]\n    img = resize(img, (IMG_HEIGHT, IMG_WIDTH), mode='constant', preserve_range=True)\n    img_mean = np.mean(img, axis=2).astype(np.uint8)\n    for c in range(IMG_CHANNELS):\n        img[:,:,c] = img_mean\n    X_test[n * 4 + 0] = img\n    X_test[n * 4 + 1] = np.fliplr(img)\n    X_test[n * 4 + 2] = np.flipud(img)\n    X_test[n * 4 + 3] = np.flipud(np.fliplr(img))    ","execution_count":14,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"ed53fcc72c59ecda2a1da85aa429f22decf83055"},"cell_type":"code","source":"with open('../data/X_train_%s.pkl' % fname, 'wb') as file:\n    pkl.dump(X_train, file, protocol=pkl.HIGHEST_PROTOCOL)\n\nwith open('../data/Z_train_%s.pkl' % fname, 'wb') as file:\n    pkl.dump(Z_train, file, protocol=pkl.HIGHEST_PROTOCOL)\n\nwith open('../data/X_test_%s.pkl' % fname, 'wb') as file:\n    pkl.dump(X_test, file, protocol=pkl.HIGHEST_PROTOCOL)","execution_count":23,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"989ef46ca8e5af22a34b162ef332ea54bdb7ad91"},"cell_type":"code","source":"with open('../data/X_train_%s.pkl' % fname, 'rb') as file:\n    X_train= pkl.load(file)\n\nwith open('../data/Z_train_%s.pkl' % fname, 'rb') as file:\n    Z_train = pkl.load(file)\n\nwith open('../data/X_test_%s.pkl' % fname, 'rb') as file:\n    X_test = pkl.load(file)","execution_count":11,"outputs":[]},{"metadata":{"_cell_guid":"c0523b03-1fc5-4505-a1b8-eb35ee617c8a","_uuid":"d4f8327802a1ec6139ce0585953986272ba62ce1"},"cell_type":"markdown","source":"Let's see if things look all right by drawing some random images and their associated masks."},{"metadata":{"trusted":false,"_uuid":"f8c93fa58855d880ceccb16a43cbe2fe67b80090"},"cell_type":"code","source":"# Check if training data looks all right\nix = 23#random.randint(0, len(train_ids))\nprint(ix, train[ix].shape)\nimshow(train[ix])\nplt.show()\n\nfor i in range(8*ix, 8*ix + 2):\n    print(i, X_train[i].shape)\n    imshow(X_train[i])\n    plt.show()\n    imshow(255*Z_train[i])\n    plt.show()","execution_count":12,"outputs":[]},{"metadata":{"_cell_guid":"c3b9f148-1dba-4b6a-981b-6cdbf394fc3c","_uuid":"986488a4c5223576be370e224426a30431911eb2"},"cell_type":"markdown","source":"# Build and train our neural network\nNext 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.\n\n![](https://lmb.informatik.uni-freiburg.de/people/ronneber/u-net/u-net-architecture.png)"},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"f99fb044ad02f3ac8f21b9bde718eed5f3b19a9a"},"cell_type":"code","source":"import keras.backend as K\n\ndef pixelwise_crossentropy(target, output):\n    _epsilon = 10e-8\n    output = tf.clip_by_value(output, _epsilon, 1. - _epsilon)\n    weight = 30 * target[:,:,:,0:1] + 3 * target[:,:,:,1:2] + 1 * target[:,:,:,2:3]\n    return - tf.reduce_sum(target * weight *  tf.log(output) +\n                           (1 - target)  *  tf.log(1 - output),\n                           len(output.get_shape()) - 1)\n","execution_count":13,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"190fd197d7647e99bc03aa7a0c46cd593487b236"},"cell_type":"code","source":"from keras.engine import Layer\nfrom keras import backend as K\n\nclass SpeckleNoise(Layer):\n    \"\"\"Apply multiplicative one-centered Gaussian noise.\n    This is useful to mitigate overfitting\n    (you could see it as a form of random data augmentation).\n    Speckle Noise (GS) is a natural choice as corruption process\n    for real valued inputs.\n    As it is a regularization layer, it is only active at training time.\n    # Arguments\n        stddev: float, standard deviation of the noise distribution.\n    # Input shape\n        Arbitrary. Use the keyword argument `input_shape`\n        (tuple of integers, does not include the samples axis)\n        when using this layer as the first layer in a model.\n    # Output shape\n        Same shape as input.\n    \"\"\"\n\n#    @interfaces.legacy_specklenoise_support\n    def __init__(self, stddev, **kwargs):\n        super(SpeckleNoise, self).__init__(**kwargs)\n        self.supports_masking = True\n        self.stddev = stddev\n\n    def call(self, inputs, training=None):\n        def noised():\n            return K.clip(inputs * K.random_normal(shape=K.shape(inputs),\n                                            mean=1.,\n                                            stddev=self.stddev), 0.0, 1.0)\n        return K.in_train_phase(noised, inputs, training=training)\n\n    def get_config(self):\n        config = {'stddev': self.stddev}\n        base_config = super(SpeckleNoise, self).get_config()\n        return dict(list(base_config.items()) + list(config.items()))\n\n    def compute_output_shape(self, input_shape):\n        return input_shape","execution_count":14,"outputs":[]},{"metadata":{"_cell_guid":"c1dbc57c-b497-4ccb-b077-2053203ab7ed","_uuid":"0aa97d66c29f45dfac9b0f45fcf74ba0e778ba5d","collapsed":true,"trusted":false},"cell_type":"code","source":"# Build U-Net model\ndef get_model(verbose=False):\n    inputs = Input((IMG_HEIGHT, IMG_WIDTH, IMG_CHANNELS))\n    s = Lambda(lambda x: x / 255) (inputs)\n    \n    s = SpeckleNoise(0.01)(s)  #skimage speckel var defaults to 0.01\n\n    c1 = Conv2D(16, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (s)\n    c1 = Dropout(0.1) (c1)\n    c1 = Conv2D(16, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c1)\n    p1 = MaxPooling2D((2, 2)) (c1)\n\n    c2 = Conv2D(32, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (p1)\n    c2 = Dropout(0.1) (c2)\n    c2 = Conv2D(32, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c2)\n    p2 = MaxPooling2D((2, 2)) (c2)\n\n    c3 = Conv2D(64, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (p2)\n    c3 = Dropout(0.2) (c3)\n    c3 = Conv2D(64, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c3)\n    p3 = MaxPooling2D((2, 2)) (c3)\n\n    c4 = Conv2D(128, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (p3)\n    c4 = Dropout(0.2) (c4)\n    c4 = Conv2D(128, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c4)\n    p4 = MaxPooling2D(pool_size=(2, 2)) (c4)\n\n    c5 = Conv2D(256, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (p4)\n    c5 = Dropout(0.3) (c5)\n    c5 = Conv2D(256, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c5)\n\n    u6 = Conv2DTranspose(128, (2, 2), strides=(2, 2), padding='same') (c5)\n    u6 = concatenate([u6, c4])\n    c6 = Conv2D(128, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (u6)\n    c6 = Dropout(0.2) (c6)\n    c6 = Conv2D(128, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c6)\n\n    u7 = Conv2DTranspose(64, (2, 2), strides=(2, 2), padding='same') (c6)\n    u7 = concatenate([u7, c3])\n    c7 = Conv2D(64, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (u7)\n    c7 = Dropout(0.2) (c7)\n    c7 = Conv2D(64, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c7)\n\n    u8 = Conv2DTranspose(32, (2, 2), strides=(2, 2), padding='same') (c7)\n    u8 = concatenate([u8, c2])\n    c8 = Conv2D(32, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (u8)\n    c8 = Dropout(0.1) (c8)\n    c8 = Conv2D(32, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c8)\n\n    u9 = Conv2DTranspose(16, (2, 2), strides=(2, 2), padding='same') (c8)\n    u9 = concatenate([u9, c1], axis=3)\n    c9 = Conv2D(16, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (u9)\n    c9 = Dropout(0.1) (c9)\n    c9 = Conv2D(16, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c9)\n\n    outputs = Conv2D(3, (1, 1), activation='sigmoid') (c9)\n\n    model = Model(inputs=[inputs], outputs=[outputs])\n    optimizer = Adam(clipvalue=5)\n    model.compile(optimizer=optimizer, loss=pixelwise_crossentropy)\n    if verbose:\n        model.summary()\n    return model","execution_count":40,"outputs":[]},{"metadata":{"_cell_guid":"72330944-6ce7-4070-b276-c3c4b20c4fe5","_uuid":"92350b6e18cc50f3fa7b6e9a02d39fcbff8238f7"},"cell_type":"markdown","source":"*Update: Changed to ELU units, added dropout.*\n\nNext 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. \n\n*Update: Added early stopping and checkpointing and increased to 30 epochs.*"},{"metadata":{"trusted":false,"_uuid":"2773a2d117917f669e9c399f3407f40c5c6fa99c"},"cell_type":"code","source":"try:\n    sess.close()\nexcept:\n    pass\n\nsess = init_seeds(0)\n\n# even number of folds because we duplicate images\nkf = KFold(6, shuffle=False)\n\nmodels = []\n\nfor fold, (train_idx, val_idx) in enumerate(kf.split(X_train)):\n    print('*' * 40)\n    print('Fold:', fold)\n    X_train_kf = X_train[train_idx]\n    X_val_kf = X_train[val_idx]\n    Z_train_kf = Z_train[train_idx]\n    Z_val_kf = Z_train[val_idx]\n\n    model = get_model(verbose=(fold==0))\n    models.append(model)\n\n    # Fit model\n    earlystopper = EarlyStopping(patience=5, verbose=1)\n    checkpointer = ModelCheckpoint('model_%s_%d.h5' % (fname, fold), \n                                   verbose=1, save_best_only=True)\n    results = model.fit(X_train_kf, Z_train_kf, \n                        validation_data = (X_val_kf, Z_val_kf),\n                        batch_size=2, epochs=40, shuffle=True,\n                        callbacks=[earlystopper, checkpointer])\n\n#sess.close()","execution_count":41,"outputs":[]},{"metadata":{"_cell_guid":"1f381f5b-1b71-4daa-a417-e02f4894540b","_uuid":"bb15226ea617cf91ed8f43179fccb5a15809e5a0"},"cell_type":"markdown","source":"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! \n\n# Make predictions\n\nLet'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."},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"ec9066a92c6e7f50717e9e5f7af221172e189e24"},"cell_type":"code","source":"try:\n    sess.close()\nexcept:\n    pass\n\nsess = init_seeds(0)\n\n# even number of folds because we duplicate images\nkf = KFold(6, shuffle=False)\n","execution_count":15,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"fb7d6958c2ee700cb7ba81e381a05c33f93e9744"},"cell_type":"code","source":"preds_train = np.zeros(Z_train.shape)\npreds_test = 0\n# Predict on train, val and test\nfor fold, (train_idx, val_idx) in enumerate(kf.split(X_train)):\n    model = load_model('model_%s_%d.h5' % (fname, fold), \n                       custom_objects={'pixelwise_crossentropy':pixelwise_crossentropy,\n                                       'SpeckleNoise':SpeckleNoise,\n                                      })\n    X_val_kf = X_train[val_idx]\n    preds_train[val_idx] = model.predict(X_val_kf, verbose=1)\n    preds_test += model.predict(X_test, verbose=1)\npreds_test /= 6","execution_count":16,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"931c6fbdee706bcb2f85e904bbd5641cc21f788b"},"cell_type":"code","source":"# Create list of upsampled train preds\npreds_train_upsampled = []\nfor i in tqdm(range(len(train_split))):\n    train_i = train_split[i]\n    shape = (train_i.shape[0], train_i.shape[1], 3)\n    pred = resize((preds_train[4*i + 0]), \n                  shape, \n                  mode='constant', preserve_range=True)\n    pred += np.fliplr(resize((preds_train[4*i + 1]), \n                  shape, \n                  mode='constant', preserve_range=True))\n    pred += np.flipud(resize((preds_train[4*i + 2]), \n                  shape, \n                  mode='constant', preserve_range=True))\n    pred += np.flipud(np.fliplr(resize((preds_train[4*i + 3]), \n                  shape, \n                  mode='constant', preserve_range=True)))\n    #pred = (pred > 4*threshold).astype(np.uint8)\n    pred /= 4\n    preds_train_upsampled.append(pred)","execution_count":17,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"74bccb0e41a3e9640d529ecdfe91ccfa8b2e33c0"},"cell_type":"code","source":"def merge(img1, img2, shape):\n    ov2 = 5\n    height = shape[0]\n    width = shape[1]\n    img = np.zeros((height, width, 3), dtype=np.float32)\n    w = np.zeros((height, width, 1), dtype=np.float32)\n    height1 = img1.shape[0]\n    width1 = img1.shape[1]\n    height2 = img2.shape[0]\n    width2 = img2.shape[1]  \n    w1 = 10*ov2*np.ones((height1, width1, 1), dtype=np.float32)\n    w2 =  10*ov2*np.ones((height2, width2, 1), dtype=np.float32)\n    for i in range(ov2, 0, -1):\n        w1[i-1,:] = 10*i\n        w1[height1 - i, :] = 10*i\n        w1[:, i-1] = 10*i\n        w1[:, width1 - i] = 10*i\n        w2[i-1,:] = 10*i\n        w2[height2 - i, :] = 10*i\n        w2[:, i-1] = 10*i\n        w2[:, width2 - i] = 10*i\n\n    if height > 2*width:\n        half = int(height//2)\n        img[:half, :, :] += w1*img1\n        img[half:, :, :] += w2*img2\n        w[:half, :] += w1\n        w[half:, :] += w2\n        img /= w\n    elif height > width:\n        img[:width, :, :] += w1*img1\n        img[height-width:, :, :] += w2*img2\n        w[:width, :] += w1\n        w[height-width:, :] += w2\n        img /= w\n    elif width > 2*height:\n        half = int(width//2)\n        img[:, :half, :] += w1*img1\n        img[:, half:, :] += w2*img2\n        w[:, :half] += w1\n        w[:, half:] += w2\n        img /= w\n    else:\n        img[:, :height, :] += w1*img1 \n        img[:, width-height:, :] += w2*img2\n        w[:, :height] += w1 \n        w[:, width-height:] += w1\n        img /= w\n    return (255*img).astype(np.uint8)","execution_count":18,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"85b237e5e58438a9b122dc19a3fce614aa13dab4"},"cell_type":"code","source":"preds_train_merged = []\nfor ix in tqdm(range(len(train))):    \n    merged = merge(preds_train_upsampled[2*ix+0], \n             preds_train_upsampled[2*ix+1],\n             train[ix].shape\n            )\n    preds_train_merged.append(merged)","execution_count":19,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"13ac9c685bd10d1f9e276d33972f160d6a05cd1f"},"cell_type":"code","source":"preds_test_upsampled = []\nfor i in tqdm(range(len(test_split))):\n    test_i = test_split[i]\n    pred = resize(np.squeeze(preds_test[4*i + 0]), \n                  (test_i.shape[0], test_i.shape[1]), \n                  mode='constant', preserve_range=True)\n    pred += np.fliplr(resize(np.squeeze(preds_test[4*i + 1]), \n                  (test_i.shape[0], test_i.shape[1]), \n                  mode='constant', preserve_range=True))\n    pred += np.flipud(resize(np.squeeze(preds_test[4*i + 2]), \n                  (test_i.shape[0], test_i.shape[1]), \n                  mode='constant', preserve_range=True))\n    pred += np.flipud(np.fliplr(resize(np.squeeze(preds_test[4*i + 3]), \n                  (test_i.shape[0], test_i.shape[1]), \n                  mode='constant', preserve_range=True)))\n    #pred = (pred > 4*threshold).astype(np.uint8)\n    pred /= 4\n    preds_test_upsampled.append(pred)\n","execution_count":20,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"afa788a0d6602d9c5c3d1fa704295311cbdc1fdf"},"cell_type":"code","source":"preds_test_merged = []\nfor ix in tqdm(range(len(test))):    \n    merged = merge(preds_test_upsampled[2*ix+0], \n                 preds_test_upsampled[2*ix+1],\n                 test[ix].shape\n            )\n    preds_test_merged.append(merged)","execution_count":21,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"f69a507052af843e59e61b4273862bb9c7dcf445"},"cell_type":"code","source":"ix = 17\nimshow(train[ix])\nplt.show()\n\nimshow(train_split[2*ix+0])\nplt.show()\n\nimshow(np.sum(train_mask_split[2*ix+0], axis=0))\nplt.show()\n\nimshow((preds_train_upsampled[2*ix+0]))\nplt.show()\nimshow(merge(preds_train_upsampled[2*ix+0], \n             preds_train_upsampled[2*ix+1],\n             train[ix].shape\n            ))\nplt.show()","execution_count":22,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"c246c7801b23ea35f86ef1145b5744c690ff4aeb"},"cell_type":"code","source":"from skimage.morphology import label\n\ndef get_labels(y):\n    labels = np.zeros((y.shape[1], y.shape[2]))\n    for i in range(y.shape[0]):\n        labels = np.where(y[i,:,:] > 0, i+1, labels)\n    return labels\n\ndef iou_score_cuk(y_true, y_pred, verbose=True, thresholds=np.arange(0.5, 1.0, 0.05)):\n    y_true = get_labels(y_true)\n    y_pred = get_labels(y_pred)\n    # Compute number of objects\n    true_objects = len(np.unique(y_true))\n    pred_objects = len(np.unique(y_pred))\n    if verbose:\n        print(\"Number of true objects:\", true_objects - 1)\n        print(\"Number of predicted objects:\", pred_objects - 1)\n    \n    intersection = np.histogram2d(y_true.flatten(), y_pred.flatten(), \n                                  bins=(true_objects, pred_objects))[0].astype('int')\n\n    area_true = np.histogram(y_true, bins = true_objects)[0]\n    area_true = np.expand_dims(area_true, -1)\n\n    area_pred = np.histogram(y_pred, bins = pred_objects)[0]\n    area_pred = np.expand_dims(area_pred, 0)\n\n    union = area_true + area_pred - intersection\n    intersection = intersection[1:,1:]\n    union = union[1:,1:]\n    union[union == 0] = 1e-9\n    \n    iou = intersection / union\n    \n    # Precision helper function\n    def precision_at(threshold, iou):\n        matches = iou > threshold\n        true_positives = np.sum(matches, axis=1) == 1   # Correct objects\n        false_positives = np.sum(matches, axis=0) == 0  # Missed objects\n        false_negatives = np.sum(matches, axis=1) == 0  # Extra objects\n        tp, fp, fn = np.sum(true_positives), np.sum(false_positives), np.sum(false_negatives)\n        return tp, fp, fn    \n    \n    prec = []\n    \n    if verbose:\n        print(\"Thresh\\tTP\\tFP\\tFN\\tPrec.\")\n    for t in thresholds:\n        tp, fp, fn = precision_at(t, iou)\n        if (tp + fp + fn) == 0:\n            p = 1\n        else:\n            p = tp / (tp + fp + fn)\n        if verbose:\n            print(\"{:1.3f}\\t{}\\t{}\\t{}\\t{:1.3f}\".format(t, tp, fp, fn, p))\n        prec.append(p)\n    if verbose:\n        print(\"AP\\t-\\t-\\t-\\t{:1.3f}\".format(np.mean(prec)))\n    return np.mean(prec)","execution_count":23,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"f2f2eaa19a8dcfc4a39f9ded63f688f23415865d"},"cell_type":"code","source":"from scipy.ndimage.morphology import binary_fill_holes\n\ndef get_pred_watershed(upsampled, area_threshold, threshold, sep_threshold, \n             spread_threshold, alpha, connectivity=2):\n    img = ((upsampled[:,:,1] > 255 * threshold) &\n                    (upsampled[:,:,0] < 255 * sep_threshold))\n    img = binary_fill_holes(img)\n    img = remove_small_objects(img, area_threshold)\n    lab_img = label(img, connectivity=connectivity)\n    distance = upsampled[:,:,1] + alpha * upsampled[:,:,0]\n    img = 1 * ((distance > 255 * spread_threshold) )\n    \n    lab_img = img * watershed(- upsampled[:,:,1], lab_img)\n\n    y_pred = np.zeros((lab_img.max(), lab_img.shape[0], lab_img.shape[1]), np.uint16)\n    i = 0\n    for lab in range(lab_img.max()):\n        tmp = (lab_img == lab+1)\n        if np.sum(tmp.ravel()) > area_threshold:\n            y_pred[i,:,:] = tmp\n            i += 1\n    return y_pred[:i]\n\nfrom scipy.ndimage.morphology import binary_fill_holes\n\ndef get_pred_random_walker(upsampled, area_threshold, threshold, sep_threshold, \n             spread_threshold, alpha, connectivity=2):\n    img = ((upsampled[:,:,1] > 255 * threshold) &\n                    (upsampled[:,:,0] < 255 * sep_threshold))\n    \n    img = binary_fill_holes(img)\n    img = remove_small_objects(img, area_threshold)\n    markers = label(img, connectivity=connectivity)\n    distance = upsampled[:,:,1] + alpha * upsampled[:,:,0]\n    mask = ((distance > 255 * spread_threshold) )\n    markers[~mask] = -1\n    \n    lab_img = random_walker(mask, markers)\n\n    y_pred = np.zeros((lab_img.max(), lab_img.shape[0], lab_img.shape[1]), np.uint16)\n    i = 0\n    for lab in range(lab_img.max()):\n        tmp = (lab_img == lab+1)\n        if np.sum(tmp.ravel()) > area_threshold:\n            y_pred[i,:,:] = tmp\n            i += 1\n    return y_pred[:i]\n\ndef get_pred(upsampled, area_threshold, threshold, sep_threshold, \n             spread_threshold, alpha, connectivity=2):\n    try:\n        return get_pred_random_walker(upsampled, area_threshold, threshold, \n                                          sep_threshold, spread_threshold, \n                                      alpha, connectivity)\n    except:\n        return get_pred_watershed(upsampled, area_threshold, threshold, \n                                          sep_threshold, spread_threshold, \n                                      alpha, connectivity)\n        ","execution_count":24,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"4ff749415bc36e4ef5eac49328a2502c003143cf"},"cell_type":"code","source":"area_threshold = 20\nthreshold = 0.75\nsep_threshold = 0.6\nspread_threshold = 0.4\nalpha=0.4\n\ndef get_pred(upsampled, area_threshold=area_threshold, \n             threshold=threshold, sep_threshold=sep_threshold, \n             spread_threshold=spread_threshold, alpha=alpha, connectivity=2):\n    try:\n        return get_pred_random_walker(upsampled, area_threshold, threshold, \n                                          sep_threshold, spread_threshold, \n                                      alpha, connectivity=2)\n    except:\n        return get_pred_watershed(upsampled, area_threshold, threshold, \n                                          sep_threshold, spread_threshold, \n                                      alpha, connectivity=2)","execution_count":25,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"96a533d9ece7ef94d6e3e6d6ceb885e08cf2a412"},"cell_type":"code","source":"# Run-length encoding stolen from https://www.kaggle.com/rakhlin/fast-run-length-encoding-python\ndef rle_encoding(x):\n    dots = np.where(x.T.flatten() == 1)[0]\n    run_lengths = []\n    prev = -2\n    for b in dots:\n        if (b>prev+1): run_lengths.extend((b + 1, 0))\n        run_lengths[-1] += 1\n        prev = b\n    return run_lengths\n            \ndef pred_to_rles(y_pred):\n    for i in range(y_pred.shape[0]):\n        tmp = y_pred[i]\n        yield rle_encoding(tmp)\n","execution_count":26,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"5f6cd77e0c7e2c62f143bac29d3c45a03cb26d7d"},"cell_type":"code","source":"score = np.zeros(len(train))\n\nnew_train_ids = []\nrles = []\n                 \nfor n, id_ in enumerate(tqdm(train_ids)):  \n    y_pred = get_pred(preds_train_merged[n])\n    score[n] = iou_score_cuk(train_mask[n], y_pred, verbose=False)\n    rle = list(pred_to_rles(y_pred))\n    rles.extend(rle)\n    new_train_ids.extend([id_] * len(rle))\n\nprint(len(rles))\n    \nsub = pd.DataFrame()\nsub['ImageId'] = new_train_ids\nsub['EncodedPixels'] = pd.Series(rles).apply(lambda x: ' '.join(str(y) for y in x))\n\nsub.to_csv('../submissions/keras_unet_67_train.csv', index=False)\ntrain_score = np.mean(score)\n\nprint('%0.5f' % train_score) ","execution_count":27,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"c13e278e6db96fb3548f77057cb082ca09bf6d04"},"cell_type":"code","source":"new_test_ids = []\nrles = []\nfor n, id_ in tqdm(enumerate(test_ids)):\n    y_pred = get_pred(preds_test_merged[n])\n    rle = list(pred_to_rles(y_pred))\n    rles.extend(rle)\n    new_test_ids.extend([id_] * len(rle))\n\nprint(len(rles))\n\n# Create submission DataFrame\nsub = pd.DataFrame()\nsub['ImageId'] = new_test_ids\nsub['EncodedPixels'] = pd.Series(rles).apply(lambda x: ' '.join(str(y) for y in x))\nsub.to_csv('../submissions/keras_unet_67_test.csv', index=False)","execution_count":28,"outputs":[]},{"metadata":{"trusted":false,"collapsed":true,"_uuid":"50bb6655eb2911b65d94e0d0ac5e88c906433d1a"},"cell_type":"code","source":"with open('../data/%s_train_pred.pkl' % fname, 'wb') as file:\n    pkl.dump(preds_train_merged, file)\n\nwith open('../data/%s_test_pred.pkl' % fname, 'wb') as file:\n    pkl.dump(preds_test_merged, file)","execution_count":29,"outputs":[]},{"metadata":{"collapsed":true,"trusted":false,"_uuid":"f2f226f0e2ba0eaa906236e43d8eae6bbb42f103"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.1"}},"nbformat":4,"nbformat_minor":1}