{"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":"markdown","source":"# CPTAC-PDA  - deepflash2 Cross Validation\n\n> Cross validation for models trained with efficient region based sampling. \n\n***\n\n## Overview\n\n1. Installation and package loading\n2. Functions and classes for prediction\n3. Configuration\n4. Prediction\n5. Submission\n\n#### Related Kernels\n\n- Train Notebook: \n- Sampling Notebook: ","metadata":{"id":"OcsetTMwKXqC"}},{"cell_type":"markdown","source":"### Installation and package loading","metadata":{}},{"cell_type":"code","source":"# Install deepflash2 and dependencies\nimport sys\nsys.path.append(\"../input/segmentation-models-pytorch-install\")\n!pip install -q ../input/deepflash2-lfs\nimport cv2, torch, gc, zarr\nimport torch.nn.functional as F\nimport deepflash2.tta as tta\nimport matplotlib.pyplot as plt\nimport pandas as pd, numpy as np\nimport segmentation_models_pytorch as smp\nfrom pathlib import Path\nfrom torch.utils.data import Dataset, DataLoader\nfrom scipy .ndimage.filters import gaussian_filter\nfrom tqdm.notebook import tqdm\nfrom sklearn.model_selection import KFold\nfrom fastcore.all import *\nfrom skimage.measure import regionprops_table\nfrom deepflash2.utils import iou\nfrom skimage.segmentation import clear_border\nfrom scipy.ndimage.morphology import binary_fill_holes\nimport warnings\nwarnings.filterwarnings(\"ignore\")","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Functions and classes for prediction","metadata":{}},{"cell_type":"code","source":"def load_model_weights(model, file, strict=True):\n    state = torch.load(file, map_location='cpu')\n    stats = state['stats']\n    model_state = state['model']\n    model.load_state_dict(model_state, strict=strict)\n    return model, stats\n\n# from https://github.com/MIC-DKFZ/nnUNet/blob/2fade8f32607220f8598544f0d5b5e5fa73768e5/nnunet/network_architecture/neural_network.py#L250\ndef _get_gaussian(patch_size, sigma_scale=1. / 8) -> np.ndarray:\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    return gaussian_importance_map\n\ndef postproc(pred):\n    pred = (pred*255).astype('uint8')\n    _ ,msk = cv2.threshold(pred,0,1,cv2.THRESH_BINARY+cv2.THRESH_OTSU)  \n    msk = clear_border(msk)\n    msk = binary_fill_holes(msk)\n    return msk\n\ndef expand_slice(slices, n_pix=100):\n    slc_list = []\n    for s in slices:\n        slc_list.append(slice(max(s.start-n_pix,0), s.stop+n_pix))\n    return tuple(slc_list)","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Some code adapted from https://www.kaggle.com/iafoss/hubmap-pytorch-fast-ai-starter-sub\nclass HubmapDataset(Dataset):\n    'HubmapDataset class that does not load the full tiff files.'\n    def __init__(self, file, stats, scale=3, shift=.8, output_shape=(512,512), s_th = 40):\n        \n        self.mean, self.std = stats\n        self.scale = scale\n        self.shift = shift\n        self.output_shape = output_shape\n        self.input_shape = tuple(int(t*scale) for t in self.output_shape)      \n        self.s_th = s_th #saturation blancking threshold\n        self.p_th = 1000*(self.output_shape[0]//256)**2 #threshold for the minimum number of pixels\n        self.data = zarr.open(file.as_posix())\n            \n        # Tiling\n        self.slices = []\n        self.out_slices = []\n        self.out_data_shape = tuple(int(x//self.scale) for x in self.data.shape[:-1])\n        start_points = [o//2 for o in self.output_shape]\n        end_points = [(s - st) for s, st in zip(self.out_data_shape, start_points)]\n        n_points = [int(s//(o*self.shift))+1 for s, o in zip(self.out_data_shape, self.output_shape)]\n        center_points = [np.linspace(st, e, num=n, endpoint=True, dtype=np.int64) for st, e, n in zip(start_points, end_points, n_points)]\n        for cx in center_points[1]:\n            for cy in center_points[0]:\n                # Calculate output slices for whole image\n                slices = tuple(slice(int((c*self.scale - o/2).clip(0, s)), int((c*self.scale + o/2).clip(max=s)))\n                                 for (c, o, s) in zip((cy, cx), self.input_shape, self.data.shape))\n                self.slices.append(slices)\n                \n                out_slices = tuple(slice(int((c - o/2).clip(0, s)), int((c + o/2).clip(max=s)))\n                                 for (c, o, s) in zip((cy, cx), self.output_shape, self.out_data_shape))\n                self.out_slices.append(out_slices)\n                \n\n    def __len__(self):\n        return len(self.slices)\n\n    def __getitem__(self, idx):\n        if torch.is_tensor(idx):\n            idx = idx.tolist()\n            \n        slices = self.slices[idx]\n        img = self.data[slices]\n\n        if self.scale!=1:\n            img = cv2.resize(img, self.output_shape, interpolation = cv2.INTER_AREA)\n        \n        #check for empty imges\n        hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)\n        h,s,v = cv2.split(hsv)\n        if (s>self.s_th).sum() <= self.p_th or img.sum() <= self.p_th:\n            # Remove if idx=-1\n            idx = -1\n        \n        img = (img/255.0 - self.mean)/self.std\n        img = img.transpose(2, 0, 1).astype('float32')\n        \n        return torch.from_numpy(img), idx\n    \nclass Model_pred:\n    'Class for prediction with multiple models'\n    def __init__(self, models, use_tta=True, batch_size=32, energy_T=1):\n        self.model = model\n        self.bs = batch_size\n        self.energy_T = energy_T\n        self.tfms = [tta.HorizontalFlip(), tta.VerticalFlip()] if use_tta else []\n        \n    def predict(self, ds):\n        dl = DataLoader(ds, self.bs, num_workers=0, shuffle=False, pin_memory=True)\n        \n        # Create zero arrays\n        softmax = np.zeros(ds.out_data_shape, dtype='float32')\n        energy = np.zeros(ds.out_data_shape, dtype='float32')\n        st_deviation = np.zeros(ds.out_data_shape, dtype='float32')\n        merge_map = np.zeros(ds.out_data_shape, dtype='float32')\n        \n        # Gaussian weights\n        gw_numpy = _get_gaussian(ds.output_shape)\n        gw = torch.from_numpy(gw_numpy).to(device)\n        gw = gw.view(1,*gw.shape)\n        \n        with torch.no_grad():\n            for images, idxs in tqdm(iter(dl), total=len(dl)):\n                if ((idxs>=0).sum() > 0): #exclude empty images\n                    images = images[idxs>=0].to(device)\n                    idxs = idxs[idxs>=0]\n                    m_softmax, m_energy = tta.Merger(), tta.Merger()\n                    for t in tta.Compose(self.tfms):\n                        aug_images = t.augment_image(images)\n                        out = self.model(aug_images)\n                        out = t.deaugment_mask(out)\n                        m_softmax.append(F.softmax(out, dim=1))\n                        # https://proceedings.neurips.cc/paper/2020/file/f5496252609c43eb8a3d147ab9b9c006-Paper.pdf\n                        e = (self.energy_T*torch.logsumexp(out/self.energy_T, dim=1)) #negative energy score\n                        m_energy.append(e)\n            \n                    # Merge preds and apply gaussian weigthing \n                    merge_ll = []\n                    batch_smx = m_softmax.result()*gw.view(1,*gw.shape)\n                    merge_ll.append([x for x in batch_smx.permute(0,2,3,1).cpu().numpy()])\n                    \n                    batch_std = torch.mean(m_softmax.result('std'), 1)*gw\n                    merge_ll.append([x for x in batch_std.cpu().numpy()])\n                    \n                    batch_energy =  m_energy.result()*gw\n                    merge_ll.append([x for x in batch_energy.cpu().numpy()])\n\n                    for idx, smx, std, eng  in zip(idxs, *merge_ll):\n                        slcs = ds.out_slices[idx]\n                        # Only using positive class here\n                        softmax[slcs] += smx[...,1]\n                        energy[slcs] += eng\n                        st_deviation[slcs] += std\n                        merge_map[slcs] += gw_numpy\n\n        softmax /= merge_map\n        st_deviation /= merge_map\n        energy /= merge_map\n        return softmax, st_deviation, energy","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Configuration","metadata":{}},{"cell_type":"code","source":"class CONFIG():\n    \n    # data paths\n    path = Path('../input/cptacpda')\n    data_path = Path('../input/cptac-pda-to-zarr/images_scale2')\n    annotations_path = Path('../input/cptac-pda-pancreas-efficient-sampling-deepflash2/masks_scale2')\n    model_path = Path('../input/cptac-pda-train-deepflash2')\n    \n    # scale factor, data is already downscaled to 2, so absulute downscale with 1.5 is 3\n    scale = 1.5 \n    # tile shift for prediction\n    shift = 0.8 \n    tile_shape = (512, 512)\n\n    # pytorch model (https://github.com/qubvel/segmentation_models.pytorch)\n    encoder_name = \"efficientnet-b2\"\n    encoder_weights = None\n    in_channels = 3\n    classes = 2\n    \n    # dataloader \n    batch_size = 32\n    \n    # test time augmentation\n    tta = True\n    # prediction threshold\n    threshold = 0.5\n    \n    # Postprocessing\n    min_pixel = 1000 # needs to be further evaluated \n    area_threshold = 10000\n    \ncfg = CONFIG()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Train annotations for ids\ndf_train = pd.read_csv(cfg.path/'dataset_information.csv')\nfiles = L([cfg.data_path/x for x in df_train.image_file])\n\ng_mask = zarr.open((cfg.annotations_path/'labels').as_posix())\nroot = zarr.group('cv_preds', overwrite=True)\ng_pred, g_smx, g_std, g_eng = root.create_groups('pred', 'smx', 'std', 'eng', overwrite=True)\n\ndevice = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n\n# Models (see https://github.com/qubvel/segmentation_models.pytorch)\nMODELS = sorted([f for f in cfg.model_path.iterdir() if f.suffix=='.pth'])\nprint(f'Found {len(MODELS)} models', *MODELS)\n\nmodels = []\nfor i, m_path in enumerate(MODELS):\n    #state_dict = torch.load(path,map_location=torch.device('cpu'))\n    model = smp.Unet(encoder_name=cfg.encoder_name, \n                     encoder_weights=cfg.encoder_weights, \n                     in_channels=cfg.in_channels, \n                     classes=cfg.classes)\n    model, stats = load_model_weights(model, m_path)\n    model.float()\n    model.eval()\n    model.to(device)\n    models.append(model)\nassert len(models)==len(MODELS)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Cross validation\n\n- Calculation metrics for predictions from k-fold cross validation\n- Plotting some expamples","metadata":{}},{"cell_type":"code","source":"def plot_selection(row):\n    fig, ax = plt.subplots(ncols=7, figsize=(30,5))\n    slc = expand_slice(row.slice, 20)\n    \n    # Image\n    ax[0].imshow(ds.data[slc])\n    ax[0].set_title(row.idx)\n    \n    # Annotation\n    ax[1].imshow(g_mask[row.idx][slc])\n    ax[1].set_title('Annotation')\n    \n    # Softmax\n    ax[2].imshow(g_smx[row.idx][slc].astype(np.single))\n    ax[2].set_title('Softmax')\n    \n    # Threshold\n    ax[3].imshow(g_smx[row.idx][slc]>0.5)\n    ax[3].set_title(f'Threshold at 0.5')\n    \n    # Postproc\n    ax[4].imshow(g_pred[row.idx][slc])\n    ax[4].set_title(f'Postprocessing \\n IoU {row.iou:.2f}')\n    \n    # Energy\n    ax[5].imshow(g_eng[row.idx][slc].astype(np.single))\n    ax[5].set_title(f'Energy \\n Mean {row.mean_intensity:.2f}')\n    \n    # Energy\n    ax[6].imshow(g_std[row.idx][slc].astype(np.single))\n    ax[6].set_title(f'TTA Standard Deviation \\n Mean {row.mean_intensity_std:.2f}')\n    \n    plt.show()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Use validation split from training\nkf = KFold(len(MODELS), shuffle=True, random_state=42)\nmetric_list = []\nrecall_list = []\nprecision_list = []\nPLOT_OVERVIEW = True\n\nfor i, (train_idx, val_idx) in enumerate(kf.split(files)):\n    print('Split', i)\n    files_val = files[val_idx]\n    mp = Model_pred([models[i]], use_tta=cfg.tta, batch_size=cfg.batch_size)\n    print('Using model', MODELS[i])\n    \n    for f in files_val:\n        print('Validating', f.name)\n        \n        ds = HubmapDataset(f, stats, scale=cfg.scale, shift=cfg.shift, output_shape=cfg.tile_shape)\n    \n        # print('Predicting...')   \n        softmax, st_deviation, energy = mp.predict(ds)\n        \n        #print('Rezising...')\n        msk = g_mask[f.name][:].astype('uint8')\n        shape = msk.shape\n        g_smx[f.name] = cv2.resize(softmax, (shape[1], shape[0])).astype(np.half)\n        del softmax\n        gc.collect()\n        \n        g_std[f.name] = cv2.resize(st_deviation, (shape[1], shape[0])).astype(np.half)\n        del st_deviation\n        gc.collect\n        \n        g_eng[f.name] = cv2.resize(energy, (shape[1], shape[0])).astype(np.half)\n        del energy\n        gc.collect()\n        \n        # Counting prediction\n        pred = (g_smx[f.name][:]>cfg.threshold).astype('uint8')\n        _, comps = cv2.connectedComponents(pred, connectivity=4)\n        props = regionprops_table(comps, properties=('area', 'slice'))\n        df_props = pd.DataFrame(props)\n\n        # Postprocessing\n        pred_pp = np.zeros_like(pred)\n        for _, row in df_props.iterrows():\n            slc = row.slice\n            if row.area>cfg.area_threshold:\n                tile = binary_fill_holes(pred[slc])\n                pred_pp[slc][tile>0] = 1\n            else:\n                if any([(s.stop-s.start)<100 for s in slc]):\n                    slc = expand_slice(slc, 10)\n                tile = 0\n                while np.sum(tile)==0:\n                    slc = expand_slice(slc, 10)\n                    tile = postproc(g_smx[f.name][slc])\n                if np.sum(tile)>2000:\n                    pred_pp[slc][tile>0] = 1\n        \n        # Recounting prediction\n        g_pred[f.name] = pred_pp\n        _, comps = cv2.connectedComponents(pred_pp, connectivity=4)\n        props_energy = regionprops_table(comps, g_eng[f.name][:],  properties=('mean_intensity', 'area', 'bbox', 'solidity', 'slice'))\n        props_std = regionprops_table(comps, g_std[f.name][:],  properties=('mean_intensity',))\n        df_props =  pd.DataFrame(props_energy).join(pd.DataFrame(props_std)[['mean_intensity']], rsuffix='_std')\n        \n        for _, row in df_props.iterrows():\n            df_props.loc[_, 'iou'] = iou(msk[row.slice], pred_pp[row.slice])\n        \n        df_props['idx'] = f.name\n        df_props['fold'] = i\n        df_props['model'] = MODELS[i]\n        \n        precision_list.append(df_props)\n        \n        # Mask vs. Prediction\n        _, comps_msk = cv2.connectedComponents(msk, connectivity=4)\n        props_msk = regionprops_table(comps_msk, properties=('area', 'bbox', 'slice'))\n        df_msk =  pd.DataFrame(props_msk)\n        \n        for _, row in df_msk.iterrows():\n            df_msk.loc[_, 'iou'] = iou(msk[row.slice], pred_pp[row.slice])\n        \n        df_msk['idx'] = f.name\n        df_msk['fold'] = i\n        df_msk['model'] = MODELS[i]\n        recall_list.append(df_props)\n        \n        # Calc prediction metrics \n        iou_score_pp = iou(pred_pp,msk)\n        dice_score_pp = 2*iou_score_pp/(iou_score_pp+1)\n        \n        iou_score = iou(pred,msk)\n        dice_score = 2*iou_score/(iou_score+1)\n        \n        precision =  len(df_props[df_props.iou>.5])/len(df_props)\n        recall = len(df_msk[df_msk.iou>.5])/len(df_msk)\n        f1_score = 2 * (precision * recall) / (precision + recall)\n        \n        tmp = pd.Series({'idx' : f.name,\n                         'fold' : i,\n                         'model': MODELS[i],\n                         'iou': iou_score,\n                         'dice': dice_score,\n                         'iou_pp': iou_score_pp,\n                         'dice_pp': dice_score_pp,\n                         'precision' : precision, \n                         'recall' : recall,\n                         'f1_score' : f1_score})\n        \n        metric_list.append(tmp)\n        \n        del msk, pred, pred_pp\n        gc.collect()\n        \n        if PLOT_OVERVIEW:\n            #print('#######################################')\n            print('Example with highest energy')\n            row = df_props.sort_values('mean_intensity', ascending=False).iloc[0]\n            plot_selection(row)\n            \n            #print('#######################################')\n            try: \n                row = df_props[df_props.iou<0.5].sort_values('mean_intensity', ascending=False).iloc[0]\n                print('Higher energy but no annotation')\n                plot_selection(row)\n            except:\n                pass\n            \n            #print('#######################################')\n            print('Difficult example:  Low energy with annotation')\n            row = df_props[df_props.iou>0.5].sort_values('mean_intensity', ascending=True).iloc[0]\n            plot_selection(row)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.DataFrame(metric_list).to_csv('cv_metrics.csv', index=False)\npd.concat(precision_list).to_csv('cv_precision.csv', index=False)\npd.concat(recall_list).to_csv('cv_recall.csv', index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}