Skip to content

Instantly share code, notes, and snippets.

@wflynny
Last active September 17, 2020 20:17
Show Gist options
  • Save wflynny/2e385a7de6361d1d471b4dfd66100a21 to your computer and use it in GitHub Desktop.
Save wflynny/2e385a7de6361d1d471b4dfd66100a21 to your computer and use it in GitHub Desktop.
Stitch together tiles of a slide exported from PE Harmony
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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment