{"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":"!pip install spams\n!pip install staintools","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pathlib\nimport json\nimport numpy as np\nimport pandas as pd\nimport cv2\nimport tifffile\nimport matplotlib.pyplot as plt\nimport staintools","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-08-07T15:49:16.392418Z","iopub.execute_input":"2022-08-07T15:49:16.392829Z","iopub.status.idle":"2022-08-07T15:49:16.398520Z","shell.execute_reply.started":"2022-08-07T15:49:16.392795Z","shell.execute_reply":"2022-08-07T15:49:16.397073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1. Introduction\n\nThere are two properties that are different across HPA and HuBMAP images and consistent for every image from the same organ type. Those are pixel size and tissue thickness of the whole slide image. In this notebook, those properties are heuristically adapted between HPA and HuBMAP images. Techniques are pretty straight-forward and experimental since there is no direct solution for this problem.","metadata":{}},{"cell_type":"code","source":"INTERNAL_DATASET = pathlib.Path('../input/hubmap-organ-segmentation')\n\ndf_train = pd.read_csv(INTERNAL_DATASET / 'train.csv')\ndf_test = pd.read_csv(INTERNAL_DATASET / 'test.csv')\n\ntrain_images = 'train_images/'\ntrain_annotations = 'train_annotations/'\ntest_images = 'test_images/'\n\ndf_train['image_filename'] = df_train['id'].apply(lambda x:  str(INTERNAL_DATASET) + '/' + train_images + str(x) + '.tiff')\ndf_train['polygon_filename'] = df_train['id'].apply(lambda x:  str(INTERNAL_DATASET) + '/' + train_annotations + str(x) + '.json')\ndf_test['image_filename'] = df_test['id'].apply(lambda x:  str(INTERNAL_DATASET) + '/' + test_images + str(x) + '.tiff')\ndf_test['age'] = np.nan\ndf_test['sex'] = np.nan\n\nprint(f'Training Set Shape: {df_train.shape} - Memory Usage: {df_train.memory_usage().sum() / 1024 ** 2:.2f} MB')\nprint(f'Test Set Shape: {df_test.shape} - Memory Usage: {df_test.memory_usage().sum() / 1024 ** 2:.2f} MB')","metadata":{"execution":{"iopub.status.busy":"2022-08-07T15:49:18.675897Z","iopub.execute_input":"2022-08-07T15:49:18.676543Z","iopub.status.idle":"2022-08-07T15:49:18.871756Z","shell.execute_reply.started":"2022-08-07T15:49:18.676509Z","shell.execute_reply":"2022-08-07T15:49:18.870564Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Pixel sizes and tissue thicknesses are provided in [data](https://www.kaggle.com/competitions/hubmap-organ-segmentation/data) tab. Pixel size is 0.4 μm and tissue thickness is 4 μm for every organ type for HPA images. Pixel sizes and tissue thicknesses are different among organs for HuBMAP images. Pixel sizes are 0.5 µm for kidney, 0.2290 µm for large intestine, 0.7562 µm for lung, 0.4945 µm for spleen, and 6.263 µm for prostate. Tissue slice thicknesses are 10 µm for kidney, 8 µm for large intestine, 4 µm for spleen, 5 µm for lung, and 5 µm for prostate.","metadata":{}},{"cell_type":"code","source":"imaging_measurements = {\n    'hpa': {\n        'pixel_size': {\n            'kidney': 0.4,\n            'prostate': 0.4,\n            'largeintestine': 0.4,\n            'spleen': 0.4,\n            'lung': 0.4\n        },\n        'tissue_thickness': {\n            'kidney': 4,\n            'prostate': 4,\n            'largeintestine': 4,\n            'spleen': 4,\n            'lung': 4\n        }\n    },\n    'hubmap': {\n        'pixel_size': {\n            'kidney': 0.5,\n            'prostate': 6.263,\n            'largeintestine': 0.229,\n            'spleen': 0.4945,\n            'lung': 0.7562\n        },\n        'tissue_thickness': {\n        'kidney': 10,\n            'prostate': 5,\n            'largeintestine': 8,\n            'spleen': 4,\n            'lung': 5\n        }\n    }\n}\n\nprint(json.dumps(imaging_measurements, indent=2))","metadata":{"execution":{"iopub.status.busy":"2022-08-07T15:20:56.855988Z","iopub.execute_input":"2022-08-07T15:20:56.856437Z","iopub.status.idle":"2022-08-07T15:20:56.865808Z","shell.execute_reply.started":"2022-08-07T15:20:56.856401Z","shell.execute_reply":"2022-08-07T15:20:56.864538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A single sample is selected from every organ type for experiments.","metadata":{}},{"cell_type":"code","source":"np.random.seed(42)\ndf_sample = df_train.groupby('organ').sample(1).reset_index(drop=True)\ndf_sample","metadata":{"execution":{"iopub.status.busy":"2022-08-07T15:34:32.579932Z","iopub.execute_input":"2022-08-07T15:34:32.580427Z","iopub.status.idle":"2022-08-07T15:34:32.608446Z","shell.execute_reply.started":"2022-08-07T15:34:32.580386Z","shell.execute_reply":"2022-08-07T15:34:32.607360Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Pixel Size\n\nAdapting the pixel size is easier. Relative scale factor can be calculated by `domain_pixel_size / target_pixel_size`. For example, HPA prostate images have 0.4 μm and HuBMAP prostate images have 6.263 μm pixel size. `0.4 / 6.263` is `0.063` which is approximately 1/16, so HPA prostate images are 16x larger than HuBMAP prostate images. Downsampling and upsampling HPA prostate images by that scale could yield similar images to HuBMAP prostate images in terms of pixel size. This is shared by [@hengck23](https://www.kaggle.com/hengck23) in [this](https://www.kaggle.com/competitions/hubmap-organ-segmentation/discussion/332941#1874042) post.\n\n### Tissue Thickness\n\nAdapting the tissue thickness is more experimental. From the figures below (Chilipala, 2020, fig. 3, 10), it can be seen that stain intensity increases when tissue is thicker. It makes a lot of sense because as the tissue gets thicker, less light can go through it.\n\nIncrease of stain intesity is very strong in H&E stained images for both hematoxylin and eosin for every organ. However, it varies for IHC methods. As same of the HuBMAP images are H&E stained, same transformation can be applied to them. \n\n![tissue_thickness1](https://i.ibb.co/jzx3p9B/Screenshot-from-2022-08-07-20-10-26.png)\n\n![tissue_thickness2](https://i.ibb.co/G9D9Bqm/Screenshot-from-2022-08-07-20-12-47.png)","metadata":{}},{"cell_type":"markdown","source":"I used this function to adapt pixel sizes and tissue thicknesses. The function is very simple at this point. \n\nFor tissue thickness, it multiplies saturation channel with `1 + alpha * (target_tissue_thickness - domain_tissue_thickness)` and multiplies value channel with `1 - alpha * (target_tissue_thickness - domain_tissue_thickness)`. It always upscales saturation and downscales value because all of HuBMAP image tissues are thicker than HPA image tissues, thus they might have more stain intensity. I added another multiplier alpha and set it to 0.15 which worked good on sample images but it might break on different images.\n\nFor pixel size, it downsamples and upsamples the same image by a factor of `domain_pixel_size / target_pixel_size`.\n\nLuminosity is standardized using `staintools` package in order reduce artifacts on the background.","metadata":{}},{"cell_type":"code","source":"def augment_image(image, domain_pixel_size, target_pixel_size, domain_tissue_thickness, target_tissue_thickness, alpha=0.15):\n    \n    \"\"\"\n    Visualize raw and augmented images \n    \n    Parameters\n    ----------\n    image (numpy.ndarray of shape (height, width, 3)): Image array\n    domain_pixel_size (float): Pixel size of the domain images in micrometers\n    target_pixel_size (float): Pixel size of the target images in micrometers\n    domain_tissue_thickness (float): Tissue thickness of the domain images in micrometers\n    target_tissue_thickness (float): Tissue thickness of the target images in micrometers\n    alpha (float): Multiplier to control saturation and value scale\n    \"\"\"\n    \n    # Augment tissue thickness\n    tissue_thickness_scale_factor = target_tissue_thickness - domain_tissue_thickness\n    image_hsv = cv2.cvtColor(image, cv2.COLOR_RGB2HSV).astype(np.float32)\n    image_hsv[:, :, 1] *= (1 + (alpha * tissue_thickness_scale_factor))\n    image_hsv[:, :, 2] *= (1 - (alpha * tissue_thickness_scale_factor))\n    image_hsv = image_hsv.astype(np.uint8)\n    image_scaled = cv2.cvtColor(image_hsv, cv2.COLOR_HSV2RGB)\n    \n    # Standardize luminosity\n    image_scaled = staintools.LuminosityStandardizer.standardize(image_scaled)\n\n    # Augment pixel size\n    pixel_size_scale_factor = domain_pixel_size / target_pixel_size\n    image_resized = cv2.resize(\n        image_scaled,\n        dsize=None,\n        fx=pixel_size_scale_factor,\n        fy=pixel_size_scale_factor,\n        interpolation=cv2.INTER_CUBIC\n    )\n    image_resized = cv2.resize(\n        image_resized,\n        dsize=(\n            image.shape[1],\n            image.shape[0]\n        ),\n        interpolation=cv2.INTER_CUBIC\n    )\n    \n    # Standardize luminosity\n    image = staintools.LuminosityStandardizer.standardize(image)\n    image_augmented = staintools.LuminosityStandardizer.standardize(image_resized)\n    \n    return image, image_augmented\n","metadata":{"execution":{"iopub.status.busy":"2022-08-07T17:48:10.882782Z","iopub.execute_input":"2022-08-07T17:48:10.883178Z","iopub.status.idle":"2022-08-07T17:48:10.893091Z","shell.execute_reply.started":"2022-08-07T17:48:10.883147Z","shell.execute_reply":"2022-08-07T17:48:10.892070Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def visualize_augmentation(image, image_augmented, metadata, path=None):\n\n    \"\"\"\n    Visualize raw and augmented images \n    \n    Parameters\n    ----------\n    image (numpy.ndarray of shape (height, width, 3)): Raw image array\n    image_augmented (numpy.ndarray of shape (height, width, 3)): Augmented image array\n    metadata (dict): Dictionary of image metadata\n    path (path-like str or None): Path of the output file or None (if path is None, plot is displayed with selected backend)\n    \"\"\"\n\n    fig, axes = plt.subplots(figsize=(36, 20), ncols=2)\n    \n    axes[0].imshow(image)\n    axes[1].imshow(image_augmented)\n\n    for i in range(2):\n        axes[i].set_xlabel('')\n        axes[i].set_ylabel('')\n        axes[i].tick_params(axis='x', labelsize=15, pad=10)\n        axes[i].tick_params(axis='y', labelsize=15, pad=10)\n\n    axes[0].set_title(f'Raw Image\\nMean: {np.mean(image):.4f} - Std: {np.std(image):.4f} - Min: {np.min(image)} - Max: {np.max(image)}', size=25, pad=15)\n    axes[1].set_title(f'Augmented Image\\nMean: {np.mean(image_augmented):.4f} - Std: {np.std(image_augmented):.4f} - Min: {np.min(image_augmented)} - Max: {np.max(image_augmented)}', size=25, pad=15)\n\n    fig.suptitle(\n        f'''\n        Image ID {metadata[\"id\"]} - {metadata[\"organ\"]} - {metadata[\"data_source\"]} - {metadata[\"age\"]} - {metadata[\"sex\"]}\n        Raw Image Shape: {metadata[\"img_height\"]}x{metadata[\"img_width\"]} - Pixel Size: {metadata[\"pixel_size\"]}µm - Tissue Thickness: {metadata[\"tissue_thickness\"]}µm\n        Augmented Image Shape: {metadata[\"img_height\"]}x{metadata[\"img_width\"]} - Pixel Size: {imaging_measurements[\"hubmap\"][\"pixel_size\"][metadata[\"organ\"]]}µm - Tissue Thickness: {imaging_measurements[\"hubmap\"][\"tissue_thickness\"][metadata[\"organ\"]]}µm\n        ''',\n        fontsize=50,\n        y=1.05\n    )\n\n    if path is None:\n        plt.show()\n    else:\n        plt.savefig(path)\n        plt.close(fig)","metadata":{"execution":{"iopub.status.busy":"2022-08-07T18:06:15.617883Z","iopub.execute_input":"2022-08-07T18:06:15.618315Z","iopub.status.idle":"2022-08-07T18:06:15.630609Z","shell.execute_reply.started":"2022-08-07T18:06:15.618276Z","shell.execute_reply":"2022-08-07T18:06:15.629629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Most intense transformations are kidney and large intestine images in terms of tissue thickness and prostate image in terms of pixel size. Other images aren't drastically transformed because their imaging measurements were close to each other. ","metadata":{}},{"cell_type":"code","source":"for idx, row in df_sample.iterrows():\n    \n    image = tifffile.imread(row['image_filename'])\n    image, image_augmented = augment_image(\n        image=image,\n        domain_pixel_size=imaging_measurements['hpa']['pixel_size'][row['organ']],\n        target_pixel_size=imaging_measurements['hubmap']['pixel_size'][row['organ']],\n        domain_tissue_thickness=imaging_measurements['hpa']['tissue_thickness'][row['organ']],\n        target_tissue_thickness=imaging_measurements['hubmap']['tissue_thickness'][row['organ']],\n        alpha=0.15\n    )\n    \n    visualize_augmentation(\n        image=image,\n        image_augmented=image_augmented,\n        metadata=row.to_dict()\n    )\n","metadata":{"execution":{"iopub.status.busy":"2022-08-07T18:05:34.044173Z","iopub.execute_input":"2022-08-07T18:05:34.044590Z","iopub.status.idle":"2022-08-07T18:05:58.277050Z","shell.execute_reply.started":"2022-08-07T18:05:34.044555Z","shell.execute_reply":"2022-08-07T18:05:58.275541Z"},"trusted":true},"execution_count":null,"outputs":[]}]}