{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Summary\n\nThe following notebook is used to separate the glomeruli tiles and visualize them in three different ways. \n* Firstly as raw tiles which can be used for further processing and analysis. \n![](http://)\n* Secondly as annotated tiles which show the glomeruli polygon on top of the raw data\n* Finally as mosaics of all the glomeruli of an image ordered by their size and labelled by their number. \n\nFurthermore, a preliminary analysis of the glomeruli sizes and shapes is presented at the final section\n\n\n# Output dataset\nThe output of the notebook is stored in the dataset [HuBMAP: glomeruli tiles and mosaics](https://www.kaggle.com/anadelta/hubmap-glomeruli-tiles)\n\n\n# Future work \n\nI am considering continuing with: \n* A more extended analysis of glomeruli geometry\n* Some type of color analysis/clustering of the glomeruli\n\nSuggestions are always welcome. "},{"metadata":{},"cell_type":"markdown","source":"# Utilities\nThe following sections contain utilities that are required for loading both the images and their annotations."},{"metadata":{},"cell_type":"markdown","source":"## Paths"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"import os\n\ninput_path = '/kaggle/input/hubmap-kidney-segmentation'\ntrain_folder = os.path.join(input_path, 'train')\ntest_folder = os.path.join(input_path, 'test')\n\noutput_path = 'output'\n\n\nclass TrainingPaths:\n\n    def __init__(self):\n\n        self._tiff_names = [name for name in os.listdir(train_folder) if name.endswith('.tiff')]\n\n        self._len = len(self._tiff_names)\n        self._index = 0\n\n    def __next__(self):\n\n        while True:\n\n            if not self._index < self._len:\n                raise StopIteration\n\n            tiff_name = self._tiff_names[self._index]\n            tiff_path = os.path.join(train_folder, tiff_name)\n            name = os.path.splitext(tiff_name)[0]\n            annotation_path = os.path.splitext(tiff_path)[0] + '.json'\n\n            self._index += 1\n\n            if (os.path.isfile(tiff_path) and os.path.isfile(annotation_path)):\n                break\n\n        return name, tiff_path, annotation_path\n\n    def __iter__(self):\n\n        return self","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Tiff loader"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"from numpy import squeeze, transpose\nfrom tifffile import imread\n\n\ndef read_tiff(file_path):\n\n    tiff_array = imread(file_path)\n    tiff_array = squeeze(tiff_array)\n    if(tiff_array.shape[0] == 3):\n        tiff_array = transpose(tiff_array, (1, 2, 0))\n\n    return tiff_array","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Annotations iterator"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"\nfrom json import load\n\nclass Annotations:\n\n    def __init__(self, file_path):\n\n        with open(file_path) as af:\n            self._glomeruli_polygons = load(af)\n        self._len = len(self._glomeruli_polygons)\n        self._index = 0\n\n    def __next__(self):\n\n        if not self._index < self._len:\n            raise StopIteration\n\n        coords = self._glomeruli_polygons[self._index]['geometry']['coordinates'][0]\n\n        i = self._index\n        x = [c[0] for c in coords]\n        y = [c[1] for c in coords]\n\n        self._index += 1\n\n        return i, x, y\n\n    def __iter__(self):\n        return self\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Geometry methods"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"\nimport numpy as np\n\n\ndef get_polygon_area(x, y):\n    return 0.5*np.abs(np.dot(x, np.roll(y, 1))-np.dot(y, np.roll(x, 1)))\n\n\ndef get_polygon_perimeter(x, y):\n    return np.sum(np.sqrt(np.power(x - np.roll(x, 1), 2) + np.power(y - np.roll(y, 1), 2)))\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Separate glomeruli tiles\n\nSeparate and store the raw and annotated glomeruli tiles"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"\nimport os\n\nfrom PIL import Image\n\n\nclass Tiles:\n\n    def __init__(self, image_name):\n        self._image_name = image_name\n        self._tiles_folder = os.path.join(output_path,\n                                          self._get_tiles_folder(),\n                                          image_name)\n        if not os.path.exists(self._tiles_folder):\n            os.makedirs(self._tiles_folder)\n\n    def _get_tiles_folder(self):\n        return 'raw_tiles'\n\n    def _get_tile_path(self, index):\n        return os.path.join(self._tiles_folder,\n                            f'{self._image_name}_{index:04}.png')\n\n    def save(self, tile_image, index):\n        tile_image.save(self._get_tile_path(index))\n\n    def open(self, index):\n        return Image.open(self._get_tile_path(index))\n\n\nclass AnnotatedTiles(Tiles):\n\n    def _get_tiles_folder(self):\n        return 'annotated_tiles'","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"import os\n\nfrom PIL import Image, ImageDraw\n\ndef separate_glomeruli_tiles():\n    print(\"Separating glomeruli tiles:\")\n\n    for name, tiff_path, annotation_path in TrainingPaths():\n        print(name)\n\n        kidney_array = read_tiff(tiff_path)\n\n        tiles = Tiles(name)\n        annotated_tiles = AnnotatedTiles(name)\n\n        for index, x, y in Annotations(annotation_path):\n\n            padd = 20  # Additional pixels outside the glomeruli bounding box\n\n            # Calculate the borders of the tile\n\n            x_min = max(min(x)-padd, 0)\n            x_max = min(max(x)+padd, kidney_array.shape[1])\n            y_min = max(min(y)-padd, 0)\n            y_max = min(max(y)+padd, kidney_array.shape[0])\n\n            # Get raw tile\n\n            tile_image = Image.fromarray(kidney_array[y_min:y_max, x_min:x_max, :])\n\n            # Save raw tile\n\n            tiles.save(tile_image, index)\n\n            # Draw the annotation polygon\n\n            draw = ImageDraw.Draw(tile_image)\n            draw.line(tuple(tuple([xd-x_min, yd-y_min]) for xd, yd in zip(x, y)),\n                      fill='#3dfe22', width=4)\n\n            # Save annotated tile\n\n            annotated_tiles.save(tile_image, index)\n\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"separate_glomeruli_tiles()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Create glomeruli mosaics\n\nSort all the glomeruli of an image according to their size and combine then into one representative mosaic.\n\nThe output of this step can be used in order to quickly separate and delete mislabelled glomeruli.\n\nFor higher resolution images please checkout the extracted [HuBMAP: glomeruli tiles and mosaics dataset](https://www.kaggle.com/anadelta/hubmap-glomeruli-tiles)\n"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"import os\n\nimport matplotlib.pyplot as plt\nfrom PIL import Image, ImageDraw, ImageFont\n\n\nclass Mosaic:\n\n    def __init__(self, width, margin):\n\n        # Initialize mosaic image\n\n        self._margin = margin\n        self._mosaic_image = Image.new('RGB', (width, 0))\n\n        # Initialize title font\n\n        # self._font = ImageFont.truetype('Pillow/Tests/fonts/FreeMono.ttf', size=self._margin)\n\n        # Initialize the row corners\n\n        self._row_ul = [self._margin, self._margin]  # Row upper left\n        self._row_lr = [self._margin, self._margin]  # Row lower right\n\n    def _extend_mosaic(self, new_width, new_height):\n\n        width, height = self._mosaic_image.size\n\n        if (new_width <= width) and (new_height <= height):\n            return\n\n        draft_mosaic = Image.new('RGB', (new_width, new_height))\n        draft_mosaic.paste(self._mosaic_image, (0, 0))\n        self._mosaic_image = draft_mosaic\n\n    def _crop_mosaic(self, new_width, new_height):\n\n        width, height = self._mosaic_image.size\n\n        if (new_width >= width) and (new_height >= height):\n            return\n\n        self._mosaic_image = self._mosaic_image.crop((0, 0, new_width, new_height))\n\n    def _get_mosaic_path(self, image_name):\n\n        mosaics_folder = os.path.join(output_path, 'annotated_mosaics')\n        if not os.path.exists(mosaics_folder):\n            os.makedirs(mosaics_folder)\n\n        return os.path.join(mosaics_folder, f'{image_name}.png')\n\n    def add_tile(self, tile, title):\n\n        width, _ = self._mosaic_image.size\n        tile_width, tile_height = tile.size\n\n        # Change row when necessary\n\n        if(self._row_lr[0]+tile_width > width):\n            self._row_ul = [self._margin, self._row_lr[1]+self._margin]\n            self._row_lr = [self._margin, self._row_lr[1]+self._margin]\n\n        # Paste tile into the mosaic\n\n        self._extend_mosaic(width, self._row_ul[1]+tile_width)\n        self._mosaic_image.paste(tile, (self._row_lr[0], self._row_ul[1]))\n\n        # Draw the title of the tile\n\n        # draw = ImageDraw.Draw(self._mosaic_image)\n        # draw.text((self._row_lr[0], self._row_ul[1]-self._margin), str(title),\n        #           font=self._font, fill='#3dfe22')\n\n        # Update the row corners\n\n        self._row_lr = [self._row_lr[0] + tile_width+self._margin,\n                        max(self._row_lr[1], self._row_ul[1]+tile_height)]\n\n    def save(self, image_name):\n\n        width, _ = self._mosaic_image.size\n\n        # Crop unused area of the mosaic\n\n        self._crop_mosaic(width, self._row_lr[1])\n\n        # Save the mosaic\n\n        self._mosaic_image.save(self._get_mosaic_path(image_name))\n\n    def show(self):\n\n        width, height = self._mosaic_image.size\n        figure_width = 20\n        figure_height = figure_width * height // width\n\n        fig, ax = plt.subplots(1, 1, sharex=True, figsize=(figure_width, figure_height))\n        ax.imshow(self._mosaic_image)\n        ax.axis('off')\n        plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"\nimport os\n\n\ndef create_glumeruli_mosaics():\n\n    for name, _, annotation_path in TrainingPaths():\n\n        annotated_tiles = AnnotatedTiles(name)\n\n        all_tiles = []\n        max_tile_size = [0, 0]\n\n        for index, x, y in Annotations(annotation_path):\n\n            area = get_polygon_area(x, y)\n\n            tile_image = annotated_tiles.open(index)\n\n            max_tile_size[0] = max(max_tile_size[0], tile_image.size[0])\n            max_tile_size[1] = max(max_tile_size[1], tile_image.size[1])\n\n            all_tiles.append({'image': tile_image, 'area': area, 'index': index})\n\n        max_width = 12000  # Maximum width of the mosaic\n        margin = 40       # Empty margin between tiles in the mosaic\n\n        mosaic = Mosaic(max_width, margin)\n\n        for tile in sorted(all_tiles,  key=lambda x: x['area']):\n\n            mosaic.add_tile(tile['image'], tile['index'])\n\n        mosaic.save(name)\n\n        mosaic.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"create_glumeruli_mosaics()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Study glomeruli"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"\nimport os\nfrom math import pi, log\n\nimport matplotlib.pyplot as plt\n\n\ndef get_glomeruli_properties():\n\n    all_props = []\n    for name, _, annotation_path in TrainingPaths():\n\n        props = []\n        for index, x, y in Annotations(annotation_path):\n\n            area = get_polygon_area(x, y)\n            perimeter = get_polygon_perimeter(x, y)\n\n            props.append({'name': name, 'index': index, 'area': area, 'perimeter': perimeter})\n\n        all_props.append(props)\n\n    return all_props\n\n\ndef hist_glomeruli_properties(properties, properties_function, x_label, output_filename):\n\n    fig, ax = plt.subplots(1, 2, sharex=True, figsize=(10, 4))\n\n    all_hist_data = []\n    for props in all_props:\n\n        hist_data = [properties_function(p) for p in props]\n        ax[0].hist(hist_data, histtype='step', label=props[0]['name'])\n\n        all_hist_data.extend(hist_data)\n\n    ax[0].legend()\n    ax[0].set_xlabel(x_label)\n    ax[0].set_ylabel('Frequency []')\n\n    ax[1].hist(all_hist_data, color='black', histtype='step', label='all')\n    ax[1].legend()\n    ax[1].set_xlabel(x_label)\n    ax[1].set_ylabel('Frequency []')\n\n    # Save plot\n\n    analysis_folder_path = os.path.join(output_path, 'analysis')\n    if not os.path.exists(analysis_folder_path):\n        os.makedirs(analysis_folder_path)\n\n    fig.savefig(os.path.join(analysis_folder_path, output_filename))\n\n\ndef hist_glomeruli_area(all_props):\n    hist_glomeruli_properties(all_props, lambda p: p['area']/10000,\n                              'Area [dpixel^2]',\n                              'hist_glomeruli_area.png')\n\n\ndef hist_glomeruli_roundness(all_props):\n    hist_glomeruli_properties(all_props, lambda p: 4*pi*p['area']/p['perimeter']/p['perimeter'],\n                              'Roundness []',\n                              'hist_glomeruli_roundness.png')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"all_props = get_glomeruli_properties()\nhist_glomeruli_area(all_props)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"hist_glomeruli_roundness(all_props)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Cleanup output"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import shutil\n\n\ndef cleanup_output():\n\n    try:\n        shutil.rmtree(output_path)\n    except OSError as e:\n        print(\"Error: %s - %s.\" % (e.filename, e.strerror))\n\n\ncleanup_output()\n","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}