Store original image downsample in LabeledImageServer (#27)

Previously all labeledservers had downsample of 1, meaning they effectively lost any concept of the original image's resolution. This lead to an inconvenient API for the most common use cases (matched use of similar labelled and original images).

Now, running

ls = LabeledImageServer(original_server.metadata, features, downsample=2)
ls.read_region(downsample=1)

will return an image of the size of original_server, rather than downsampled by 2.

It also means that region requests are resized to the labeledserver's downsampled coordinate space, without needing to run region.downsample_region() on every request.

notebooks/working_with_objects.ipynb

@@ -199,12 +199,12 @@
     "\n",
     "To begin, let's use the GeoJSON representation to create masks and labeled images.\n",
     "\n",
-    "To create masks and labeled images, Qubalab has a `LabeledImageServer` class. This class is an implementation of the Qubalab `ImageServer` which is described in the *opening_images.ipynb* notebook, so it is recommended that you go through this notebook first. In short, `ImageServer` is a class to access metadata and pixel values of images.\n",
+    "To create masks and labeled images, Qubalab has a `LabeledImageServer` class. This class is an implementation of the Qubalab `ImageServer` class described in the *opening_images.ipynb* notebook, so it is recommended that you go through this notebook first. In short, `ImageServer` is a class to access metadata and pixel values of images.\n",
     "\n",
     "This server needs:\n",
     "- Some metadata representing the image containing the objects. Since we are working with the image that is opened in QuPath, we can read the metadata of the `QuPathServer`, as described in *communicating_with_qupath.ipynb*.\n",
     "- The objects to represent. We will give the annotations we've been working with.\n",
-    "- A downsample to apply to the image.\n",
+    "- A downsample to apply to the image features.\n",
     "\n",
     "Once the server is created, all functions described in *opening_images.ipynb* (such as `read_region()` to read the image) are also available."
    ]
@@ -225,10 +225,12 @@
     "# Set a downsample. The labeled image will be 20 times smaller than the image currently opened in QuPath\n",
     "downsample = 20\n",
     "\n",
-    "# Create the LabeledImageServer. This doesn't create labeled image yet\n",
+    "# Create the LabeledImageServer. This doesn't create labeled image yet, it just creates a downsampled version of the image features\n",
     "labeled_server = LabeledImageServer(qupath_server.metadata, annotations, downsample=downsample)\n",
     "\n",
     "# Request the pixel values of the entire labeled image. Pixel values will be created as they are requested \n",
+    "# note that when reading regions, LabeledImageServer considers the downsample relative to the original image, even if a downsample != 1 is provided on creation of the server\n",
+    "# this means that here, we receive a labeled image the size of the currently-open image, relating to downsampled image features (annotations)\n",
     "label_image = labeled_server.read_region()\n",
     "\n",
     "\n",

qubalab/images/image_server.py

