Napari custom lazy-classification layer

Hi,

I’m trying to use Napari as a frontend for a deep learning classifier by creating a custom lazy-classification layer. My goal is to only run a region through a classifier when a user zooms into it. What’s the best way to go about doing this?

I originally thought the best way to do this would be to subclass Layer, and update the output of data() whenever it was accessed. However, I tried to create a basic subclass of Layer (just outputs a static array), but I can’t seem to even get Napari to display the layer:

class MyLayer(Layer):
    def __init__(self):
        self.array = np.arange(512**2, dtype=np.uint8).reshape((512, 512))
        self.array_extent = np.array([[0, 0], [512, 512]])
    @property
    def data(self):
        return self.array
    @data.setter
    def data(self, array):
        self.array = array
    @property
    def _extent_data(self):
        return array_extent
    @_extent_data.setter
    def _extent_data(self, array_extent):
        self.array_extent = array_extent
    def _set_view_slice(self, slice):
        pass
    def _basename(self):
        return 'MyLayer'
    def _get_ndim(self):
        return 2
    def _get_state(self):
        return {}
    def _get_value(self):
        return (12, 34)
    def _update_thumbnail(self):
        self.thumbnail = np.ones((1,1), dtype=np.uint8)*128

Is the best approach to create a subclass of Layer, or is there a better way to repurpose the existing ImageLayer? If I should be subclassing Layer, what methods do I have to implement to get it to work?

Thanks,
Arjun

I think you could probably do this without creating a custom layer.

It should be possible to write a function that takes in an image, applies the pixel classifier, and returns a label image. Since you don’t want to run the classifier when there are too many pixels (presumably that would be too slow), probably we want something in there that checks the zoom/extent of the input image before running the classifier.

@talley has an example here that lazily computes something (in this case an edge filter) from a dask image, and shows both the input and output of this in napari. I think this functional approach is going to be simpler than creating new custom classes.

3 Likes

I’ll add one more suggestion: depending on how zoomed in “zoomed in” is, and how big your data is, you can sort of hack this with a multiscale layer: add_image([classifier_output_array, lower_resolution_source_data]) will display the source data when zoomed out, and then the classifier data once the data pixels match your screen pixels.

2 Likes

I’d go with two multiscale layers in this case, one with the source data, one with the classification. I’ve done something like that before here, where pyramid was a multiscale set of rgb images as dask array and theshold_block was a simple thresholding function.

    viewer.add_image(pyramid, name='source data', multiscale=True)
    thresholded = [ar.map_blocks(lambda block: theshold_block(block, value), drop_axis=2, chunks=ar.chunksize[:-1], dtype=np.int)[:ar.shape[0], :ar.shape[1]] for ar in pyramid]
    viewer.add_labels(thresholded, name='thresholded', multiscale=True)

I could put together a longer example if it would help too

2 Likes

@sofroniewn Multiple layers is great but note that @Arjun.Vikram asked to only show the classifier output when zoomed in.

Anyway @Arjun.Vikram when you do have something working it would be great to see a gif here, seems like a cool application! :blush: In the meantime, keep the questions coming!

2 Likes

I haven’t had a chance to try these out yet, but I will certainly post a GIF and the code I used once I get it working!

@GenevieveBuckley @jni @sofroniewn thanks for your help so far!

I’ve successfully gotten lazy classification to work using Dask’s map_blocks and map_overlap!

Unforunately, Napari eagerly computes the entire lazy dask array (which is time-consuming) when the viewer window opens because the full image is in view. I attempted to use a fake pyramid with a blank lower-resolution array, but I am not able to control when napari switches to the full-resolution classifier output.

Here is a minimal example:

# %%
import napari
import numpy as np
import dask.array as da
from skimage.data import astronaut
from time import sleep

# %%
img = np.tile(astronaut()[:,:,0], (4, 4))
for i in range(0, img.shape[0], 512):
    for j in range(0, img.shape[1], 512):
        img[i,j]=10*i//512+j//512 # tag upper left pixel with 10*r+c
imgd = da.from_array(img, chunks=(512, 512))
blank = np.zeros(
            (img.shape[0]-1,  img.shape[1]-1)
            # (img.shape[0]//2, img.shape[1]//2)
            # (img.shape[0]//4, img.shape[1]//4)
        , dtype=np.bool)

# %%
def transform(arr, log_label="LAZY "):
    """ This is normally where the classifier would be run """
    sleep(1)
    print(f"[{log_label}] transform(shape={arr.shape}, tag=(row: {arr[0,0]//10}, col: {arr[0,0]%10}))\n", end='', flush=True)
    return arr>128

# %%
with napari.gui_qt():
    viewer = napari.view_image(img, name="Original")
    # viewer.add_labels(transform(img, log_label="EAGER"), name="Eager")
    viewer.add_labels([imgd.map_overlap(transform, depth=0), blank], name="Lazy")

This is what happens when I run it and slowly zoom in:


I was zooming in continuously and slowly, the lag seen is because of the slow computation of the entire lazy array. The log shows that all 16 chunks were computed, starting after I zoomed in the second time.

Edit: without the sleep(1) adding artificial delay, this is what happens:

