{"cells":[{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"# Importing all the needed libraries\nimport os\nimport openslide\nimport skimage.io\nimport pandas as pd\nimport numpy as np\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport PIL\nfrom IPython.display import Image, display\nimport plotly.graph_objs as go","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Train and Testing data\nimport pandas as pd\ntrain = pd.read_csv('../input/prostate-cancer-grade-assessment/train.csv')\ntrain.columns.to_list()\ntest = pd.read_csv('../input/prostate-cancer-grade-assessment/test.csv')\ntest.columns.to_list()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\n# Plot histogram\ntrain.data_provider.hist(bins=30, alpha=0.5)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Location of the training images\ndata_dir = '/kaggle/input/prostate-cancer-grade-assessment/train_images'\nmask_dir = '/kaggle/input/prostate-cancer-grade-assessment/train_label_masks'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def print_slide_details(slide, show_thumbnail=True, max_size=(600,400)):\n    \"\"\"Print some basic information about a slide\"\"\"\n    # Generate a small image thumbnail\n    if show_thumbnail:\n        display(slide.get_thumbnail(size=max_size))\n\n    # Here we compute the \"pixel spacing\": the physical size of a pixel in the image.\n    # OpenSlide gives the resolution in centimeters so we convert this to microns.\n    spacing = 1 / (float(slide.properties['tiff.XResolution']) / 10000)\n    \n    print(f\"File id: {slide}\")\n    print(f\"Dimensions: {slide.dimensions}\")\n    print(f\"Microns per pixel / pixel spacing: {spacing:.3f}\")\n    print(f\"Number of levels in the image: {slide.level_count}\")\n    print(f\"Downsample factor per level: {slide.level_downsamples}\")\n    print(f\"Dimensions of levels: {slide.level_dimensions}\")\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Open the image (does not yet read the image into memory)\nimage = openslide.OpenSlide(os.path.join(data_dir, '005e66f06bce9c2e49142536caf2f6ee.tiff'))\n\n# Read a specific region of the image starting at upper left coordinate (x=17800, y=19500) on level 0 and extracting a 256*256 pixel patch.\n# At this point image data is read from the file and loaded into memory.\npatch = image.read_region((17800,19500), 0, (256, 256))\n\n# Display the image\ndisplay(patch)\n\n# Close the opened slide after use\nimage.close()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"example_slides = [\n    '005e66f06bce9c2e49142536caf2f6ee',\n    '00928370e2dfeb8a507667ef1d4efcbb',\n    '007433133235efc27a39f11df6940829',\n    '024ed1244a6d817358cedaea3783bbde',\n]\n\nfor case_id in example_slides:\n    biopsy = openslide.OpenSlide(os.path.join(data_dir, f'{case_id}.tiff'))\n    print_slide_details(biopsy)\n    biopsy.close()\n    \n    # Print the case-level label\n    print(f\"ISUP grade: {train_labels.loc[case_id, 'isup_grade']}\")\n    print(f\"Gleason score: {train_labels.loc[case_id, 'gleason_score']}\\n\\n\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"biopsy = openslide.OpenSlide(os.path.join(data_dir, '00928370e2dfeb8a507667ef1d4efcbb.tiff'))\n\nx = 5150\ny = 21000\nlevel = 0\nwidth = 512\nheight = 512\n\nregion = biopsy.read_region((x,y), level, (width, height))\ndisplay(region)\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"x = 5140\ny = 21000\nlevel = 1\nwidth = 512\nheight = 512\n\nregion = biopsy.read_region((x,y), level, (width, height))\ndisplay(region)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"biopsy.close()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def print_mask_details(slide, center='radboud', show_thumbnail=True, max_size=(400,400)):\n    \"\"\"Print some basic information about a slide\"\"\"\n\n    if center not in ['radboud', 'karolinska']:\n        raise Exception(\"Unsupported palette, should be one of [radboud, karolinska].\")\n\n    # Generate a small image thumbnail\n    if show_thumbnail:\n        # Read in the mask data from the highest level\n        # We cannot use thumbnail() here because we need to load the raw label data.\n        mask_data = slide.read_region((0,0), slide.level_count - 1, slide.level_dimensions[-1])\n        # Mask data is present in the R channel\n        mask_data = mask_data.split()[0]\n\n        # To show the masks we map the raw label values to RGB values\n        preview_palette = np.zeros(shape=768, dtype=int)\n        if center == 'radboud':\n            # Mapping: {0: background, 1: stroma, 2: benign epithelium, 3: Gleason 3, 4: Gleason 4, 5: Gleason 5}\n            preview_palette[0:18] = (np.array([0, 0, 0, 0.5, 0.5, 0.5, 0, 1, 0, 1, 1, 0.7, 1, 0.5, 0, 1, 0, 0]) * 255).astype(int)\n        elif center == 'karolinska':\n            # Mapping: {0: background, 1: benign, 2: cancer}\n            preview_palette[0:9] = (np.array([0, 0, 0, 0.5, 0.5, 0.5, 1, 0, 0]) * 255).astype(int)\n        mask_data.putpalette(data=preview_palette.tolist())\n        mask_data = mask_data.convert(mode='RGB')\n        mask_data.thumbnail(size=max_size, resample=0)\n        display(mask_data)\n\n    # Compute microns per pixel (openslide gives resolution in centimeters)\n    spacing = 1 / (float(slide.properties['tiff.XResolution']) / 10000)\n    \n    print(f\"File id: {slide}\")\n    print(f\"Dimensions: {slide.dimensions}\")\n    print(f\"Microns per pixel / pixel spacing: {spacing:.3f}\")\n    print(f\"Number of levels in the image: {slide.level_count}\")\n    print(f\"Downsample factor per level: {slide.level_downsamples}\")\n    print(f\"Dimensions of levels: {slide.level_dimensions}\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mask = openslide.OpenSlide(os.path.join(mask_dir, '08ab45297bfe652cc0397f4b37719ba1_mask.tiff'))\nprint_mask_details(mask, center='radboud')\nmask.close()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mask = openslide.OpenSlide(os.path.join(mask_dir, '090a77c517a7a2caa23e443a77a78bc7_mask.tiff'))\nprint_mask_details(mask, center='karolinska')\nmask.close()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mask = openslide.OpenSlide(os.path.join(mask_dir, '08ab45297bfe652cc0397f4b37719ba1_mask.tiff'))\nmask_data = mask.read_region((0,0), mask.level_count - 1, mask.level_dimensions[-1])\n\nplt.figure()\nplt.title(\"Mask with default cmap\")\nplt.imshow(np.asarray(mask_data)[:,:,0], interpolation='nearest')\nplt.show()\n\nplt.figure()\nplt.title(\"Mask with custom cmap\")\n# Optional: create a custom color map\ncmap = matplotlib.colors.ListedColormap(['black', 'gray', 'green', 'yellow', 'orange', 'red'])\nplt.imshow(np.asarray(mask_data)[:,:,0], cmap=cmap, interpolation='nearest', vmin=0, vmax=5)\nplt.show()\n\nmask.close()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def overlay_mask_on_slide(slide, mask, center='radboud', alpha=0.8, max_size=(800, 800)):\n    \"\"\"Show a mask overlayed on a slide.\"\"\"\n\n    if center not in ['radboud', 'karolinska']:\n        raise Exception(\"Unsupported palette, should be one of [radboud, karolinska].\")\n\n    # Load data from the highest level\n    slide_data = slide.read_region((0,0), slide.level_count - 1, slide.level_dimensions[-1])\n    mask_data = mask.read_region((0,0), mask.level_count - 1, mask.level_dimensions[-1])\n\n    # Mask data is present in the R channel\n    mask_data = mask_data.split()[0]\n\n    # Create alpha mask\n    alpha_int = int(round(255*alpha))\n    if center == 'radboud':\n        alpha_content = np.less(mask_data.split()[0], 2).astype('uint8') * alpha_int + (255 - alpha_int)\n    elif center == 'karolinska':\n        alpha_content = np.less(mask_data.split()[0], 1).astype('uint8') * alpha_int + (255 - alpha_int)\n    \n    alpha_content = PIL.Image.fromarray(alpha_content)\n    preview_palette = np.zeros(shape=768, dtype=int)\n    \n    if center == 'radboud':\n        # Mapping: {0: background, 1: stroma, 2: benign epithelium, 3: Gleason 3, 4: Gleason 4, 5: Gleason 5}\n        preview_palette[0:18] = (np.array([0, 0, 0, 0.5, 0.5, 0.5, 0, 1, 0, 1, 1, 0.7, 1, 0.5, 0, 1, 0, 0]) * 255).astype(int)\n    elif center == 'karolinska':\n        # Mapping: {0: background, 1: benign, 2: cancer}\n        preview_palette[0:9] = (np.array([0, 0, 0, 0, 1, 0, 1, 0, 0]) * 255).astype(int)\n    \n    mask_data.putpalette(data=preview_palette.tolist())\n    mask_rgb = mask_data.convert(mode='RGB')\n\n    overlayed_image = PIL.Image.composite(image1=slide_data, image2=mask_rgb, mask=alpha_content)\n    overlayed_image.thumbnail(size=max_size, resample=0)\n\n    display(overlayed_image)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"slide = openslide.OpenSlide(os.path.join(data_dir, '08ab45297bfe652cc0397f4b37719ba1.tiff'))\nmask = openslide.OpenSlide(os.path.join(mask_dir, '08ab45297bfe652cc0397f4b37719ba1_mask.tiff'))\noverlay_mask_on_slide(slide, mask, center='radboud')\nslide.close()\nmask.close()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"slide = openslide.OpenSlide(os.path.join(data_dir, '090a77c517a7a2caa23e443a77a78bc7.tiff'))\nmask = openslide.OpenSlide(os.path.join(mask_dir, '090a77c517a7a2caa23e443a77a78bc7_mask.tiff'))\noverlay_mask_on_slide(slide, mask, center='karolinska', alpha=0.6)\nslide.close()\nmask.close()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"biopsy = skimage.io.MultiImage(os.path.join(data_dir, '0b373388b189bee3ef6e320b841264dd.tiff'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"for i,level in enumerate(biopsy):\n    print(f\"Biopsy level {i} dimensions: {level.shape}\")\n    print(f\"Biopsy level {i} memory size: {level.nbytes / 1024**2:.1f}mb\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"display(PIL.Image.fromarray(biopsy[-1]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Deleting the object frees up memory\ndel biopsy","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"If you are only interested in the lowest level (highest magnification), you can also load level 0 using [imread](https://scikit-image.org/docs/0.16.x/api/skimage.io.html?highlight=imread#skimage.io.imread):","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"biopsy_level_0 = skimage.io.imread(os.path.join(data_dir, '0b373388b189bee3ef6e320b841264dd.tiff'))\nprint(biopsy_level_0.shape)\ndel biopsy_level_0","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"biopsy = skimage.io.MultiImage(os.path.join(data_dir, '00928370e2dfeb8a507667ef1d4efcbb.tiff'))\n\nx = 5150\ny = 21000\nlevel = 0\nwidth = 512\nheight = 512\n\npatch = biopsy[0][y:y+width, x:x+height]\n\n# You can also visualize patches with matplotlib\nplt.figure()\nplt.imshow(patch)\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"x = 5150 // 4\ny = 21000 // 4\nwidth = 512\nheight = 512\n\npatch = biopsy[1][y:y+width, x:x+height]\n\nplt.figure()\nplt.imshow(patch)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"x = 5150 // (4*4)\ny = 21000 // (4*4)\nwidth = 512\nheight = 512\n\npatch = biopsy[2][y:y+width, x:x+height]\n\nplt.figure()\nplt.imshow(patch)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Free up memory\ndel biopsy","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"maskfile = skimage.io.MultiImage(os.path.join(mask_dir, '090a77c517a7a2caa23e443a77a78bc7_mask.tiff'))\nmask_level_2 = maskfile[-1][:,:,0]\n\nplt.figure()\nplt.imshow(mask_level_2)\nplt.colorbar()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del maskfile","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class WSIViewer(object):\n    def __init__(self, plot_size = 1000):\n        self._plot_size = plot_size\n        \n    def set_slide(self, slide_path):      \n        self._slide = openslide.open_slide(slide_path)\n        self._base_dims = self._slide.level_dimensions[-1]\n        self._base_ds = self._slide.level_downsamples[-1]\n        img_arr = self._slide.read_region((0,0), len(self._slide.level_dimensions[-1]), (self._base_dims[0], self._base_dims[1]))\n        \n        self._fig = go.FigureWidget(data=[{'x': [0, self._base_dims[0]], \n                                           'y': [0, self._base_dims[1]], \n                                           'mode': 'markers',\n                                           'marker': {'opacity': 0}}], # invisible trace to init axes and to support autoresize\n                                    layout={'width': self._plot_size, 'height': self._plot_size, 'yaxis' : dict(scaleanchor = \"x\", scaleratio = 1)})  \n        # Set background image\n        self._fig.layout.images = [go.layout.Image(\n            source = img_arr,  # plotly now performs auto conversion of PIL image to png data URI\n            xref = \"x\",\n            yref = \"y\",\n            x = 0,\n            y = 0,\n            sizex = self._base_dims[0],\n            sizey = self._base_dims[1],\n            sizing = \"stretch\",\n            layer = \"below\")]\n        self._fig.update_layout(plot_bgcolor='rgba(0,0,0,0)',xaxis_showgrid=False, yaxis_showgrid=False, xaxis_zeroline=False, yaxis_zeroline=False);        \n        self._fig.layout.on_change(self._update_image, 'xaxis.range', 'yaxis.range', 'width', 'height')          \n\n    def _gen_zoomed_image(self, x_range, y_range):\n        # Below is a workaround which rounds image requests to multiples of 4, once the libpixman fix is in place these can be removed\n        #xstart = x_range[0] * self._base_ds\n        #ystart = (self._base_dims[1] - y_range[1]) * self._base_ds \n        xstart = 4 * round(x_range[0] * self._base_ds / 4)\n        ystart = 4 * round((self._base_dims[1] - y_range[1]) * self._base_ds / 4)\n        xsize0 = (x_range[1] - x_range[0]) * self._base_ds\n        ysize0 = (y_range[1] - y_range[0]) * self._base_ds\n        if (xsize0 > ysize0):\n            req_downs = xsize0 / self._plot_size\n        else:\n            req_downs = ysize0 / self._plot_size\n        req_level = self._slide.get_best_level_for_downsample(req_downs)\n        level_downs = self._slide.level_downsamples[req_level]\n        # Nasty workaround for buggy container\n        level_size_x = int(xsize0 / level_downs)\n        level_size_y = int(ysize0 / level_downs)\n        new_img = self._slide.read_region((int(xstart), int(ystart)), req_level, (level_size_x, level_size_y)).resize((1000,1000)) # Letting PIL do the resize is faster than plotly\n        return new_img\n    \n    def _update_image(self, layout, x_range, y_range, plot_width, plot_height):\n        img = self._fig.layout.images[0]\n        # Update with batch_update so all updates happen simultaneously\n        with self._fig.batch_update():\n            new_img = self._gen_zoomed_image(x_range, y_range)\n            img.x = x_range[0]\n            img.y = y_range[1]\n            img.sizex = x_range[1] - x_range[0]\n            img.sizey = y_range[1] - y_range[0]\n            img.source = new_img\n\n    def show(self):\n        return self._fig","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"viewer = WSIViewer()\nviewer.set_slide(os.path.join(data_dir, '08ab45297bfe652cc0397f4b37719ba1.tiff'))\nviewer.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_sub = pd.read_csv('../input/prostate-cancer-grade-assessment/sample_submission.csv')","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}