Renders RGB JPEG pyramid tiles from OME-TIFF pyramids for the melanoma triplet dataset

composit_color_image.py

import numpy as np
from matplotlib import colors
from matplotlib import pyplot as plt
from PIL import Image


def composite_channel(target, image, color, range_min, range_max):
    ''' Render _image_ in pseudocolor and composite into _target_
    Args:
        target: Numpy float32 array containing composition target image
        image: Numpy uint16 array of image to render and composite
        color: Color as r, g, b float array, 0-1
        range_min: Threshhold range minimum, 0-65535
        range_max: Threshhold range maximum, 0-65535
    '''
    f_image = (image.astype('float32') - range_min) / (range_max - range_min)
    f_image = f_image.clip(0,1, out=f_image)
    for i, component in enumerate(color):
        target[:, :, i] += f_image * component

# Create single channel images with numpy, users load their single-channel
# images using method of preferences
channel_1 = np.roll(np.pad(np.ones((5, 5)), 4), 2, axis=0)
channel_2 = np.roll(np.pad(np.ones((5, 5)), 4), -2, axis=0)
channel_3 = np.roll(np.pad(np.ones((5, 5)), 4), 2, axis=1)

# Create array to hold the output image
target = np.zeros(channel_1.shape + (3, ), np.float32)

# Composite channels
composite_channel(target, channel_1, colors.to_rgb('#ff0000'), 0, 1)
composite_channel(target, channel_2, colors.to_rgb('lime'), 0, 1)
composite_channel(target, channel_3, colors.to_rgb('blue'), 0, 1)

# Visual verification
plt.figure(); plt.imshow(target)

# Clear out of range numbers
np.clip(target, 0, 1, out=target)

# Next line does gamma correction which generate more perceptually accurate 
# color gradient. Comment it out if you need to match the result image with 
# OMERO/Pathviewer.
np.power(target, 1/2.2, out=target)

# Convert to 8-bit/channel RGB image and save
target_u8 = (target * 255).astype(np.uint8)
img = Image.frombytes('RGB', target.T.shape[1:], target_u8.tobytes())
img.save('out.png')

melanoma_triplet_render.py

import re
import csv

share_link = 'https://omero.hms.harvard.edu/webgateway/img_detail/1203770/?c=-1|3000:40000$0000FF,-2|2000:30000$00FF00,-3|2000:30000$FFFFFF,-4|2000:30000$FF0000,-5|3000:45000$0000FF,-6|2000:30000$00FF00,-7|2000:30000$FFFFFF,-8|2000:30000$FF0000,-9|3000:50000$0000FF,-10|4000:16000$00FF00,-11|2500:15000$FFFFFF,-12|4000:16000$FF0000,13|3000:45000$0000FF,14|3000:12000$00FF00,15|4000:50000$FFFFFF,16|3000:20000$FF0000,-17|3000:50000$0000FF,-18|3000:12000$00FF00,-19|5000:50000$FFFFFF,-20|3000:12000$FF0000,-21|3000:50000$0000FF,-22|4000:50000$00FF00,-23|2500:12000$FFFFFF,-24|2500:12000$FF0000,-25|3000:45000$0000FF,-26|3000:10000$00FF00,-27|2000:10000$FFFFFF,-28|1500:8000$FF0000,-29|3000:45000$0000FF,-30|3000:25000$00FF00,-31|2000:12000$FFFFFF,-32|4000:20000$FF0000,-33|3000:35000$0000FF,-34|2000:12000$00FF00,-35|2000:15000$FFFFFF,-36|2500:12500$FF0000,-37|2500:25000$0000FF,-38|2000:12000$00FF00,-39|1500:9000$FFFFFF,-40|5000:50000$FF0000&m=c&p=normal&ia=0&q=0.9&t=1&z=1&zm=3.126150550590755&x=14938.5&y=14090.81209850107&maps=[{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}},{"inverted":{"enabled":false}}]'
marker_names = 'Hoechst0, AF488, AF555, AF647, Hoechst1, A488, A555, A647, Hoechst2, CD3, NaKATPase, CD45RO, Hoechst3, Ki67, panKeratin, aSMA, Hoechst4, CD4, CD45, PD-1, Hoechst5, CD20, CD68, CD8a, Hoechst6, CD163, FOXP3, PD-L1, Hoechst7, E-Cadherin, Vimentin, CDX2, Hoechst8, LaminABC, Desmin, CD31, Hoechst9, PCNA, Ki67, Collagen'