I have tried many dimensions for blank, but I can’t find one that delays the lazy computation for long enough. A larger blank appears to delay the computation for longer, but I can’t make blank larger than the image itself.

I was hoping to have finer control over how small an area must be visible before the lazy array in that region is computed. I don’t need exact control, but I would like to be able to tune it to the point where local classification over that area is almost instantaneous on our hardware (so napari doesn’t get frozen while computing the slide). How can I best achieve this?

Could this maybe be because the blank array isn’t chunked?

Playing around with this a bit more, I think my earlier suggestion is not correct. I think our mental model of when dask computes stuff we see in napari is incorrect.

What I thought happened: Dask only computes results for the image chunks visible on the napari canvas.
What actually happens: Dask computes the whole array, regardless of which chunks are visible.

Perhaps the fancy stuff is only happening in vispy, what do you think @sofroniewn?

Here’s what I did:

%gui qt

# %%
import napari
import numpy as np
import dask.array as da
from skimage.data import astronaut
from time import sleep

# %%
img = np.tile(astronaut()[:,:,0], (4, 4))
for i in range(0, img.shape[0], 512):
    for j in range(0, img.shape[1], 512):
        img[i,j]=10*i//512+j//512 # tag upper left pixel with 10*r+c
imgd = da.from_array(img, chunks=(512, 512))

# %%
def transform(arr, log_label="LAZY "):
    """ This is normally where the classifier would be run """
    sleep(1)
    print(f"[{log_label}] transform(shape={arr.shape}, tag=(row: {arr[0,0]//10}, col: {arr[0,0]%10}))\n", end='', flush=True)
    return arr>128

# %%
viewer = napari.view_image(img, name="Original")
# ZOOM IN MANUALLY TO JUST ONE OF THE TILES

# %%
lbls = imgd.map_overlap(transform, depth=0)
viewer.add_labels(lbls)
# ZOOM BACK OUT MANUALLY, OBSERVE 
# WHAT HAPPENED: We see the labels are already computed for the whole tiled image (otherwise there would be one seoncd of delay)
# WHAT I EXPECTED: I expected the map_overlap function to only compute image regions visible on the napari canvas

Did you mean to add a viewer.add_labels(lbls) in the last cell?

Anyway, I believe the entire dask array is being computed when Napari creates a pyramid out of the provided image.

I’ll try playing around with manually providing a pyramid to see if that can be changed

Did you mean to add a viewer.add_labels(lbls) in the last cell?

Oops yes I did, thank you!

Anyway, I believe the entire dask array is being computed when Napari creates a pyramid out of the provided image.

There’s nothing multiscale happening in the examples above, btw. If you’d like to try that next, I suggest taking a look at the example here: https://github.com/napari/napari/blob/2c3e696bd33558ea5ca50b936dfadece8c87a433/examples/add_multiscale_image.py

From the napari docs,

This example had precomputed image pyramids stored in a zarr file, which is best for performance. If, however you don’t have a precomputed pyramid but try and show a exceptionally large image napari will try and compute pyramids for you unless you tell it not too.

Do you know what the threshold for “exceptionally large image” is? I originally misread this and assumed that all images were pyramidized by Napari.

I don’t know, but I did find this guess_multiscale function that looks like it’s the bit that does that.

Without multiscale, napari does indeed compute the whole 2D image. When displaying multiscale, though, it only computes what’s on the canvas.

Regarding controlling when napari loads the lowest-res tiles, the answer is no, you can’t actually control this, except through hacks Update: I’m an idiot, and forgot about @sofroniewn’s awesome camera model PR. See below!

napari loads the highest resolution tile exactly when the number of pixels on the canvas (ie on your monitor) exceeds the number of pixels displayed in the currently loaded layer. Having said this, napari does no work if a layer is not visible! So we can actually tie layer visibility to a function that responds to events on viewer.camera.events.zoom. Code:

import numpy as np
import toolz as tz
import napari
from skimage import data, filters
import dask.array as da


@tz.curry
def toggle_edges_vis_on_zoom(
        viewer, event, layer_name='edges', zoom_threshold=9
        ):
    layer = viewer.layers[layer_name]
    layer.visible = viewer.camera.zoom > zoom_threshold


coins = data.coins()
edges_high_res = filters.farid(coins)
# make a dummy lower-res array to trigger multi-scale rendering
dummy_edges = da.zeros(tuple(np.array(edges_high_res.shape) // 2))
edges = [edges_high_res, dummy_edges]

with napari.gui_qt():
    viewer = napari.view_image(coins)
    viewer.add_image(
            edges,
            colormap='magenta',
            blending='additive',
            visible=False,
            )
    viewer.camera.events.zoom.connect(toggle_edges_vis_on_zoom(viewer))

Result:

Pretty darn cool if I do say so myself! :smiley:

3 Likes

This is great, thank you @jni

Also, I think we should add a version of this example (maybe with a dask array of tiled/chunked coins) to the napari examples/ folder

2 Likes

These docs are out of data too I think and we don’t do that computation anymore. Needs to be fixed. Sorry.

This works really well, thank you so much!

I’ll try to make a GIF of this in action when I figure out dependency issues and get it applied to our actual slides with real CNN classifiers.

3 Likes