Writing multi-channel pyramids

[petroslk]

Hey there,

I had a question regarding the writing of pyramids with more than 3 channels. In my case, these are multiplex immunofluorescence images with ~20 channels per magnification layer.

When reading them I can see that all the pages are correctly read:

image = pyvips.Image.new_from_file("/home/storage/co_registration_IF_HE/slides_to_reg/test_slide.er.qptiff")
image.get_n_pages()
133

But when I write this to a tiff file, it only seems to write one channel per pyramid layer:

image.tiffsave(f'test_multi_channel.tif', tile=True, pyramid=True)
image2 = pyvips.Image.new_from_file("test_multi_channel.tif")
image2.get_n_pages()
9

Can I actually write such images directly?

Thanks a lot!

Petros

[jcupitt]

Hello @petroslk,

TIFF doesn't support multichannel images that well -- for example:

$ vips black x.tif 100 100 --bands 100
$ tiffinfo x.tif 
=== TIFF directory 0 ===
TIFF Directory at offset 0xf4248 (1000008)
  Image Width: 100 Image Length: 100
  Resolution: 25.4, 25.4 pixels/inch
  Bits/Sample: 8
  Sample Format: unsigned integer
  Compression Scheme: None
  Photometric Interpretation: RGB color
  Extra Samples: 97<unassoc-alpha, unassoc-alpha, ...

So a 100 channel image is saved as RGB, with 97 extra unassociated alpha channels (!!). This works, but not many systems will be able to read files like this.

The usual trick with TIFF is to save these as multi-page images instead, so this would be a 100 page image, with each page holding one channel. OME-TIFF adds some extra XML metadata to the IMAGEDESCRIPTION tag to explain what these pages mean.

However, this means you can't easily save pyramids, since TIFF usually uses the "pages" dimension to keep pyramid levels. OME-TIFF uses an extra, hidden, sub-IFD dimension hanging off each page to hold the pyramid. The image looks like:

page0 (channel or band 0)
   |
   +-> page0, subifd 0 (band 0 at half size)
         |
         +-> page0, subifd0, subifd0 (band 0 at quarter size)
page1 (band 1)
   |
   ... etc.

libvips can read and write this layout, but you need to make the many-page image yourself, and write the XML yourself. Something like:

#!/usr/bin/python3

import sys
import pyvips

im = pyvips.Image.new_from_file(sys.argv[1])

image_height = im.height

# split to separate image planes and stack vertically ready for OME 
im = pyvips.Image.arrayjoin(im.bandsplit(), across=1)

# set minimal OME metadata
# before we can modify an image (set metadata in this case), we must take a 
# private copy
im = im.copy()
im.set_type(pyvips.GValue.gint_type, "page-height", image_height)
im.set_type(pyvips.GValue.gstr_type, "image-description",
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="{im.bands}"
                SizeT="1"
                SizeX="{im.width}"
                SizeY="{image_height}"
                SizeZ="1"
                Type="uint8">
        </Pixels>
    </Image>
</OME>""")

im.tiffsave(sys.argv[2], compression="lzw", tile=True,
            tile_width=512, tile_height=512,
            pyramid=True, subifd=True)

I'd use something like QuPath for testing. Perhaps you are using this already.

[petroslk]

Hey @jcupitt

Thanks a lot for your feedback! Makes a lot more sense now. I did try your short snippet but still it only wrote one of the channels.

Tried something like this aswell but it only writes 2 channels?

#!/usr/bin/python3

import sys
import pyvips

# Load the input image
im = pyvips.Image.new_from_file(sys.argv[1])

# Check the number of channels
num_channels = im.bands

# Create an empty list to hold the separate channel images
channel_images = []

# Split the image into separate channels and add each channel to the list
for i in range(num_channels):
    # Extract the i-th channel
    channel = im.bandjoin(im[i])
    channel_images.append(channel)

# Stack the separate channels vertically to prepare for OME-TIFF
combined_image = pyvips.Image.arrayjoin(channel_images, across=1)

# Set minimal OME metadata
combined_image = combined_image.copy()
combined_image.set_type(pyvips.GValue.gint_type, "page-height", im.height)
combined_image.set_type(pyvips.GValue.gstr_type, "image-description",
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="{num_channels}"
                SizeT="1"
                SizeX="{im.width}"
                SizeY="{im.height}"
                SizeZ="1"
                Type="uint8">
        </Pixels>
    </Image>
</OME>""")

# Save the image as an OME-TIFF file
combined_image.tiffsave(sys.argv[2], compression="lzw", tile=True,
                        tile_width=512, tile_height=512,
                        pyramid=True, subifd=True)

Would it be easier to write each channel as a separate tif and then merge them?

