Python writing and appending numpy on ome.tiff files

Hi.

I have a problem saving multi-dimensional data (and the metadata) to a bigtiff ome-tif.

I’m doing GPU deconvolution with flowdec from the hammerlab on my live cell imaging files are 1.5TB acquired over two days in my lab’s Nikon microscope.
I open the file with pims (either ND2 reader or bio-formats) - to have sequential reading from disk as it’s impossible to have the whole set on RAM. (I either have Nikon’s metadata, or bioformats metada format).

frames =  pims.ND2_Reader('./antibody4.nd2')
#frames.iter_axes = 't'  # 't' is the default already
frames.bundle_axes = 'zyx'  # when 'z' is available, this will be default
meta = frames.metadata
print(meta)
{'width': 2048, 'width_bytes': 20480, 'height': 2044, 'components': 5, 'bitsize_memory': 16, 'bitsize_significant': 16, 'sequence_count': 88, 'tile_width': 2048, 'tile_height': 2044, 'compression': None, 'compression_quality': 4294967197, 'plane_count': 5, 'angle': -1.5867835237991448, 'calibration_um': 0.108333333333333, 'time_start_jdn': 2459012.426468472, 'time_start': datetime.datetime(2020, 6, 11, 17, 14, 6, 875980), 'time_start_utc': datetime.datetime(2020, 6, 11, 22, 14, 6, 875980), 'objective': 'Plan Apo λ 60x Oil', 'magnification': -1.0, 'NA': 1.4, 'refractive_index1': 1.515, 'refractive_index2': 1.515, 'pinhole': 0.0, 'zoom': 1.0, 'projective_mag': -1.0, 'image_type': 'normal', 'z_home': 5, 'plane_0': {'components': 1, 'rgb_value': (0.0, 0.0, 1.0), 'name': 'DAPI', 'oc': 'DAPI', 'emission_nm': 535.0}, 'plane_1': {'components': 1, 'rgb_value': (0.0, 1.0, 0.0), 'name': 'GFP', 'oc': '', 'emission_nm': 535.0}, 'plane_2': {'components': 1, 'rgb_value': (1.0, 0.3333333333333333, 0.0), 'name': 'mRuby', 'oc': '', 'emission_nm': 620.0}, 'plane_3': {'components': 1, 'rgb_value': (0.9058823529411765, 0.0, 0.0), 'name': 'a647', 'oc': '', 'emission_nm': 670.0}, 'plane_4': {'components': 1, 'rgb_value': (0.3568627450980392, 1.0, 0.0), 'name': 'mneongreen', 'oc': '', 'emission_nm': 535.0}}

Now this is the problem:
I do deconvolution at an amazing speed and quality with the GPU. But then I need to save the data back to an OME.tiff for every channel and basically appending to the file as I go along(and when I get this going, I plan to do it for every time point and every visit point)
(data in example is 11 Z’s, 2 channels, 2044x2048 XY.)
The example below, is the attempt I did with python-bioformats, and tried with ‘TCZYX’ or to transpose to XYZCT. The file saves, adds and appends data (I see the file growing). But ImageJ/Fiji or Icy can’t pick up the different channels. only the first 11z’s

for ch in channels:
    for i in range(meta['components']):
        if (meta['plane_'+str(i)]['name']) == ch:
            print("Channel{} i:",ch,i)
            if ch == 'DAPI': 
                frames.default_coords['c'] = i
                psf = psfdapi

                res =algo.run(
                fd_data.Acquisition(
                data=frames[0],
                kernel=psf)
                , niter=10)
                #print(res.data.shape)


                for j in range(11):
                    #print(res.data[j,:,:].shape)
                    bioformats.write_image(
                        '/tmp/psf.tif', res.data[j,:,:].transpose().astype(np.int16)
                        ,bioformats.PT_UINT16
                        ,c=0
                        ,z=j
                        ,size_c=2
                        ,size_z=12
                        #,channel_names='DAPI'
                    )

            else: 
                frames.default_coords['c'] = i
                psf = psfruby
                res =algo.run(
                fd_data.Acquisition(
                data=frames[0],
                kernel=psf)
                , niter=10)
                #print(res.data.shape)
                for j in range(11):
                    bioformats.write_image(
                        '/tmp/psf.tif', res.data[j,:,:].transpose().astype(np.int16)
                        ,bioformats.PT_UINT16
                        ,c=1
                        ,z=j
                        ,size_c=2
                        ,size_z=12
                        #,channel_names='mRuby3'
                    )

