Napari benchmark with custom object layer

Dear Napari team (and everyone else!),

I am working with @joeljo on using Napari using APR data (GitHub - AdaptiveParticles/LibAPR: Library for producing and processing on the Adaptive Particle Representation (APR).). Briefly, APR allows to represent data more efficiently than using pixels. The sparser the data, the higher the speed up in computation and memory footprint.

As a first step we implemented a simple hack consisting in giving Napari an object with a getitem method so that the pixel image is reconstructed on the fly. Doing that we can benefit from all Napari features out of the box for 3D, blending, etc - which is great.

Now, we are using Napari to display the stitching results of many tiles and it gets really slow. The conversion from APR to pixel should be fast (It takes Napari more than 1 second to display 16 tiles even though we are able to load, stitch and globally optimise the whole 4x4 tiles volume in 2 s.) and I tried to profile Napari to understand where is the bottleneck. I first tried to print some elapsed time in the getitem method but it does not print in PyCharm terminal nor in Napari terminal. I tried to return the object used to create the layer after using calling Napari but the object is the same as the input object, as if Napari deep copy the data object under the hood. I quickly looked into the source code but couldn’t find anything useful. What’s the best way we could profile (time and memory) Napari?

I am also curious on the best way to implement the support of APR in Napari. Should we create a custom layer or an entire plugin? I’d happy to give you more detail on APR or what we are trying to achieve.

Cheers,
Jules

Here is some info on performance monitoring in napari Performance Monitoring - napari using some custom tooling we’ve added and profiling Add profiling documentation by Czaki · Pull Request #1998 · napari/napari · GitHub.

My initial hunch is that it would be best as an io plugin, see Plugins - napari.

You might take some inspiration from something like GitHub - manzt/napari-lazy-openslide: An experimental plugin to lazily load multiscale whole-slide tiffs with openslide and dask. which is a zarr interface to a multiresolution tiff format that @manzt made. In your case though you’d try and create a valid Zarr store from your APR. If slicing into that zarr and working with tiles within that zarr is fast then napari will be able to use it. We’re in general working on improving our multi-tiled rendering, see Asynchronous Rendering - napari, but at least we’ll know that things are good outside napari, and you APR will be usable now with tools like dask and anyone else who wants a numpy like representation of your data.

Let me know if i’m way off with that idea, happy to chat more. I’ve read one of the APR papers but not used it myself, so having napari support for it would be very exciting!

Thanks for the suggestions @sofroniewn !

An io plugin would be awesome, but I think that can wait a little. The first priority would be to enable napari to actually visualize APRs in an efficient way. I’ll try to provide some more details on our current approach/hack, and things that we would like to do.

Brief explanation of our approach

Our situation is the following: we have a sparse multi-resolution data structure representing a pixel image, with the ability to reconstruct the original pixel values with a bounded error. So as a quick hack to visualize APRs by z-slice in napari, we created a class that mimics the behavior of a uniform pixel array. Basically, it holds the sparse data structure, and reconstructs pixels on the fly via the __getitem__ method:

import numpy as np
import pyapr
import napari

class FakeArray:
    def __init__(self, apr, particle_values):
        self.apr = apr
        self.parts = particle_values
        # preallocate pixel buffer (assuming only z-slices are requested)
        self.slice = np.empty(shape=(self.shape[1], self.shape[2]))

    @property
    def shape(self):
        return self.apr.org_dims()
    
    def __getitem__(self, item):
        # reconstruct pixels in 'item' slice region
        pyapr.reconstruct(self.apr, self.parts, self.slice, ...)
        return self.slice

We can then wrap a FakeArray object in a napari Image layer for visualization. This works pretty well now, at least for viewing a single APR. Previously (at the time of the original post) there was a large overhead when requesting the first slice, as if the entire volume was reconstructed. This was on our side though, no issues with napari. I’m not sure if this was also the culprit behind @Jules issues.

What we would like

If we can tap a bit deeper into napari we could implement quite a lot of features for APR visualization. For example:

  1. 3D rendering (maximum intensity ray cast) directly on the APR. That is, producing 2D projections as pixels for display, without ever reconstructing the 3D pixel volume.
  2. Slice reconstruction and ray casting on the GPU.
  3. There is a fairly large overlap between the APR and other concepts you support (e.g. multi-res pyramids and octrees). It may be possible to come up with some clever solutions in the intersection there. For example, from a single APR object one could reconstruct/project at different resolution levels depending on the zoom.

It would also be cool to have some additional widgets in the GUI to manually control e.g. the resolution of the reconstruction and the interpolation mode from APR to pixels.