Thanks a lot!

Petros

[jcupitt]

Ooop! You're right, it should get num_channels before splitting the image, sorry about that. I tried this version:

#!/usr/bin/env python3

import sys
import pyvips

im = pyvips.Image.new_from_file(sys.argv[1])
num_channels = im.bands

image_height = im.height

# split to separate image planes and stack vertically ready for OME 
im = pyvips.Image.arrayjoin(im.bandsplit(), across=1)

# set minimal OME metadata
# before we can modify an image (set metadata in this case), we must take a 
# private copy
im = im.copy()
im.set_type(pyvips.GValue.gint_type, "page-height", image_height)
im.set_type(pyvips.GValue.gstr_type, "image-description",
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="{num_channels}"
                SizeT="1"
                SizeX="{im.width}"
                SizeY="{image_height}"
                SizeZ="1"
                Type="uint8">
        </Pixels>
    </Image>
</OME>""")

im.tiffsave(sys.argv[2], compression="lzw", tile=True,
            tile_width=512, tile_height=512,
            pyramid=True, subifd=True)

I made a six-band test image:

$ vips bandjoin "k2.jpg k2.jpg" x.v
$ vipsheader x.v
x.v: 1450x2048 uchar, 6 bands, srgb, jpegload

Then ran that code and examined the TIFF structure:

$ ./vips2ome.py x.v x.tif
$ tiffinfo x.tif
=== TIFF directory 0 ===
TIFF Directory at offset 0x184b58 (1592152)
  Subfile Type: multi-page document (2 = 0x2)
  Image Width: 1450 Image Length: 2048
  Tile Width: 512 Tile Length: 512
  Resolution: 72.009, 72.009 pixels/inch
  Bits/Sample: 8
  Sample Format: unsigned integer
  Compression Scheme: LZW
  Photometric Interpretation: min-is-black
  Orientation: row 0 top, col 0 lhs
  Samples/Pixel: 1
  Planar Configuration: single image plane
  Page Number: 0-6
  SubIFD Offsets: 2278676 2461530
  ImageDescription: <?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="6"
                SizeT="1"
                SizeX="1450"
                SizeY="2048"
                SizeZ="1"
                Type="uint8">
        </Pixels>
    </Image>
</OME>
  Predictor: horizontal differencing 2 (0x2)
--- SubIFD image descriptor tag within TIFF directory 0 with array of 2 SubIFD chains ---
--- SubIFD 0 of chain 0 at offset 0x22c514 (2278676):
TIFF Directory at offset 0x22c514 (2278676)
  Subfile Type: reduced-resolution image (1 = 0x1)
... snip ...
=== TIFF directory 1 ===
TIFF Directory at offset 0x3da3ec (4039660)
  Subfile Type: multi-page document (2 = 0x2)
  Image Width: 1450 Image Length: 2048
  Tile Width: 512 Tile Length: 512
  Resolution: 72.009, 72.009 pixels/inch
  Bits/Sample: 8
  Sample Format: unsigned integer
  Compression Scheme: LZW
  Photometric Interpretation: min-is-black
  Orientation: row 0 top, col 0 lhs
  Samples/Pixel: 1
  Planar Configuration: single image plane
  Page Number: 1-6
  SubIFD Offsets: 4725898 4907982
... snip ...

=== TIFF directory 5 ===
TIFF Directory at offset 0xd75e14 (14114324)
  Subfile Type: multi-page document (2 = 0x2)
  Image Width: 1450 Image Length: 2048
  Tile Width: 512 Tile Length: 512
... snip ...

You can see it's written a six page TIFF from the six-band image, and each page has a pyramid on it.

It loads into QuPath as a 6-channel image:

You'd need to adjust the XML if you want QuPath to label your fluorescence images correctly, of course.

[jcupitt]

... I meant to say, vipsdisp is useful for testing as well. With that TIFF I see:

You can see the six bands in the info bar at the bottom. That's the default "Pages as bands", but if you select "Toilet roll" in the menu on the bottom left you see:

So you can examine every band separately.

https://github.com/jcupitt/vipsdisp

[petroslk]

Dear John,

To experiment a bit more I experted the qptiff image as an ome.tif via qupath. When opening the ome.tif file on vipsdisp, it looks good.

Following conversion with the script, the resulting image looks like this:

I then checked the header of the ome.tif and, it looks different from the image you generated:

doron_test_lzw.ome.tif: 26880x30240 ushort, 1 band, grey16, tiffload

Could this be the reason?

[jcupitt]

Did you update your vips2ome.py script? The first version I posted had a bug that could cause this.

[petroslk]

I think I did yes:

#!/usr/bin/env python3

import sys
import pyvips

im = pyvips.Image.new_from_file(sys.argv[1])
num_channels = im.bands

image_height = im.height

# split to separate image planes and stack vertically ready for OME 
im = pyvips.Image.arrayjoin(im.bandsplit(), across=1)

# set minimal OME metadata
# before we can modify an image (set metadata in this case), we must take a 
# private copy
im = im.copy()
im.set_type(pyvips.GValue.gint_type, "page-height", image_height)
im.set_type(pyvips.GValue.gstr_type, "image-description",
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="{num_channels}"
                SizeT="1"
                SizeX="{im.width}"
                SizeY="{image_height}"
                SizeZ="1"
                Type="uint8">
        </Pixels>
    </Image>
</OME>""")

