{"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 pandas as pd\nimport numpy as np\nimport cv2\nimport matplotlib.pyplot as plt\nfrom matplotlib.patches import Rectangle\n%matplotlib inline\nfrom pathlib import Path\nimport rasterio\nfrom rasterio.windows import Window\nimport warnings\nimport gc","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-09-27T18:48:03.814929Z","iopub.execute_input":"2022-09-27T18:48:03.815856Z","iopub.status.idle":"2022-09-27T18:48:03.825197Z","shell.execute_reply.started":"2022-09-27T18:48:03.81579Z","shell.execute_reply":"2022-09-27T18:48:03.823635Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Patch:\n    \"\"\"Holds data of a sampled patch.\n    \n    Patches are sampled from a grid of larger tiles of an image.\n    \n    Attributes:\n        crop (array): patch cropped from a tile.\n        image_id (str): image id.\n        yrange, xrange (tuple: int, int): upper, lower (y) and left, right (x) border of patch\n            relative to the tile the patch originates from.\n        yoff, xoff (int): coordinates of upper left corner of tile used as offset to\n            compute patch coordinates relative to the whole image.\n    \"\"\"\n    \n    def __init__(self, crop, image_id, yrange, xrange, yoff=0, xoff=0, label=\"unknown\"):\n        self._crop = crop\n        self.image_id = image_id\n        self.label = label\n        y1, y2 = yrange\n        x1, x2 = xrange\n        self.y1, self.y2 = y1 + yoff, y2 + yoff\n        self.x1, self.x2 = x1 + xoff, x2 + xoff\n        \n    def get_fname(self):\n        \"\"\"Return representative string used as file name.\"\"\"\n        return f\"{self.image_id}_{self.y1}-{self.y2}_{self.x1}-{self.x2}_{self.label}.tif\"\n    \n    def get_crop(self):\n        return self._crop\n\n    def get_mask(self):\n        return self._mask\n    \n    def get_x(self):\n        return self.x1, self.x2\n    \n    def get_y(self):\n        return self.y1, self.y2\n    \n    def __repr__(self):\n        return self.get_fname()","metadata":{"execution":{"iopub.status.busy":"2022-09-27T18:48:04.067651Z","iopub.execute_input":"2022-09-27T18:48:04.068068Z","iopub.status.idle":"2022-09-27T18:48:04.078981Z","shell.execute_reply.started":"2022-09-27T18:48:04.068033Z","shell.execute_reply":"2022-09-27T18:48:04.077357Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class PatchSampler:\n    \"\"\"Sample patches from an image.\n    \n    To control memory consumption, a grid of tiles generated and patches are sampled\n    from each tile at a time.\n    \n    Attributes:\n        image_path (Path): Path to the image file.\n        patch_size, step_size (tuple: int, int): Patch size and step size to sample subsequent\n            patches. Step size controls the extent of overlapping of sampled patches.\n        tile_size (int): imaga grid consists of square tiles with size tile_size x tile_size.\n        tile_foreground_ratio, patch_foreground_ratio (float): only tiles/patches with a\n            foreground ratio greater or equal to this cutoffs are sampled.\n        hsv_ratio (float): Threshold used to distinguish mask background artefacts. It is\n            applied over a HSV converted representation of a tile.\n        masking (bool): apply a mask to overlay background?\n    \"\"\"\n    \n    def __init__(self, image_path, patch_size=(224, 224), step_size=(112, 112), tile_size=2240,\n                 tile_foreground_ratio=0.05, patch_foreground_ratio=0.5, hsv_ratio=0.10, masking=False):\n        self.image_path = image_path\n        self.image_id = image_path.stem\n        self.patch_h, self.patch_w = patch_size\n        self.step_h, self.step_w = step_size\n        self.tile_size = tile_size\n        self.tile_foreground_ratio = tile_foreground_ratio\n        self.patch_foreground_ratio = patch_foreground_ratio\n        self.hsv_ratio = hsv_ratio\n        self.masking = masking\n        self.patch_list = []\n        with rasterio.open(self.image_path) as src:\n            width = src.width\n            height = src.height\n        self.height = height\n        self.width = width\n        tile_gen = self._tile_generator(height, width)\n        self.tile_coords = list(tile_gen)\n    \n    \n    def _get_tile_mask(self, im_arr):\n        \"\"\"Generate a background mask of a crop.\n        \n        Args:\n            im_arr (array): Image array.\n        \n        Returns:\n            array: background mask.\n        \"\"\"\n        \n        laplace = np.abs(cv2.Laplacian(cv2.cvtColor(im_arr, cv2.COLOR_RGB2GRAY), ddepth=cv2.CV_8U))\n        blurred = cv2.GaussianBlur(laplace, (21,21), 5)\n        _, mask = cv2.threshold(blurred, blurred.mean(), 255, cv2.THRESH_BINARY)\n        mask = cv2.morphologyEx(mask, cv2.MORPH_OPEN, np.ones((3,3), np.uint8), iterations=2)\n        mask = cv2.morphologyEx(mask, cv2.MORPH_CLOSE, np.ones((3,3), np.uint8), iterations=5)\n        return mask\n        \n    \n    def _tile_generator(self, height, width):\n        \"\"\"Divide an image into tiles.\n        \n        Square tiles with size tile_size x tile_size are generated. On each side\n        of the grid a margin is kept so that a grid of same sized tiles fits over an image.\n        \n        Args:\n            height, width (int): Height & width of the whole image.\n        \n        Yields:\n            array: A tile of size tile_size x tile_size of the image.\n        \"\"\"\n        \n        tile_size = self.tile_size\n        xp = width // tile_size\n        yp = height // tile_size\n        x_margin = (width - xp*tile_size) // 2\n        y_margin = (height - yp*tile_size) // 2\n        if xp > 0:\n            xcoords = [(x_margin+tile_size*i, x_margin+tile_size*(i+1)) for i in range(0, xp)]\n        else:\n            xcoords = [(0, width)]\n        if yp > 0:\n            ycoords = [(y_margin+tile_size*i, y_margin+tile_size*(i+1)) for i in range(0, yp)]\n        else:\n            ycoords = [(0, height)]\n        for iy in ycoords:\n            for ix in xcoords:\n                yield iy, ix\n    \n    \n    def _crop_mask_generator(self, im_arr, tile_mask, patch_h, patch_w, step_h, step_w):\n        \"\"\"Generate patches from an image tile.\n        \n        Args:\n            im_arr (array): Array representing a tile.\n            tile_mask (array): Tile background mask.\n            patch_h, patch_w, step_h, step_w (int): patch & step size\n        \n        Yields:\n            array: A patch sampled from a tile.\n        \"\"\"\n        \n        img_h, img_w, _ = im_arr.shape\n        for x1 in range(0, img_w - patch_w, step_w):\n            for y1 in range(0, img_h - patch_h, step_h):\n                # y-coordinate, i.e. height, goes first in array\n                # x-coordinate, i.e. width, comes 2nd\n                x2 = x1 + patch_w\n                y2 = y1 + patch_h\n                im_crop = im_arr[y1:y2, x1:x2, :]\n                mask_crop = tile_mask[y1:y2, x1:x2, np.newaxis]\n                yield im_crop, mask_crop, y1, y2, x1, x2\n    \n    \n    def _make_patches(self, im_arr, tile_mask, tile_y0=0, tile_x0=0):\n        \"\"\"Sample overlapping patches from a grid tile.\n        \n        A mask is generated for each generated patch. If the foreground constitutes\n        >= foreground_ratio of the patch it is sampled.\n        \n        Args:\n            im_arr (array): Array representing a tile.\n            tile_mask (array): Tile background mask.\n            tile_y0, tile_x0 (int): coordinates of the tile's upper left corner used\n                as offset to compute patch coordinates relative to the whole image.\n        Returns:\n            an array of sampled patches.\n        \"\"\"\n        \n        patch_h, patch_w = self.patch_h, self.patch_w\n        step_h, step_w = self.step_h, self.step_w\n        img_h, img_w, _ = im_arr.shape\n        patch_list = []\n        crop_gen = self._crop_mask_generator(im_arr, tile_mask, patch_h, patch_w, step_h, step_w)\n        for crop, mask, y1, y2, x1, x2 in crop_gen:\n            patch_fg_ratio = np.sum(mask > 0)/(patch_h * patch_w)\n            if patch_fg_ratio >= self.patch_foreground_ratio:\n                if self.masking:\n                    crop = crop & mask\n                p = Patch(crop, self.image_id, (y1,y2), (x1, x2), tile_y0, tile_x0)\n                patch_list.append(p)\n        return patch_list\n    \n    \n    def get_tile(self, idx, tile_y0=0, tile_x0=0):\n        \"\"\"Crop a tile from the image.\n        \n        Tiles are of size tile_size x tile_size.\n        Arguments:\n            idx (int): index of the tile into the tile coordinates array.\n            tile_y0, tile_x0 (int): coordinates of the tile's upper left corner used\n                as offset to compute patch coordinates relative to the whole image.\n        Returns:\n            A tuple of:\n                array: A grid tile.\n                float: tile foreground ratio.\n                array: array of sampled patches.\n            \n        \"\"\"\n        \n        if idx >= len(self.tile_coords):\n            return None\n        tile_coord = self.tile_coords[idx]\n        y1, y2 = tile_coord[0]\n        x1, x2 = tile_coord[1]\n        with rasterio.open(self.image_path) as src:\n            arr = src.read(window=Window(x1, y1, x2-x1, y2-y1))\n            im_arr = np.swapaxes(np.swapaxes(arr,0,2),0,1)\n            tile_mask = self._get_tile_mask(im_arr)\n        tile_fg_ratio = np.sum(tile_mask > 0)/np.prod(tile_mask.shape)\n        hsv = cv2.cvtColor(im_arr, cv2.COLOR_RGB2HSV)\n        # gets rid of bubble artefacts in the background\n        hsv_mask = cv2.inRange(hsv, (0, 11, 0), (180, 255, 255))\n        hsv_ratio = np.sum(hsv_mask == 255)/np.prod(hsv_mask.shape)\n        if tile_fg_ratio >= self.tile_foreground_ratio and hsv_ratio >= self.hsv_ratio:\n            patch_list = self._make_patches(im_arr, tile_mask, tile_y0, tile_x0)\n        else:\n            patch_list = []\n        if self.masking:\n            res = im_arr & tile_mask[:,:,np.newaxis]\n        else:\n            res = im_arr\n        return res, tile_fg_ratio, patch_list\n        \n        \n    def sample_patches(self):\n        \"\"\"Samples overlapping patches over all tiles of an image.\"\"\"\n        \n        for i, coords in enumerate(self.tile_coords):\n            y1, y2 = coords[0]\n            x1, x2 = coords[1]\n            tile, tile_fg_ratio, patch_list = self.get_tile(i, tile_y0=y1, tile_x0=x1)\n            if tile_fg_ratio >= self.tile_foreground_ratio:\n                self.patch_list += patch_list","metadata":{"execution":{"iopub.status.busy":"2022-09-27T18:48:04.282521Z","iopub.execute_input":"2022-09-27T18:48:04.282979Z","iopub.status.idle":"2022-09-27T18:48:04.318814Z","shell.execute_reply.started":"2022-09-27T18:48:04.28294Z","shell.execute_reply":"2022-09-27T18:48:04.316837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"warnings.filterwarnings(\"ignore\", category=rasterio.errors.NotGeoreferencedWarning)\nimage_path = Path(\"../input/mayo-clinic-strip-ai/test/00c058_0.tif\")\nps = PatchSampler(image_path, tile_foreground_ratio=0.05, masking=True)\n#ps.sample_patches()","metadata":{"execution":{"iopub.status.busy":"2022-09-27T18:48:04.549189Z","iopub.execute_input":"2022-09-27T18:48:04.549646Z","iopub.status.idle":"2022-09-27T18:48:04.573453Z","shell.execute_reply.started":"2022-09-27T18:48:04.549609Z","shell.execute_reply":"2022-09-27T18:48:04.572157Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tile, fg_ratio, patch_list = ps.get_tile(44)\nfig, ax = plt.subplots(figsize=(10, 10))\nax.imshow(tile)\nfor i, patch in enumerate(patch_list):\n    y1, y2 = patch.get_y()\n    x1, x2 = patch.get_x()\n    orig = (x1, y1)\n    h, w = y2-y1, x2-x1\n    rect = Rectangle(orig, w, h, linewidth=2, edgecolor='yellow', facecolor='none')\n    ax.add_patch(rect)","metadata":{"execution":{"iopub.status.busy":"2022-09-27T18:48:04.972661Z","iopub.execute_input":"2022-09-27T18:48:04.97324Z","iopub.status.idle":"2022-09-27T18:48:06.487892Z","shell.execute_reply.started":"2022-09-27T18:48:04.973172Z","shell.execute_reply":"2022-09-27T18:48:06.486591Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# May exhaust available memory depending on the size of the tif image\nps.sample_patches()\nim_arr = cv2.cvtColor(cv2.imread(str(image_path)), cv2.COLOR_BGR2RGB)\nfig, ax = plt.subplots(figsize=(50, 50))\nax.imshow(im_arr)\nfor i, patch in enumerate(ps.patch_list):\n    y1, y2 = patch.get_y()\n    x1, x2 = patch.get_x()\n    orig = (x1, y1)\n    h, w = y2-y1, x2-x1\n    rect = Rectangle(orig, w, h, linewidth=2, edgecolor='yellow', facecolor='none')\n    ax.add_patch(rect)","metadata":{"execution":{"iopub.status.busy":"2022-09-27T18:48:06.490429Z","iopub.execute_input":"2022-09-27T18:48:06.490973Z","iopub.status.idle":"2022-09-27T18:51:28.601992Z","shell.execute_reply.started":"2022-09-27T18:48:06.490921Z","shell.execute_reply":"2022-09-27T18:51:28.600768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}