LabeledImageServer coordinate problems

[alanocallaghan]

• if you specify downsample on init then the coords of all geometries are transformed, which makes the whole affair very confusing, since now the server coordinates don't correspond to anything in the base image.. We don't need to allow this, since we can just generate images at a suitable downsample on request without needing to resize the geometries on init.
  

  qubalab/qubalab/images/labeled_server.py

  Line 57 in 2bc2027

  self._geometries = [shapely.affinity.scale(shapely.geometry.shape(f.geometry), 1/self._downsample, 1/self._downsample, origin=(0, 0, 0)) for f in self._features]

• if you specify a region and/or downsample later on, then we try to draw geometries on an image that doesn't share the same coordinate space

  qubalab/qubalab/images/labeled_server.py

  Lines 117 to 124 in 2bc2027

  	drawing_context = PIL.ImageDraw.Draw(image)
  	for i in self._tree.query(region.geometry):
  	draw_geometry(
  	image.size,
  	drawing_context,
  	self._geometries[i],
  	self._feature_index_to_label[i]
  	)

  
  and

  qubalab/qubalab/images/labeled_server.py

  Lines 96 to 112 in 2bc2027

  	full_image = np.zeros((self.metadata.n_channels, region.height, region.width), dtype=self.metadata.dtype)
  	feature_indices = self._tree.query(region.geometry)
  	labels = set(self._feature_index_to_label.values())
  	for label in labels:
  	image = PIL.Image.new('1', (region.width, region.height))
  	drawing_context = PIL.ImageDraw.Draw(image)
  	for i in feature_indices:
  	if label == self._feature_index_to_label[i]:
  	draw_geometry(
  	image.size,
  	drawing_context,
  	self._geometries[i],
  	1
  	)
  	full_image[label, :, :] = np.asarray(image, dtype=self.metadata.dtype)

  
  if we fix the downsample on init, then we transform the geometries on request based on the downsample and region, which seems alright

We could allow both, which makes the code a bit more complex and arguably more confusing.

[petebankhead]

I see in the docs for labeled_server.py:

# The image will only have one resolution level.

This is also the case for QuPath’s Java LabeledImageServer.

I don’t know / remember the exact reason, but I suspect it is related to the fact that downsampling labeled images can have weird effects because of interpolation - and we’d likely get very different results if we downsample a rasterized version vs. downsampling the geometries, and a standard ImageServer uses rasterized resizing between pyramid levels.

The thinking was that you’d create a labeled image at the exact downsample you’d want to export, then discard it afterwards. This means ‘later on’ is not expected to be relevant, and export at the downsample specified from the beginning should be free of any rasterization artifacts.

Our use of arbitrary downsamples in requests later has posed problems in the past within QuPath. If (for example) we want to export 512x512 tiles with downsample 2.3235245 then this can be very awkward to do from external code because we can’t specify the output tile size. Rather, we have to figure out the coordinates to request in the full-resolution image and trust that QuPath will make the same internal resizing/rounding decisions that we expect.

Furthermore, we can have trouble on tile boundaries, where some tiles will have coordinates rounded up and some rounded down. For most ‘natural’ images this isn’t a big deal, but could be more obvious for labeled or binary images if it can result in one-pixel gaps or line thickening. These problems are avoided/reduced if the size is fixed at the point the ImageServer is created, and so tiles are defined only at that resolution.

For a ‘natural’ image in QuPath, we can use PyramidGeneratingImageServer to fix the downsample value and compute more reliable tiles - but it will resize raster images, so shouldn’t be used with binary or labeled images.

Lastly, the QuPath version supports outlining regions - and then defining one downsample is crucial.

Changing this may well may sense, but I think we need to proceed with caution and lots of example use-cases - and we should consider both the QuPath and QuBaLab behavior together.

[alanocallaghan]

I don’t know / remember the exact reason, but I suspect it is related to the fact that downsampling labeled images can have weird effects because of interpolation - and we’d likely get very different results if we downsample a rasterized version vs. downsampling the geometries, and a standard ImageServer uses rasterized resizing between pyramid levels.

