Using pyvips to stitch component data images generated by inform

UnInform.py

# -*- coding: utf-8 -*-
"""
This script is used for stitching the component_data.tif files generated by un-mixing workflow in inForm. This
script uses tifffile to read the metadata (Channel names & colors, spatial calibration and stage cordinates) from
input files. pyvips (libvips) is used to read the pixel data from image files and create a stitched image. The final stitched image is saved
in ome-tiff compatible format with minimal ome-xml.
Created on Wed Aug 22 15:53:08 2021
Modified on 2023-07-18: add resolution, channel names and channel colors to ome-xml. some code cleanup.
@author: Ajay Zalavadia (zalavadia.ajay@gmail.com)
"""
import math
import os
import time
import glob
from xml.etree import ElementTree as ET
import tifffile

# Path to libvips binaries
vipshome = 'c:\\Apps\\vips-dev-8.14\\bin'
os.environ['PATH'] = vipshome + ';' + os.environ['PATH']
import pyvips
print("vips version: " + str(pyvips.version(0))+"."+str(pyvips.version(1))+"."+str(pyvips.version(2)))

imagedir= "C:\\Component-Data"
outFile = os.path.join(imagedir, 'stitched.ome.tif')

class region:
    def __init__(self, XPosition, YPosition, Height, Width, filePath, XCord, YCord):
        self.XPosition = XPosition # Stage position
        self.YPosition = YPosition # Stage position
        self.Height = Height
        self.Width = Width
        self.filePath = filePath
        self.XCord = XCord # pixel position on the stitched image
        self.YCord = YCord # pixel position on the stitched image

def get_region(filePath):
    with tifffile.TiffFile(filePath) as tif:
        tag_xpos = tif.pages[0].tags['XPosition']
        tag_ypos = tif.pages[0].tags['YPosition']
        tag_xres = tif.pages[0].tags['XResolution']
        tag_yres = tif.pages[0].tags['YResolution']
        
        xpos = 10000*tag_xpos.value[0]/tag_xpos.value[1]
        xres = tag_xres.value[0]/(tag_xres.value[1]*10000)
        ypos = 10000*tag_ypos.value[0]/tag_ypos.value[1]
        yres = tag_yres.value[0]/(tag_yres.value[1]*10000)
        height = tif.pages[0].tags['ImageLength'].value
        width = tif.pages[0].tags['ImageWidth'].value

    rg = region(xpos*xres, ypos*yres, height, width, filePath, xpos*xres, ypos*yres)
    return rg 

#collect regions 
def collect_regions(imagedir):
    fileList = glob.glob(imagedir + '/*component_data.tif')
    region_collection=[];
    for idx, file in enumerate(fileList):
        region_collection.append(get_region(file))
    return region_collection

#Calculate image cordinates from the stage cordinates
def zero_center_regions(region_collection):
    minX = min(region.XPosition for region in region_collection)
    minY = min(region.YPosition for region in region_collection)
    for idx, region in enumerate(region_collection):
        region.XCord = region.XPosition-minX
        region.YCord = region.YPosition-minY
        
    return region_collection

#Return the number of channels from first image
def count_channels(region_collection):
    temp = pyvips.Image.new_from_file(region_collection[0].filePath, access='sequential')
    channel_count = int(temp.get("n-pages"))-1 # get n-pages from 1st image, -1 for thumbnail
    return channel_count

def interleaved_tile(filePath, channels):
    tile = pyvips.Image.new_from_file(filePath, n=channels)
    page_height = tile.get("page-height")
    # chop into pages
    pages = [tile.crop(0, y, tile.width, page_height) for y in range(0, tile.height, page_height)]
    # join pages band-wise to make an interleaved image
    tile = pages[0].bandjoin(pages[1:])
    tile = tile.copy(interpretation="multiband")
    return tile

