JPEG2000 codec support (attempt)

[xinaesthete]

Attempting to implement #774

Background

I've been trying to implement this and I'm running into some issues.

As suggested, I've registered a codec for 34712 with https://www.npmjs.com/package/@cornerstonejs/codec-openjpeg and it's been able to correctly process one sample I passed (which was very small, with only one channel & a few Z slices, fwiw), but not others.

I'm trying to get Xenium images to render (from here, for example - I should track down where my smaller successful test image was from), but would also like to use JP2K more generally. This could be via support added to Viv, or by an app using it by registering the decoder itself - but it's useful for experimentation to have a version running in Avivator.

I've also been experimenting with using a local build of openjpegjs - and may also consider publishing a slightly more ergonomic/typed version as an npm module, but as of now I am not getting any different results in terms of actually properly parsing and rendering the data. I don't think the issue is within the codec, but wouldn't rule it out.

One thing I would note as well is that I have found JP2K images I have converted with bftools don't render correctly in QuPath either - I would appreciate some good reference parameters if such are known.

The issue I currently have when trying xenium images is that geotiff.js throws an error

geotiff.js:461 Uncaught (in promise) GeoTIFFImageIndexError: No image at index 1 at geotiff.js:461:15 
at async geotiff.js:459:27
at async geotiff.js:459:27

I've started debugging this - the part of Viv which calls it is in packages/loaders/src/tiff/lib/indexers.ts:

