Is this conversion to TIFF correct, it's very slow

[jonasteuwen]

Hi,

Thanks for relicensing the library! I wanted to write an isyntax -> tiff converter using libtiff. I came up with this:

#include "libisyntax.h"

#include <stdint.h>
#include <assert.h>

#include "tiffio.h"
#include <string.h>

#define LOG_VAR(fmt, var) printf("%s: %s=" fmt "\n", __FUNCTION__, #var, var)



uint32_t bgra_to_rgba(uint32_t val) {
    return ((val & 0xff) << 16) | (val & 0x00ff00) | ((val & 0xff0000) >> 16) | (val & 0xff000000);
}


void write_pyramid_level(isyntax_t* isyntax, isyntax_cache_t *isyntax_cache, isyntax_level_t* level, TIFF* output_tiff, int32_t tile_size, int32_t num_channels) {
    int32_t num_tiles_x = libisyntax_level_get_width_in_tiles(level);
    int32_t num_tiles_y = libisyntax_level_get_height_in_tiles(level);

    int32_t width = tile_size * libisyntax_level_get_width_in_tiles(level);
    int32_t height = tile_size * libisyntax_level_get_height_in_tiles(level);

    int32_t scale = libisyntax_level_get_scale(level);
    // TODO: Compute the spacing

    TIFFSetField(output_tiff, TIFFTAG_IMAGEWIDTH, width);
    TIFFSetField(output_tiff, TIFFTAG_IMAGELENGTH, height);

    for (uint32_t tile_y = 0; tile_y < num_tiles_y; ++tile_y) {
        for (uint32_t tile_x = 0; tile_x < num_tiles_x; ++tile_x) {
            // Let's print some progress for each percentage of the image that we have processed
            int32_t progress = (int32_t) (100.0 * (tile_y * num_tiles_x + tile_x) / (num_tiles_x * num_tiles_y));
            if (progress % 5 == 0) {
                printf("Progress: %d%%\r", progress);
                fflush(stdout);
            }

            uint32_t *pixels = NULL;
            assert(libisyntax_tile_read(isyntax, isyntax_cache, scale, tile_x, tile_y, &pixels) == LIBISYNTAX_OK);

            // convert data to the correct pixel format (bgra->rgba).
            for (int i = 0; i < tile_size * tile_size; ++i) {
                pixels[i] = bgra_to_rgba(pixels[i]);
            }

            // Copy the tile data to the output buffer
            uint32_t tile_width = (tile_x < num_tiles_x - 1) ? tile_size : (width % tile_size);
            uint32_t tile_height = (tile_y < num_tiles_y - 1) ? tile_size : (height % tile_size);

            if (tile_x == num_tiles_x - 1 || tile_y == num_tiles_y - 1) {
                // TODO: Not sure we need to adjust here for the border tiles?
                // Create a new buffer of the correct size based on the adjusted tile_width and tile_height
                unsigned char *buf = (unsigned char *)_TIFFmalloc(tile_size * tile_size * num_channels);

                // Copy the relevant data from the original buffer (pixels) to the new buffer
                for (uint32_t y = 0; y < tile_height; ++y) {
                    memcpy(buf + y * tile_size * num_channels, pixels + y * tile_width * num_channels, tile_width * num_channels);
                }

                TIFFWriteTile(output_tiff, buf, tile_x * tile_size, tile_y * tile_size, 0, 0);
                // Free the adjusted buffer
                _TIFFfree(buf);
            } else {
                // Pass the original buffer (pixels) to TIFFWriteTile
                TIFFWriteTile(output_tiff, pixels, tile_x * tile_size, tile_y * tile_size, 0, 0);
            }

            libisyntax_tile_free_pixels(pixels);

        }
    }
}