def parseChannelsInfos(url_link, marker_list):
    markersNamesList = marker_list.split(', ')
    channelsInfos = []
    for text in url_link.split('|'):
        m = re.match(r"(\d+):(\d+)\$(\w+)", text)
        if m:
            ms = re.split(r":|\$", m.group(0))
            channelsInfos.append({'Low': ms[0], 'High': ms[1], 'Color': '#{}'.format(ms[2])})
    for idx, marker in enumerate(markersNamesList):
        channelsInfos[idx].update({ 'Marker Name': marker, 'Channel Number': idx })
    return channelsInfos
channels_infos = parseChannelsInfos(share_link, marker_names)

with open('render_setting.csv', 'w') as f:
    writer = csv.DictWriter(f, channels_infos[0].keys())
    writer.writeheader()
    writer.writerows(channels_infos)

import sys
import itertools
try:
    import pathlib
except ImportError:
    import pathlib2 as pathlib
import json
import numpy as np
import pytiff
from PIL import Image

def composite_channel(target, image, color, range_min, range_max):
    ''' Render _image_ in pseudocolor and composite into _target_
    Args:
        target: Numpy float32 array containing composition target image
        image: Numpy uint16 array of image to render and composite
        color: Color as r, g, b float array, 0-1
        range_min: Threshhold range minimum, 0-65535
        range_max: Threshhold range maximum, 0-65535
    '''
    f_image = (image.astype('float32') - range_min) / (range_max - range_min)
    f_image = f_image.clip(0,1, out=f_image)
    for i, component in enumerate(color):
        target[:, :, i] += f_image * component

        
TILE_SIZE = 1024
EXT = 'jpg'



# List markers as: green, white, red
RENDER_GROUPS = [list(range(36))[i*4:(i+1)*4] for i in range(9)]
COLORS = [
    [0, 0, 1], # blue
    [0, 1, 0], # green
    [1, 1, 1], # white
    [1, 0, 0], # red
]

input_file_path = pathlib.Path('TMA0809.ome.tif')
output_path = pathlib.Path('./')

tiff = pytiff.Tiff(str(input_file_path))

num_channels = len(channels_infos)
assert tiff.number_of_pages % num_channels == 0, "Pyramid/channel mismatch"

num_levels = tiff.number_of_pages // num_channels

for level in range(num_levels):
# for level in range(num_levels)[4:]:

    page_base = level * num_channels
    tiff.set_page(page_base)
    ny = int(np.ceil(tiff.shape[0] / TILE_SIZE))
    nx = int(np.ceil(tiff.shape[1] / TILE_SIZE))
    print(f"level {level} ({ny} x {nx})")

    for ty, tx in itertools.product(range(0, ny), range(0, nx)):

        iy = ty * TILE_SIZE
        ix = tx * TILE_SIZE
        filename = f'{level}_{tx}_{ty}.{EXT}'
#         print(f"    {filename}")

        for idx, marker_list in enumerate(RENDER_GROUPS):
            group_dir = pathlib.Path('---'.join(
                str(x) + '___' + channels_infos[x]['label'] for x in marker_list
            ))
            (output_path / group_dir).mkdir(exist_ok=True)
#             print(f"        {group_dir}")
            for i, (marker, color) in enumerate(zip(marker_list, COLORS)):
                ch_info = channels_infos[marker]
                tiff.set_page(page_base + marker)
                tile = tiff[iy:iy+TILE_SIZE, ix:ix+TILE_SIZE]
                if i == 0:
                    target = np.zeros(tile.shape + (3,), np.float32)
                composite_channel(
                    target, tile, color, float(ch_info['start']), float(ch_info['end'])
                )
            np.clip(target, 0, 1, out=target)
#             np.power(target, 1/2.2, out=target)
            target_u8 = (target * 255).astype(np.uint8)
            img = Image.frombytes('RGB', target.T.shape[1:], target_u8.tobytes())
            img.save(str(output_path / group_dir / filename), quality=85)