Fix bioformats processing of MONOCHROME1 dicom data

CHANGELOG.md

@@ -15,6 +15,7 @@
 ### Changes
 
 - Require a minimum version of Pillow to avoid a CVE ([#1828](../../pull/1828))
+- Work around how bioformats handles DICOMs in Monochrome1 ([#1834](../../pull/1834))
 
 ## 1.31.1
 

sources/bioformats/large_image_source_bioformats/__init__.py

@@ -332,6 +332,7 @@
             delattr(self, '_lastGetTileException')
         except Exception as exc:
             raise TileSourceError('Bioformats cannot read a tile: %r' % exc)
+        self._checkForOffset()
         self._populatedLevels = len([
             v for v in self._metadata['frameSeries'][0]['series'] if v is not None])
 
@@ -348,6 +349,40 @@
                 if javabridge.get_env():
                     javabridge.detach()
 
+    def _checkForOffset(self):
+        """
+        The bioformats DICOM reader does unfortunate things to MONOCHROME1
+        16-bit images.  Store an offset to undo it, if appropriate.
+        """
+        if self._metadata.get('readerClassName') != 'loci.formats.in.DicomReader':
+            return
+        if self._metadata.get('seriesMetadata', {}).get(
+                '0028,0004 Photometric Interpretation') != 'MONOCHROME1':
+            return
+        if np.issubdtype(self.dtype, np.uint8):
+            self._fix_offset = 255
+            return
+        if not np.issubdtype(self.dtype, np.int16) and not np.issubdtype(self.dtype, '>i2'):
+            return
+        # This is bioformats behavior
+        try:
+            maxPixelRange = int(self._metadata['seriesMetadata'].get(
+                '0028,1051 Window Width', 0))
+        except Exception:
+            maxPixelRange = -1
+        try:
+            centerPixelValue = int(self._metadata['seriesMetadata'].get(
+                '0028,1050 Window Center', 0))
+        except Exception:
+            centerPixelValue = -1
+        maxPixelValue = maxPixelRange + (centerPixelValue // 2)
+        maxAllowRange = 2 ** int(self._metadata['seriesMetadata'].get(
+            '0028,0101 Bits Stored', 16)) - 1
+        if maxPixelRange == -1 or centerPixelValue < maxPixelRange // 2:
+            maxPixelValue = maxAllowRange
+        if maxPixelValue:
+            self._fix_offset = maxPixelValue
+
     def _metadataForCurrentSeries(self, rdr):
         self._metadata = getattr(self, '_metadata', {})
         self._metadata.update({
@@ -710,6 +745,8 @@
             retile[0:min(tile.shape[0], finalHeight), 0:min(tile.shape[1], finalWidth)] = tile[
                 0:min(tile.shape[0], finalHeight), 0:min(tile.shape[1], finalWidth)]
             tile = retile
+        if hasattr(self, '_fix_offset') and format == TILE_FORMAT_NUMPY:
+            tile = self._fix_offset - tile
         return self._outputTile(tile, format, x, y, z, pilImageAllowed, numpyAllowed, **kwargs)
 
     def getAssociatedImagesList(self):

test/datastore.py

@@ -129,6 +129,10 @@
     # --rgb --quality=0.015 --compression='JPEG-2000 Lossy' parameters to make
     # the file small
     'TCGA-55-8207-01Z-00-DX1.ome.tiff': 'sha512:50cf63f0e8bfa3054d3532b7dd0237b66aeb4c7609da874639a28bc068dbd157f786e84d3eb76a3b0e6636a042c56c3b96d3be2ad66f7589d0542a5d20cecdb4',  # noqa
+    # Extracted from a sample dicom file from issue #1823 on github
+    # This is a dicom with monochrome1 format data that bioformats incorrectly
+    # inverts and offsets
+    'monochrome1.dcm': 'sha512:a67c38a5e26aba31b68e40eec0f260acf0ba638f0bf7fd99d41c3e16ddd8fd43b2737aeb3759ce92f2c4a10a81995bdae3b93b101ebcaa59543ea6eff7a4c8f2',  # noqa
 }
 
 

test/test_source_base.py

@@ -33,7 +33,7 @@
     'deepzoom': {},
     'dicom': {
         'read': r'\.dcm$',
-        'noread': r'tcia.*\.dcm$',
+        'noread': r'(tcia.*|monochrome1)\.dcm$',
     },
     'dummy': {'any': True, 'skipTiles': r''},
     'gdal': {
@@ -63,7 +63,7 @@
     'openjpeg': {'read': r'\.(jp2)$'},
     'openslide': {
         'read': r'\.(ptif|svs|ndpi|tif.*|qptiff|dcm)$',
-        'noread': r'(oahu|DDX58_AXL|huron\.image2_jpeg2k|landcover_sample|d042-353\.crop|US_Geo\.|extraoverview|imagej|bad_axes|synthetic_untiled|indica|tcia.*dcm|multiplane.*ndpi)',  # noqa
+        'noread': r'(oahu|DDX58_AXL|huron\.image2_jpeg2k|landcover_sample|d042-353\.crop|US_Geo\.|extraoverview|imagej|bad_axes|synthetic_untiled|indica|tcia.*dcm|monochrome1.dcm|multiplane.*ndpi)',  # noqa
         'skip': r'nokeyframe\.ome\.tiff|TCGA-55.*\.ome\.tiff$',
         'skipTiles': r'one_layer_missing',
     },
@@ -83,14 +83,14 @@
         'skipTiles': r'(sample_image\.ptif|one_layer_missing_tiles)'},
     'tifffile': {
         'read': r'',
-        'noread': r'((\.(nc|nd2|yml|yaml|json|czi|png|jpg|jpeg|jp2|zarr\.db|zarr\.zip)|(nokeyframe\.ome\.tiff|XY01\.ome\.tif|level.*\.dcm|tcia.*dcm)$)' +  # noqa
+        'noread': r'((\.(nc|nd2|yml|yaml|json|czi|png|jpg|jpeg|jp2|zarr\.db|zarr\.zip)|(nokeyframe\.ome\.tiff|XY01\.ome\.tif|level.*\.dcm|tcia.*dcm|monochrome1.dcm)$)' +  # noqa
                   (r'|bad_axes' if sys.version_info < (3, 9) else '') +
                   r')',
         'skip': r'indica' if sys.version_info < (3, 9) else '^$',
     },
     'vips': {
         'read': r'',
-        'noread': r'(\.(nc|nd2|yml|yaml|json|czi|png|svs|scn|zarr\.db|zarr\.zip)|tcia.*dcm)$',
+        'noread': r'(\.(nc|nd2|yml|yaml|json|czi|png|svs|scn|zarr\.db|zarr\.zip)|tcia.*dcm|monochrome1.dcm)$',  # noqa
         'skipTiles': r'(sample_image\.ptif|one_layer_missing_tiles|JK-kidney_B-gal_H3_4C_1-500sec\.jp2|extraoverview|synthetic_untiled)',  # noqa
     },
     'zarr': {'read': r'\.(zarr|zgroup|zattrs|db|zarr\.zip)$'},

test/test_source_bioformats.py

@@ -1,3 +1,5 @@
+import large_image
+
 from . import utilities
 from .datastore import datastore
 
@@ -31,3 +33,12 @@ def testBioformatsJarVersion():
     import large_image_source_bioformats
 
     assert '.' in large_image_source_bioformats._getBioformatsVersion()
+
+
+def testBioformatsDicomMonochome1():
+    import large_image_source_bioformats
+
+    imagePath = datastore.fetch('monochrome1.dcm')
+    source = large_image_source_bioformats.open(imagePath)
+    img, _ = source.getRegion(format=large_image.constants.TILE_FORMAT_NUMPY)
+    assert img[255, 191, 0] == 618