Performance comparison between tifffile and iohub's custom OME-TIFF implementation

[ziw-liu]

A custom OME-TIFF reader (iohub.multipagetiff.MicromanagerOmeTiffReader) was implemented because historically tifffile and AICSImageIO was slow when reading large OME-TIFF series generated by Micro-Manager acquisitions.

While debugging #65 I found that this implementation does not guarantee data integrity during reading. Before investing more time in fixing it, I think it is worth revisiting the topic of whether it is worth maintaining a custom OME-TIFF reader, given that the more widely adopted solutions have evolved since waveorder.io's designation.

Here is a simple read speed benchmark of tifffile and iohub's custom reader:

The test was done on a 123 GB dataset with TCZYX=(8, 9, 3, 81, 2048, 2048) dimensions. Voxels from 2 non-sequential positions was read into RAM in each iteration (N=5).

Test script (click to expand):

Environment: Python 3.10.8, Linux 4.18 (x86_64, AMD EPYC 7H12@2.6GHz)

# %%
import os
from timeit import timeit
import zarr
import pandas as pd

# readers tested
from tifffile import TiffSequence  # 2023.2.3
from iohub.multipagetiff import MicromanagerOmeTiffReader  # 0.1.dev368+g3d62e6f


# %%
# 123 GB total
DATASET = (
  "/hpc/projects/comp_micro/rawdata/hummingbird/Soorya/"
  "2022_06_27_A549cellMembraneStained/"
  "A549_CellMaskDye_Well1_deltz0.25_63X_30s_2framemin/"
  "A549_CellMaskdye_Well1_30s_2framemin_1"
)

POSITIONS = (2, 0)


# %%
def read_tifffile():
  sequence = TiffSequence(os.scandir(DATASET))
  data = zarr.open(sequence.aszarr(), mode="r")
  for p in POSITIONS:
      _ = data[p]
  sequence.close()


# %%
def read_custom():
  reader = MicromanagerOmeTiffReader(DATASET)
  for p in POSITIONS:
      _ = reader.get_array(p)


# %%
def repeat(n=5):
  tf_times = []
  wo_times = []
  for _ in range(n):
      tf_times.append(
          timeit(
              "read_tifffile()", number=1, setup="from __main__ import read_tifffile"
          )
      )
      wo_times.append(
          timeit("read_custom()", number=1, setup="from __main__ import read_custom")
      )
  return pd.DataFrame({"tifffile": tf_times, "waveorder": wo_times})


# %%
timings = repeat()