export function createOmeImageIndexerFromResolver(
  resolveBaseResolutionImageLocation: OmeTiffResolver,
  image: {
    size: { z: number; t: number; c: number };
  }
) {
  const ifdCache: ImageFileDirectory[] = [];
  return async (sel: OmeTiffSelection, pyramidLevel: number) => {
    const { tiff, ifdIndex } = await resolveBaseResolutionImageLocation(sel);
    // this is throwing in geotiff.js `requestIFD(1)` with test xenium jp2k image.
    const baseImage = await tiff.getImage(ifdIndex);

Checklist

• Update JSdoc types if there is any API change.
• Make sure Avivator works as expected with your change.

[xinaesthete]

Just noticed that the test task is timing out in CI, which is also reproducible locally.

[xinaesthete]

ok, so I think the issue is the way xenium stores these in a multi-file pyramid structure - so this is liable to entail a bit more work, and I think the way I've been testing by dragging a particular image into Avivator won't be able to work as the browser wouldn't allow accessing these siblings like that. Maybe if the directory was dragged it'd be ok.

[xinaesthete]

https://www.10xgenomics.com/support/software/xenium-onboard-analysis/latest/analysis/xoa-output-understanding-outputs#images

The morphology_focus/ directory contains the 2D autofocus projection images for the nuclei DAPI stain image, as well as three additional stain images for Xenium outputs generated with the multimodal cell segmentation assay workflow. These files are in multi-file OME-TIFF format. They each contain a pyramid of images including full resolution and downsampled images. The image order is specified in OME-XML metadata.

I'm having a look at whether I can patch loadMultiTiff to understand this better (with multiple files selected in the avivator file chooser), but might move on to other things soon.

[keller-mark]

You should be able to pull the latest changes from main, where I recently changed the test framework. I think the hanging previously indicated test failures under the old framework

[xinaesthete]

Thanks @keller-mark, those tests are now passing (not including any new ones for these changes).

I could do with checking on some other JP2K test data - the problems I have now seem to be to do with the way Xenium do multi-file TIFF, which is technically orthogonal to JP2K as I understand it, but it might be useful to combine those things in a PR since I think easy reading of Xenium is relevant to our common interests.

When I've converted things to JP2K OME-TIFF with bfconvert I've ended up with images that also render incorrectly elsewhere (QuPath)... Maybe user error, but I'd like to have at least one sample where that does work, but really if images are being converted it's probably better for the output to be zarr - I don't really want to be running any new mass-conversions to OME-TIFF. The decoder used in this PR is only associated with TIFF as of now. I guess that support in zarr would look quite different, and I'm not entirely clear on the codec support etc currently - https://www.blosc.org/posts/blosc2-grok-release/ is of interest...

I guess maybe some of this discussion could be on image.sc or something.

[ilan-gold]

@xinaesthete If we're creating a xenium reader, I think we can make it a separate package. The current state of loaders is a bit of a footgun when it comes to avivator and I wouldn't want to add to it unless we implemented something like #875 (comment) in which case then I could be open to having something here, although I still would lean towards making it its own package

[keller-mark]

I think Xenium data is important enough that it would be beneficial to have the support in Viv, especially if the decoding adapter code is non-trivial. Would it be possible to not activate the Decoder up-front, but instead either:

• just export the Decoder, allowing a consumer to globally activate it
• or, allow passing an array of additional decoders to loadOmeTiff, which Viv can activate prior to loading the OME-TIFF

Ideally, then with tree-shaking, a consumer who does not need the JP2K support does not receive an increased bundle size.

[keller-mark]

the problems I have now seem to be to do with the way Xenium do multi-file TIFF, which is technically orthogonal to JP2K as I understand it

I don't follow. Why is the Xenium usage of this compression special? Can you point to any references or docs regarding this?

Edit: I understand now what you are saying. Then yes, I think testing with a single-file OME-TIFF that uses JP2K compression would be fine to validate that the Decoder works. Then the multi-file tiff issues can be addressed independently/orthogonally.

[keller-mark]

Regarding multi-file tiff, I think that functionality is independent of the Decoder/decompression functionality. If there are issues with the Viv multi-file tiff support (or the issue that it is not being exposed via the Avivator demo), that could be its own PR.

[xinaesthete]

@xinaesthete If we're creating a xenium reader, I think we can make it a separate package. The current state of loaders is a bit of a footgun when it comes to avivator and I wouldn't want to add to it unless we implemented something like #875 (comment) in which case then I could be open to having something here, although I still would lean towards making it its own package

Yes - I don't disagree (although interesting to hear your phrasing on the state of loaders). This particular PR isn't necessarily meant as something that should be merged in anything very close to the current form - I suppose hypothetically it could be as a relatively simple "jp2k tiff codec support"... but it's clear that xenium support is a wider scope.

I started out by just adding a similar codec handler in MDV, but then when it didn't quite go as easily as I hoped it was useful to test it on a simpler Avivator base, and as I say, it seems like loading xenium is in our common interest... but I started this PR more as a place to experiment with the implementation and hopefully get some more feedback on how it should work.

Just noticed @keller-mark's comments as I was writing this - it seems like the xenium use of multi-tiff may apparently be somewhat standard, but the support viv has for this is designed for something different - and I agree that seems like a separate issue. I have slowly been starting to decipher this a bit better, but I think following through from the xenium spec I linked above is the best reference I have currently.

I haven't thought a huge amount about the design or tested whether it is properly doing so, but the decoder potentially only lazily imports any heavy dependencies when it is actually invoked. Still - I don't know if it definitely wants to be there adding extra bloat to viv in general.

Hope that's somewhat coherent - only intended to unplug my laptop after a couple of work social beers...

[xinaesthete]

ok, so now I'm using raw2ometiff to generate test images, which I seem to have more luck with than bfconvert (although with lossy compression it seems to have a ceiling of ~23mb for the output, from an ~800mb input, the lossless version is 500mb).

The viv rendering has an issue with byte-order or something getting mangled, so the results look very glitchy. I'll see if I can try a few more variations of precision. The image I've been testing with is uint16.

[xinaesthete]

It looks as though the decoder has a block for when componentCount !== 1 && bitsPerSample > 8 which is commented out in the source:

      // Convert from int32 to native size
      int comp_num;
      for (int y = 0; y < sizeAtDecompositionLevel.height; y++)
      {
        size_t lineStartPixel = y * sizeAtDecompositionLevel.width;
        size_t lineStart = lineStartPixel * frameInfo_.componentCount * bytesPerPixel;
        if(frameInfo_.componentCount == 1) {
          int* pIn = (int*)&(image->comps[0].data[y * sizeAtDecompositionLevel.width]);
          if(frameInfo_.bitsPerSample <= 8) {
              unsigned char* pOut = (unsigned char*)&decoded_[lineStart];
              for (size_t x = 0; x < sizeAtDecompositionLevel.width; x++) {
                int val = pIn[x];;
                pOut[x] = std::max(0, std::min(val, UCHAR_MAX));
              }
          } else {
            if(frameInfo_.isSigned) {
              short* pOut = (short*)&decoded_[lineStart];
              for (size_t x = 0; x < sizeAtDecompositionLevel.width; x++) {
                int val = pIn[x];;
                pOut[x] = std::max(SHRT_MIN, std::min(val, SHRT_MAX));
              }
            } else {
              unsigned short* pOut = (unsigned short*)&decoded_[lineStart];
              for (size_t x = 0; x < sizeAtDecompositionLevel.width; x++) {
                int val = pIn[x];;
                pOut[x] = std::max(0, std::min(val, USHRT_MAX));
              }
            }
          }
        } else {
            if(frameInfo_.bitsPerSample <= 8) {
              uint8_t* pOut = &decoded_[lineStart];
              for (size_t x = 0; x < sizeAtDecompositionLevel.width; x++) {
                pOut[x*3+0] = image->comps[0].data[lineStartPixel + x];
                pOut[x*3+1] = image->comps[1].data[lineStartPixel + x];
                pOut[x*3+2] = image->comps[2].data[lineStartPixel + x];
              }
            } /*else {
              // This should work but has not been tested yet
              if(frameInfo.isSigned) {
                short* pOut = (short*)&decoded_[lineStart] + c;
                for (size_t x = 0; x < sizeAtDecompositionLevel.width; x++) {
                  int val = line->i32[x];
                  pOut[x * frameInfo.componentCount] = std::max(SHRT_MIN, std::min(val, SHRT_MAX));
                }
              } else {
                unsigned short* pOut = (unsigned short*)&decoded_[lineStart] + c;
                for (size_t x = 0; x < sizeAtDecompositionLevel.width; x++) {
                    int val = line->i32[x];
                    pOut[x * frameInfo.componentCount] = std::max(0, std::min(val, USHRT_MAX));
                }
              }
            }*/
        }
      }

So I'm going to try a build which includes that, fingers crossed this may help.

edit: I get build errors when I uncomment that block, will investigate further / raise an issue upstream.

edit2: componentCount is 1 there in the cases I'm testing.

[xinaesthete]

ok, so it needs more testing with a few different dtypes but the main issue I was seeing with uint16 images was endianness which I've patched in the C++ now locally. I'll try to get some better tests in place and probably raise a PR / possibly publish a different version of that decoder module with some other updates to export types etc as well.

[xinaesthete]

I'd like to be able to process float32 data as well - but this seems not to be implemented in OpenJpeg. https://github.com/jasper-software/jasper/ may be an alternative - the releases include WebAssembly editions, but target WASI & I think it'll need a bit of a tweak to work.

I also noticed a pure-js jpeg2000 in npm, first pass on that is not working, likely user error, but I think really WASM (or maybe WebGPU in an ideal dream-world) is likely to be necessary for decent performance anyway.

[keller-mark]

Hi @xinaesthete,
Is this PR ready for review? I would be fine to merge the functionality if it is only supporting a limited set of dtypes such as uint16, especially if it would be enough for the Xenium use-case. Two other questions

• Can you share the raw2ometiff commands you are using to generate test images?
• What is the status of the code in packages/loaders/src/tiff/lib/jpeg2000? Would prefer to depend on an NPM package rather than vendor the minified/compiled JS+WASM code.

[xinaesthete]

Hi @xinaesthete,
Is this PR ready for review? I would be fine to merge the functionality if it is only supporting a limited set of dtypes such as uint16, especially if it would be enough for the Xenium use-case.

Unfortunately in the current form I probably wouldn't merge it (I did merge equivalent changes in MDV, but with the intention of future revisions). It'd probably be worth checking out and looking at a few test images.

I'm on vacation until September, but thought I'd respond.

I don't think we necessarily need to aim for supporting float32 etc, but right now the way it only looks right with big endian is a bit of an issue. The xenium images which are able to load directly (morphology.ome.tiff in root iirc) look bad because they use little endian (they were ok before I altered the build of the decoder).

It seems like it should be simple to read the first bytes of the file to determine which order to use, but I'm not sure how to do that in the context of a geotiff decoder. I looked at the lzw decoder implementation and it seems to have a bit-order flag hard coded, but I also believe I've tried lzw images with both orders and it is correct either way. So I'm not sure if some byte switching is happening elsewhere or... I'm not sure what the issue is.

Two other questions

• Can you share the raw2ometiff commands you are using to generate test images?

I don't think there are any other flags than setting the compression type and quality. Certainly couldn't see a way to specify endianness.

I also raised another issue there that the quality seems to flatten out at a relatively low level, but it does get remarkably good image quality at pretty massive reductions in size. Quite a remarkable thing to behold, so I would recommend giving it a try, also maybe if you can find that issue and let them know if you do or don't reproduce that'd be good.

• What is the status of the code in packages/loaders/src/tiff/lib/jpeg2000? Would prefer to depend on an NPM package rather than vendor the minified/compiled JS+WASM code.

Yes, I agree, I would prefer it to be an npm package... I was in a bit of a rush before my holiday. Currently this messy state is what is making biome complain in CI for this PR, which having it installed from npm should fix.

May be worth reaching out to Chris chafey or the cornerstone people who published the existing npm module I started out using.

The state is that I made a few changes to update some dependencies etc so I could get it to build on arm64 Mac, with emitted TypeScript and configured as ESM. I also enabled a SIMD flag which was commented in a build script somewhere as not making a difference at the time. I don't think it makes a big difference, but it seemed maybe somewhat better. Aside from that, I made a somewhat dubious change to the byte order handling, but as commented this part is not really satisfactorily resolved as yet. Also notable that it doesn't change their test results either way iirc, so that'd be another thing to look at.

Another comment while I'm here: I worked a bit on some hobby GIS stuff a few years ago with a version of openjpegjs that handles HTJ2K variant of the codec... It'd be potentially good if that was well supported by ome etc at some point. Faster encoding/decoding. Way out of scope for this PR of course, but worth bearing in mind for future evolution of tools. I'd like to be using that.