Do you have any recommendations for implementing functionality like this? How much control could one get via, for example, a custom layer implementation? I tried to wrap my head around that, but it seems like custom rendering will require going deeper (perhaps a custom vispy layer?)

I hope this made things somewhat clear. If not, please let me know - I’d be happy to discuss further and provide more details.

Thanks for the additional details here! I think the FakeArray approach will be of general utility in the python ecosystem and being able to get the dense nD representation for either the whole data or parts of the data easily will be beneficial. We’ve recently been doing some work with the sparse package Sparse — sparse 0.11.2+0.g900c9f9.dirty documentation to try and get sparse arrays working for our labels layer.

That being said I do see how for some things like 3D rendering/ max intensity ray casting it would be really nice not to move into the dense representation and do the projection directly on the APR. That would be a huge time save I imagine.

Similarly doing the slice reconstruction and ray casting on the GPU sounds very nice too.

For both of these things we’d probably need to go the custom layer and vispy visual route.

For intersecting the APR with our multires pyramid/ octree work I’d again go back to zarr. @joshmoore has been working on a multiscale zarr spec see some discussion here Multiscale arrays v0.1.

Stepping back I see two important things:

(1) would be to try and get a multiscale zarr store backed by an APR, curious @joshmoore if you have any advice here too, but I think that would be a really nice way to get the APR into the hands of the scientific python community in an API that many other tools will be able to handle right away. Given that we should be able to do some nice visualization in napari with no additional modifications. napari wouldn’t actually know about the APR, it just knows about how to work with the zarr spec, and you could feed this dense representation now into any algorithm that wanted one.

(2) would be to pass the APR representation directly to napari via a new APRLayer which has a new vispy APRVisual associated with it. This visual would do the 2D and 3D slice reconstruction on it (we’d probably keep the nD slicing on napari and any tiled stuff we could defer till later, or use some of our experimental TiledVisuals). You’d need both the custom image layer and custom visual.

Ultimately option (2) is going to be much higher performance, and much nicer to work with, but it will be a bigger lift. Even then you’ll probably still want option (1) as some people will still want to get to the dense numpy-like array interface for analysis until there are APR optimized algos for everything.

So I guess I’d encourage giving (1) some investigation. For option (2) I can start trying to guide you through the code, but we might want to hop on a video call some time to go into details. Let me know if you want to come to one of our developer meetings and what timezone you are in, we have a couple.

The last time I discussed APR with the community (I imagine I2K 2018), I think the consensus was: to make full use of it, e.g. in Bio-Formats and OMERO, there would need to be specialized APIs. Otherwise, it amounts to a codec which is always converting into the dense representation. (option (1))

The same holds for OME-Zarr. Independent of APR, I think there’s already the need to allow returning something other than [DaskArray,...] from ome-zarr-py. I guess I would turn the question around and say what would that look like in the case of an APR array? Happy to work through that as a first extension of the binary return type in ome-zarr-py. (option (2))

~Josh

1 Like

I agree! I will try to expand it with a setter, which can interpolate a (sub-) image back to the particles. Then one can use pixel algorithms on reconstructed chunks and either a) interpolate the result back to the particles (if the particle/sampling locations still make sense) or b) store the result in a uniform array and possibly compute a new APR from that.

So my idea for enabling multiscale functionality in napari for the APR was quite simple. Basically, the FakeArray class can return reconstructions at different resolution levels. So instead of emulating a pixel array we could emulate a pyramid / list of arrays. This is basically possible already by giving napari a list of FakeArray objects, but I figured it could be made more elegant.

Perhaps it is possible to take this one step further and emulate an OME-Zarr reader? That is, the underlying data is an APR (in memory) and we’d create virtual nodes corresponding to different image regions and resolutions. When requesting data from a node the pixels are reconstructed and returned as an array. Is this also how you imagined it @sofroniewn ?

On a related note, we are just now starting to use the APR for tiled datasets (by converting each tile into an APR), but we don’t have any IO support for this yet. So in the near future we’ll looking into storing groups of APRs and enabling something similar to OME-Zarr in a pure-APR way (i.e. return type is APR data structures). I would also be happy to discuss this further, to see where there may be points of overlap and potential collaborations.

Option 2 sounds like exactly what we want for APR visualization. Any help to make it happen would be most appreciated! I’d be happy to chat 1on1, or join a developer meeting, or both. I’m in Central European Time, but have no issue with late evening meetings if you’re in the US.

1 Like

That sounds very nice!

Ok great, we’ve also recently made some improvements to our single scale tiled rendering that will be helpful here too.

Nice, I will send you a DM with details. - Also anyone reading that is interested in attending one of our developer meetings please just reachout, we coordinate them using our napari zulip - napari