{"cells":[{"metadata":{},"cell_type":"markdown","source":"# PanDa tiles preprocessing\nThis notebook is about the [PanDa Dataset 2020](https://www.kaggle.com/c/prostate-cancer-grade-assessment) - Prostate Cancer Grade Assessment.\nI propose here an approach to generate tile information without having to save them locally but just by keeping any relevant information like coordinates and details Gleason score (from Radboud mask). This notebook only concern image preprocessing before any model training.\nI will particularly focus on runtime issue because it can be very frustating to wait too long between each iteration. So I mention everytime is possible\n\nYou can follow my details reflexion in ***part I - Explorating code*** or jump and directly to the final code in ***part II - Final code***. \nA stand-alone script final version is also available: [PanDa Final script](https://www.kaggle.com/huynhdoo/panda-wsi-tiles-preprocessing-script). Feel free to copy, reuse and fork anything that you find usefull!\n\nAuthor : Do Huynh / Date : June 1st 2020","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# I. Explorating code","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Import common libraries\nimport os\nimport glob\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport numpy as np\nfrom tqdm.notebook import tqdm\nimport openslide\nfrom openslide import OpenSlideError\nfrom IPython.display import Image\nimport seaborn as sns\nimport multiprocessing\nimport datetime\n\nclass Time:\n  \"\"\"\n  Class for displaying elapsed time.\n  \"\"\"\n\n  def __init__(self):\n    self.start = datetime.datetime.now()\n\n  def elapsed_display(self):\n    time_elapsed = self.elapsed()\n    print(\"Time elapsed: \" + str(time_elapsed))\n\n  def elapsed(self):\n    self.end = datetime.datetime.now()\n    time_elapsed = self.end - self.start\n    return time_elapsed\n\n# Open a slide\ndef open_slide(filename):\n    \"\"\"\n    Open a whole-slide image (*.svs, etc).\n    :filename : Name of the slide file.\n    return: an OpenSlide object representing a whole-slide image.\n    \"\"\"\n    try:\n        slide = openslide.open_slide(filename)\n    except OpenSlideError:\n        slide = None\n    except FileNotFoundError:\n        slide = None\n    return slide    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# PARAMETERS\nBASE_DIR = '/kaggle/input/prostate-cancer-grade-assessment/'\nOUTPUT_DIR = './'\nTRAIN_DIR = os.path.join(BASE_DIR, \"train_images\")\nTRAIN_EXT = \".tiff\"\nMASK_DIR = os.path.join(BASE_DIR, \"train_label_masks\")\nMASK_EXT = \"_mask.tiff\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Number slide > Mask so we take the mask as the minimum value\ntrain = glob.glob1(TRAIN_DIR, \"*\" + TRAIN_EXT)\nlabel = glob.glob1(MASK_DIR, \"*\" + MASK_EXT)\n\n# Keep only image_id\ntrain = [x[:-len(TRAIN_EXT)] for x in train]\nlabel = [y[:-len(MASK_EXT)] for y in label]\n\nlen(train), len(label)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Add filenames to dataframe\ntrain_df = pd.read_csv(BASE_DIR + 'train.csv')\n\n# Add train file column for each existing file in train folder\ntrain_df['train_file'] = list(map(lambda x : x + TRAIN_EXT if x in set(train) else '', \n                              train_df['image_id']))\n# Add label file column for each existing file in mask folder\ntrain_df['label_file'] = list(map(lambda y : y + MASK_EXT if y in set(label) else '', \n                              train_df['image_id']))\ntrain_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Split dataframe by provider / we keep radboud scoring because their mask labels are more details\nprint('Dataframe original:', len(train_df))\ntrain_radboud = train_df[train_df['data_provider'] == 'radboud'].copy()\nprint('Dataframe after split:', len(train_radboud))\n\n# Keep only row with both train and label file\ntrain_radboud = train_radboud[train_radboud['train_file'] != '']\nprint('Dataframe after file select:', len(train_radboud))\ntrain_radboud = train_radboud[train_radboud['label_file'] != '']\nprint('Dataframe after label select:', len(train_radboud))\n\ntrain_radboud.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Check time to open and close a slide\nfile = train_radboud['train_file'].values[0]\nfilepath = os.path.join(TRAIN_DIR, file)\n\n# Open\nt = Time()\nprint('Open slide')\nbiopsy = open_slide(filepath)\nt.elapsed_display()\n\n# Close\nt = Time()\nprint('Close slide')\nbiopsy.close()\nt.elapsed_display()\n\n# OBSERVATION : <20ms to open and close a slide","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Check on all the dataset\nt = Time()\nfiles = train_radboud['train_file'].values[:1]\nfor file in tqdm(files):\n    filepath = os.path.join(TRAIN_DIR, file)    \n    # Open\n    biopsy = open_slide(filepath)\n    # Do something\n    biopsy.close()\nprint('Open and close', len(files),'slides')\nt.elapsed_display()\n\n# OBSERVATION : ~1-2 min to open and close more 5000 slides => ~20ms by slide","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Keep only ISUP grade = 5\ntrain_radboud5 = train_radboud[train_radboud['isup_grade'] == 5].copy()\nprint('Dataframe after grade select:', len(train_radboud5))\ntrain_radboud5.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Check times on this dataset\nt = Time()\nfiles = train_radboud5['train_file'].values\nfor file in tqdm(files):\n    filepath = os.path.join(TRAIN_DIR, file)    \n    # Open\n    biopsy = open_slide(filepath)\n    # Do something\n    biopsy.close()\nprint('Open and close', len(files),'slides')\nt.elapsed_display()\n\n# OBSERVATION : ~10s to open and close 964 slides","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Check time to open and close the lowest level image\nt = Time()\nfiles = train_radboud5['train_file'].values[:1]\nfor file in tqdm(files):\n    filepath = os.path.join(TRAIN_DIR, file)    \n    # Open\n    biopsy = open_slide(filepath)\n    # Read lowest definition image\n    level = biopsy.level_count - 1\n    dimensions = biopsy.level_dimensions[level]\n    sample = biopsy.read_region((0, 0), level, dimensions)\n    # Close\n    biopsy.close()\n    sample = None\nprint('Open, read image and close', len(files),'slides')\nt.elapsed_display()\n\n# OBSERVATION 2: <1s to open, read and close 964 slides\n# OBSERVATION 2: ~40s to open, read and close 964 slides","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Check time to open, save and close the lowest level image\nDEST_TRAIN_DIR = 'train_png'\nDEST_TRAIN_EXT = '.png'\n\nt = Time()\nfiles = train_radboud5['train_file'].values[:1]\nfor file in tqdm(files):\n    filepath = os.path.join(TRAIN_DIR, file)    \n    # Open\n    biopsy = open_slide(filepath)\n    # Read lowest definition image\n    level = biopsy.level_count - 1\n    dimensions = biopsy.level_dimensions[level]\n    sample = biopsy.read_region((0, 0), level, dimensions)\n    # Save\n    if not os.path.exists(DEST_TRAIN_DIR):\n        os.makedirs(DEST_TRAIN_DIR)\n    sample.save(os.path.join(DEST_TRAIN_DIR, file + DEST_TRAIN_EXT))\n    # Close\n    biopsy.close()\n    sample = None\nprint('Open, read, save and close', len(files),'slides')\nt.elapsed_display()\n\n# OBSERVATION 1: ~400ms to open, read, save and close 1 slide\n# OBSERVATION 2: < 3min to open, read, save and close 964 slides","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Observations about I/O times :\nAs usually, read/write operations on the filesystem tend to increase a lot the global operations times. In the same conditions, it will take more than 30 min. to proceed all the given dataset. Not so long in absolute time but for the purpose to iterate fastly on the dataset, we should avoid any unusefull I/O operation on the filesystem. Note that with a multiprocess split, we could drastically reduce the operations times.\n\nIn the ideal, on each file of the dataset, we may follow  these steps:\n    1. Open once a file\n    2. Make any operation (read, extract, ...)\n    3. Save once any export file\n    3. Close once a file","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Count the proportion of blank white pixel (= 255) by slide\nt = Time()\nfiles = train_radboud5['train_file'].values\nwhite_pixel = []\n\nfor file in tqdm(files):\n    filepath = os.path.join(TRAIN_DIR, file)    \n    # Open\n    biopsy = open_slide(filepath)\n    # Read lowest definition image\n    level = biopsy.level_count - 1\n    dimensions = biopsy.level_dimensions[level]\n    sample = biopsy.read_region((0, 0), level, dimensions)\n    num_pixels = dimensions[0]*dimensions[1]\n    sample = sample.convert(\"1\") #Convert to black and white\n    white_pixel.append(np.count_nonzero(sample)/num_pixels)  \n    # Close\n    biopsy.close()\n    sample = None\nprint('Open, read, save and close', len(files),'slides')\nt.elapsed_display()\n\ntrain_radboud5['white_proportion'] = white_pixel\nwhite_pixel = None\n\n# OBSERVATION 1: varying the level definition as no impact on the white proportion so we can keep the lowest level\n# OBSERVATION 2: 50s to count white pixel on 964 slides","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_radboud5.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Observation about unsignifiant pixel\nAs we can observe there are a lot of unsignifiant regions on each slide with many empty-white pixel (from 80% to 98% of a slide !). One of the first step will be to drop all these unsignifiant pixels by cropping our slide.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Displaying a very 'blank' slide\nt = Time()\nfiles = train_radboud5.loc[train_radboud5['white_proportion']>0.98]['train_file'].values[:1]\n\nfor file in tqdm(files):\n    filepath = os.path.join(TRAIN_DIR, file)    \n    # Open\n    biopsy = open_slide(filepath)\n    # Read lowest definition image\n    level = biopsy.level_count - 1\n    dimensions = biopsy.level_dimensions[level]\n    sample = biopsy.read_region((0, 0), level, dimensions)\n    display(sample)\n    # Close\n    biopsy.close()\n    sample = None\nprint('Open, show, save and close', len(files),'slides')\nt.elapsed_display()\n\n# OBSERVATION : 250ms to display a slide as low definition","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Grid tile selection\nAs we have seen previously, opening and saving an important number of big files is not a great idea. First, if you want to benefit to the best definition, the slide is to huge for the local memory and you will fastly face with not enough ressource. Secondly, in the specific case of prostate biopsy, there are at least 80% unusefull informations in the slide so you spend 80% of the time waiting for nothing relevant! Thirdly, saving any part of the original files is data redundant.\n\nBut how can we deal with this so big images to send it to any learning algorithm in an efficient way? Some Kaggle players have propose a tile extraction that keep the most relevant part of the image (= with most appearing tissues and less blank) and save them to disk (see [Iafoss](https://www.kaggle.com/iafoss/panda-16x128x128-tiles) and [PAB97](https://www.kaggle.com/rftexas/better-image-tiles-removing-white-spaces) kernels). Then you can load this tiles selection to train your model. But, after testing some of this solutions, I was disappointed by the time that it need (several hours...). Moreover, the tiles are often extract from the lower definition and if you want to keep more information or a different tile size, you must put a new coin in the time machine :(\n\nAlternatively, my assumption is that a \"reading on the fly\" of the slide is more efficient without local saving step. Like layers in Photoshop, you can create a tile grid on top of a slide that reveal only usefull informations like this proposition from [Deron Eriksson, Fei Hu](https://developer.ibm.com/articles/an-automatic-method-to-identify-tissues-from-big-whole-slide-images-pt4/). Let see how we can code this idea.\n","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Draw a grid on an image\n# https://stackoverflow.com/questions/20368413/draw-grid-lines-over-an-image-in-matplotlib\nimport matplotlib.ticker as plticker\n\n# Open image file\nt = Time()\nfiles = train_radboud5['train_file'].values[:1]\n\nfor file in tqdm(files):\n    filepath = os.path.join(TRAIN_DIR, file)    \n    # Open\n    biopsy = open_slide(filepath)\n    level = biopsy.level_count - 1\n    dimensions = biopsy.level_dimensions[level]\n    sample = biopsy.read_region((0, 0), level, dimensions)\n\n    # Resolution\n    dpi=100.\n\n    # Set up figure\n    fig=plt.figure(figsize=(float(sample.size[0])/dpi,float(sample.size[1])/dpi),dpi=dpi)\n    ax=fig.add_subplot(111)\n\n    # Set the gridding interval: here we use the major tick interval\n    interval=32\n    loc = plticker.MultipleLocator(base=interval)\n    ax.xaxis.set_major_locator(loc)\n    ax.yaxis.set_major_locator(loc)\n    ax.xaxis.set_ticklabels([])\n    ax.yaxis.set_ticklabels([])\n\n    # Add the grid\n    ax.grid(which='major', axis='both', linestyle='-')\n\n    # Add the image\n    ax.imshow(sample)\n\n    # Find number of gridsquares in x and y direction\n    nx=abs(int(float(ax.get_xlim()[1]-ax.get_xlim()[0])/float(interval)))\n    ny=abs(int(float(ax.get_ylim()[1]-ax.get_ylim()[0])/float(interval)))\n\n    # Add some labels to the gridsquares\n    for i in range(nx):\n        x=interval/2.+float(i)*interval\n        ax.text(x,interval/2,i,color='black',ha='center',va='center')\n    for j in range(ny):\n        y=interval/2+j*interval\n        ax.text(interval/2,y,j,color='black',ha='center',va='center')\n\n\n    # Save the figure\n    #fig.savefig('myImageGrid.tiff',dpi=my_dpi)\n    \n    # Close\n    biopsy.close()\n    sample = None\nprint('Open, draw grid and close', len(files),'slides')\nt.elapsed_display()\n\n# OBSERVATION : 210ms to draw a grid on slide at low definition","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"With this grid, we can select and display at different level any specific slide zone (x,y):","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# See in detail the (19,5) tile on different level size\nt = Time()\nfiles = train_radboud5['train_file'].values[:1]\nindex = (12, 47)\ninterval = 32\n\nfor file in tqdm(files):\n    filepath = os.path.join(TRAIN_DIR, file)    \n    # Open\n    biopsy = open_slide(filepath)\n    for level in range(biopsy.level_count):\n        scale = int(16/biopsy.level_downsamples[level]) # Scale factor of the given level\n        size = interval*scale # Tile size depend on scale factor\n        dimensions = (size, size)\n        x, y = index[0]*interval*16, index[1]*interval*16 #Localisation from the level 0 => * max scale interval to get coordinate\n        sample = biopsy.read_region((x, y), level, dimensions)\n\n        # Display\n        print('tile:', index, '- level:', level, '- scale:', scale,'- size:', size)\n        display(sample)\n\n    # Close\n    biopsy.close()\n    sample = None\nprint('Open, show tiles and close', len(files),'slides')\nt.elapsed_display()\n\n# OBSERVATION : 330ms to display tiles on each level definition","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We can see here that a tile size of 128x128px on scale 4 can be a good trade-off between size and quality details. It is equivalent of a grid of 32x32 tile on the low level 3.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Grid tile score\nKnowing how to extract a specific tile from a slide, we want now to score each tile for keeping only relevant information. We can apply the same previous process of white pixel counting on each tile.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Display a heatmap of tile color intensities\nt = Time()\nfiles = train_radboud5['train_file'].values[:1]\nwhite_pixel = []\ninterval = 32\n\nfor file in tqdm(files):\n    filepath = os.path.join(TRAIN_DIR, file)    \n    # Open\n    biopsy = open_slide(filepath)\n    # Read lowest definition image\n    level = biopsy.level_count - 1\n    dimensions = biopsy.level_dimensions[level]\n\n    # Get number of gridsquares in x and y direction\n    nx=int(dimensions[0]/interval)\n    ny=int(dimensions[1]/interval)\n    tiles = np.zeros((nx, ny))\n\n    # Browse each tiles\n    level = 1\n    scale = 4\n    size = interval*scale # Tile size depend on scale factor\n    dimensions = (size, size)\n    num_pixels = dimensions[0]*dimensions[1]\n    \n    for i in range(nx):\n        for j in range(ny):  \n            x, y = i*interval*16, j*interval*16 #Localisation from the level 0 => * max scale interval to get coordinate\n            sample = biopsy.read_region((x, y), level, dimensions)\n            sample = sample.convert(\"1\") #Convert to black and white\n            tiles[i][j] = 1-np.count_nonzero(sample)/num_pixels\n    white_pixel.append(tiles)\n    \n    # Close\n    biopsy.close()\n    sample = None\n    \nprint('Open, score tiles and close', len(files),'slides')\n\n# Generate a heatmap\ngrid = white_pixel[0]\nsns.set(style=\"white\")\nplt.subplots(figsize=(grid.shape[0]/5, grid.shape[1]/5))\n\nmask = np.zeros_like(grid)\nmask[np.where(grid < 0.1)] = True #Mask blank tiles\n\nsns.heatmap(grid.T, square=True, linewidths=.5, mask=mask.T, cbar=False, vmin=0, vmax=1, cmap=\"Reds\")\nplt.show()\nprint('Not-blank tiles:', np.count_nonzero(grid), 'on', grid.size, 'total tiles')\ngrid = None\nwhite_pixel = None\n\nt.elapsed_display()\n# OBSERVATION: 1.5s to count white pixel on 1 slide","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Tiles generation\nWe can extract tile informations (slide region coordinate) for only usefull tiles (= not empty)","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Generate tiles for one slide\nt = Time()\nfiles = train_radboud5['train_file'].values[:1]\ntiles = []\ninterval = 32\n\nfor file in tqdm(files):\n    filepath = os.path.join(TRAIN_DIR, file)\n    \n    # Open\n    biopsy = open_slide(filepath)\n    # Read lowest definition image\n    level = biopsy.level_count - 1\n    dimensions = biopsy.level_dimensions[level]\n\n    # Get number of gridsquares in x and y direction\n    nx=int(dimensions[0]/interval)\n    ny=int(dimensions[1]/interval)\n    #tiles = np.zeros((nx, ny))\n\n    # Browse each tiles\n    level = 1\n    scale = 4\n    size = interval*scale # Tile size depend on scale factor\n    dimensions = (size, size)\n    num_pixels = dimensions[0]*dimensions[1]\n    \n    for i in range(nx):\n        for j in range(ny):  \n            x, y = i*interval*16, j*interval*16 #Localisation from the level 0 => * max scale interval to get coordinate\n            sample = biopsy.read_region((x, y), level, dimensions)\n            sample = sample.convert(\"1\") #Convert to black and white\n            score = 1-np.count_nonzero(sample)/num_pixels\n            # Keep only not blank tiles\n            if score > 0.1:\n                tiles.append({'image_id':file[:-len('.tiff')], 'tile':(i,j), 'x':x, 'y':y, 'level':level, 'size':size})\n    \n    # Close\n    biopsy.close()\n    sample = None\n    \nprint('Extract',len(tiles),'tiles from', len(files),'slides')\nt.elapsed_display()\n\n# OBSERVATION 1: ~1,5s to extract and score tiles on 1 slides\n# OBSERVATION 2: to keep an equivalent time, we could use multiprocessing on a larger dataset of 5000 slides.","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Generate tiles from files in a given range\ndef generate_tiles(file):\n    interval = 32\n    tiles = []\n    filepath = os.path.join(TRAIN_DIR, file)\n\n    # Open\n    biopsy = open_slide(filepath)\n    \n    # Read lowest definition image\n    level = biopsy.level_count - 1\n    dimensions = biopsy.level_dimensions[level]\n\n    # Get number of gridsquares in x and y direction\n    nx=int(dimensions[0]/interval)\n    ny=int(dimensions[1]/interval)\n    #tiles = np.zeros((nx, ny))\n\n    # Browse each tiles\n    level = 1\n    scale = 4\n    size = interval*scale # Tile size depend on scale factor\n    dimensions = (size, size)\n    num_pixels = dimensions[0]*dimensions[1]\n\n    for i in range(nx):\n        for j in range(ny):  \n            x, y = i*interval*16, j*interval*16 #Localisation from the level 0 => * max scale interval to get coordinate\n            sample = biopsy.read_region((x, y), level, dimensions)\n            sample = sample.convert(\"1\") #Convert to black and white\n            score = 1-np.count_nonzero(sample)/num_pixels\n            # Keep only not blank tiles\n            if score > 0.1:\n                tiles.append({'image_id':file[:-len('.tiff')], 'tile':(i,j), 'x':x, 'y':y, 'level':level, 'size':size})\n\n    # Close\n    biopsy.close()\n    sample = None\n    return tiles","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Extract tiles for all slides (with multiprocessing)\nt = Time()\nfiles = train_radboud5['train_file'].values\n\n# Processes available\nnum_processes = multiprocessing.cpu_count()\npool = multiprocessing.Pool(num_processes)\n\n# Image per process split\nnum_files = len(files)\nif num_processes > num_files:\n    num_processes = num_files\nfiles_per_process = num_files / num_processes\n\nprint(\"Number of processes: \" + str(num_processes))\nprint(\"Number of files: \" + str(num_files))\n\n# start tasks pooling\ntiles = []\n\ndef get_tiles(result):\n    tiles.append(result)\n    \nfor file in tqdm(files):\n    result = pool.apply_async(generate_tiles, args = (file,), callback = get_tiles)\n\npool.close()\npool.join()\n\ntiles = np.concatenate(tiles)\nprint('Extract',len(tiles),'tiles from', len(files),'slides')\nt.elapsed_display()\n\n# OBSERVATION 1: ~1,5s to extract and score tiles from 1 slides\n# OBSERVATION 2: ~10m to extract and score tiles of 961 slides from 1 process\n# OBSERVATION 3: ~5m to extract and score tiles of 961 slides from 4 process","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Export tiles for further usage\nt = Time()\ntiles_df = pd.DataFrame(tiles.tolist())\ntiles_df.to_csv('PANDA_tiles_EDA.csv')\nt.elapsed_display()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Import tiles dataframe for checkpoint purpose\ntiles_df = None # Free memory\n# tiles_df = pd.read_csv('/kaggle/input/panda-tiles-EDA/PANDA_tiles_EDA.csv', index_col=0)\n# tiles_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Improve tiles generation\nThat seem a good way to decompose each in tiles that could be then send to any model. But can we add more usefull information during this generation process in order to save time for further exploration? It is time to retrieve information from mask. As we work with Radboud slides, we will find details mask evaluated pixel by pixel. Like for white pixels, we can use the corresponding masks to add 6 scores proportion column.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Display a mask with tiles scoring\nt = Time()\nfiles = train_radboud5['label_file'].values[:1]\ntiles = []\ninterval = 32\n\nfor file in tqdm(files):\n    filepath = os.path.join(MASK_DIR, file)    \n    # Open\n    gleason = open_slide(filepath)\n    \n    # Read lowest definition image\n    level = gleason.level_count - 1\n    dimensions = gleason.level_dimensions[level]\n\n    # Get number of tiles in x and y direction\n    nx=int(dimensions[0]/interval)\n    ny=int(dimensions[1]/interval)\n    labels = np.zeros((nx, ny, 6)) #tiles with score dimension\n\n    # Browse each tiles\n    level = 1\n    scale = 4\n    size = interval*scale # Tile size depend on scale factor\n    dimensions = (size, size)\n    num_pixels = dimensions[0]*dimensions[1]\n    \n    for i in range(nx):\n        for j in range(ny):  \n            x, y = i*interval*16, j*interval*16 #Localization from the level 0 => * max scale interval to get coordinate\n            sample = gleason.read_region((x, y), level, dimensions)\n            sample = np.array(sample.convert('RGB'))\n            labels[i][j] = np.zeros(6, dtype='uint') #Create an empty score list\n            key, value = np.unique(sample[:,:,0], return_counts=True) # Count by pixel score present on the first color channel\n            scores = dict(zip(key, value)) #Create a score dict\n            for k in scores.keys():\n                labels[i][j][k] = scores[k]/num_pixels #Update score list\n    tiles.append(labels)\n    \n    # Close\n    gleason.close()\n    #sample = None\n    \nprint('Open, score tiles and close', len(files),'slides')\n\n# Generate a heatmap\ngrade = 5\ngrid = tiles[0][:,:,grade] # Get negative zone score\nsns.set(style=\"white\")\nplt.subplots(figsize=(grid.shape[0]/5, grid.shape[1]/5))\n\nmask = np.zeros_like(grid)\nmask[np.where(grid < 0.1)] = True #Mask threshold\n\nsns.heatmap(grid.T, square=True, linewidths=.5, mask=mask.T, cbar=False, vmin=0, vmax=1, cmap=\"Reds\")\nplt.show()\n\nprint('Tiles with Gleason grade', grade, '> 10% on', grid.size, 'tiles')\ngrid = None\n\nt.elapsed_display()\n# OBSERVATION: 2s to score tiles from 1 slide","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Sounds good! We can now add proportional gleason score on each tiles:","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# Generate score tiles from files and masks\ndef generate_tiles_labels(file):\n    interval = 32\n    tiles = []\n    filepath = os.path.join(TRAIN_DIR, file)\n    image_id = file[:-len(TRAIN_EXT)]\n    maskpath = os.path.join(MASK_DIR, image_id + MASK_EXT)\n\n    # Open files\n    biopsy = open_slide(filepath)\n    mask = open_slide(maskpath)\n    \n    # Read lowest definition image\n    level = biopsy.level_count - 1\n    dimensions = biopsy.level_dimensions[level]\n\n    # Get number of gridsquares in x and y direction\n    nx=int(dimensions[0]/interval)\n    ny=int(dimensions[1]/interval)\n    #tiles = np.zeros((nx, ny))\n\n    # Browse each tiles\n    level = 1\n    scale = 4\n    size = interval*scale # Tile size depend on scale factor\n    dimensions = (size, size)\n    num_pixels = dimensions[0]*dimensions[1]\n\n    for i in range(nx):\n        for j in range(ny):  \n            x, y = i*interval*16, j*interval*16 #Localization from the level 0 => * max scale interval to get coordinate\n            \n            # Read biopsy file\n            sample = biopsy.read_region((x, y), level, dimensions)\n            sample = sample.convert(\"1\") #Convert to black and white\n            score = 1-np.count_nonzero(sample)/num_pixels\n\n            # Keep only not blank tiles\n            if score > 0.1:\n                # Read mask file\n                sample = mask.read_region((x, y), level, dimensions)\n                sample = np.array(sample.convert('RGB'))\n                \n                key, value = np.unique(sample[:,:,0], return_counts=True) # Count by pixel score present on the first color channel\n                scores = dict(zip(key, value)) #Create a score dict\n                \n                PREFIX = 'gleason_'\n                labels = {PREFIX+str(k) : 0 for k in range(6)} #Create an empty score list from 0 to 5\n                for k in scores.keys():\n                    labels[PREFIX+str(k)] = scores[k]/num_pixels #Update score list\n                \n                # Add tile\n                tile = {'image_id':file[:-len('.tiff')], 'tile':(i,j), 'x':x, 'y':y, 'level':level, 'size':size,}   \n                tiles.append({**tile, **labels})\n\n    # Close\n    biopsy.close()\n    mask.close()\n    sample = None\n    return tiles","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Testing generate_tiles_labels\nfile = train_radboud5['train_file'].values[0]\ntiles = generate_tiles_labels(file)\ntest_df = pd.DataFrame(tiles)\ntest_df.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df = None # Free memory","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Extract label tiles for all slides (with multiprocessing)\nt = Time()\nfiles = train_radboud5['train_file'].values\n\n# Processes available\nnum_processes = multiprocessing.cpu_count()\npool = multiprocessing.Pool(num_processes)\n\n# Image per process split\nnum_files = len(files)\nif num_processes > num_files:\n    num_processes = num_files\nfiles_per_process = num_files / num_processes\n\nprint(\"Number of processes: \" + str(num_processes))\nprint(\"Number of files: \" + str(num_files))\n\n# start tasks pooling\ntiles = []\n\ndef get_tiles(result):\n    tiles.append(result)\n    \nfor file in tqdm(files):\n    result = pool.apply_async(generate_tiles_labels, args = (file,), callback = get_tiles)\n\npool.close()\npool.join()\n\ntiles = np.concatenate(tiles)\nprint('Extract',len(tiles),'tiles from', len(files),'slides')\nt.elapsed_display()\n\n# OBSERVATION 1: ~2s to extract and score tiles from 1 slides\n# OBSERVATION 2: ~6m to extract and score tiles from 964 slides on 4 process","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Use this to interrupt the multiprocess pool and free your CPUs after a 'Cancel run' command\npool.terminate()\npool.close()\npool.join()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Export tiles for further usage\nt = Time()\ntiles_labels_df = pd.DataFrame(tiles.tolist())\ntiles_labels_df.to_csv('PANDA_tiles_labels_EDA.csv')\nprint('Export dataframe')\nt.elapsed_display()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Import tiles with labels dataframe for checkpoint purpose\ntiles_labels_df = None # Free memory\n# tiles_labels_df = pd.read_csv('/kaggle/input/panda-tiles-labels_EDA/PANDA_tiles_labels_EDA.csv', index_col=0)\n# tiles_labels_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# II. Final code\nBellow, here is our final code to generate a tile+label dataframe on all 5060 Radboud slides. It can be run alone without all previous codes.\nThe overall runtime is about 30min. with 4 CPU process. So you can have a break, take a coffe or play with your children! \nYou may also run & reuse this separate [script](https://www.kaggle.com/huynhdoo/panda-wsi-tiles-preprocessing-script).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"# DEPENDANCIES ###########################################################################\nimport os\nimport glob\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport numpy as np\nfrom tqdm.notebook import tqdm\nimport openslide\nfrom openslide import OpenSlideError\nfrom IPython.display import Image\nimport seaborn as sns\nimport multiprocessing\nimport datetime\n# / DEPENDANCIES #########################################################################\nprint(\"Dependencies loaded\")\n\n# UTILITIES ##############################################################################\nclass Time:\n  \"\"\"\n  Class for displaying elapsed time.\n  \"\"\"\n\n  def __init__(self):\n    self.start = datetime.datetime.now()\n\n  def elapsed_display(self):\n    time_elapsed = self.elapsed()\n    print(\"Time elapsed: \" + str(time_elapsed))\n\n  def elapsed(self):\n    self.end = datetime.datetime.now()\n    time_elapsed = self.end - self.start\n    return time_elapsed\n# / UTILITIES ############################################################################\nprint('Utilities loaded')\n\n# PARAMETERS #############################################################################\nBASE_DIR = '/kaggle/input/prostate-cancer-grade-assessment/'\nOUTPUT_DIR = './'\nTRAIN_DIR = os.path.join(BASE_DIR, \"train_images\")\nTRAIN_EXT = \".tiff\"\nMASK_DIR = os.path.join(BASE_DIR, \"train_label_masks\")\nMASK_EXT = \"_mask.tiff\"\n\nprint(\"Parameters loaded\")\n# /PARAMETERS ############################################################################\n\n# DATASET ################################################################################\n# Get train/label slides ID\ntrain = glob.glob1(TRAIN_DIR, \"*\" + TRAIN_EXT)\nlabel = glob.glob1(MASK_DIR, \"*\" + MASK_EXT)\n\n# Keep only image_id\ntrain = [x[:-len(TRAIN_EXT)] for x in train]\nlabel = [y[:-len(MASK_EXT)] for y in label]\n\n# Add filenames to dataframe\ntrain_df = pd.read_csv(BASE_DIR + 'train.csv')\n\n# Add train file column for each existing file in train folder\ntrain_df['train_file'] = list(map(lambda x : x + TRAIN_EXT if x in set(train) else '', \n                              train_df['image_id']))\n# Add label file column for each existing file in mask folder\ntrain_df['label_file'] = list(map(lambda y : y + MASK_EXT if y in set(label) else '', \n                              train_df['image_id']))\n\n# Split dataframe by provider / we keep radboud scoring because their mask labels are more details\nprint('Dataframe original:', len(train_df))\ntrain_radboud = train_df[train_df['data_provider'] == 'radboud'].copy()\nprint('Dataframe after provider select:', len(train_radboud))\n# Keep only row with both train and label file\ntrain_radboud = train_radboud[train_radboud['train_file'] != '']\nprint('Dataframe after file select:', len(train_radboud))\ntrain_radboud = train_radboud[train_radboud['label_file'] != '']\nprint('Dataframe after label select:', len(train_radboud))\n\n# Release memory\ntrain_df = None\n# / DATASET ##############################################################################\n\n\n# FUNCTIONS ##############################################################################\n# Open a slide\ndef open_slide(filename):\n    \"\"\"\n    Open a whole-slide image (*.svs, etc).\n    :filename : Name of the slide file.\n    return: an OpenSlide object representing a whole-slide image.\n    \"\"\"\n    try:\n        slide = openslide.open_slide(filename)\n    except OpenSlideError:\n        slide = None\n    except FileNotFoundError:\n        slide = None\n    return slide    \n\n# Generate score tiles from files and masks\ndef generate_tiles_labels(file):\n    \"\"\"\n    Generate a list of tiles with coordonnate and label from file/mask whole-slide image .tiff\n    :file : Name of the slide file (must start and end with define directory and extension)\n    return: a list of dictionnary tiles with gleason labels\n    \"\"\"\n    interval = 32\n    tiles = []\n    filepath = os.path.join(TRAIN_DIR, file)\n    image_id = file[:-len(TRAIN_EXT)]\n    maskpath = os.path.join(MASK_DIR, image_id + MASK_EXT)\n\n    # Open files\n    biopsy = open_slide(filepath)\n    mask = open_slide(maskpath)\n    \n    # Read lowest definition image\n    level = biopsy.level_count - 1\n    dimensions = biopsy.level_dimensions[level]\n\n    # Get number of gridsquares in x and y direction\n    nx=int(dimensions[0]/interval)\n    ny=int(dimensions[1]/interval)\n    #tiles = np.zeros((nx, ny))\n\n    # Browse each tiles\n    level = 1\n    scale = 4\n    size = interval*scale # Tile size depend on scale factor\n    dimensions = (size, size)\n    num_pixels = dimensions[0]*dimensions[1]\n\n    for i in range(nx):\n        for j in range(ny):  \n            x, y = i*interval*16, j*interval*16 #Localization from the level 0 => * max scale interval to get coordinate\n            \n            # Read biopsy file\n            sample = biopsy.read_region((x, y), level, dimensions)\n            sample = sample.convert(\"1\") #Convert to black and white\n            score = 1-np.count_nonzero(sample)/num_pixels #Normalize the value between 0 and 1 (0=white, 1=black)\n\n            # Keep only not empty tiles\n            if score > 0.1:\n                # Read mask file\n                sample = mask.read_region((x, y), level, dimensions)\n                sample = np.array(sample.convert('RGB'))\n                \n                key, value = np.unique(sample[:,:,0], return_counts=True) # Count by pixel score present on the first color channel\n                scores = dict(zip(key, value)) #Create a score dict\n                \n                PREFIX = 'gleason_'\n                labels = {PREFIX+str(k) : 0 for k in range(6)} #Create an empty score list from 0 to 5\n                for k in scores.keys():\n                    labels[PREFIX+str(k)] = scores[k]/num_pixels #Update score list\n                \n                # Add tile\n                tile = {'image_id':file[:-len('.tiff')], 'tile':(i,j), 'x':x, 'y':y, 'level':level, 'size':size,}   \n                tiles.append({**tile, **labels})\n\n    # Close\n    biopsy.close()\n    mask.close()\n    sample = None\n    return tiles\nprint('Functions loaded')\n# / FUNCTIONS ############################################################################\n\nt = Time() # Launch timer\n\n# EXTRACTION #############################################################################\n# Extract label tiles for all slides (with multiprocessing)\nprint('Start tiles generation...')\nfiles = train_radboud['train_file'].values\n\n# Processes available\nnum_processes = multiprocessing.cpu_count()\npool = multiprocessing.Pool(num_processes)\n\n# Image per process split\nnum_files = len(files)\nif num_processes > num_files:\n    num_processes = num_files\nfiles_per_process = num_files / num_processes\n\nprint(\"Number of processes: \" + str(num_processes))\nprint(\"Number of files: \" + str(num_files))\n\n# start tasks pooling\ntiles = []\n\ndef get_tiles(result):\n    tiles.append(result)\n    \nfor file in files:\n    result = pool.apply_async(generate_tiles_labels, args = (file,), callback = get_tiles)\n\npool.close()\npool.join()\n\ntiles = np.concatenate(tiles)\nprint('Extract',len(tiles),'tiles from', len(files),'slides')\n\n# OBSERVATION: ~30m to extract and score 603602 tiles from 5060 slides on 4 process\n# / EXTRACTION ###########################################################################\n\n# OUTPUT #################################################################################\n# Export tiles for further usage\ntiles_final_df = pd.DataFrame(tiles.tolist())\ntiles_final_df.to_csv(OUTPUT_DIR + 'PANDA_tiles_labels_final.csv')\nprint('Tiles exported')\n\n# OBSERVATION: output csv file of size 72,4Mb\n# / OUTPUT ###############################################################################\n\nt.elapsed_display() # Print timer","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Use this to interrupt the multiprocess pool and free your CPUs after a 'Cancel run' command\npool.terminate()\npool.close()\npool.join()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tiles_final_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tiles_final_df.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"raw","source":"","execution_count":null}],"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}