{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Load modules"},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"%matplotlib inline\nimport os\nimport glob\nimport tqdm\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport multiprocessing\nfrom functools import partial\n\nimport pydicom as dicom\nimport nibabel as nib","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Set path"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"# TODO: modify the following paths\ntrain_image_folder = \"../input/train-images/image/\"\ntrain_label_folder = \"../input/train-labels/label/\"\ntest_image_folder = \"../input/test-images/image/\"\n\ntrain_list = os.listdir(train_image_folder)\n\n# Ignore this data\nif 'hmvsa0loxh3ek2y8rzmcyb6zrrh9mwyp' in train_list:\n    train_list.remove('hmvsa0loxh3ek2y8rzmcyb6zrrh9mwyp')\n    \nprint('Train data:', len(train_list))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_dicom_volume(src_dir, suffix='*.dcm'):\n    \"\"\"Load DICOM volume and get meta data.\n    \"\"\"\n    encode_name = src_dir.split('/')[-1]\n    # Read dicom files from the source directory\n    # Sort the dicom slices in their respective order by slice location\n    dicom_scans = [dicom.read_file(sp) \\\n                   for sp in glob.glob(os.path.join(src_dir, suffix))]\n    # dicom_scans.sort(key=lambda s: float(s.SliceLocation))\n    dicom_scans.sort(key=lambda s: float(s[(0x0020, 0x0032)][2]))\n\n    # Convert to int16, should be possible as values should always be low enough\n    # Volume image is in z, y, x order\n    volume_image = np.stack([ds.pixel_array \\\n                             for ds in dicom_scans]).astype(np.int16)\n    return encode_name, volume_image\n\ndef load_label(label_fpath, transpose=False):\n    encode_name = label_fpath[-39: -7]\n    label_data = nib.load(label_fpath)\n    label_array = label_data.get_fdata()\n    if transpose:\n        label_array = np.transpose(label_array, axes=(2, 1, 0))\n    return encode_name, label_array","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## DICOM to npz"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_image_npz_folder = './npz/train_images/'\n\nif not os.path.exists(train_image_npz_folder):\n    os.makedirs(train_image_npz_folder)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for encode in tqdm.tqdm(train_list):\n    _, volume_image = load_dicom_volume(os.path.join(train_image_folder, encode))\n    npz_folder = os.path.join(train_image_npz_folder, encode)\n    if not os.path.exists(npz_folder):\n        os.mkdir(npz_folder) \n        \n    num_slice = volume_image.shape[0]\n    for _z in range(0, num_slice):\n        npz_path = os.path.join(npz_folder, \"%03d.npz\"%(_z))\n        np.savez_compressed(npz_path, image=volume_image[_z])\n        \n    del volume_image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!ls './npz/train_images/'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## NIFTI to npz"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_label_npz_folder = './npz/train_labels/'\n\nif not os.path.exists(train_label_npz_folder):\n    os.makedirs(train_label_npz_folder)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for encode in tqdm.tqdm(train_list):\n    _, label_array = load_label(os.path.join(train_label_folder, encode + '.nii.gz'), transpose=True)\n    npz_folder = os.path.join(train_label_npz_folder, encode)\n    if not os.path.exists(npz_folder):\n        os.mkdir(npz_folder) \n        \n    num_slice = label_array.shape[0]\n    for _z in range(0, num_slice):\n        npz_path = os.path.join(npz_folder, \"%03d.npz\"%(_z))\n        np.savez_compressed(npz_path, label=label_array[_z])\n        \n    del label_array","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Try to apply UNet "},{"metadata":{"trusted":true},"cell_type":"code","source":"from keras.models import Model, load_model\nfrom keras import layers as klayers\nfrom keras.optimizers import Adam\nfrom keras import utils as kutils\nfrom keras import backend as K\nfrom keras.callbacks import ModelCheckpoint\n\n# Make sure keras running on GPU\nK.tensorflow_backend._get_available_gpus()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Data generator"},{"metadata":{"trusted":true},"cell_type":"code","source":"map_image_list = sorted(glob.glob(os.path.join(train_image_npz_folder, '*/*.npz')))\nmap_label_list = sorted(glob.glob(os.path.join(train_label_npz_folder, '*/*.npz')))\n\nmap_df = pd.DataFrame(data={'image': map_image_list, 'label': map_label_list})\nmap_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class LungSliceModelGenerator(kutils.Sequence):\n    'Generates data for Keras'\n    def __init__(self, mapping_df, batch_size, shuffle=True):\n        'Initialization'\n        self.mapping_df = mapping_df\n        self.data_num   = mapping_df.shape[0]\n        self.batch_size = batch_size\n        self.shuffle = shuffle\n        self.on_epoch_end()\n\n    def __len__(self):\n        'Denotes the number of batches per epoch'\n        return int(np.floor(self.data_num / self.batch_size))\n\n\n    def __getitem__(self, index):\n        'Generate one batch of data'\n        # Generate indexes of the batch\n        batch_mapping_df = \\\n            self.mapping_df.iloc[index*self.batch_size: (index+1)*self.batch_size]\n\n        # Generate data\n        X, y = self.__data_generation(batch_mapping_df)\n        return X, y\n\n\n    def on_epoch_end(self):\n        'Updates indexes after each epoch'\n        if self.shuffle:\n            self.mapping_df = self.mapping_df.sample(frac=1).reset_index(drop=True)\n            \n    def __data_generation(self, batch_mapping_df):\n        'Generates data containing batch_size samples' \n        # Initialization\n        X = np.zeros((self.batch_size, 512, 512, 1))\n        y = np.zeros((self.batch_size, 512, 512, 1))\n\n        # Generate data\n        cnt = 0\n        for i, row in batch_mapping_df.iterrows():\n            X[cnt, :, :, 0] = np.load(row['image'])['image']\n            y[cnt, :, :, 0] = np.load(row['label'])['label']\n            cnt += 1\n        return X, y","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"batch_size = 16\nslice_generator = LungSliceModelGenerator(map_df, batch_size=batch_size)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Loss function"},{"metadata":{"trusted":true},"cell_type":"code","source":"def _dice_coefficient(threshold = 0.3):\n    def hard_dice_coefficient(y_true, y_pred, smooth=1.0):\n        y_true_f = K.flatten(K.cast(y_true > threshold, dtype=float))\n        y_pred_f = K.flatten(K.cast(y_pred > threshold, dtype=float))\n        intersection = K.sum(y_true_f * y_pred_f)\n        return (2. * intersection + smooth) / (K.sum(y_true_f) + K.sum(y_pred_f) + smooth)\n    return hard_dice_coefficient\n\ndef dice_coefficient_loss(y_true, y_pred):\n    return 1 - _dice_coefficient()(y_true, y_pred)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Model architecture"},{"metadata":{"trusted":true},"cell_type":"code","source":"def unet(pretrained_weights=None, input_size=[512, 512, 1], depth=3, init_filter=8, \n         filter_size=3, padding='same', pool_size=[2, 2], strides=[2, 2]):\n    \n    inputs = klayers.Input(input_size)\n    \n    current_layer = inputs\n    encoding_layers = []\n    \n    # Encoder path\n    for d in range(depth + 1):\n        num_filters = init_filter * 2 ** d\n        \n        conv = klayers.Conv2D(num_filters, filter_size, padding=padding, kernel_initializer='he_normal')(current_layer)\n        conv = klayers.BatchNormalization()(conv)\n        conv = klayers.Activation('relu')(conv)\n        conv = klayers.Conv2D(num_filters * 2, filter_size, padding=padding, kernel_initializer='he_normal')(conv)\n        conv = klayers.BatchNormalization()(conv)\n        conv = klayers.Activation('relu')(conv)\n        encoding_layers.append(conv)\n    \n        pool = klayers.MaxPooling2D(pool_size=pool_size)(conv)\n        \n        if d == depth:\n            # Bridge\n            current_layer = conv\n        else:\n            current_layer = pool\n\n        \n    # Decoder path\n    for d in range(depth, 0, -1):\n        num_filters = init_filter * 2 ** d\n        up = klayers.Deconvolution2D(num_filters * 2, pool_size, strides=strides)(current_layer)\n\n        crop_layer = encoding_layers[d - 1]\n        # Calculate two layers shape\n        up_shape = np.array(up._keras_shape[1:-1])\n        conv_shape = np.array(crop_layer._keras_shape[1:-1])\n\n        # Calculate crop size of left and right\n        crop_left = (conv_shape - up_shape) // 2\n\n        crop_right = (conv_shape - up_shape) // 2 + (conv_shape - up_shape) % 2\n        crop_sizes = tuple(zip(crop_left, crop_right))\n\n        crop = klayers.Cropping2D(cropping=crop_sizes)(crop_layer)\n\n        # Concatenate\n        up = klayers.Concatenate(axis=-1)([crop, up])\n        conv = klayers.Conv2D(num_filters, filter_size, padding=padding, kernel_initializer='he_normal')(up)\n        conv = klayers.BatchNormalization()(conv)\n        conv = klayers.Activation('relu')(conv)\n        conv = klayers.Conv2D(num_filters, filter_size, padding=padding, kernel_initializer='he_normal')(conv)\n        conv = klayers.BatchNormalization()(conv)\n        conv = klayers.Activation('relu')(conv)\n        \n        current_layer = conv\n    \n    \n    outputs = klayers.Conv2D(1, 1, padding=padding, kernel_initializer='he_normal')(current_layer)\n    outputs = klayers.Activation('sigmoid')(outputs)\n    model = Model(inputs=inputs, outputs=outputs)\n\n    if(pretrained_weights):\n        model.load_weights(pretrained_weights)\n\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = unet(depth=3)\nmodel.compile(optimizer=Adam(lr=1e-3), loss='binary_crossentropy', metrics=[_dice_coefficient(0.5)])\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model_folder = os.path.join('./model', 'sample-code')\n\nif not os.path.exists(model_folder):\n    os.makedirs(model_folder)\n\ncallbacks = []\ncallbacks.append(ModelCheckpoint(os.path.join(model_folder, 'model-{epoch:03d}.h5'), \n                                 save_best_only=False, \n                                 period=5))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"history = model.fit_generator(slice_generator,\n                              epochs=15,\n                              verbose=1, \n                              callbacks=callbacks)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Predict test images"},{"metadata":{"trusted":true},"cell_type":"code","source":"def retrieve_pred_str(src_dir, model, threshold=0.4):\n    encode_name = src_dir.split('/')[-1]\n    \n    _, test_volume = load_dicom_volume(src_dir, suffix='*.dcm')\n    \n    pred_label = model.predict(np.expand_dims(test_volume, axis=-1))\n    pred_label = np.transpose(pred_label[:, :, :, 0], axes=(2, 1, 0))\n    \n    pred_label = (pred_label > threshold).astype(np.int)\n\n    label_flatten = pred_label.flatten()\n\n    label_flatten_idx = np.where(label_flatten == 1)[0]\n\n    label_str = ''\n    \n    if label_flatten_idx.size > 0:\n        prev_idx = label_flatten_idx[0]\n        idx_start = label_flatten_idx[0]\n        cnt = 1\n        for _idx in label_flatten_idx[1:]:\n            if _idx == prev_idx+1:\n                cnt += 1\n            else:\n                label_str += str(idx_start) + ' ' + str(cnt) + ' '\n\n                cnt = 1\n                idx_start = _idx\n            prev_idx = _idx\n\n        label_str = label_str.rstrip(' ')\n    return (encode_name, label_str)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"sample_submission = np.genfromtxt('../input/sample_submission.csv', \n                                  delimiter=',', \n                                  dtype='str',\n                                  skip_header = 1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false},"cell_type":"code","source":"test_encode_list = sample_submission[:, 0]","execution_count":null,"outputs":[]},{"metadata":{"scrolled":true,"trusted":false},"cell_type":"code","source":"pred_pair_list = []\n\nfor encode_name in tqdm.tqdm(test_encode_list, total=len(test_encode_list)):\n    (encode, label_str) = retrieve_pred_str(os.path.join(test_image_folder, encode_name), model, threshold=0.4)\n    pred_pair_list.append((encode, label_str))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"solution_path = './sample-code_pred.csv'\nwith open(solution_path, 'w') as f:\n    f.write('encode,pixel_value\\n')\n    for _pair in pred_pair_list:\n        encode = _pair[0]\n        label_str = _pair[1]\n        f.write(encode + ',' + label_str + '\\n')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Remove all npz files"},{"metadata":{"trusted":true},"cell_type":"code","source":"!rm -r ./npz","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.4"}},"nbformat":4,"nbformat_minor":1}