At least in this test, the latest tifffile consistently out-performs the iohub implementation. While a comprehensive benchmark will take more time (#57), I think as long as a widely used library is not significantly slower, the reduction of maintenance overhead and increased user testing can make a strong case for us to reconsider maintaining the custom code in iohub.

[ziw-liu]

A note on interpreting the results (why no error bars) from timeit docs:

Note It’s tempting to calculate mean and standard deviation from the result vector and report these. However, this is not very useful. In a typical case, the lowest value gives a lower bound for how fast your machine can run the given code snippet; higher values in the result vector are typically not caused by variability in Python’s speed, but by other processes interfering with your timing accuracy. So the min() of the result is probably the only number you should be interested in. After that, you should look at the entire vector and apply common sense rather than statistics.

In our case if the file is small FS caching might decrease that lower bound, but the general idea holds true.

[ziw-liu]

Repeated the test with a much larger (1036 GB) dataset, and this time iohub's custom reader is faster by a large margin:

Test script (click to expand):

Environment: Python 3.10.8, Linux 4.18 (x86_64, AMD EPYC 7H12@2.6GHz)

# %%
import os
from timeit import timeit
from tqdm import tqdm
import zarr
import pandas as pd

# readers tested
from tifffile import TiffSequence  # 2023.2.3
from iohub.multipagetiff import MicromanagerOmeTiffReader  # 0.1.dev368+g3d62e6f


# %%
# 1036 GB total
DATASET = (
    "/hpc/projects/comp_micro/rawdata/hummingbird/Janie/"
    "2022_03_15_orgs_nuc_mem_63x_04NA/all_21_3"
)
POSITIONS = (50, 100, 150, 200, 250)


# %%
def read_tifffile():
    sequence = TiffSequence(os.scandir(DATASET))
    data = zarr.open(sequence.aszarr(), mode="r")
    for p in POSITIONS:
        _ = data[p]
    sequence.close()


# %%
def read_custom():
    reader = MicromanagerOmeTiffReader(DATASET)
    for p in POSITIONS:
        _ = reader.get_array(p)


# %%
def repeat(n=5):
    tf_times = []
    wo_times = []
    for _ in tqdm(range(n)):
        tf_times.append(
            timeit(
                "read_tifffile()", number=1, setup="from __main__ import read_tifffile"
            )
        )
        wo_times.append(
            timeit("read_custom()", number=1, setup="from __main__ import read_custom")
        )
    return pd.DataFrame({"tifffile": tf_times, "waveorder": wo_times})


# %%
def main():
    timings = repeat()
    print(timings)
    timings.to_csv("large_tiff_time.csv")


# %%
if __name__ == "__main__":
    main()

[ziw-liu]

Given these tests it seems that the inflection point of tifffile/iohub performance is in the vicinity of 500 GB dataset size. Since our users do expect ~TB OME-TIFF datasets frequently, it makes sense to keep investing in a custom solution.

One thing of concern, however, is that 7a6675a makes the custom reader a lot slower too. Edit: fixed.

[cgohlke]

Allow me to chime in. 350 vs 50 seconds to read 25 files seems excessive. What is the number of files, size of files, and number of pages in each file? If this is a single, multi-file MicroManager dataset, the use of TiffSequence is counterproductive. Tifffile uses OME-XML metadata in the first file, which describes the whole multi-file dataset and triggers indexing all TIFF pages in all files of the dataset...

[ziw-liu]

@cgohlke Thanks for letting us know of the improvement we can make in benchmarking! We are very much looking forward to switching to a more widely used and tested implementation with comparable performance.

What is the number of files, size of files, and number of pages in each file? If this is a single, multi-file MicroManager dataset, the use of TiffSequence is counterproductive.

There's 260 .ome.tif files, each 4 GB (except for those written at the end of the acquisition). It is a single MM acquisition of shape PTCZYX=(336, 1, 4, 97, 2048, 2048) and type uint16.

Can you point us to the recommended entry point for such a dataset with tifffile? I tried using tifs = imread(FIRST_FILE, aszarr=True) to get the Zarr store but it only seems to take even longer time.

[cgohlke]

Can you point us to the recommended entry point for such a dataset with tifffile? I tried using tifs = imread(FIRST_FILE, aszarr=True) to get the Zarr store but it only seems to take even longer time.

Unfortunately it is currently not possible to get a single Zarr store of a multi-position MicroManager dataset. Tifffile uses the OME-XML metadata, which describes positions as distinct OME series. imread(FIRST_FILE, ...) returns the zarr store for the first position only. Instead, use the lower level interface, e.g.:

with TiffFile(FIRST_FILE) as tif:
    series = tif.series  # this parses OME-XML and indexes all pages in all files
    for position in POSITIONS:
        # read series as numpy array ...
        im = series[position].asarray()
        # ... or via zarr store
        store = series[position].aszarr()
        z = zarr.open(store, mode='r')
        im = z[:]
        store.close()

A major bottleneck when only accessing a few positions is that tifffile needs to index all TIFF pages/IFDs in the dataset when creating the series. That requires a lot of seeks and little reads over many files. I have a patch that uses the MicroManager IndexMap, which speeds up indexing about 10x for a smaller 8 GB dataset.

It would require changes in tifffile to return a single series for the whole MicroManager datastore. Either special casing the OME-XML parser, or adding a dedicated MicroManager parser similar to NDTiff. A special parser would be preferable, but (1) the format seems poorly documented compared to OME-TIFF, (2) there are variations between MicroManager versions, (3) there is no public set of test files, and (4) MicroManager frequently produces corrupted files (at least on my system).

[ziw-liu]

Tifffile uses the OME-XML metadata, which describes positions as distinct OME series. imread(FIRST_FILE, ...) returns the zarr store for the first position only. Instead, use the lower level interface, e.g.:

Thanks for showing this! This way is much faster in most runs. Changing the tifffile reading function to:

def read_tifffile():
    with TiffFile(os.path.join(DATASET, FIRST_FILE)) as tif:
        series = tif.series
        for p in POSITIONS:
            _ = series[p].asarray()

while the rest of the script is kept the same gives:

The 1000+ second outlier happened in the first run:

      tifffile  iohub
0  1456.934305  28.823575
1    47.703576  28.290186
2    47.394221  27.825434
3    47.895452  27.775385
4    46.929777  27.824112

I'll have to run it a couple more times to see if that's a consistent pattern or just a network glitch.

[ziw-liu]

I'll have to run it a couple more times to see if that's a consistent pattern or just a network glitch.

Tested again with iohub running first in each loop:

    tifffile  iohub
0  47.075930  27.999947
1  47.212810  27.868444
2  47.222090  27.899368
3  47.039293  27.784547
4  47.265484  27.865451

[cgohlke]

Interesting. I guess that "glitch" is due to caching and abysmal performance of parsing the TIFF IFD structures on a network drive.

I have just released tifffile 2023.2.27, which uses the Micro-Manager indexmap instead of parsing the IFD chain. That should improve timing, especially when files are on a network drive.

Another thought: are the image data returned by iohub and tifffile the same, i.e. (4, 97, 2048, 2048) uint16 arrays?

[mattersoflight]

A special parser would be preferable, but (1) the format seems poorly documented compared to OME-TIFF, (2) there are variations between MicroManager versions, (3) there is no public set of test files, and (4) MicroManager frequently produces corrupted files (at least on my system).

The special parser in iohub was developed primarily with MicroManager gamma datasets. We will need to ask @nicost how data formats have evolved. We can provide example datasets and collaborate on the parser for MM2 gamma maintained with tifffile package.

cc: @edyoshikun, @talonchandler.

[ziw-liu]

Another thought: are the image data returned by iohub and tifffile the same, i.e. (4, 97, 2048, 2048) uint16 arrays?

Yes, they are the same.

I have just released tifffile 2023.2.27, which uses the Micro-Manager indexmap instead of parsing the IFD chain. That should improve timing, especially when files are on a network drive.

Great work! The timing did improve a lot. I just tested the 1 TB dataset on other compute nodes linked to the same storage server, and the initial delay seems to be a pattern of the network storage caching.

Tifffile running first:

     tifffile  iohub
0  133.257984  31.177287
1   17.239525  29.169270
2   17.044940  29.175739
3   17.088925  29.221248
4   17.046381  29.282445

iohub running first (2 different nodes):

    tifffile  iohub
0  15.708963  62.636006
1  13.058120  20.977946
2  12.750226  20.705275
3  12.794518  20.277593
4  12.784348  20.557426

    tifffile   iohub
0  29.763436  188.030651
1  18.895259   32.523990
2  18.771151   31.426238
3  18.842250   31.213109
4  19.288997   31.242613

In later iterations tifffile is now consistently faster, and its performance on 'fresh' nodes is now usable. A remaining question is the first shot performance between the two.