Python script to create a Max Intensity Projection for an OMERO Image

max_intensity_projection.py

# see https://forum.image.sc/t/maximum-intensity-projection-feature-in-omero-web/54392/5

import argparse
import sys
import omero.clients
from omero.cli import cli_login
from omero.rtypes import wrap
from omero.gateway import BlitzGateway
import numpy as np


def main(argv):
    parser = argparse.ArgumentParser()
    parser.add_argument('image_id', type=int, help='Image ID')
    args = parser.parse_args(argv)

    with cli_login() as cli:
        conn = BlitzGateway(client_obj=cli._client)
        conn.SERVICE_OPTS.setOmeroGroup(-1)

        image_id = args.image_id
        image = conn.getObject("Image", image_id)
        print("Image", image_id, image.name)

        group_id = image.getDetails().group.id.val
        print('group_id', group_id)
        conn.SERVICE_OPTS.setOmeroGroup(group_id)

        dataset = None
        parents = list(image.listParents())
        if len(parents) > 0:
            dataset = parents[0]._obj
            print("Dataset", dataset.id.val)

        sizeZ = image.getSizeZ()
        sizeC = image.getSizeC()
        sizeT = image.getSizeT()

        clist = range(sizeC)
        pixels = image.getPrimaryPixels()

        def planeGen():
            for c in clist:
                for t in range(sizeT):
                    # get Z-stack...
                    print(f"C: {c}, T: {t}")
                    zctList = [(z,c,t) for z in range(sizeZ)]
                    # planes is a generator - no data loaded yet...
                    planes = pixels.getPlanes(zctList)
                    data = []
                    for p in planes:
                        # could add a sleep here if you want to reduce load on server?
                        data.append(p)
                    z_stack = np.stack(data, axis=0)
                    # return the Max-Intensity projection for C and T
                    yield np.amax(z_stack, axis=0)

        # Use sourceImageId to copy channels metadata etc.
        new_image = conn.createImageFromNumpySeq(
            planeGen(), f'{image.name}_proj', sizeZ=1, sizeC=sizeC, sizeT=sizeT,
            sourceImageId=image_id, channelList=clist, dataset=dataset)

        print("Projected Image", new_image)

if __name__ == '__main__':
    main(sys.argv[1:])