5

I am trying to generate heat map, or probability map, for Whole Slide Images (WSIs) using probability values. I have coordinate points (which determine areas on the WSIs) and corresponding probability values.

Basic Introduction on WSI: WSIs are large is size (almost 100000 x 100000 pixels). Hence, can't open these images using normal image viewer. The WSIs are processed using OpenSlide software.

I have seen previous posts in Stack Overflow on related to heat map, but as WSIs are processed in a different way, I am unable to figure out how to apply these solutions. Some examples that I followed: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, etc.

Js541
  • 100
  • 1
  • 13
  • Have you tried reading in the image in tiles, creating the heat map for each tile independently, and saving the tiles back to file? – Cris Luengo May 14 '19 at 06:15
  • @CrisLuengo Yes, I can create heat map for each tile independently. But I don't have idea how can I stitch all the tiles together or, saving the tiles back to file. Please suggest. – Js541 May 14 '19 at 09:41
  • I was going to suggest you write a tiled TIFF file out, which is not too complicated, but you are likely best off using the vips library. It has a Python binding: https://jcupitt.github.io/libvips/2017/09/01/libvips-for-Python.html — I don’t know how to use it, so won’t write a full answer, but I’m sure you can figure it out, I hear it’s easy to use. – Cris Luengo May 14 '19 at 12:47
  • Yes, I know. but I am looking for a specific answer. Thanks – Js541 May 14 '19 at 21:03

2 Answers2

4

To generate heat map on WSIs, follow below instructions:

First of all Extract image patches and save the coordinates. Use below code for patch extraction. The code require some changes as per the requirements. The code has been copied from: patch extraction code link

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import argparse
import logging
try:
    import Image
except:
    from PIL import Image
import math
import numpy as np
import openslide
import os
from time import strftime,gmtime

parser = argparse.ArgumentParser(description='Extract a series of patches from a whole slide image')
parser.add_argument("-i", "--image", dest='wsi',  nargs='+', required=True, help="path to a whole slide image")
parser.add_argument("-p", "--patch_size", dest='patch_size', default=299, type=int, help="pixel width and height for patches")
parser.add_argument("-b", "--grey_limit", dest='grey_limit', default=0.8, type=float, help="greyscale value to determine if there is sufficient tissue present [default: `0.8`]")
parser.add_argument("-o", "--output", dest='output_name', default="output", help="Name of the output file directory [default: `output/`]")

parser.add_argument("-v", "--verbose",
            dest="logLevel",
            choices=['DEBUG', 'INFO', 'WARNING', 'ERROR', 'CRITICAL'],
            default="INFO",
            help="Set the logging level")
args = parser.parse_args()

if args.logLevel:
        logging.basicConfig(level=getattr(logging, args.logLevel))

wsi=' '.join(args.wsi)


""" Set global variables """
mean_grey_values = args.grey_limit * 255
number_of_useful_regions = 0
wsi=os.path.abspath(wsi)
outname=os.path.abspath(args.output_name)
basename = os.path.basename(wsi)
level = 0

def main():
    img,num_x_patches,num_y_patches = open_slide()
    logging.debug('img: {}, num_x_patches = {}, num_y_patches: {}'.format(img,num_x_patches,num_y_patches))

    for x in range(num_x_patches):
        for y in range(num_y_patches):
            img_data = img.read_region((x*args.patch_size,y*args.patch_size),level, (args.patch_size, args.patch_size))
            print_pics(x*args.patch_size,y*args.patch_size,img_data,img)

    pc_uninformative = number_of_useful_regions/(num_x_patches*num_y_patches)*100
    pc_uninformative = round(pc_uninformative,2)
    logging.info('Completed patch extraction of {} images.'.format(number_of_useful_regions))
    logging.info('{}% of the image is uninformative\n'.format(pc_uninformative))


