Stitch together tiles of a slide exported from PE Harmony

stitch_harmony_export.py

import re
from pathlib import Path
from xml.etree import ElementTree as ET

import numpy as np
from skimage import exposure
from skimage.io import imread, imsave, imshow
from skimage.filters import threshold_otsu, gaussian


NS = {"harmony": "http://www.perkinelmer.com/PEHH/HarmonyV5"}
def parse_image_positions(xml_file):
    rootpath = Path(xml_file).parent
    
    root = ET.parse(xml_file).getroot()
    
    wells = root.findall("./*/harmony:Well", NS)
    positions = {}
    for well in wells:
        r = int(well.find("harmony:Row", NS).text)
        c = int(well.find("harmony:Col", NS).text)
        positions[(r, c)] = {}        
    
    images = root.findall("./*/harmony:Image", NS)
    for image in images:
        x = float(image.find("harmony:PositionX", NS).text)
        y = float(image.find("harmony:PositionY", NS).text)
        r = int(image.find("harmony:Row", NS).text)
        c = int(image.find("harmony:Col", NS).text)
        url = rootpath / (image.find("harmony:URL", NS).text)
        img_w = int(image.find("harmony:ImageSizeX", NS).text)
        img_h = int(image.find("harmony:ImageSizeY", NS).text)

        positions[(r, c)][(x, y)] = url
    
    sorted_positions = [
        dict(((x, y), p) for (x, y), p in sorted(d.items(), key=lambda t: (-t[0][1], t[0][0])))
        for d in positions.values()
    ]
    
    widths = []
    for d in sorted_positions:
        xs = [x for x,y in d.keys()]
        N = np.argmax(xs) + 1
        widths += [N]
    return widths, sorted_positions


def create_img_stack(width, p_dict):
    imgs = []
    sub = []
    for (_, path) in p_dict.items():
        sub.append(imread(str(path)))
        if len(sub) == width:
            imgs.append(sub)
            sub = []
    return np.block(imgs)


def fix_fiducials(img, q1=40, q2=51):
    """
    I wrote this to bring out *extremely* faint features in the images (these were 10X Visium fiducials).
    """
    c1, c2 = img.copy(), img.copy()
    c1 = exposure.adjust_gamma(c1, gamma=0.5, gain=10)
    p1 = np.percentile(c1, q1)
    p2 = np.percentile(c1, q2)
    c1[c1 < p1] = p1
    c2 = exposure.rescale_intensity(c2, out_range=(0, 2**16))
    g = gaussian(c2, sigma=100)
    inds = c1 > p2
    inds &= g > threshold_otsu(g)
    c2[inds] = exposure.rescale_intensity(c1[inds], out_range=(0, 2**16))
    return c2.astype(np.uint16)


def save_image(path, img):
    imsave(path, img.astype(np.uint16), bigtiff=False)
   
    
if __name__ == "__main__":
    slide_widths, slide_positions = parse_image_positions(
        "path/to/harmony_export/Images/Index.idx.xml"
    )
    slide_img = create_img_stack(slide_widths[0], slide_positions[0])
    
    save_image(slide_img, "path/to/slide.tiff")
    save_image(fix_fiducials(slide_img, q1=33, q2=37), "path/to/slide_with-fid.tiff")