@@ -294,6 +294,8 @@ def _resize(image: Union[np.ndarray, Image.Image], target_size: tuple[int, int],
                     pilImage = Image.fromarray(image)
                 elif np.issubdtype(image.dtype, np.integer):
                     pilImage = Image.fromarray(image.astype(np.int32), mode='I')
+                elif np.issubdtype(image.dtype, np.bool_):
+                    pilImage = Image.fromarray(image, "1")
                 else:
                     pilImage = Image.fromarray(image.astype(np.float32), mode='F')
                 pilImage = ImageServer._resize(pilImage, target_size=target_size, resample=resample)

qubalab/images/labeled_server.py

@@ -18,8 +18,9 @@ class LabeledImageServer(ImageServer):
     present on an image.
 
     The returned image will have one timepoint and one z-stack. The size of the remaining dimensions depend
-    on the parameters given during the server creation.
-    The image will only have one resolution level.
+    on the metadata provided when creating the server --- usually, the same as the ImageServer that the labeled image corresponds to.
+
+    The image will only have one resolution level; the downsample for this level may be greater than or less than 1, and consequently region requests and downsamples should be considered relative to the metadata provided at server creation, **not** relative to the downsampled (or upsampled) LabeledImageServer coordinates.
     """
 
     def __init__(
@@ -29,97 +30,127 @@ def __init__(
         label_map: dict[Classification, int] = None,
         downsample: float = None,
         multichannel: bool = False,
-        **kwargs
+        resize_method=PIL.Image.Resampling.NEAREST,
+        **kwargs,
     ):
         """
         :param base_image_metadata: the metadata of the image containing the image features
         :param features: the image features to draw
         :param label_map: a dictionary mapping a classification to a label. The value of pixels where an image feature with
                           a certain classification is present will be taken from this dictionnary. If not provided, each feature
-                          will be assigned a unique integer. All labels must be greater than 0 
+                          will be assigned a unique integer. All labels must be greater than 0
         :param downsample: the downsample to apply to the image. Can be omitted to use the full resolution image
         :param multichannel: if False, the image returned by this server will have a single channel where pixel values will be unsigned
                              integers representing a label (see the label_map parameter). If True, the number of channels will be
                              equal to the highest label value + 1, and the pixel located at (c, y, x) is a boolean indicating if an annotation
                              with label c is present on the pixel located at (x, y)
-        :param resize_method: the resampling method to use when resizing the image for downsampling. Bicubic by default
+        :param resize_method: the resampling method to use when resizing the image for downsampling. Nearest neighbour by default for labeled images.
         :raises ValueError: when a label in label_map is less than or equal to 0
         """
-        super().__init__(**kwargs)
+        super().__init__(resize_method=resize_method, **kwargs)
 
         if label_map is not None and any(label <= 0 for label in label_map.values()):
-            raise ValueError('A label in label_map is less than or equal to 0: ' + str(label_map))
+            raise ValueError(
+                "A label in label_map is less than or equal to 0: " + str(label_map)
+            )
 
         self._base_image_metadata = base_image_metadata
         self._downsample = 1 if downsample is None else downsample
         self._multichannel = multichannel
-        self._features = [f for f in features if label_map is None or f.classification in label_map]
-        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]
+        self._features = [
+            f for f in features if label_map is None or f.classification in label_map
+        ]
+        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
+        ]
         self._tree = shapely.STRtree(self._geometries)
 
         if label_map is None:
-            self._feature_index_to_label = {i: i+1 for i in range(len(self._features))}
+            self._feature_index_to_label = {
+                i: i + 1 for i in range(len(self._features))
+            }
         else:
-            self._feature_index_to_label = {i: label_map[self._features[i].classification] for i in range(len(self._features))}
+            self._feature_index_to_label = {
+                i: label_map[self._features[i].classification]
+                for i in range(len(self._features))
+            }
 
     def close(self):
         pass
 
     def _build_metadata(self) -> ImageMetadata:
         return ImageMetadata(
             self._base_image_metadata.path,
-            f'{self._base_image_metadata.name} - labels',
-            (ImageShape(
-                int(self._base_image_metadata.width / self._downsample),
-                int(self._base_image_metadata.height / self._downsample),
-                1,
-                max(self._feature_index_to_label.values(), default=0)+1 if self._multichannel else 1,
-                1,
-            ),),
+            f"{self._base_image_metadata.name} - labels",
+            (
+                ImageShape(
+                    int(self._base_image_metadata.width),
+                    int(self._base_image_metadata.height),
+                    1,
+                    max(self._feature_index_to_label.values(), default=0) + 1
+                    if self._multichannel
+                    else 1,
+                    1,
+                ),
+            ),
             PixelCalibration(
                 PixelLength(
-                    self._base_image_metadata.pixel_calibration.length_x.length * self._downsample,
-                    self._base_image_metadata.pixel_calibration.length_x.unit
+                    self._base_image_metadata.pixel_calibration.length_x.length,
+                    self._base_image_metadata.pixel_calibration.length_x.unit,
                 ),
                 PixelLength(
-                    self._base_image_metadata.pixel_calibration.length_y.length * self._downsample,
-                    self._base_image_metadata.pixel_calibration.length_y.unit
+                    self._base_image_metadata.pixel_calibration.length_y.length,
+                    self._base_image_metadata.pixel_calibration.length_y.unit,
                 ),
-                self._base_image_metadata.pixel_calibration.length_z
+                self._base_image_metadata.pixel_calibration.length_z,
             ),
             False,
-            bool if self._multichannel else np.uint32
+            bool if self._multichannel else np.uint32,
+            downsamples=[self._downsample],
         )
 
     def _read_block(self, level: int, region: Region2D) -> np.ndarray:
         if self._multichannel:
-            full_image = np.zeros((self.metadata.n_channels, region.height, region.width), dtype=self.metadata.dtype)
+            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))
+                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,
-                            shapely.affinity.translate(self._geometries[i], -region.x, -region.y),
-                            1
+                            shapely.affinity.translate(
+                                self._geometries[i], -region.x, -region.y
+                            ),
+                            1,
                         )
                 full_image[label, :, :] = np.asarray(image, dtype=self.metadata.dtype)
 
             return full_image
         else:
-            image = PIL.Image.new('I', (region.width, region.height))
+            image = PIL.Image.new("I", (region.width, region.height))
             drawing_context = PIL.ImageDraw.Draw(image)
             for i in self._tree.query(region.geometry):
                 draw_geometry(
                     image.size,
                     drawing_context,
-                    shapely.affinity.translate(self._geometries[i], -region.x, -region.y),
-                    self._feature_index_to_label[i]
+                    shapely.affinity.translate(
+                        self._geometries[i], -region.x, -region.y
+                    ),
+                    self._feature_index_to_label[i],
                 )
             return np.expand_dims(np.asarray(image, dtype=self.metadata.dtype), axis=0)

tests/images/test_labeled_server.py

@@ -37,7 +37,7 @@
 
 def test_image_width_with_downsample():
     downsample = 1.5
-    expected_width = int(sample_metadata.shape.x / downsample)
+    expected_width = sample_metadata.shape.x
     labeled_server = LabeledImageServer(sample_metadata, [], downsample=downsample)
 
     width = labeled_server.metadata.width
@@ -49,7 +49,7 @@ def test_image_width_with_downsample():
 
 def test_image_height_with_downsample():
     downsample = 1.5
-    expected_height = int(sample_metadata.shape.y / downsample)
+    expected_height = sample_metadata.shape.y
     labeled_server = LabeledImageServer(sample_metadata, [], downsample=downsample)
 
     height = labeled_server.metadata.height
@@ -174,7 +174,7 @@ def test_image_n_resolutions():
 
 def test_x_pixel_length_with_downsample():
     downsample = 1.5
-    expected_length_x = sample_metadata.pixel_calibration.length_x.length * downsample
+    expected_length_x = sample_metadata.pixel_calibration.length_x.length
     labeled_server = LabeledImageServer(sample_metadata, [], downsample=downsample)
 
     length_x = labeled_server.metadata.pixel_calibration.length_x.length
@@ -327,7 +327,7 @@ def test_read_points_in_single_channel_image_without_label_map_with_downsample()
     )
     labeled_server = LabeledImageServer(sample_metadata, features, multichannel=False, downsample=downsample)
 
-    image = labeled_server.read_region(1, Region2D(0, 0, labeled_server.metadata.width, labeled_server.metadata.height))
+    image = labeled_server.read_region(downsample, Region2D(0, 0, labeled_server.metadata.width, labeled_server.metadata.height))
 
     np.testing.assert_array_equal(image, expected_image)
 
@@ -359,7 +359,7 @@ def test_read_line_in_single_channel_image_without_label_map_with_downsample():
     )
     labeled_server = LabeledImageServer(sample_metadata, features, multichannel=False, downsample=downsample)
 
-    image = labeled_server.read_region(1, Region2D(0, 0, labeled_server.metadata.width, labeled_server.metadata.height))
+    image = labeled_server.read_region(downsample, Region2D(0, 0, labeled_server.metadata.width, labeled_server.metadata.height))
 
     np.testing.assert_array_equal(image, expected_image)
 
@@ -391,7 +391,7 @@ def test_read_polygon_in_single_channel_image_without_label_map_with_downsample(
     )
     labeled_server = LabeledImageServer(sample_metadata, features, multichannel=False, downsample=downsample)
 
-    image = labeled_server.read_region(1, Region2D(0, 0, labeled_server.metadata.width, labeled_server.metadata.height))
+    image = labeled_server.read_region(downsample, Region2D(0, 0, labeled_server.metadata.width, labeled_server.metadata.height))
 
     np.testing.assert_array_equal(image, expected_image)
 
@@ -409,10 +409,7 @@ def rands():
             (x, y + 1)
         )
 
-    coords = [rands() for i in range(max_objects)]
-
-    n_objects = len(coords)
-    features = [ImageFeature(geojson.Polygon([coords[i]]), Classification("Some classification")) for i in range(n_objects)]
+    features = [ImageFeature(geojson.Polygon([rands()])) for i in range(max_objects)]
     labeled_server = LabeledImageServer(large_metadata, features, multichannel=False, downsample=downsample)
 
     image = labeled_server.read_region(1, Region2D(0, 0, labeled_server.metadata.width, labeled_server.metadata.height))
@@ -434,8 +431,41 @@ def test_single_channel_labeled_image_with_region_request():
 
     np.testing.assert_array_equal(image, expected_image)
 
+
+
+def test_single_channel_labeled_image_with_starting_downsample():
+    features = [ImageFeature(geojson.LineString([(6, 5), (9, 5)]))]
+    # when resizing, we lose the labels with bicubic
+    expected_image = np.array(
+        [[[0, 0, 0, 0, 0],
+          [0, 0, 0, 0, 0],
+          [0, 0, 0, 1, 1]]]
+    )
+    labeled_server = LabeledImageServer(sample_metadata, features, multichannel=False, downsample=1)
+    downsample = 2
+    region = Region2D(0, 0, labeled_server.metadata.width, labeled_server.metadata.height)
+    image = labeled_server.read_region(downsample, region)
+
+    np.testing.assert_array_equal(image, expected_image)
+
+
+def test_single_channel_labeled_image_with_request_downsample():
+    # we downsample 
+    features = [ImageFeature(geojson.LineString([(6, 5), (9, 5)]))]
+    expected_image = np.array(
+        [[[0, 0, 0, 0, 0],
+          [0, 0, 0, 0, 0],
+          [0, 0, 0, 1, 1]]]
+    )
+    labeled_server = LabeledImageServer(sample_metadata, features, multichannel=False, downsample=1)
+    region = Region2D(0, 0, labeled_server.metadata.width, labeled_server.metadata.height)
+    image = labeled_server.read_region(2, region)
+
+    np.testing.assert_array_equal(image, expected_image)
+
+
+
 def test_multi_channel_labeled_image_with_region_request():
-    downsample = 1
     features = [ImageFeature(geojson.LineString([(7, 5), (9, 5)]))]
     expected_image = np.array(
         [[[False, False, False, False, False],
@@ -445,7 +475,66 @@ def test_multi_channel_labeled_image_with_region_request():
           [False, False, False, False, False],
           [False, False, True, True, True]]]
     )
-    labeled_server = LabeledImageServer(sample_metadata, features, multichannel=True, downsample=downsample)
+    labeled_server = LabeledImageServer(sample_metadata, features, multichannel=True, downsample=1)
     region = Region2D(5, 3, labeled_server.metadata.width-5, labeled_server.metadata.height-3)
     image = labeled_server.read_region(1, region)
+
+    np.testing.assert_array_equal(image, expected_image)
+
+
+
+
+def test_multi_channel_labeled_image_with_starting_downsample():
+    # we downsample the feature, then request at the same downsample
+    features = [ImageFeature(geojson.LineString([(6, 5), (9, 5)]))]
+    expected_image = np.array(
+        [[[False, False, False, False, False],
+          [False, False, False, False, False],
+          [False, False, False, False, False]],
+          [[False, False, False, False, False],
+          [False, False, False, False, False],
+          [False, False, False, True, True]]]
+    )
+    downsample = 2
+    labeled_server = LabeledImageServer(sample_metadata, features, multichannel=True, downsample=downsample)
+    region = Region2D(0, 0, sample_metadata.width, sample_metadata.height)
+    image = labeled_server.read_region(2, region)
+
+    np.testing.assert_array_equal(image, expected_image)
+
+def test_multi_channel_labeled_image_with_request_downsample():
+    features = [ImageFeature(geojson.LineString([(6, 5), (9, 5)]))]
+    ## because we resize the image after reading, we lose the small region
+    expected_image = np.array(
+        [[[False, False, False, False, False],
+          [False, False, False, False, False],
+          [False, False, False, False, False]],
+          [[False, False, False, False, False],
+          [False, False, False, False, False],
+          [False, False, False, False, False]]]
+    )
+    labeled_server = LabeledImageServer(sample_metadata, features, multichannel=True, downsample=1)
+    downsample = 2
+    region = Region2D(0, 0, labeled_server.metadata.width, labeled_server.metadata.height)
+    image = labeled_server.read_region(downsample, region)
+
+    np.testing.assert_array_equal(image, expected_image)
+
+
+def test_multi_channel_labeled_image_with_starting_downsample_upsampled():
+    # we downsample the feature, then request at a downsample of 1, so upsampled!
+    # therefore the feature gets much bigger
+    features = [ImageFeature(geojson.LineString([(5, 5), (9, 5)]))]
+    expected_image = np.array(
+        [[[False, False, False, False, False],
+        [False, False, False, False, False],
+        [False, False, False, False, False]],
+
+       [[False, False, False, False, False],
+        [False, False, False, False, False],
+        [False, False,  True,  True,  True]]]
+    )
+    labeled_server = LabeledImageServer(sample_metadata, features, multichannel=True, downsample=2)
+    image = labeled_server.read_region(2)
+
     np.testing.assert_array_equal(image, expected_image)