I also tried with Tifffile (append=True). It does append all channels (but they are just recognized as extra Zs). I tried to append the metadata in different ways (even using TifffileWriter and appending the metadata at such moment but it always says Json cannot be concatenated (I don’t recall the exact message)

import tifffile as tf

for ch in channels:
    for i in range(meta['components']):
        if (meta['plane_'+str(i)]['name']) == ch:
            print("Channel{} i:",ch,i)
            if ch == 'DAPI': 
                frames.default_coords['c'] = i
                psf = psfdapi
                
                res =algo.run(
                fd_data.Acquisition(
                data=frames[0],
                kernel=psf)
                , niter=10)
                
                tf.imwrite('/tmp/test.tiff', res.data 
                #,resolution= 0.108333333333333
                #,metadata=meta
                ,append=True
                , bigtiff=True
                )

            else: 
                frames.default_coords['c'] = i
                psf = psfruby
                res =algo.run(
                fd_data.Acquisition(
                data=frames[0],
                kernel=psf)
                , niter=10)
                
                tf.imwrite('/tmp/test.tiff', res.data
                #,resolution= 0.108333333333333
                #,metadata=meta
                ,append=True
                , bigtiff=True
                )

Thanks for any help.

1 Like

Hi @jmamede,

This might be a distraction/side note, as I am not very familiar with tiff writing or metadata writing. My question, however, is whether you have considered writing to zarr/n5 instead of tiff? I think it will be a more suitable format for almost any downstream visualisation/analysis you may want, and there are readers in both Fiji and Python, see e.g. https://github.com/ome/ome-zarr-py/

The caveat is that the metadata specification for zarr is still evolving (search “metadata zarr” in the forum for discussions). But the core format is stable and metadata should be easy enough to migrate. (And of course if you keep your raw data around somewhere you can always run the pipeline again with new versions!)

Hope this helps!

I will try. I also thought about it.
I think I can use navapi with zarr.

However I also wanted to import the tiffs to NIS-elements.

I managed to have something usefull with tifffile and set ImageJ=True.
Metadata will be the next problem to solve.

@jmamede

I am not sure if it directly helps, but I use this code when writing numpy array back to OME-TIFF (If i really hav too). It relies on the AICSImageIO python libary. And it allows to save also the basic metadata. You certainly have to adapt it to your needs and it might have some bugs still … :slight_smile:

def write_ometiff_aicsimageio(savepath, imgarray, metadata,
                              reader='aicsimageio',
                              overwrite=False):
    """Write an OME-TIFF file from an image array based on the metadata.

    :param filepath: savepath of the OME-TIFF stack
    :type filepath: str
    :param imgarray: multi-dimensional image array
    :type imgarray: NumPy.Array
    :param metadata: metadata dictionary with the required information
    to create an correct OME-TIFF file
    :type metadata: dict
    :param reader: string (aicsimagio or czifile) specifying
    the used reader, defaults to aicsimageio
    :type metadata: str
    :param overwrite: option to overwrite an existing OME-TIFF, defaults to False
    :type overwrite: bool, optional
    """

    # define scaling from metadata or use defualt scaling
    try:
        pixels_physical_size = [metadata['XScale'],
                                metadata['YScale'],
                                metadata['ZScale']]
    except KeyError as e:
        print('Key not found:', e)
        print('Use default scaling XYZ=1')
        pixels_physical_size = [1, 1, 1]

    # define channel names list from metadata
    try:
        channel_names = []
        for ch in metadata['Channels']:
            channel_names.append(ch)
    except KeyError as e:
        print('Key not found:', e)
        channel_names = None

    # get the dimensions and their position inside the dimension string
    if reader == 'aicsimageio':

        dims_dict, dimindex_list, numvalid_dims = get_dimorder(metadata['Axes_aics'])

        # if the array has more than 5 dimensions then remove the S dimension
        # because it is not supported by OME-TIFF
        if len(imgarray.shape) > 5:
            try:
                imgarray = np.squeeze(imgarray, axis=dims_dict['S'])
            except Exception:
                print('Could not remover S Dimension from string.)')

        # remove the S character from the dimension string
        new_dimorder = metadata['Axes_aics'].replace('S', '')

    if reader == 'czifile':

        new_dimorder = metadata['Axes']
        dims_dict, dimindex_list, numvalid_dims = get_dimorder(metadata['Axes'])
        """
        '0': 'Sample',  # e.g. RGBA
        'X': 'Width',
        'Y': 'Height',
        'C': 'Channel',
        'Z': 'Slice',  # depth
        'T': 'Time',
        'R': 'Rotation',
        'S': 'Scene',  # contiguous regions of interest in a mosaic image
        'I': 'Illumination',  # direction
        'B': 'Block',  # acquisition
        'M': 'Mosaic',  # index of tile for compositing a scene
        'H': 'Phase',  # e.g. Airy detector fibers
        'V': 'View',  # e.g. for SPIM
        """

        to_remove = []

        # list of unspupported dims for writing an OME-TIFF
        dims = ['R', 'I', 'M', 'H', 'V', 'B', 'S', '0']

        for dim in dims:
            if dims_dict[dim] >= 0:
                # remove the CZI DIMENSION character from the dimension string
                new_dimorder = new_dimorder.replace(dim, '')
                # add dimension index to the list of axis to be removed
                to_remove.append(dims_dict[dim])
                print('Remove Dimension:', dim)

        # create tuple with dimensions to be removed
        dims2remove = tuple(to_remove)
        # remove dimensions from array
        imgarray = np.squeeze(imgarray, axis=dims2remove)

    # write the array as an OME-TIFF incl. the metadata
    try:
        with ome_tiff_writer.OmeTiffWriter(savepath, overwrite_file=overwrite) as writer:
            writer.save(imgarray,
                        channel_names=channel_names,
                        image_name=os.path.basename((savepath)),
                        pixels_physical_size=pixels_physical_size,
                        channel_colors=None,
                        dimension_order=new_dimorder)
            writer.close()
    except Exception as error:
        print(error.__class__.__name__ + ": " + error.msg)
        print('Could not write OME-TIFF')
        savepath = None

    return savepath
1 Like

Hi,

this is also not exactly a OME solution, but at least a working piece of code that saves sequentially a 4D Numpy array which is then recognized by Fiji as having time, channel and Z dimensions. I only used it for TCXY data (a long time ago), not TCXYZ data so this is provided as an example you might have to fix. Also it uses skimage.external.tifffile which I think is going to be dropped (or has been already). So you might want to adapt it to the original tifffile package. If I remembre correctly, key things were: 1) to set the correct metadata, and 2) to use append=force when saving each frame:

