{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport os\nimport cv2\nimport gc\nfrom scipy.ndimage.filters import gaussian_filter\n\nimport tensorflow as tf\nfrom tensorflow.keras.utils import *\nfrom tensorflow.keras.layers import *\nfrom tensorflow.keras.models import *\nfrom tensorflow.keras.optimizers import *\nfrom tensorflow.keras.callbacks import *\nfrom tensorflow.keras.losses import *\nfrom tensorflow.keras.metrics import *\n\nimport tensorflow.keras.backend as K\nimport glob","metadata":{"execution":{"iopub.status.busy":"2022-07-14T07:07:38.497442Z","iopub.execute_input":"2022-07-14T07:07:38.498780Z","iopub.status.idle":"2022-07-14T07:07:45.826178Z","shell.execute_reply.started":"2022-07-14T07:07:38.498619Z","shell.execute_reply":"2022-07-14T07:07:45.824954Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def _get_gaussian(patch_size, sigma_scale = 1 / 8.):\n    tmp = np.zeros(patch_size)\n    center_coords = [i // 2 for i in patch_size]\n    sigmas = [i * sigma_scale for i in patch_size]\n    tmp[tuple(center_coords)] = 1\n    gaussian_importance_map = gaussian_filter(tmp, sigmas, 0, mode='constant', cval=0)\n    gaussian_importance_map = gaussian_importance_map / np.max(gaussian_importance_map) * 1\n    gaussian_importance_map = gaussian_importance_map.astype(np.float32)\n\n    # gaussian_importance_map cannot be 0, otherwise we may end up with nans!\n    gaussian_importance_map[gaussian_importance_map == 0] = np.min(\n        gaussian_importance_map[gaussian_importance_map != 0])\n    \n    result = np.empty((*gaussian_importance_map.shape, 3))\n    \n    for i in range(3):\n        result[..., i] = gaussian_importance_map\n\n    return result","metadata":{"execution":{"iopub.status.busy":"2022-07-14T07:07:45.828820Z","iopub.execute_input":"2022-07-14T07:07:45.829171Z","iopub.status.idle":"2022-07-14T07:07:45.839985Z","shell.execute_reply.started":"2022-07-14T07:07:45.829125Z","shell.execute_reply":"2022-07-14T07:07:45.838760Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\n    CONFIGS\n'''\nscan_lists = glob.glob(f\"../input/uw-madison-gi-tract-image-segmentation/test/**/*.png\", recursive=True)\nDebug = len(scan_lists) == 0\nclass_map = [\"large_bowel\", \"small_bowel\", \"stomach\"]\n\n'''\n    UNET\n'''\nNUM_CLASSES = 3\nMODEL_NAME = \"Unet3D_DS\"\nMODEL_WEIGHT = glob.glob(\"../input/uw-model-weights/residual_ds_unet3d_deep.h5\") + glob.glob(\"../input/uw-model-weights/unet3d_ds_deep.h5\") + glob.glob(\"../input/uw-model-weights/residual_ds_unet3d.h5\") + glob.glob(\"../input/uw-model-weights/unet3d_ds_1000.h5\")\nprint(MODEL_WEIGHT)\n\n'''\n    Inference mode\n'''\nWINDOW_SIZE = (80, 160, 160)\nSTRIDE = tuple([x // 2 for x in WINDOW_SIZE]) # Half overlap\ngaussian_importance_map = _get_gaussian(WINDOW_SIZE, sigma_scale = 1. / 8)\n\n'''\n    Post-preprocessing params\n'''\nmin_mask_area = [11, 29, 36]\n\ntest_data_dir = \"../input/uw-madison-gi-tract-image-segmentation/test\" if not Debug else \"../input/uw-madison-gi-tract-image-segmentation/train\"\nsubmission = pd.read_csv(\"../input/uw-madison-gi-tract-image-segmentation/sample_submission.csv\") if not Debug else pd.read_csv(\"../input/uw-madison-gi-tract-image-segmentation/train.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-07-14T07:24:53.936921Z","iopub.execute_input":"2022-07-14T07:24:53.937512Z","iopub.status.idle":"2022-07-14T07:24:54.487780Z","shell.execute_reply.started":"2022-07-14T07:24:53.937468Z","shell.execute_reply":"2022-07-14T07:24:54.486720Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\n    Utilities functions\n'''\ndef rle_encoding(mask):\n    \n    pixels = mask.flatten()\n    pixels = np.concatenate([[0], pixels,[0]])\n    runs = np.where(pixels[1:] != pixels[:-1])[0] + 1\n    runs[1::2] -= runs[0::2]\n    \n    return ' '.join(str(x) for x in runs)\n\ndef rle_decode(mask_rle, shape):\n    s = mask_rle.split()\n    starts, lengths = [np.asarray(x, dtype=int) for x in (s[0:][::2], s[1:][::2])]\n    starts -= 1\n    ends = starts + lengths\n    img = np.zeros(shape[0]*shape[1], dtype=np.uint8)\n    for lo, hi in zip(starts, ends):\n        img[lo:hi] = 1\n    return img.reshape(shape)\n\ndef get_image_size(image_path):\n    image = cv2.imread(image_path, -1)\n    return image.shape[0:2]\n\ndef get_id_from_path(image_path):\n    result = \"\"\n    filename = os.path.basename(image_path)\n    scans = filename.split(\"_\")\n    result = scans[0] + \"_\" + scans[1]\n    image_path = os.path.dirname(image_path)\n    image_path = os.path.dirname(image_path)\n    result = os.path.basename(image_path) + \"_\" + result\n    return result\n\ndef get_path_from_id(scan_ids):\n    # case#_day#_slice_#\n    paths = []\n    for scan_id in scan_ids:\n        case, day, _, slice_num = scan_id.split(\"_\")\n        path = glob.glob(f\"{test_data_dir}/{case}/{case}_{day}/scans/slice_{slice_num}*.png\")[0]\n        paths.append(path)\n    return paths","metadata":{"execution":{"iopub.status.busy":"2022-07-14T07:07:46.944483Z","iopub.execute_input":"2022-07-14T07:07:46.944868Z","iopub.status.idle":"2022-07-14T07:07:46.960369Z","shell.execute_reply.started":"2022-07-14T07:07:46.944820Z","shell.execute_reply":"2022-07-14T07:07:46.958998Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\n    Add file path\n'''\nsubmission[\"file_path\"] = get_path_from_id(submission[\"id\"])\nsubmission.head()","metadata":{"execution":{"iopub.status.busy":"2022-07-14T07:07:46.962194Z","iopub.execute_input":"2022-07-14T07:07:46.963426Z","iopub.status.idle":"2022-07-14T07:09:15.410339Z","shell.execute_reply.started":"2022-07-14T07:07:46.963363Z","shell.execute_reply":"2022-07-14T07:09:15.409285Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\n    TEST DATA LOADER\n'''\nclass DataLoader(Sequence):\n    def __init__(self, scan_paths):\n        self.scan_paths = np.array(scan_paths)\n        case_index = [0]\n        index = 0\n        while index < len(scan_paths):\n            index += len(os.listdir(os.path.dirname(scan_paths[index])))\n            case_index.append(index)\n        self.case_index = case_index\n        \n    def load_data(self, scan_path):\n        X = cv2.imread(scan_path, -1)\n        return X\n    \n    def preprocess(self, X):\n        # Clipping to remove noise\n        top = np.percentile(X.flatten(), 95)\n        bottom = np.percentile(X.flatten(), 5)\n        X[X > top] = top\n        X[X < bottom] = bottom\n        \n        # Minmax normalization\n        X = X.astype('float32')\n        min_val = X.min()\n        max_val = X.max()\n        return (X - min_val) / (max_val - min_val)\n\n    def __len__(self):\n        return len(self.case_index) - 1\n\n    def __getitem__(self, index):\n        indices = self.scan_paths[self.case_index[index]:self.case_index[index + 1]]\n        X = []\n        for i in range(len(indices)):\n            scan_path = indices[i]\n            X.append(self.load_data(scan_path))\n        X = np.array(X, dtype='float32')\n        X = self.preprocess(X)\n        original_shape = X.shape\n        X = np.expand_dims(X,axis=-1)\n        X_out = np.array([X], dtype='float32')\n        return np.array([X], dtype='float32'), original_shape, indices","metadata":{"execution":{"iopub.status.busy":"2022-07-14T07:09:15.412471Z","iopub.execute_input":"2022-07-14T07:09:15.412806Z","iopub.status.idle":"2022-07-14T07:09:15.427039Z","shell.execute_reply.started":"2022-07-14T07:09:15.412760Z","shell.execute_reply":"2022-07-14T07:09:15.425961Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Unet3D:\n    def __init__(self, num_classes = 3, \n                 input_shape = (None, None, None, 1),\n                 conv_settings = [32, 64, 128, 256, 320], \n                 deep_supervision = True,\n                 activation = 'sigmoid'):\n        self.input_shape = input_shape\n        self.num_classes = num_classes\n        self.conv_settings = conv_settings\n        self.deep_supervision = deep_supervision\n        self.activation = activation\n    \n    def conv_in_relu(self, inp, kernels, kernel_size = 3, stride = 1, bn_relu = True):\n        x = inp\n        x = Conv3D(kernels, \n                    kernel_size = kernel_size,\n                    padding = 'same',\n                    strides = stride)(x)\n        if bn_relu:\n            x = BatchNormalization()(x)\n            x = ReLU()(x)\n        return x\n    \n    def conv_block(self, inp, kernels, downsample = False):\n        '''\n            Basic convolution block vgg styles\n        '''\n        x = inp\n        if downsample:\n            x = MaxPooling3D(pool_size=(2, 2, 2))(x)\n        x = self.conv_in_relu(x, kernels)\n        x = self.conv_in_relu(x, kernels)\n        return x\n\n    def up_conv_block(self, inp, kernels, connect):\n        x = UpSampling3D()(inp)\n        x = Concatenate(axis=-1)([x, connect]) # Skip connection\n        x = self.conv_in_relu(x, kernels)\n        x = self.conv_in_relu(x, kernels)\n        return x\n\n    def __call__(self):\n        conv_settings = self.conv_settings\n        \n        num_blocks = len(conv_settings)\n        \n        inp = Input(self.input_shape)\n\n        outputs = []\n\n        encoder_blocks = []\n\n        # Encoder\n        conv = BatchNormalization()(inp)\n        \n        for i in range(0, num_blocks - 1):\n            if i == 0:\n                conv = self.conv_block(conv, conv_settings[i])\n                encoder_blocks.append(conv)\n            else:\n                conv = self.conv_block(conv, conv_settings[i], downsample = True)\n                encoder_blocks.append(conv)\n                \n        out = self.conv_block(conv, conv_settings[-1], downsample = True)\n\n        # Decoder\n        for i in range(num_blocks - 1, 0, -1):\n            out = self.up_conv_block(out, conv_settings[i - 1], encoder_blocks[i - 1])\n            if self.deep_supervision and 1 < i < 4:\n                pool_size = 2**(i - 1)\n                pred = Conv3D(self.num_classes, \n                            kernel_size = (1, 1, 1), \n                            padding = 'same')(out)\n                pred = UpSampling3D(pool_size)(pred)\n                pred = Activation(self.activation, name = f'output_{i}')(pred)\n                outputs.append(pred)\n\n        out = Conv3D(self.num_classes, \n                    kernel_size = (1, 1, 1), \n                    padding = 'same')(out)\n        out = Activation(self.activation)(out)\n        outputs.append(out)\n\n        model = Model(inputs = inp, outputs = outputs)\n        return model","metadata":{"execution":{"iopub.status.busy":"2022-07-14T07:09:15.429084Z","iopub.execute_input":"2022-07-14T07:09:15.429865Z","iopub.status.idle":"2022-07-14T07:09:15.454457Z","shell.execute_reply.started":"2022-07-14T07:09:15.429816Z","shell.execute_reply":"2022-07-14T07:09:15.453117Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class ResUnet3D:\n    def __init__(self, num_classes = 3, \n                 input_shape = (None, None, None, 1),\n                 conv_settings = [32, 64, 128, 256, 320], \n                 repeat = [2, 2, 2, 2, 2],\n                 deep_supervision = False,\n                 activation = 'sigmoid'):\n        self.input_shape = input_shape\n        self.num_classes = num_classes\n        self.conv_settings = conv_settings\n        self.repeat = repeat\n        self.deep_supervision = deep_supervision\n        self.activation = activation\n    \n    def conv_bn_relu(self, inp, kernels, kernel_size = 3, stride = 1, bn_relu = True):\n        x = inp\n        x = Conv3D(kernels, \n                    kernel_size = kernel_size,\n                    padding = 'same',\n                    strides = stride)(x)\n        if bn_relu:\n            x = BatchNormalization()(x)\n            x = ReLU()(x)\n        return x\n    \n    def residual_conv_block(self, inp, kernels, repeat = 2, downsample = False):\n        '''\n            Residual convolution block\n        '''\n        x = inp\n        for i in range(repeat):\n            stride = 1\n            if i == 0 and downsample:\n                stride = 2\n\n            skip_conn = self.conv_bn_relu(x, kernels, 1, stride, False)\n\n            x = self.conv_bn_relu(x, kernels, 3, stride)\n            x = self.conv_bn_relu(x, kernels, bn_relu=False)\n            \n            # Residual connection\n            x = Add()([x, skip_conn])\n            x = BatchNormalization()(x)\n            x = ReLU()(x)\n        return x\n\n\n    def up_conv_block(self, inp, kernels, connect):\n        '''\n            Decoder block\n        '''\n        x = UpSampling3D()(inp)\n        x = Concatenate(axis=-1)([x, connect]) # Skip connection\n        x = self.residual_conv_block(x, kernels)\n        return x\n\n    def __call__(self):\n        conv_settings = self.conv_settings\n        \n        num_blocks = len(conv_settings)\n        \n        inp = Input(self.input_shape)\n\n        outputs = []\n\n        encoder_blocks = []\n\n        # Encoder\n        conv = BatchNormalization()(inp)\n        \n        for i in range(0, num_blocks - 1):\n            if i == 0:\n                conv = self.residual_conv_block(conv, conv_settings[i], self.repeat[i])\n                encoder_blocks.append(conv)\n            else:\n                conv = self.residual_conv_block(conv, conv_settings[i], self.repeat[i], downsample = True)\n                encoder_blocks.append(conv)\n                \n        out = self.residual_conv_block(conv, conv_settings[-1], self.repeat[-1], downsample = True)\n\n        # Decoder\n        for i in range(num_blocks - 1, 0, -1):\n            out = self.up_conv_block(out, conv_settings[i - 1], encoder_blocks[i - 1])\n            if self.deep_supervision and 1 < i < 4:\n                pool_size = 2**(i - 1)\n                pred = Conv3D(self.num_classes, \n                            kernel_size = (1, 1, 1), \n                            padding = 'same')(out)\n                pred = UpSampling3D(pool_size)(pred)\n                pred = Activation(self.activation, name = f'output_{i}')(pred)\n                outputs.append(pred)\n\n        out = Conv3D(self.num_classes, \n                     kernel_size = (1, 1, 1), \n                     padding = 'same')(out)\n        out = Activation(self.activation, name=f'output_final')(out)\n        outputs.append(out)\n        \n        model = Model(inputs = inp, outputs = outputs)\n        return model","metadata":{"execution":{"iopub.status.busy":"2022-07-14T07:09:15.457812Z","iopub.execute_input":"2022-07-14T07:09:15.458207Z","iopub.status.idle":"2022-07-14T07:09:15.485241Z","shell.execute_reply.started":"2022-07-14T07:09:15.458154Z","shell.execute_reply":"2022-07-14T07:09:15.484104Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\n    MODEL\n'''\nmodels = [ResUnet3D(conv_settings = [64, 128, 256, 512, 1024], input_shape = (None, None, None, 1), deep_supervision=True)(),\n          Unet3D(conv_settings = [64, 128, 256, 512, 1024], input_shape = (None, None, None, 1), deep_supervision=True)(),\n          ResUnet3D(input_shape = (None, None, None, 1), deep_supervision=True)(),\n          Unet3D(input_shape = (None, None, None, 1), deep_supervision=True)(),]\nfor i in range(len(models)):\n    print(MODEL_WEIGHT[i])\n    models[i].load_weights(MODEL_WEIGHT[i])\n    models[i].trainable = False","metadata":{"execution":{"iopub.status.busy":"2022-07-14T07:25:20.193004Z","iopub.execute_input":"2022-07-14T07:25:20.193354Z","iopub.status.idle":"2022-07-14T07:25:37.154690Z","shell.execute_reply.started":"2022-07-14T07:25:20.193321Z","shell.execute_reply":"2022-07-14T07:25:37.153541Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\n    SLIDING WINDOW INFERENCE\n'''\ndef sliding_window_inference(volume, models, window_size=WINDOW_SIZE, stride=STRIDE, deep_supervision=True):\n    '''\n        Sliding window inference\n        --------\n        volume : numpy.ndarray\n            Input 3D volume with shape = (1, Depth, Height, Width, 1)\n        model : tf.keras.Model\n            Inference model output with shape = (1, Depth, Height, Width, 3)\n\n        --------\n        return : numpy.ndarray\n            Output segmentation\n    '''\n    d, h, w = volume.shape[1:4]\n    w_d, w_h, w_w = window_size\n    s_d, s_h, s_w = stride\n    result = np.zeros((*volume.shape[:4], NUM_CLASSES), dtype='float32')\n    overlap = np.zeros((*volume.shape[:4], NUM_CLASSES), dtype='float32')\n    starting_points = [(x, y, z) for x in set( list(range(0, d - w_d, s_d)) + [d - w_d] ) \n                                 for y in set( list(range(0, h - w_h, s_h)) + [h - w_h] ) \n                                 for z in set( list(range(0, w - w_w, s_w)) + [w - w_w] )]\n\n    patches = np.empty((1, *window_size, 1), dtype='float32')\n    for i, (x, y, z) in enumerate(starting_points):\n        patches[0] = volume[0, x:x + w_d, y:y + w_h, z:z + w_w, :]\n\n        for model in models:\n            y_pred = model.predict(patches, batch_size = 1)\n            if deep_supervision:\n                y_pred = y_pred[-1]\n\n            result[:, x:x + w_d, y:y + w_h, z:z + w_w] += y_pred[0] * gaussian_importance_map\n            overlap[:, x:x + w_d, y:y + w_h, z:z + w_w] += gaussian_importance_map\n\n    assert np.sum(overlap == 0.) == 0, \"Sliding window does not cover all volume\"\n\n    return result / overlap","metadata":{"execution":{"iopub.status.busy":"2022-07-14T07:26:03.680281Z","iopub.execute_input":"2022-07-14T07:26:03.680615Z","iopub.status.idle":"2022-07-14T07:26:03.696865Z","shell.execute_reply.started":"2022-07-14T07:26:03.680581Z","shell.execute_reply":"2022-07-14T07:26:03.695849Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\n    INFERENCE\n'''\nfrom timeit import default_timer as timer\nimport tqdm\n# Get scans list\nif Debug:\n    submission = submission[submission.columns][:432]\n    submission[\"predicted\"] = \"\"\n\ndf_idx = 0\npaths = [submission[\"file_path\"][i] for i in range(0, len(submission[\"file_path\"]), 3)]\ndata_loader = DataLoader(paths)\nfor image, original_shape, image_path in tqdm.tqdm(data_loader):\n    gc.collect()\n    predictions = sliding_window_inference(image, models)\n    predictions = np.squeeze(predictions, axis=0)\n    predictions = (predictions > 0.5).astype('uint8')\n    \n    # Write prediction to csv\n    for i, prediction in enumerate(predictions):\n        for j in range(3):\n            scan_id = get_id_from_path(image_path[i])\n            class_id = class_map[j]\n            submission.loc[df_idx + j, \"class\"] = class_id\n            if np.sum(prediction[:,:,j]) < min_mask_area[j]:\n                continue\n            rle = rle_encoding(prediction[:, :, j])\n            submission.loc[df_idx + j, \"predicted\"] = rle \n            \n        df_idx += 3\n        \n    gc.collect()","metadata":{"execution":{"iopub.status.busy":"2022-07-14T07:26:06.608795Z","iopub.execute_input":"2022-07-14T07:26:06.609101Z","iopub.status.idle":"2022-07-14T07:28:12.592278Z","shell.execute_reply.started":"2022-07-14T07:26:06.609068Z","shell.execute_reply":"2022-07-14T07:28:12.591117Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_idx)\nprint(submission)","metadata":{"execution":{"iopub.status.busy":"2022-07-14T07:21:30.832802Z","iopub.execute_input":"2022-07-14T07:21:30.833130Z","iopub.status.idle":"2022-07-14T07:21:30.847712Z","shell.execute_reply.started":"2022-07-14T07:21:30.833095Z","shell.execute_reply":"2022-07-14T07:21:30.846372Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"'''\n    CREATE SUBMISSION FILE\n'''\nsubmission = submission[[\"id\", \"class\", \"predicted\"]]\nsubmission.to_csv(\"submission.csv\", index=False)\nsubmission.head()","metadata":{"execution":{"iopub.status.busy":"2022-07-14T07:10:26.272194Z","iopub.status.idle":"2022-07-14T07:10:26.273148Z","shell.execute_reply.started":"2022-07-14T07:10:26.272793Z","shell.execute_reply":"2022-07-14T07:10:26.272829Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}