int main(int argc, char** argv) {

    if (argc <= 1) {
        printf("Usage: %s <isyntax_file> <output.tiff> - convert an isyntax image to a tiff.", argv[0]);
        return 0;
    }

    char* filename = argv[1];

    libisyntax_init();

    isyntax_t* isyntax;
    if (libisyntax_open(filename, /*is_init_allocators=*/0, &isyntax) != LIBISYNTAX_OK) {
        fprintf(stderr, "Failed to open %s\n", filename);
        return -1;
    }
    printf("Successfully opened %s\n", filename);

    TIFF *output_tiff = TIFFOpen(argv[2], "w8");
    if (!output_tiff) {
        fprintf(stderr, "Failed to open output TIFF file\n");
        libisyntax_close(isyntax);
        return 1;
    }
    LOG_VAR("%s", argv[2]);

    int32_t tile_width = libisyntax_get_tile_width(isyntax);
    int32_t tile_height = libisyntax_get_tile_height(isyntax);
    LOG_VAR("%d", tile_width);
    LOG_VAR("%d", tile_height);

    isyntax_cache_t *isyntax_cache = NULL;
    assert(libisyntax_cache_create("tiff conversion cache", 2000, &isyntax_cache) == LIBISYNTAX_OK);
    assert(libisyntax_cache_inject(isyntax_cache, isyntax) == LIBISYNTAX_OK);

    int wsi_image_idx = libisyntax_get_wsi_image_index(isyntax);
    const isyntax_image_t* wsi_image = libisyntax_get_image(isyntax, wsi_image_idx);
    int32_t num_levels = libisyntax_image_get_level_count(wsi_image);
    int32_t num_channels = 4;
    int32_t bits_per_sample = 8;

    int32_t quality = 70;
    TIFFSetField(output_tiff, TIFFTAG_SAMPLESPERPIXEL, num_channels);
    TIFFSetField(output_tiff, TIFFTAG_BITSPERSAMPLE, bits_per_sample);
    TIFFSetField(output_tiff, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
    TIFFSetField(output_tiff, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
    TIFFSetField(output_tiff, TIFFTAG_TILEWIDTH, tile_width);
    TIFFSetField(output_tiff, TIFFTAG_TILELENGTH, tile_height);
    TIFFSetField(output_tiff, TIFFTAG_COMPRESSION, COMPRESSION_JPEG);
    TIFFSetField(output_tiff, TIFFTAG_JPEGQUALITY, quality);

    assert(tile_width == tile_height);

    for (int32_t level_idx = 0; level_idx < num_levels; ++level_idx) {
        LOG_VAR("%d", level_idx);
        isyntax_level_t* current_level = libisyntax_image_get_level(wsi_image, level_idx);
        write_pyramid_level(isyntax, isyntax_cache,current_level, output_tiff, tile_width, num_channels);
        if (level_idx < num_levels - 1) {
            TIFFWriteDirectory(output_tiff);
        }
    }

    // Close the iSyntax file and the TIFF file
    TIFFClose(output_tiff);
    libisyntax_cache_destroy(isyntax_cache);
    libisyntax_close(isyntax);
    return 0;
}

However, this is very slow, and my profiler tells me a major part is actually in the read_tile. Am I using the library as intended here?

[Falcury]

I haven't looked into this in detail yet, but my suspicion is that each each invocation of read_tile is currently repeating a lot of work, not caching the intermediate results of the progressive decoding properly for re-use. Maybe @alex-virodov could shine some light as to what might be going on.

In general, for iSyntax decoding to be performant, a smart strategy is needed for (1) I/O of compressed data, i.e. a strategy for reading from disk and (pre)caching the compressed codeblocks/codeblock clusters (2) a caching scheme for the coefficients stored in the codeblocks, including the successive intermediate results of the progressive decoding. We should be able to optimize this further.

Some general stuff that might help a bit:

• Set the cache size as large as possible, maybe multiple gigabytes or as large as RAM allows. (might need to do some tests on this)
• Decode the levels in reverse order, this better matches the progressive decoding scheme.

Also note that the origin offset for each level shifts due to the way the wavelet transform is set up (this is kind of annoying, and one of the reasons I'd like for libisyntax to supply its own read_region implementation to be able to correct for this). For the tiles to line up you need to take into account an offset:

libisyntax/src/isyntax/isyntax.c

Lines 2970 to 2983 in f142dce

	// When recursively decoding the tiles, at each iteration the image is slightly offset
	// to the top left.
	// The amount of shift seems correspond to the level padding: ((3 >> (scale-1)) - 2).
	// (I guess this is related to the way the wavelet transform works.)
	// Put another way: the highest (zoomed out levels) are shifted the to the bottom right
	// (this is also reflected in the x and y coordinates of the codeblocks in the iSyntax header).
	for (i32 scale = 1; scale < wsi_image->level_count; ++scale) {
	isyntax_level_t* level = wsi_image->levels + scale;
	float offset_in_pixels = (float) (get_first_valid_coef_pixel(scale - 1));
	level->origin_offset_in_pixels = offset_in_pixels;
	float offset_in_um_x = offset_in_pixels * wsi_image->levels[0].um_per_pixel_x;
	float offset_in_um_y = offset_in_pixels * wsi_image->levels[0].um_per_pixel_y;
	level->origin_offset = (v2f){offset_in_um_x, offset_in_um_y};
	}

[jonasteuwen]

I've indeed seen this offset. But do I understand correctly that we would need to read more tiles (taking into account the offset) and subsequently crop?

Perhaps for this use-case it would then be best to just read level 0, and use something like libvips to create a dense tiff pyramid.

[Falcury]

I've indeed seen this offset. But do I understand correctly that we would need to read more tiles (taking into account the offset) and subsequently crop?

For now, yes I think so.

Perhaps for this use-case it would then be best to just read level 0, and use something like libvips to create a dense tiff pyramid.

Yes I think this would be a good way to do it. Especially since you would also minimize artifacts due to recompression of levels > 0.

[alex-virodov]

I haven't looked into this in detail yet, but my suspicion is that each each invocation of read_tile is currently repeating a lot of work, not caching the intermediate results of the progressive decoding properly for re-use. Maybe @alex-virodov could shine some light as to what might be going on.

There is some caching. The default cache is small though, tuned for random or sequential access for multiple slides (ML usecase). For one-slide-at-a-time usecase, you probably want to crank it up so more work is reused.

In general, for iSyntax decoding to be performant, a smart strategy is needed for (1) I/O of compressed data, i.e. a strategy for reading from disk and (pre)caching the compressed codeblocks/codeblock clusters (2) a caching scheme for the coefficients stored in the codeblocks, including the successive intermediate results of the progressive decoding. We should be able to optimize this further.

I hope so too. Not currently optimized in libisyntax_tile_read.

Set the cache size as large as possible, maybe multiple gigabytes or as large as RAM allows. (might need to do some tests on this)

+1 :)
I think I see some recent work doing that.

Decode the levels in reverse order, this better matches the progressive decoding scheme.

I'm afraid with small cache it won't matter - if you decode a whole level, it will be flushed out from a small cache later on anyway.
In any case, top level tiles should stay in cache, as it is LRU (least recent use), and when you read a tile at low level, all tiles above it are bumped (marked used) in cache as well.

[alex-virodov]

Also what would help speed is internal threading. There is quite some space for parallelization withing a single tile read - we have to decode all the tiles above and around, if not cached. This is not too hard, but not done yet.

We can also allow invocation of libisyntax_tile_read from multiple threads to allow even more parallelization. This one is harder, as it requires tracking data dependencies across libisyntax_tile_read requests. Also not done.

[jonasteuwen]

I would think having a thread-safe read_region/tile would be very useful, as that also can translate better to ML contexts.

I did some preliminary experiments using memory mapping rather than a standard file pointer to check the effect on performance. Would you expect it to help in our use-cases too? The OS might help with caching too in such a case.

@alex-virodov I will profile the TIFF writing. Is there a preliminary setup for OpenSlide? I can also profile and actual ML setup, I have many annotated slides + model.

[alex-virodov]

Good question, I was just about to ask how I can do unit-tests :)

I do have some unit-tests in OpenSlide that should probably be in libisyntax (and OpenSlide too, as integration tests). I will make a PR with my proposal for libisyntax/test folder, and maybe we can further discuss there.

After we have some canned tests with known good outputs, we can more freely do optimization (having baselines) and other code refactoring (having a way to check the code doesn't change output).

One last note - if you don't see a PR by tomorrow for libisyntax/test, or an update by me, that means I got pulled into other work. Please don't wait and proceed as you see fit.

[alex-virodov]

memory mapping rather than a standard file pointer to check the effect on performance. Would you expect it to help in our use-cases too?

Measurements will tell. I don't see why not, if it results in faster reading (due to OS caching or more efficient OS data move), it helps all usecases, because at the end all uscases need to read the required blocks from disk in some way.

Another thing to keep in mind is that we may want to allow filesystem virtualization - it was discussed in other threads, but I raise it here just to keep in mind that we may want to keep these optimizations separate enough to allow other access patterns to be injected.

[jonasteuwen]

A proper speed test for ML would be perhaps this: Take a list of about 100 files, use a bunch of threads to open/close them and read a random region. Pass it to a function that waits N ms (can figure that out by measuring the time it takes to pass a tile to say a u-net). Then we can collect some profiler information, right? Should give a clear idea where the bottleneck is.

[jonasteuwen]

This is after 3 minutes, only converting level 0 on my Apple M1:

• isyntax-to-tiff`signed_magnitude_to_twos_complement_16_block: 7% of all operations
• isyntax-to-tiff`isyntax_hulsken_decompress: 23% of all (above function is called here)
• ycocggr: 7%

So I don't know if we can avoid many calls to isyntax_hulsken_decompress, but probably not. There is some space to improve ycocggr and signed_magnitude_to_twos_complement_16_block with SIMD instructions. I'll see if I can do that in a separate PR

[alex-virodov]

I think there's 2 separate questions:

1. How many tiles should have been read at minimum for a given user-requested tile sequence (or scenario if we're multithreading).
2. How long each tile read takes in isolation.

Note that in finite cache scenario, it is expected that for 1 the number of actual tiles read would be larger than minimum, due to cache eviction. I would be curious to know how much larger, I do not have a measurement. It will of course depend on cache size and scenario. With cache_size=MAX_INT, there should be no evictions, but there will be cache management overhead.

For 2 we can optimize separately, possibly profiling 1 tile read at a time.

With the profiler picture above, I don't think we can separate 1 and 2 - so we won't know if it takes long because of repeated reads or some inefficiency in decoding/file access.

Take a list of about 100 files, use a bunch of threads to open/close them and read a random region

Exactly, I can set up this test and we can measure that as we optimize.

[jonasteuwen]

Great, that test would be absolutely great! I can try to scan some mouse tissue on the scanner so we have an image available that can be shared.

[alex-virodov]

Having a test slide or two we all consistently use would be good.

In PR #12 I use a test slide posted by OpenPhi:
https://github.com/amspath/libisyntax/pull/12/files#diff-1e7de1ae2d059d21e1dd75d5812d5a34b0222cef273b7c3a2af62eb747f9d20aR74

file(DOWNLOAD https://zenodo.org/record/5037046/files/testslide.isyntax?download=1 ${CMAKE_BINARY_DIR}/testslide.isyntax SHOW_PROGRESS)

It is a large slide, ~435Mb. If you want to scan a smaller slide, that would probably be nice too. I looked at our isyntax library, and the smallest slide we have is 30Mb, probably something we can just drop into the GitHub repo. But our slides are deidentified, so they are not good tests for macro/label images for example.

Also note that OpenSlide seems to host their own test slides (https://openslide.cs.cmu.edu/download/openslide-testdata/)
Maybe @bgilbert can tell us if they are interested in hosting an isyntax test slide, large and/or small.