The problem I have with the current design is that we allow both. First we do this conversion of geometries to the downsampled scale, but then read_region calls the version of _read_block defined in the labeled server class, and possibly does rasterized resizing on the result, which is computed on the fly from the geometries anyway! Why not just downsample the geometries to the correct level on request?

The thinking was that you’d create a labeled image at the exact downsample you’d want to export, then discard it afterwards. This means ‘later on’ is not expected to be relevant, and export at the downsample specified from the beginning should be free of any rasterization artifacts.

Then we should at least disallow people from setting the downsample twice by overriding read_region. Furthermore, as far as I can tell we can accomplish the same goal of avoiding rasterisation effects by just computing the images at the correct scale every time based on whatever downsample is requested --- this is what I was suggesting.

Anyway for the time being I've just added failing unit tests to demonstrate that it's currently broken.

[petebankhead]

If that’s all #26 does, can you give it a more meaningful title and description please?

Admittedly I haven’t tried to run any of this since the major refactoring, but based on the description here, I don’t understand why the current behavior is ‘broken’.

And if I understand the alternatives correctly, I fear they might easily be broken in worse ways.

If I remember correctly _read_block should be always be called at a fixed downsample level where tile sizes can be known and computed without rounding errors. Any call to read_region with a downsample that doesn’t exactly match a ‘fixed’ level will always risk issues with rasterization (no matter what ImageServer implementation we have) - because it will call _read_block at the most appropriate fixed level, then resample the pixels.

Therefore it’s perfectly valid to call read_region at any downsample you like, but there has to be a way to ensure that a LabeledImageServer can be created that doesn’t have rasterization effects at the specific downsample you know that you’re going to want to use.

I can think of two alternatives, but both are problematic:

1. We support multiple ‘levels’, i.e. fixed downsamples in a LabeledImageServer, rendering directly by rescaling the geometries at only these downsamples. However, there could be very distinct ‘jumps’ in how the image appears when requesting at specific downsamples, depending upon which level is supplying the pixels and if rasterization is needed.
2. We avoid _read_block and always rasterize at the downsample requested by read_region. Then we either can’t cache tiles, or we risk filling up the cache quickly by caching all tiles for any arbitrary downsample. Also, if we support outlines, then the line thickness effectively ‘grows’ as the downsample increases (e.g. it’s 1 pixel at the requested resolution)… or else we get horrible artifacts by decreasing the thickness to something much less than 1 pixel when using a high downsample.

Think how it might work in a viewer. There, pixels can be displayed at any downsample. If we always request at the display downsample, computing this every time would be extremely computationally intensive - we’d need to make an enormous number of requests when zooming in or out.

What we probably want are that the labels are consistent and trustworthy at one resolution, but can be requested (with raster-based downsampling) at specific lower resolutions for a fast display. Also, the viewer should be able to rendered at any necessary resolution in between, again with raster-based scaling.

[alanocallaghan]

If that’s all #26 does, can you give it a more meaningful title and description please?

At the moment it's just failing unit tests. If the implementation does not change, then I'll change the tests so they pass and document the current behaviour more clearly. I'll add a proper description either way, but there's not content to document yet.

Admittedly I haven’t tried to run any of this since the major refactoring, but based on the description here, I don’t understand why the current behavior is ‘broken’.

Maybe I'm misunderstanding, but it's certainly not intuitive at the moment. That's why I started with failing tests. Hopefully the tests demonstrate what is currently unintuitive (it arose from a user testing the notebooks).

If I remember correctly _read_block should be always be called at a fixed downsample level where tile sizes can be known and computed without rounding errors.

As far as I can tell the relevant calls are in the base ImageServer class, where for a LabeledImageServer with only one level, the blocks are requested with no downsample, then resized with rasterisation, or served as a dask array with no resizing. It's not clear to me why:

labeled_server = LabeledImageServer(qupath_server.metadata, annotations, downsample=2)
labeled_server.read_region(downsample=1, x, y, width, height)

should be substantially different to

labeled_server = LabeledImageServer(qupath_server.metadata, annotations, downsample=1)
labeled_server.read_region(downsample=2, x, y, width, height)

in terms of how we process annotations, or in terms of the expected output, or in terms of the x,y,width,height coordinates we supply. Currently, the two calls won't return the same region of the parent server, and one will use rasterised resizing. We should definitely be letting users know that one of these is "right" and one is "wrong", if we keep the current implementation.

