{"cells":[{"metadata":{},"cell_type":"markdown","source":"It aggregates several public tile strategies and provides pipeline for finding optimal SIZE, COLS, ROWS params form the human perspective. It also shows how to save the results.\n\nderives from\n* https://www.kaggle.com/iafoss/panda-16x128x128-tiles\n* https://www.kaggle.com/akensert/panda-optimized-tiling-tf-data-dataset\n* https://www.kaggle.com/debanga/let-s-enhance-the-images\n* https://www.kaggle.com/harupy/visualization-panda-16x128x128-tiles\n* https://www.kaggle.com/raghaw/panda-medium-resolution-dataset-25x256x256\n\noutput dataset examples\n* 6x6 256 lafoss_tiles https://www.kaggle.com/dlarionov/prostate-cancer-tiles-medium\n* 4x4 312 akensert_tiles https://www.kaggle.com/dlarionov/akensert-1-312-4x4 (this kernel result)","execution_count":null},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport warnings\nimport math\nfrom multiprocessing import Pool\nfrom tqdm.notebook import tqdm\nimport matplotlib.pyplot as plt\nfrom PIL import Image\nfrom openslide import OpenSlide\nimport cv2\n\nIMG_DIR = '../input/prostate-cancer-grade-assessment/train_images/'\n\nLAYER = 1 # medium\nSIZE = 312\nCOLS = 4\nROWS = 4\nN = COLS*ROWS","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.read_csv('../input/prostate-cancer-grade-assessment/train.csv')\ndf","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def split_tiles(img:np.ndarray)->np.ndarray:\n    reshaped = img.reshape(\n        img.shape[0] // SIZE,\n        SIZE,\n        img.shape[1] // SIZE,\n        SIZE,\n        3,\n    )\n    transposed = reshaped.transpose(0, 2, 1, 3, 4)\n    return transposed.reshape(-1, SIZE, SIZE, 3)\n\ndef join_tiles(img:np.ndarray)->np.ndarray:\n    reshaped = img.reshape(\n        COLS,\n        ROWS,    \n        img.shape[1],\n        img.shape[2],\n        3\n    )\n    transposed = reshaped.transpose(0, 2, 1, 3, 4)\n    return transposed.reshape(COLS * SIZE, ROWS * SIZE, 3)\n\ndef lafoss_tiles(img:np.ndarray)->np.ndarray:\n    \n    # calculate paddings\n    H, W, _ = img.shape\n    pad_w = (SIZE - W % SIZE) % SIZE\n    pad_h = (SIZE - H % SIZE) % SIZE\n    \n    # implement padding\n    padded = np.pad(\n        img,\n        [[pad_h // 2, pad_h - pad_h // 2],\n         [pad_w // 2, pad_w - pad_w // 2],\n         [0, 0]],\n        constant_values=255, # 255 - white\n    )\n    \n    # split image into tiles\n    tiles = split_tiles(padded)\n    \n    # calculate sums of all pixsels for each tile\n    sums = tiles.reshape(tiles.shape[0], -1).sum(axis=-1)\n    \n    # take top N tiles by minimum sum value\n    idxs_selected = np.argsort(sums)[:N]\n    selected = tiles[idxs_selected]\n    \n    # append white tiles if necessary\n    if len(selected)<N:\n        selected = np.pad(\n            selected,\n            [[0,N-len(selected)],[0,0],[0,0],[0,0]],\n            constant_values=255\n        )\n    \n    # join selected tiles into one image\n    merged = join_tiles(selected)\n\n    return merged\n\ndef akensert_tiles(img:np.ndarray, debug=False)->np.ndarray:    \n    \n    # get tile coords\n    img, coords = compute_coords(\n        img,\n        patch_size=SIZE,\n        precompute=False, # returns new padded img\n        min_patch_info=0.35,\n        min_axis_info=0.35,\n        min_consec_axis_info=0.35,\n        min_decimal_keep=0.7)\n    \n    # sort coords (high info -> low info)\n    coords = sorted(coords, key= lambda x: x[0], reverse=False)\n    \n    # select top N tiles\n    tiles = []\n    for i in range(len(coords)):\n        if i == N:\n            break;\n        _, x, y = coords[i]\n        tiles.append(img[x:x+SIZE,y:y+SIZE])\n    \n    # append white tiles if necessary\n    selected = np.array(tiles)\n    if len(selected)<N:\n        selected = np.pad(\n            selected,\n            [[0,N-len(selected)],[0,0],[0,0],[0,0]],\n            constant_values=255\n        )\n    \n    # merge tiles to one image\n    merged = join_tiles(selected)\n    \n    if debug:\n        for (v, y, x) in coords:\n            img = cv2.rectangle(img, (x, y), (x+SIZE, y+SIZE), color=(0, 0, 0), thickness=5)\n            img = cv2.circle(img, (x, y), radius=5, color=(255, 0, 0), thickness=-1)\n            img = cv2.circle(img, (x+SIZE, y+SIZE), radius=5, color=(0, 255, 0), thickness=-1)\n        return merged, img\n    else:\n        return merged","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# copypaste https://www.kaggle.com/akensert/panda-optimized-tiling-tf-data-dataset\n# copypaste https://www.kaggle.com/debanga/let-s-enhance-the-images\n\ndef enhance_image(image, contrast=1, brightness=15):\n    \"\"\"\n    Enhance constrast and brightness of images\n    \"\"\"\n    img_enhanced = cv2.addWeighted(image, contrast, image, 0, brightness)\n    return img_enhanced\n\ndef unsharp_masking(img):\n    \"\"\" Unsharp masking of an RGB image\"\"\"\n    img_gaussian = cv2.GaussianBlur(img, (21,21), 10.0)\n    return cv2.addWeighted(img, 1.8, img_gaussian, -0.8, 0, img)\n\ndef _mask_tissue(image, kernel_size=(7, 7), gray_threshold=220):\n    \"\"\"Masks tissue in image. Uses gray-scaled image, as well as\n    dilation kernels and 'gap filling'\n    \"\"\"\n    # Define elliptic kernel\n    kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, kernel_size)\n    # Convert rgb to gray scale for easier masking\n    gray = cv2.cvtColor(image, cv2.COLOR_RGB2GRAY)\n    # Now mask the gray-scaled image (capturing tissue in biopsy)\n    mask = np.where(gray < gray_threshold, 1, 0).astype(np.uint8)\n    # Use dilation and findContours to fill in gaps/holes in masked tissue\n    mask = cv2.dilate(mask, kernel, iterations=1)\n    contour, _ = cv2.findContours(mask, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_SIMPLE)\n    for cnt in contour:\n        cv2.drawContours(mask, [cnt], 0, 1, -1)\n    return mask\n\ndef _pad_image(image, pad_len, pad_val):\n    \"\"\"Pads inputted image, accepts both \n    2-d (mask) and 3-d (rgb image) arrays\n    \"\"\"\n    if image is None:\n        return None\n    elif image.ndim == 2:\n        return np.pad(\n            image, ((pad_len, pad_len), (pad_len, pad_len)), pad_val)\n    elif image.ndim == 3:\n        return np.pad(\n            image, ((pad_len, pad_len), (pad_len, pad_len), (0, 0)), pad_val)\n    return None\n\ndef _transpose_image(image):\n    \"\"\"Inputs an image and transposes it, accepts \n    both 2-d (mask) and 3-d (rgb image) arrays\n    \"\"\"\n    if image is None:\n        return None\n    elif image.ndim == 2:\n        return np.transpose(image, (1, 0)).copy()\n    elif image.ndim == 3:\n        return np.transpose(image, (1, 0, 2)).copy()\n    return None\n\ndef _get_tissue_parts_indices(tissue, min_consec_info):\n    \"\"\"If there are multiple tissue parts in 'tissue', 'tissue' will be \n    split. Each tissue part will be taken care of separately (later on), \n    and if the tissue part is less than min_consec_info, it's considered \n    to small and won't be returned.\n    \"\"\"\n    split_points = np.where(np.diff(tissue) != 1)[0]+1\n    tissue_parts = np.split(tissue, split_points)\n    return [\n        tp for tp in tissue_parts if len(tp) >= min_consec_info\n    ]\n\ndef _get_tissue_subparts_coords(subtissue, patch_size, min_decimal_keep):\n    \"\"\"Inputs a tissue part resulting from '_get_tissue_parts_indices'.\n    This tissue part is divided into N subparts and returned.\n    Argument min_decimal_keep basically decides if we should reduce the\n    N subparts to N-1 subparts, due to overflow.\n    \"\"\"\n    start, end = subtissue[0], subtissue[-1]\n    num_subparts = (end-start)/patch_size\n    if num_subparts % 1 < min_decimal_keep and num_subparts >= 1:\n        num_subparts = math.floor(num_subparts)\n    else:\n        num_subparts = math.ceil(num_subparts)\n\n    excess = (num_subparts*patch_size) - (end-start)\n    shift = excess // 2\n\n    return [\n        i * patch_size + start - shift \n        for i in range(num_subparts)\n    ]\n\ndef _eval_and_append_xy_coords(coords,\n                               image, \n                               mask, \n                               patch_size, \n                               x, y, \n                               min_patch_info,\n                               transposed,\n                               precompute):\n    \"\"\"Based on computed x and y coordinates of patch: \n    slices out patch from original image, flattens it,\n    preprocesses it, and finally evaluates its mask.\n    If patch contains more info than min_patch_info,\n    the patch coordinates are kept, along with a value \n    'val1' that estimates how much information there \n    is in the patch. Smaller 'val1' assumes more info.\n    \"\"\"\n    patch_1d = (\n        image[y: y+patch_size, x:x+patch_size, :]\n        .mean(axis=2)\n        .reshape(-1)\n    )\n    idx_tissue = np.where(patch_1d <= 210)[0]\n    idx_black = np.where(patch_1d < 5)[0]\n    idx_background = np.where(patch_1d > 210)[0]\n\n    if len(idx_tissue) > 0:\n        patch_1d[idx_black] = 210\n        patch_1d[idx_background] = 210\n        val1 = int(patch_1d.mean())\n        val2 = mask[y:y+patch_size, x:x+patch_size].mean()\n        if val2 > min_patch_info:\n            if precompute:\n                if transposed:\n                    coords = np.concatenate([\n                        coords, [[val1, x-patch_size, y-patch_size]]\n                    ])\n                else:\n                    coords = np.concatenate([\n                        coords, [[val1, y-patch_size, x-patch_size]]\n                    ])\n            else:\n                coords = np.concatenate([\n                    coords, [[val1, y, x]]\n                ])\n               \n    return coords\n\ndef compute_coords(image,\n                   patch_size=256,\n                   precompute=False,\n                   min_patch_info=0.35,\n                   min_axis_info=0.35,\n                   min_consec_axis_info=0.35,\n                   min_decimal_keep=0.7):\n\n    \"\"\"\n    Input:\n        image : 3-d np.ndarray\n        patch_size : size of patches/tiles, will be of \n            size (patch_size x patch_size x 3)\n        precompute : If True, only coordinates will be returned,\n            these coordinates match the inputted 'original' image.\n            If False, both an image and coordinates will be returned,\n            the coordinates does not match the inputted image but the\n            image that it is returned with.\n        min_patch_info : Minimum required information in patch\n            (see '_eval_and_append_xy_coords')\n        min_axis_info : Minimum fraction of on-bits in x/y dimension to be \n            considered enough information. For x, this would be fraction of \n            on-bits in x-dimension of a y:y+patch_size slice. For y, this would \n            be the fraction of on-bits for the whole image in y-dimension\n        min_consec_axis_info : Minimum consecutive x/y on-bits\n            (see '_get_tissue_parts_indices')\n        min_decimal_keep : Threshold for decimal point for removing \"excessive\" patch\n            (see '_get_tissue_subparts_coords')\n    \n    Output:\n        image [only if precompute is False] : similar to input image, but fits \n            to the computed coordinates\n        coords : the coordinates that will be used to compute the patches later on\n    \"\"\"\n    \n    \n    if type(image) != np.ndarray:\n        # if image is a Tensor\n        image = image.numpy()\n    \n    # masked tissue will be used to compute the coordinates\n    mask = _mask_tissue(image)\n\n    # initialize coordinate accumulator\n    coords = np.zeros([0, 3], dtype=int)\n\n    # pad image and mask to make sure no tissue is potentially missed out\n    image = _pad_image(image, patch_size, 'maximum')\n    mask = _pad_image(mask, patch_size, 'minimum')\n    \n    y_sum = mask.sum(axis=1)\n    x_sum = mask.sum(axis=0)\n    # if on bits in x_sum is greater than in y_sum, the tissue is\n    # likely aligned horizontally. The algorithm works better if\n    # the image is aligned vertically, thus the image will be transposed\n    if len(np.where(x_sum > 0)[0]) > len(np.where(y_sum > 0)[0]):\n        image = _transpose_image(image)\n        mask = _transpose_image(mask)\n        y_sum, _ = x_sum, y_sum\n        transposed = True\n    else:\n        transposed = False\n    \n    # where y_sum is more than the minimum number of on-bits\n    y_tissue = np.where(y_sum >= (patch_size*min_axis_info))[0]\n    \n    if len(y_tissue) < 1:\n        warnings.warn(\"Not enough tissue in image (y-dim)\", RuntimeWarning)\n        if precompute: return [(0, 0, 0)]\n        else: return image, [(0, 0, 0)]\n    \n    y_tissue_parts_indices = _get_tissue_parts_indices(\n        y_tissue, patch_size*min_consec_axis_info)\n    \n    if len(y_tissue_parts_indices) < 1: \n        warnings.warn(\"Not enough tissue in image (y-dim)\", RuntimeWarning)\n        if precompute: return [(0, 0, 0)]\n        else: return image, [(0, 0, 0)]\n    \n    # loop over the tissues in y-dimension\n    for yidx in y_tissue_parts_indices:\n        y_tissue_subparts_coords = _get_tissue_subparts_coords(\n            yidx, patch_size, min_decimal_keep)\n        \n        for y in y_tissue_subparts_coords:\n            # in y_slice, where x_slice_sum is more than the minimum number of on-bits\n            x_slice_sum = mask[y:y+patch_size, :].sum(axis=0)\n            x_tissue = np.where(x_slice_sum >= (patch_size*min_axis_info))[0]\n            \n            x_tissue_parts_indices = _get_tissue_parts_indices(\n                x_tissue, patch_size*min_consec_axis_info)\n            \n            # loop over tissues in x-dimension (inside y_slice 'y:y+patch_size')\n            for xidx in x_tissue_parts_indices:\n                x_tissue_subparts_coords = _get_tissue_subparts_coords(\n                    xidx, patch_size, min_decimal_keep)\n                \n                for x in x_tissue_subparts_coords:\n                    coords = _eval_and_append_xy_coords(\n                        coords, image, mask, patch_size, x, y, \n                        min_patch_info, transposed, precompute\n                    )     \n    \n    if len(coords) < 1:\n        warnings.warn(\"Not enough tissue in image (x-dim)\", RuntimeWarning)\n        if precompute: return [(0, 0, 0)]\n        else: return image, [(0, 0, 0)]\n    \n    if precompute: return coords\n    else: return image, coords","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def imread(path:str, layer:int)->Image:\n    if not os.path.exists(path):\n        return None\n\n    with OpenSlide(path) as slide:\n        im = slide.read_region((0,0), layer, slide.level_dimensions[layer])\n        im = im.convert('RGB') # drops A\n        return im","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for i in range(12):\n    f, axarr = plt.subplots(1, 3, figsize=(20,20))\n    idx = np.random.randint(0, len(df))        \n    row = df.loc[idx]\n    im = imread(f\"{IMG_DIR}{row['image_id']}.tiff\", layer=LAYER)        \n    arr = np.asarray(im)\n    arr = enhance_image(arr)\n    arr = unsharp_masking(arr)\n    tiles1, img_map = akensert_tiles(arr, debug=True)\n    tiles2 = lafoss_tiles(arr)    \n    axarr[0].set_title(f'id:{idx}')\n    axarr[0].imshow(img_map.squeeze())\n    axarr[1].set_title('akensert')\n    axarr[1].imshow(tiles1.squeeze())\n    axarr[2].set_title('lafoss')\n    axarr[2].imshow(tiles2.squeeze())        ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Any kernel provides its own file system where you can save up to 4Tb while kernel is alive. 5Gb restriction applies onlty to /kaggle/working directory.\nThen you can use https://github.com/Kaggle/kaggle-api to create new kaggle dataset.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"TILES_DIR = '/tmp/tiles'\n!mkdir -p {TILES_DIR}\n\ndef process_image(idx):    \n    row = df.loc[idx]\n    im = imread(os.path.join(IMG_DIR, f\"{row.image_id}.tiff\"), layer=LAYER)\n    im = np.asarray(im)\n    im = akensert_tiles(im)\n    im = Image.fromarray(im)\n    im.save(os.path.join(TILES_DIR, f\"{row.image_id}.jpg\"), format='JPEG', quality=90)\n\nbatch = df#.head(100)\nwith Pool(processes=4) as pool:\n    res = list(\n        tqdm(pool.imap(process_image, list(batch.index)), total = len(batch))\n    )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"DS_DIR = '/tmp/dataset/'\n!mkdir -p {DS_DIR}/tiles\n!tar czf {DS_DIR}/tiles/tiles.tar.gz -C {TILES_DIR} .","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"mykaggleapi is a private dataset. https://github.com/Kaggle/kaggle-api for more.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"!mkdir -p /root/.kaggle/\n!cp ../input/mykaggleapi/kaggle.json /root/.kaggle/\n!chmod 600 /root/.kaggle/kaggle.json","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!kaggle datasets init -p {DS_DIR}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import json\n\nwith open(f'{DS_DIR}/dataset-metadata.json', 'r+') as f:\n    data = json.load(f)\n    data['title'] = f'akensert {LAYER} {SIZE} {COLS}x{ROWS} '\n    data['id'] = f'dlarionov/akensert-{LAYER}-{SIZE}-{COLS}x{ROWS}'\n    f.seek(0)\n    json.dump(data, f, indent=4)\n    f.truncate()\n\n!cat {DS_DIR}/dataset-metadata.json","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!kaggle datasets create -p {DS_DIR} -q -r tar\n!rm -rf {TILES_DIR}\n!rm -rf {DS_DIR}","execution_count":null,"outputs":[]}],"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":4,"nbformat_minor":4}