im.tiffsave(sys.argv[2], compression="lzw", tile=True,
            tile_width=512, tile_height=512,
            pyramid=True, subifd=True)

Then tried to convert the ome.tif to essentially another ome.tif just to test:

python script.py doron_test_lzw.ome.tif test2.tif

The vipsheader result for the test2.tif file looks like this, identical to the original ome.tif, but only has a single channel when oppened on vipsdisp or QuPath:

test2.tif: 26880x30240 ushort, 1 band, grey16, tiffload

Once again, thank you for your help! If you're not sure what the source of the problem might be, I will try some other solutions, or possibly transfer the test file! Thank you for your patience.

[jcupitt]

Wait, is your source image already an OME TIFF? But then why are you trying to convert it? Are you loading it, doing some processing, and trying to save again?

Just as you have to do that stuff with bandsplit to save as an OME TIFF, you have to do the opposite conversion when you load an OME TIFF (turn the Y dimension into colour). Or if you are copying one OME TIFF to another, you can skip the bandsplit / pagesplit and just copy the pixels over.

Yes, I think a sample image would make everything clearer.

[petroslk]

Yes, pretty much, at this point the only thing I need to do is take the one tif per IF channel that I have performed the modifications on and merge all of those TIFs together into an OME.TIF. How would you go about merging a series of tif files into a multi channel ome.tif file with pyvips? Thanks a lot!

[jcupitt]

Could you post a sample image? It'd help me understand exactly what you're working with.

[petroslk]

I have just sent you an email with the example file and the use case! Hope that helps! Thanks a lot!!

[jcupitt]

