PANDA_create_adv_tiles.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import cv2\n",
    "import skimage.io\n",
    "from tqdm.notebook import tqdm\n",
    "import numpy as np\n",
    "import math\n",
    "%matplotlib inline\n",
    "from matplotlib import pyplot as plt\n",
    "import warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "TRAIN = 'D:/ML/PANDA/data/train_images'\n",
    "OUT_TRAIN = 'data/train_adv_256x256'\n",
    "sz = 128\n",
    "\n",
    "os.makedirs(OUT_TRAIN, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#the functions below are borrowed from https://www.kaggle.com/akensert/panda-optimized-tiling-tf-data-dataset\n",
    "def _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",
    "\n",
    "def _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",
    "\n",
    "def _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",
    "\n",
    "def _get_tissue_parts_indices(tissue, min_consec_info):\n",
    "    \"\"\"If there are multiple tissue parts in 'tissue', 'tissue' will be \n",
    "    split it. Each tissue part will be taken care of separately, and if \n",
    "    the tissue part is less than min_consec_info, it's considered to small\n",
    "    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",
    "\n",
    "def _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 squeeze in the \n",
    "    N subparts in an area bigger than the sum of the N subparts or not.\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",
    "\n",
    "def _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 \n",
    "    value 'val1' that estimates how much information\n",
    "    there is in the patch.\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",
    "\n",
    "def 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 from 'min_axis_info'\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 True] : 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\n",
    "    \n",
    "def get_tiles(img,img0=None):\n",
    "    coords = compute_coords(img,\n",
    "                            patch_size=sz,\n",
    "                            precompute=True,\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",
    "    tiles = []\n",
    "    if img0 is None:\n",
    "        for (v, y, x) in coords:\n",
    "            p = img[y:y+sz,x:x+sz]\n",
    "            if p.shape[0] != sz or p.shape[1] != sz: continue\n",
    "            tiles.append(p)\n",
    "    else:\n",
    "        for (v, y, x) in coords:\n",
    "            p = img0[(4*y):(4*y+4*sz),(4*x):(4*x+4*sz)]\n",
    "            if p.shape[0] != 4*sz or p.shape[1] != 4*sz: continue\n",
    "            p = cv2.resize(p, (2*sz,2*sz), interpolation = cv2.INTER_AREA)\n",
    "            tiles.append(p)\n",
    "    if len(tiles) <= 1: return tiles\n",
    "    return [x for y, x in sorted(zip([t.sum() for t in tiles], tiles), key=lambda pair: pair[0])]    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c47c11c6c80f4ca8ac399942497458e5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=10616.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "from joblib import Parallel, delayed\n",
    "import multiprocessing\n",
    "\n",
    "names = [name[:-5] for name in os.listdir(TRAIN)]\n",
    "def convert(name):\n",
    "    x_tot,x2_tot = [],[]\n",
    "    img0,img = skimage.io.MultiImage(os.path.join(TRAIN,name+'.tiff'))[:2]\n",
    "    tiles = get_tiles(img,img0)\n",
    "    for idx,img in enumerate(tiles):\n",
    "        x_tot.append((img/255.0).reshape(-1,3).mean(0))\n",
    "        x2_tot.append(((img/255.0)**2).reshape(-1,3).mean(0)) \n",
    "        #if read with PIL RGB turns into BGR\n",
    "        cv2.imwrite(os.path.join(OUT_TRAIN,f'{name}_{idx}.png'),cv2.cvtColor(img, cv2.COLOR_RGB2BGR))\n",
    "    return x_tot,x2_tot\n",
    "\n",
    "num_cores = 6#multiprocessing.cpu_count()\n",
    "results = Parallel(n_jobs=num_cores)(delayed(convert)(name) for name in tqdm(names))\n",
    "x_tot,x2_tot = [],[]\n",
    "for r in results:\n",
    "    x_tot += r[0]\n",
    "    x2_tot += r[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('mean:',img_avr, ', std:', np.sqrt(img_std))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}