We avoid _read_block and always rasterize at the downsample requested by read_region. Then we either can’t cache tiles, or we risk filling up the cache quickly by caching all tiles for any arbitrary downsample.

Is there caching currently happening in qubalab? I can't identify any. If you mean on the qupath side, I understand the desire to match behaviour exactly, but I also fear writing a bad implementation here just to match behaviour in QuPath.

Therefore it’s perfectly valid to call read_region at any downsample you like

If it's something that will produce unintuitive results we should at least be warning users that they're probably not doing things the "right" way?

[alanocallaghan]

There is also one very concrete way in which the current implementation is definitely broken. We try to draw on the pillow based on the coordinates of the entire LabeledImageServer, without considering that we may actually be in a region which is offset from the origin.

[petebankhead]

No caching in QuBaLab currently I think, but there is in QuPath (and potentially could be in QuBaLab in the future).

I think working off specific use cases would be useful. To me, the ImageServer stuff here should act in pretty much the same way as it does in QuPath but most of the time it shouldn’t be used directly, because dask arrays are more intuitive in Python.

Potential use cases:

1. We have an ImageServer with the the original pixels, and a companion LabeledImageServer. We want to export corresponding single-resolution regions at a fixed downsample for training a model.
2. We want to view the labeled image in napari as an overlay on the original image.
3. We want to export a pyramidal labeled image with fixed resolution levels. These may or may not correspond to the resolution levels of any companion image.

In all cases, we can use dask arrays.

I think that, ideally, we’d hide as much ImageServer logic as we can, and provide alternatives - which may well by backed by ImageServer logic, but most users should never have to worry about that. The main exceptions are users who are already more familiar with QuPath scripting, for whom the least confusing behavior is whatever is consistent across both implementations.

It's not clear to me why:

labeled_server = LabeledImageServer(qupath_server.metadata, annotations, downsample=2)
labeled_server.read_region(downsample=1, x, y, width, height)

should be substantially different to

labeled_server = LabeledImageServer(qupath_server.metadata, annotations, downsample=1)
labeled_server.read_region(downsample=2, x, y, width, height)

in terms of how we process annotations, or in terms of the expected output, or in terms of the x,y,width,height coordinates we supply.

Because in the first you’re requesting that the labels are generated at a downsample of 2, then you are (weirdly) requesting them to be upsampled by a factor of 2. It’s like you’ve exported the image with downsample=2 and then zoomed in. The upsampling cannot create new information or detail that is not present in the original (downsample=2) version.

In the second, you’re requesting the labels to be generated at the full image resolution, then you are requesting a quick downsampled version of them. This will lose information and have rasterization artifacts, as raster-based downsampling usually has. It’s like you’ve exported the image at full resolution and zoomed out.

The general idea of an ImageServer is that it should act like it is reading pixels from a pyramidal image - even if it is actually generating them on-the-fly. Any normal pyramidal image server doesn’t recompute pixels for any arbitrary resolution: the pixels are available only at fixed levels, determined when the pyramid is first written.

If a user simply wants a labeled image as a dask arrays, then they shouldn’t need to worry about ImageServer internals… but they will need to specify the resolution for the dask arrays to be generated. Internally, we can do that using a LabeledImageServer.

With regard to coordinates that need to be passed, I don’t know… I haven’t used QuBaLab recently enough to know what it’s doing or to confirm it’s the same as QuPath. But as long as I can request my original image at one downsample and a labeled image at the same downsample, I’d rather get them both as dask arrays of the same size and just use indexing - not have to deal with lots of read_region requests anyway.

There is also one very concrete way in which the current implementation is definitely broken. We try to draw on the pillow based on the coordinates of the entire LabeledImageServer, without considering that we may actually be in a region which is offset from the origin.

Fair enough, that does indeed sound broken.

[alanocallaghan]

Because in the first you’re requesting that the labels are generated at a downsample of 2, then you are (weirdly) requesting them to be upsampled by a factor of 2. It’s like you’ve exported the image with downsample=2 and then zoomed in.