#Stitch regions and return stitched image
def stitch_regions(region_collection):
    channels = count_channels(region_collection)
    stitched_img_size_x = int(math.ceil(max(region.XCord + region.Width for region in region_collection)))
    stitched_img_size_y = int(math.ceil(max(region.YCord + region.Height for region in region_collection)))
    print('Final image dimension Width: ', stitched_img_size_x, ' Height: ', stitched_img_size_y)
    # create a reference image of the final size
    stitched_img = pyvips.Image.black(stitched_img_size_x, stitched_img_size_y, bands=channels)

    for idx, region in enumerate(region_collection):
        tile = interleaved_tile(region.filePath,stitched_img.bands)
        # insert tile to the reference image
        stitched_img = stitched_img.insert(tile, region.XCord, region.YCord, expand=1, background=[0])

    return stitched_img

def get_colors(colors):
    new_colors = []
    for cl in colors:
        cl = str.split(cl, ',')
        RGBint = int.from_bytes([int(cl[0]), int(cl[1]), int(cl[2]), 1], byteorder="big", signed=True)
        new_colors.append(RGBint)
    return new_colors

def channel_xml(region_collection):
    img = tifffile.TiffFile(region_collection[0].filePath)
    imgData_s0 = img.series[0]
    names = [(ET.fromstring(page.description).find('Name').text) for page in imgData_s0]
    colors = [(ET.fromstring(page.description).find('Color').text) for page in imgData_s0]
    colors = get_colors(colors)
    channels_xml = ('\n').join(
        [f"""<Channel ID="Channel:0:{i}" Name="{channel_name}" Color="{colors[i]}" SamplesPerPixel="1"/>"""
         for i, channel_name in enumerate(names)]
    )
    return channels_xml

def generate_ome_xml(region_collection,final_width,final_height,bands):
    img = tifffile.TiffFile(region_collection[0].filePath)
    tag_xres = img.series[0].pages[0].tags['XResolution']
    micron_per_pixel = 1 / (tag_xres.value[0] / (tag_xres.value[1] * 10000))
    CalibrationUnits = 'µm'
    channel_ome_xml = channel_xml(region_collection)

    ome_xml=f"""<?xml version="1.0" encoding="UTF-8"?>
    <OME xmlns="http://www.openmicroscopy.org/Schemas/OME/2016-06"
        xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
        xsi:schemaLocation="http://www.openmicroscopy.org/Schemas/OME/2016-06 http://www.openmicroscopy.org/Schemas/OME/2016-06/ome.xsd">
        <Image ID="Image:0">
            <!-- Minimum required fields about image dimensions -->
            <Pixels DimensionOrder="XYCZT"
                    ID="Pixels:0"
                    SizeC="{bands}"
                    SizeT="1"
                    SizeX="{final_width}"
                    SizeY="{final_height}"
                    SizeZ="1"
                    Type="float"
                    PhysicalSizeX="{micron_per_pixel}"
                    PhysicalSizeXUnit="{CalibrationUnits}"
                    PhysicalSizeY="{micron_per_pixel}"
                    PhysicalSizeYUnit="{CalibrationUnits}">
                    {channel_ome_xml}
            </Pixels>
        </Image>
    </OME>"""
    return ome_xml

def save_ome_tiff(stitched_img, outFile, ome_xml, compression):
    height = stitched_img.height

    # OME compatible image
    stitched_img = pyvips.Image.arrayjoin(stitched_img.bandsplit(), across=1)
    # Make a copy and set tiff tags necessary for OME-TIFF
    stitched_img = stitched_img.copy()
    stitched_img.set_type(pyvips.GValue.gint_type, "page-height", height)
    # Include metadata as OME-XML for OME-TIFF
    stitched_img.set_type(pyvips.GValue.gstr_type, "image-description", ome_xml)
    stitched_img.write_to_file(outFile, compression=compression, tile=True, tile_width=512, tile_height=512,  pyramid=True, subifd=True, bigtiff=True)
    return 0

# Collect regions from the list of component_data.tif files
t0 = time.time()
regions = zero_center_regions(collect_regions(imagedir))
t1 = time.time()

#Stitch regions
stitched_img = stitch_regions(regions)
t2 = time.time()

# Create minimal ome-xml
ome_xml =generate_ome_xml(regions,stitched_img.width,stitched_img.height,stitched_img.bands)

# Save as ome-tiff
save_ome_tiff(stitched_img, outFile, ome_xml, "lzw")
t3= time.time()

print('Time to collect regions: ', t1-t0)
print('Time to stitch: ', t2-t1)
print('Time to save (sec): ', t3-t2)