filename = 'myimage.tif'

#define metadata for imagej
maxframe = 50
slices = 4
colors = ['GFP','YFP','RFP']
metadata = {'channels':len(colors),'slices':slices,
            'frames':maxframe,'hyperstack':True,'loop':False}

for t in range(maxframe):
    
    ## code to process images. The final output is a 3D numpy array with dimensions XYC
    ## called image_stack. Here it's just a random image with increasing intensity over time
    ## as a control
    image_stack = np.random.randint(0, int(t)+1, (slices,20,20,3))
    
    for j in range(slices):
        for i in range(len(colors)):
            imtosave = image_stack[j,:,:,i]

            # I think you need to flip to have correct orientation in imageJ
            imtosave_flip = np.flipud(imtosave.T)

            skimage.external.tifffile.imsave(
                filename,
                imtosave_flip.astype(np.uint16),
                append = 'force',
                imagej = True, 
                metadata = metadata)

Good luck!

2 Likes

thank you Both of you.

I managed to implement ome.tiff solution.
I still have one little warning message that I want to solve though.

This example works (however it’s not standard ome.tiff and imageJ tries to parse the file before hand to fix something).
All other metadata is conserved though.

import numpy as np
import bioformats.omexml as ome
import tifffile as tf
import sys


def writeplanes(pixel, SizeT=1, SizeZ=1, SizeC=1, order='TZCYX'
            , verbose=False):

    if order == 'TZCYX':

        p.DimensionOrder = ome.DO_XYCZT
        counter = 0
        for t in range(SizeT):
            for z in range(SizeZ):
                for c in range(SizeC):

                    if verbose:
                        print('Write PlaneTable: ', t, z, c),
                        sys.stdout.flush()

                    pixel.Plane(counter).TheT = t
                    pixel.Plane(counter).TheZ = z
                    pixel.Plane(counter).TheC = c
                    counter = counter + 1

    return pixel