Okay I've got you backed into a corner now Pete 😆 because that's not what the code will do. The labeled server has one downsample level, which is 1, so I'm actually retrieving the region at a downsample of 2 with respect to the original image. There'll be no transform happening in read_region, just requesting the labeled image directly. Analogously to a regular ImageServer, I'd be just reading the pixels at the one available level --- the one that matches a downsample of 1.

In [2]: labeled_server.metadata.downsamples
   ...: 
Out[2]: (1,)

[petebankhead]

Fair, that sounds like a bug :)

I haven’t actually run any of the code, I’m really just describing what I think/hope it does. I maintain that the two different scenarios should do something different… even if it’s not the different thing they currently do. And that the original proposed fix could open up a whole new set of problems.

I remember encountering this in QuPath with PyramidGeneratingImageServer.

I think the approach I eventually took was to always assume downsample=1 in readRegion requests, because that was the only way to easily match regions across images that correspond but are at different resolutions - which is what you normally have to do when exporting an image and its labels.

Can’t 100% guarantee that’s what happens though, as I’m not currently using a proper computer and haven’t checked anything.

As mentioned at #26 (comment) the entire readRegion / read_region idea has some major flaws… so I really hope we’ll be able to avoid users having to deal with it directly very much. And after QuPath v1.0 we can replace it entirely…

[alanocallaghan]

Yes, the design sounds sensible; the idea behind adding the failing tests was to show that our design isn't really working as stated, or at least not in an intuitive/transparent way. Ideally we resolve any confusion like this behind the scenes, properly test it, and present a nice streamlined user interface.

Also worth noting that qubalab will still be here on Monday regardless of whether or not you spend annual leave time responding to github threads :-)

[Rylern]

labeled_server = LabeledImageServer(qupath_server.metadata, annotations, downsample=2)
labeled_server.read_region(downsample=1, x, y, width, height)

Because in the first you’re requesting that the labels are generated at a downsample of 2, then you are (weirdly) requesting them to be upsampled by a factor of 2. It’s like you’ve exported the image with downsample=2 and then zoomed in.

I think I see things differently. When you create labeled_server, you get a new ImageServer that doesn't have to exactly correspond to qupath_server. They are two differents servers, so labeled_server can represent an image smaller than qupath_server if a downsample > 1 is specified (in that case, labeled_server.metadata.width < qupath_server.metadata.width).

So we have now a new image. When you call labeled_server.read_region with a specific downsample, this downsample corresponds to labeled_server:

• If the downsample is 1, then we get the full resolution of labeled_server (which has the size of qupath_server divided by 2 here).
• If the downsample is 2, then we get a downsampled view of labeled_server (which has the size of qupath_server divided by 4 here).

In short, I'm not bothered by the current behaviour, but it seems I'm the only one. Or maybe I don't clearly see the problem.

[petebankhead]

Hi @Rylern I think I used to see it the same way as you, and I expect that at some point that's also what QuPath did.... until it caused too much trouble.

The problem is that the main use case for a LabeledImageServer is for generating labels that correspond to an existing image. And then it helps if we can

1. ensure that the same x,y,width,height values can be used for the same region
2. encode the downsample in the LabeledImageServer - while returning the same values for getWidth() and getHeight() as in the original image

This gist shows how that works in QuPath. Basically, ImageServer.getDownsampleForLevel(0) does not have to return 1.0. It always does for any 'standalone' image, but does not necessarily if we generate a 'companion' image at a different resolution.

One place where this is alluded to in the QuPath docs is for RegionRequest.createInstance(server)

https://github.com/qupath/qupath/blob/5c4911628e9792fd1f7102a17f468e7b70de36ca/qupath-core/src/main/java/qupath/lib/regions/RegionRequest.java#L80-L88

This is the main reason I think we should make API changes here based on the limitations we find when we write sample notebooks illustrating specific use-cases. It's much easier to export images and labeled images that match if we can use the same RegionRequest for both, and we don't need to rescale coordinates and worry about rounding errors.

[Rylern]

OK. I always thought that the full resolution image had a downsample of 1, so it may be needed to update some things in qubalab (or at least document this). Should I do it or @alanocallaghan are you on it?

[alanocallaghan]

I am writing a short notebook or two to document things

[alanocallaghan]

There are some examples here:

demo-resolution.zip

I think #27 seems to resolve any confusing cases here