Wow that's a pretty crazy pyramid, I've not seen one like that before. I wish manufactures would follow conventions :(

$ for i in {0..44}; do echo -n "page $i: "; vipsheader Tonsil_Scan1.er.qptiff[page=$i]; done
page 0: Tonsil_Scan1.er.qptiff: 26880x36000 ushort, 1 band, grey16, tiffload
page 1: Tonsil_Scan1.er.qptiff: 26880x36000 ushort, 1 band, grey16, tiffload
page 2: Tonsil_Scan1.er.qptiff: 26880x36000 ushort, 1 band, grey16, tiffload
page 3: Tonsil_Scan1.er.qptiff: 26880x36000 ushort, 1 band, grey16, tiffload
page 4: Tonsil_Scan1.er.qptiff: 26880x36000 ushort, 1 band, grey16, tiffload
page 5: Tonsil_Scan1.er.qptiff: 26880x36000 ushort, 1 band, grey16, tiffload
page 6: Tonsil_Scan1.er.qptiff: 26880x36000 ushort, 1 band, grey16, tiffload
page 7: Tonsil_Scan1.er.qptiff: 210x281 uchar, 3 bands, srgb, tiffload
page 8: Tonsil_Scan1.er.qptiff: 13440x18000 ushort, 1 band, grey16, tiffload
page 9: Tonsil_Scan1.er.qptiff: 13440x18000 ushort, 1 band, grey16, tiffload
page 10: Tonsil_Scan1.er.qptiff: 13440x18000 ushort, 1 band, grey16, tiffload
page 11: Tonsil_Scan1.er.qptiff: 13440x18000 ushort, 1 band, grey16, tiffload
page 12: Tonsil_Scan1.er.qptiff: 13440x18000 ushort, 1 band, grey16, tiffload
page 13: Tonsil_Scan1.er.qptiff: 13440x18000 ushort, 1 band, grey16, tiffload
page 14: Tonsil_Scan1.er.qptiff: 13440x18000 ushort, 1 band, grey16, tiffload
page 15: Tonsil_Scan1.er.qptiff: 6720x9000 ushort, 1 band, grey16, tiffload
page 16: Tonsil_Scan1.er.qptiff: 6720x9000 ushort, 1 band, grey16, tiffload
page 17: Tonsil_Scan1.er.qptiff: 6720x9000 ushort, 1 band, grey16, tiffload
page 18: Tonsil_Scan1.er.qptiff: 6720x9000 ushort, 1 band, grey16, tiffload
page 19: Tonsil_Scan1.er.qptiff: 6720x9000 ushort, 1 band, grey16, tiffload
page 20: Tonsil_Scan1.er.qptiff: 6720x9000 ushort, 1 band, grey16, tiffload
page 21: Tonsil_Scan1.er.qptiff: 6720x9000 ushort, 1 band, grey16, tiffload
page 22: Tonsil_Scan1.er.qptiff: 3360x4500 ushort, 1 band, grey16, tiffload
page 23: Tonsil_Scan1.er.qptiff: 3360x4500 ushort, 1 band, grey16, tiffload
page 24: Tonsil_Scan1.er.qptiff: 3360x4500 ushort, 1 band, grey16, tiffload
page 25: Tonsil_Scan1.er.qptiff: 3360x4500 ushort, 1 band, grey16, tiffload
page 26: Tonsil_Scan1.er.qptiff: 3360x4500 ushort, 1 band, grey16, tiffload
page 27: Tonsil_Scan1.er.qptiff: 3360x4500 ushort, 1 band, grey16, tiffload
page 28: Tonsil_Scan1.er.qptiff: 3360x4500 ushort, 1 band, grey16, tiffload
page 29: Tonsil_Scan1.er.qptiff: 1680x2250 ushort, 1 band, grey16, tiffload
page 30: Tonsil_Scan1.er.qptiff: 1680x2250 ushort, 1 band, grey16, tiffload
page 31: Tonsil_Scan1.er.qptiff: 1680x2250 ushort, 1 band, grey16, tiffload
page 32: Tonsil_Scan1.er.qptiff: 1680x2250 ushort, 1 band, grey16, tiffload
page 33: Tonsil_Scan1.er.qptiff: 1680x2250 ushort, 1 band, grey16, tiffload
page 34: Tonsil_Scan1.er.qptiff: 1680x2250 ushort, 1 band, grey16, tiffload
page 35: Tonsil_Scan1.er.qptiff: 1680x2250 ushort, 1 band, grey16, tiffload
page 36: Tonsil_Scan1.er.qptiff: 840x1125 ushort, 1 band, grey16, tiffload
page 37: Tonsil_Scan1.er.qptiff: 840x1125 ushort, 1 band, grey16, tiffload
page 38: Tonsil_Scan1.er.qptiff: 840x1125 ushort, 1 band, grey16, tiffload
page 39: Tonsil_Scan1.er.qptiff: 840x1125 ushort, 1 band, grey16, tiffload
page 40: Tonsil_Scan1.er.qptiff: 840x1125 ushort, 1 band, grey16, tiffload
page 41: Tonsil_Scan1.er.qptiff: 840x1125 ushort, 1 band, grey16, tiffload
page 42: Tonsil_Scan1.er.qptiff: 840x1125 ushort, 1 band, grey16, tiffload
page 43: Tonsil_Scan1.er.qptiff: 1536x2736 uchar, 3 bands, srgb, tiffload
page 44: Tonsil_Scan1.er.qptiff: 960x720 uchar, 3 bands, srgb, tiffload

Assuming OME-TIFF is OK, you need to extract the first 7 pages (looks like those are the full-res scans), process them, then save. Something like:

#!/usr/bin/env python3

import sys
import pyvips

filename = "Tonsil_Scan1.er.qptiff"

# you should find this number from examining the XML stored in the
# IMAGEDESCRIPTION
n_pages = 7

# extract full res bands
bands = [pyvips.Image.new_from_file(filename, page=i)
         for i in range(n_pages)]

# process each band in some way ... here I just flip the greyscale
bands = [(65535 - band).cast("ushort") for band in bands] 

# reassemble into a tall, thin strip
image = pyvips.Image.arrayjoin(bands, across=1)

# set the OME-TIFF metainfo
image = image.copy()
image.set_type(pyvips.GValue.gint_type, "page-height", bands[0].height)
image.set_type(pyvips.GValue.gstr_type, "image-description",
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="{n_pages}"
                SizeT="1"
                SizeX="{bands[0].width}"
                SizeY="{bands[0].height}"
                SizeZ="1"
                Type="uint16">
        </Pixels>
    </Image>
</OME>""")

print(f"saving to x.tif ...")
image.tiffsave("x.tif", compression="lzw", tile=True,
               tile_width=512, tile_height=512,
               pyramid=True, subifd=True, bigtiff=True)

It seems to load into QuPath OK:

[petroslk]

Dear @jcupitt

This works perfectly! Thank you so much for your help! Indeed this file was quite tricky to work with, but at least we figured it out now!

Cheers,

Petros