# Dimension TZCXY
SizeT = 1
SizeZ = 11
SizeC = 2
SizeX = 2044
SizeY = 2044
Series = 0


scalex = 0.10833
scaley = scalex
scalez = 0.3
pixeltype = 'uint16'
dimorder = 'TZCYX'
output_file = r'/tmp/stack.ome.tiff' #this does nothing in this example

# create numpy array with correct order

# Getting metadata info
omexml = ome.OMEXML()
omexml.image(Series).Name = output_file
p = omexml.image(Series).Pixels
#p.ID = 0
p.SizeX = SizeX
p.SizeY = SizeY
p.SizeC = SizeC
p.SizeT = SizeT
p.SizeZ = SizeZ
p.PhysicalSizeX = np.float(scalex)
p.PhysicalSizeY = np.float(scaley)
p.PhysicalSizeZ = np.float(scalez)
p.PixelType = pixeltype
p.channel_count = SizeC
p.plane_count = SizeZ * SizeT * SizeC
p = writeplanes(p, SizeT=SizeT, SizeZ=SizeZ, SizeC=SizeC, order=dimorder)

for c in range(SizeC):
    if pixeltype == 'unit8':
        p.Channel(c).SamplesPerPixel = 1
    if pixeltype == 'unit16':
        p.Channel(c).SamplesPerPixel = 2


omexml.structured_annotations.add_original_metadata(
            ome.OM_SAMPLES_PER_PIXEL, str(SizeC))

# Converting to omexml
xml = omexml.to_xml()


img5d = np.random.randn(
        SizeT, SizeZ, SizeC, SizeY, SizeX).astype(np.uint16)

# ~ write file and save OME-XML as description
tf.imwrite(r'/tmp/test1.ome.tiff', img5d#,
    , description=xml)


with tf.TiffWriter('/tmp/test2.ome.tiff'
               #, bigtiff=True
               #, imagej=True
                      ) as tif:
    for t in range(SizeT):
        for z in range(SizeZ):
            for c in range(SizeC):      
                # ~ print(img5d[t,z,c,:,:].shape)   # -> (2044, 2044)
                tif.save(img5d[t,z,c,:,:]
            #                     ,shape=res.shape

                        #,resolution= (.1083,0.1083,3)
                         , description = xml
                        , photometric='minisblack'
                        #, datetime= True
                        , metadata={'axes': 'TZCYX'
                            , 'DimensionOrder' : 'TZCYX'
                            , 'Resolution': 0.10833}
                            )

with tf.TiffWriter('/tmp/test3.ome.tiff'
               #, bigtiff=True
               #, imagej=True
                      ) as tif:
    for t in range(SizeT):
        # ~ print(img5d[t,z,c,:,:].shape)   # -> (2044, 2044)
        tif.save(img5d[t,:,:,:,:]
    #                     ,shape=res.shape

                #,resolution= (.1083,0.1083,3)
                 , description = xml
                , photometric='minisblack'
                #, datetime= True
                , metadata={'axes': 'TZCYX'
                    , 'DimensionOrder' : 'TZCYX'
                    , 'Resolution': 0.10833}
                    )

with tf.TiffWriter('/tmp/test4.ome.tiff'
               #, bigtiff=True
               , imagej=True
                      ) as tif:
    for t in range(SizeT):
        # ~ print(img5d[t,z,c,:,:].shape)   # -> (2044, 2044)
        tif.save(img5d[t,:,:,:,:]
    #                     ,shape=res.shape

                #,resolution= (.1083,0.1083,3)
                 , description = xml
                , photometric='minisblack'
                #, datetime= True
                , metadata={'axes': 'TZCYX'
                    , 'DimensionOrder' : 'TZCYX'
                    , 'Resolution': 0.10833}
                    )



I had a XML dump from pims that gave me a string from the original metadata when I open with bioformats.
That eventually failed and I created my library to make a correct OME xml.

Are you able to read your metadata and write the same one to the exit file?

I had tried with append = ‘force’ as well. I settled with using the writer function to write frame by frame.
I read somewhere that the append = ‘force’ parses the whole file each time it needs to write and that would make it slower.