def print_pics(x_top_left,y_top_left,img_data,img):
    if x_top_left % 100 == 0 and y_top_left % 100 == 0 and x_top_left != 0:
        pc_complete = round(x_top_left /img.level_dimensions[0][0],2) * 100
        logging.info('{:.2f}% Complete at {}'.format(pc_complete,strftime("%a, %d %b %Y %H:%M:%S +0000", gmtime())))
        exit()

    img_data_np = np.array(img_data)
    """ Convert to grayscale"""
    grey_img = rgb2gray(img_data_np)
    if np.mean(grey_img) < mean_grey_values:
        logging.debug('Image grayscale = {} compared to threshold {}'.format(np.mean(grey_img),mean_grey_values))
        global number_of_useful_regions
        number_of_useful_regions += 1
                wsi_base = os.path.basename(wsi)
                wsi_base = wsi_base.split('.')[0]
        img_name = wsi_base + "_" + str(x_top_left) + "_" + str(y_top_left) + "_" + str(args.patch_size)
        #write_img_rotations(img_data_np,img_name)
                logging.debug('Saving {} {} {}'.format(x_top_left,y_top_left,np.mean(grey_img)))
        save_image(img_data_np,1,img_name)

def gen_x_and_y(xlist,ylist,img):
    for x in xlist:
        for y in ylist:
            img_data = img.read_region((x*args.patch_size,y*args.patch_size),level, (args.patch_size, args.patch_size))
            yield (x, y,img_data)

def open_slide():
    """
    The first level is always the main image
    Get width and height tuple for the first level
    """

    logging.debug('img: {}'.format(wsi))

    img = openslide.OpenSlide(wsi)
    img_dim = img.level_dimensions[0]

    """
    Determine what the patch size should be, and how many iterations it will take to get through the WSI
    """
    num_x_patches = int(math.floor(img_dim[0] / args.patch_size))
    num_y_patches = int(math.floor(img_dim[1] / args.patch_size))

    remainder_x = img_dim[0] % num_x_patches
    remainder_y = img_dim[1] % num_y_patches

    logging.debug('The WSI shape is {}'.format(img_dim))
    logging.debug('There are {} x-patches and {} y-patches to iterate through'.format(num_x_patches,num_y_patches))
    return img,num_x_patches,num_y_patches

def validate_dir_exists():
    if os.path.isdir(outname) == False:
        os.mkdir(outname)
    logging.debug('Validated {} directory exists'.format(outname))
    if os.path.exists(wsi):
        logging.debug('Found the file {}'.format(wsi))
    else:
        logging.debug('Could not find the file {}'.format(wsi))
        exit()

def rgb2gray(rgb):
    """Converts an RGB image into grayscale """
    r, g, b = rgb[:,:,0], rgb[:,:,1], rgb[:,:,2]
    gray = 0.2989 * r + 0.5870 * g + 0.1140 * b
    return gray

def save_image(img,j,img_name):
        tmp = os.path.join(outname,img_name+"_"+str(j)+".png")
    try:
        im = Image.fromarray(img)
        im.save(tmp)
    except:
        print('Could not print {}'.format(tmp))
                exit()
if __name__ == '__main__':
    validate_dir_exists()
    main()

Secondly, generate the probability values of each patches. Finally, replace all the pixel values within a coordinates with the corresponding probability values and display the results using color maps.

This is the basic idea of generating heat map on WSIs. You can modify the code and concept to get a heat map as per your wish.

user1410665
  • 719
  • 7
  • 23
1

We have developed a python package for processing whole-slide-images: https://github.com/amirakbarnejad/PyDmed Here is a tutorial for getting heatmaps for whole-slide-images: https://amirakbarnejad.github.io/Tutorial/tutorial_section5.html. Also here is a sample notebook that gets heatmaps for WSIs using PyDmed: Link to the sample notebook.

The benefit of PyDmed is that it is multi-processed. The dataloader sends a stream of patches to GPU(s), and the StreamWriter writes to disk in a separate process. Therefore, it is highly efficient. The running time of course depends on the machine, the size of WSIs, etc. On a good machine with a good GPU, PyDmed can generate heatmaps for ~120 WSIs in one day.

akbarnejad
  • 29
  • 2