Merge pull request #1427 from girder/tiff-optimizations

Speed up reading tiff tiles slightly

CHANGELOG.md

@@ -6,6 +6,7 @@
 - Read config values from the environment variables ([#1422](../../pull/1422))
 - Optimizing when reading arrays rather than images from tiff files ([#1423](../../pull/1423))
 - Better filter DICOM adjacent files to ensure they share series instance IDs ([#1424](../../pull/1424))
+- Optimizing small getRegion calls and some tiff tile fetches ([#1427](../../pull/1427)
 
 ## 1.27.0
 

large_image/tilesource/base.py

@@ -24,11 +24,12 @@
                          TileOutputPILFormat, dtypeToGValue)
 from .jupyter import IPyLeafletMixin
 from .tiledict import LazyTileDict
-from .utilities import (JSONDict, _encodeImage,  # noqa: F401
-                        _encodeImageBinary, _gdalParameters, _imageToNumpy,
-                        _imageToPIL, _letterboxImage, _makeSameChannelDepth,
-                        _vipsCast, _vipsParameters, dictToEtree, etreeToDict,
-                        getPaletteColors, histogramThreshold, nearPowerOfTwo)
+from .utilities import (JSONDict, _addSubimageToImage,  # noqa: F401
+                        _encodeImage, _encodeImageBinary, _gdalParameters,
+                        _imageToNumpy, _imageToPIL, _letterboxImage,
+                        _makeSameChannelDepth, _vipsCast, _vipsParameters,
+                        dictToEtree, etreeToDict, getPaletteColors,
+                        histogramThreshold, nearPowerOfTwo)
 
 
 class TileSource(IPyLeafletMixin):
@@ -513,7 +514,7 @@ def _getRegionBounds(self, metadata, left=None, top=None, right=None,
 
         return left, top, right, bottom
 
-    def _tileIteratorInfo(self, **kwargs):
+    def _tileIteratorInfo(self, **kwargs):  # noqa
         """
         Get information necessary to construct a tile iterator.
           If one of width or height is specified, the other is determined by
@@ -589,6 +590,14 @@ def _tileIteratorInfo(self, **kwargs):
             :edges: if True, then the edge tiles will exclude the overlap
                 distance.  If unset or False, the edge tiles are full size.
 
+        :param tile_offset: if present, adjust tile positions so that the
+            corner of one tile is at the specified location.
+
+            :left: the left offset in pixels.
+            :top: the top offset in pixels.
+            :auto: a boolean, if True, automatically set the offset to align
+                with the region's left and top.
+
         :param kwargs: optional arguments.  Some options are encoding,
             jpegQuality, jpegSubsampling, tiffCompression, frame.
         :returns: a dictionary of information needed for the tile iterator.
@@ -734,17 +743,26 @@ def _tileIteratorInfo(self, **kwargs):
             tile_overlap['x'] = int(math.ceil(tile_overlap['x'] * requestedScale))
             tile_overlap['y'] = int(math.ceil(tile_overlap['y'] * requestedScale))
 
+        offset_x = kwargs.get('tile_offset', {}).get('left', 0)
+        offset_y = kwargs.get('tile_offset', {}).get('top', 0)
+        if kwargs.get('tile_offset', {}).get('auto'):
+            offset_x = left
+            offset_y = top
+        offset_x = (left - left % tile_size['width']) if offset_x > left else offset_x
+        offset_y = (top - top % tile_size['height']) if offset_y > top else offset_y
         # If the overlapped tiles don't run over the edge, then the functional
         # size of the region is reduced by the overlap.  This factor is stored
         # in the overlap offset_*.
-        xmin = int(left / tile_size['width'])
-        xmax = max(int(math.ceil((float(right) - tile_overlap['range_x']) /
+        xmin = int((left - offset_x) / tile_size['width'])
+        xmax = max(int(math.ceil((float(right - offset_x) - tile_overlap['range_x']) /
                                  tile_size['width'])), xmin + 1)
-        ymin = int(top / tile_size['height'])
-        ymax = max(int(math.ceil((float(bottom) - tile_overlap['range_y']) /
+        ymin = int((top - offset_y) / tile_size['height'])
+        ymax = max(int(math.ceil((float(bottom - offset_y) - tile_overlap['range_y']) /
                                  tile_size['height'])), ymin + 1)
         tile_overlap.update({'xmin': xmin, 'xmax': xmax,
                              'ymin': ymin, 'ymax': ymax})
+        tile_overlap['offset_x'] += offset_x
+        tile_overlap['offset_y'] += offset_y
 
         # Use RGB for JPEG, RGBA for PNG
         mode = 'RGBA' if kwargs.get('encoding') in {'PNG', 'TIFF', 'TILED'} else 'RGB'
@@ -2162,6 +2180,13 @@ def getRegion(self, format=(TILE_FORMAT_IMAGE, ), **kwargs):
         if 'tile_position' in kwargs:
             kwargs = kwargs.copy()
             kwargs.pop('tile_position', None)
+        tiled = TILE_FORMAT_IMAGE in format and kwargs.get('encoding') == 'TILED'
+        if not tiled and 'tile_offset' not in kwargs and 'tile_size' not in kwargs:
+            kwargs = kwargs.copy()
+            kwargs['tile_size'] = {
+                'width': max(self.tileWidth, 4096),
+                'height': max(self.tileHeight, 4096)}
+            kwargs['tile_offset'] = {'auto': True}
         iterInfo = self._tileIteratorInfo(**kwargs)
         if iterInfo is None:
             image = PIL.Image.new('RGB', (0, 0))
@@ -2173,26 +2198,20 @@ def getRegion(self, format=(TILE_FORMAT_IMAGE, ), **kwargs):
         mode = None if TILE_FORMAT_NUMPY in format else iterInfo['mode']
         outWidth = iterInfo['output']['width']
         outHeight = iterInfo['output']['height']
-        tiled = TILE_FORMAT_IMAGE in format and kwargs.get('encoding') == 'TILED'
         image = None
         for tile in self._tileIterator(iterInfo):
             # Add each tile to the image
             subimage, _ = _imageToNumpy(tile['tile'])
             x0, y0 = tile['x'] - left, tile['y'] - top
-            if x0 < 0:
-                subimage = subimage[:, -x0:]
-                x0 = 0
-            if y0 < 0:
-                subimage = subimage[-y0:, :]
-                y0 = 0
-            subimage = subimage[:min(subimage.shape[0], regionHeight - y0),
-                                :min(subimage.shape[1], regionWidth - x0)]
             if tiled:
                 image = self._addRegionTileToTiled(
                     image, subimage, x0, y0, regionWidth, regionHeight, tile, **kwargs)
             else:
-                image = self._addRegionTileToImage(
-                    image, subimage, x0, y0, regionWidth, regionHeight, **kwargs)
+                image = _addSubimageToImage(
+                    image, subimage, x0, y0, regionWidth, regionHeight)
+            # Somehow discarding the tile here speeds things up.
+            del tile
+            del subimage
         # Scale if we need to
         outWidth = int(math.floor(outWidth))
         outHeight = int(math.floor(outHeight))
@@ -2213,35 +2232,6 @@ def getRegion(self, format=(TILE_FORMAT_IMAGE, ), **kwargs):
             image = _letterboxImage(_imageToPIL(image, mode), maxWidth, maxHeight, kwargs['fill'])
         return _encodeImage(image, format=format, **kwargs)
 
-    def _addRegionTileToImage(self, image, subimage, x, y, width, height, **kwargs):
-        """
-        Add a subtile to a larger image.
-
-        :param image: the output image record.  None for not created yet.
-        :param subimage: a numpy array with the sub-image to add.
-        :param x: the location of the upper left point of the sub-image within
-            the output image.
-        :param y: the location of the upper left point of the sub-image within
-            the output image.
-        :param width: the output image size.
-        :param height: the output image size.
-        :returns: the output image record.
-        """
-        if image is None:
-            if (x, y, width, height) == (0, 0, subimage.shape[1], subimage.shape[0]):
-                return subimage
-            try:
-                image = np.zeros(
-                    (height, width, subimage.shape[2]),
-                    dtype=subimage.dtype)
-            except MemoryError:
-                raise exceptions.TileSourceError(
-                    'Insufficient memory to get region of %d x %d pixels.' % (
-                        width, height))
-        image, subimage = _makeSameChannelDepth(image, subimage)
-        image[y:y + subimage.shape[0], x:x + subimage.shape[1], :] = subimage
-        return image
-
     def _vipsAddAlphaBand(self, vimg, *otherImages):
         """
         Add an alpha band to a vips image.  The alpha value is either 1, 255,
@@ -2477,8 +2467,8 @@ def tileFrames(self, format=(TILE_FORMAT_IMAGE, ), frameList=None,
                 image = self._addRegionTileToTiled(
                     image, subimage, offsetX, offsetY, outWidth, outHeight, tile, **kwargs)
             else:
-                image = self._addRegionTileToImage(
-                    image, subimage, offsetX, offsetY, outWidth, outHeight, **kwargs)
+                image = _addSubimageToImage(
+                    image, subimage, offsetX, offsetY, outWidth, outHeight)
         if tiled:
             return self._encodeTiledImage(image, outWidth, outHeight, iterInfo, **kwargs)
         return _encodeImage(image, format=format, **kwargs)
@@ -2751,6 +2741,14 @@ def tileIterator(self, format=(TILE_FORMAT_NUMPY, ), resample=True,
                 returned are: 012, 0123, 01234, 12345, 23456, 34567, 4567, 567,
                 with the non-overlapped area of each as 0, 1, 2, 3, 4, 5, 6, 7.
 
+        :param tile_offset: if present, adjust tile positions so that the
+            corner of one tile is at the specified location.
+
+            :left: the left offset in pixels.
+            :top: the top offset in pixels.
+            :auto: a boolean, if True, automatically set the offset to align
+                with the region's left and top.
+
         :param encoding: if format includes TILE_FORMAT_IMAGE, a valid PIL
             encoding (typically 'PNG', 'JPEG', or 'TIFF') or 'TILED' (identical
             to TIFF).  Must also be in the TileOutputMimeTypes map.

large_image/tilesource/tiledict.py

@@ -122,15 +122,14 @@ def _retileTile(self):
         xmax = int((self['x'] + self.width - 1) // self.metadata['tileWidth'] + 1)
         ymin = int(max(0, self['y'] // self.metadata['tileHeight']))
         ymax = int((self['y'] + self.height - 1) // self.metadata['tileHeight'] + 1)
-        for x in range(xmin, xmax):
-            for y in range(ymin, ymax):
+        for y in range(ymin, ymax):
+            for x in range(xmin, xmax):
                 tileData = self.source.getTile(
                     x, y, self.level,
                     numpyAllowed='always', sparseFallback=True, frame=self.frame)
                 tileData, _ = _imageToNumpy(tileData)
                 if retile is None:
-                    retile = np.zeros(
-                        (self.height, self.width) if len(tileData.shape) == 2 else
+                    retile = np.empty(
                         (self.height, self.width, tileData.shape[2]),
                         dtype=tileData.dtype)
                 x0 = int(x * self.metadata['tileWidth'] - self['x'])

large_image/tilesource/utilities.py

@@ -763,6 +763,32 @@ def _makeSameChannelDepth(arr1, arr2):
     return arrays['arr1'], arrays['arr2']
 
 
+def _addSubimageToImage(image, subimage, x, y, width, height):
+    """
+    Add a subimage to a larger image as numpy arrays.
+
+    :param image: the output image record.  None for not created yet.
+    :param subimage: a numpy array with the sub-image to add.
+    :param x: the location of the upper left point of the sub-image within
+        the output image.
+    :param y: the location of the upper left point of the sub-image within
+        the output image.
+    :param width: the output image size.
+    :param height: the output image size.
+    :returns: the output image record.
+    """
+    if image is None:
+        if (x, y, width, height) == (0, 0, subimage.shape[1], subimage.shape[0]):
+            return subimage
+        image = np.empty(
+            (height, width, subimage.shape[2]),
+            dtype=subimage.dtype)
+    elif len(image.shape) != len(subimage.shape) or image.shape[-1] != subimage.shape[-1]:
+        image, subimage = _makeSameChannelDepth(image, subimage)
+    image[y:y + subimage.shape[0], x:x + subimage.shape[1]] = subimage
+    return image
+
+
 def _computeFramesPerTexture(opts, numFrames, sizeX, sizeY):
     """
     Compute the number of frames for each tile_frames texture.

sources/tiff/large_image_source_tiff/tiff_reader.py

@@ -53,6 +53,21 @@
 libtiff_ctypes.suppress_errors()
 
 
+_ctypesFormattbl = {
+    (8, libtiff_ctypes.SAMPLEFORMAT_UINT): np.uint8,
+    (8, libtiff_ctypes.SAMPLEFORMAT_INT): np.int8,
+    (16, libtiff_ctypes.SAMPLEFORMAT_UINT): np.uint16,
+    (16, libtiff_ctypes.SAMPLEFORMAT_INT): np.int16,
+    (16, libtiff_ctypes.SAMPLEFORMAT_IEEEFP): np.float16,
+    (32, libtiff_ctypes.SAMPLEFORMAT_UINT): np.uint32,
+    (32, libtiff_ctypes.SAMPLEFORMAT_INT): np.int32,
+    (32, libtiff_ctypes.SAMPLEFORMAT_IEEEFP): np.float32,
+    (64, libtiff_ctypes.SAMPLEFORMAT_UINT): np.uint64,
+    (64, libtiff_ctypes.SAMPLEFORMAT_INT): np.int64,
+    (64, libtiff_ctypes.SAMPLEFORMAT_IEEEFP): np.float64,
+}
+
+
 def patchLibtiff():
     libtiff_ctypes.libtiff.TIFFFieldWithTag.restype = \
         ctypes.POINTER(libtiff_ctypes.TIFFFieldInfo)
@@ -286,6 +301,7 @@ def _loadMetadata(self):
         self._tileHeight = info.get('tilelength') or info.get('rowsperstrip')
         self._imageWidth = info.get('imagewidth')
         self._imageHeight = info.get('imagelength')
+        self._tilesAcross = (self._imageWidth + self._tileWidth - 1) // self._tileWidth
         if not info.get('tilelength'):
             self._stripsPerTile = int(max(1, math.ceil(256.0 / self._tileHeight)))
             self._stripHeight = self._tileHeight
@@ -404,12 +420,11 @@ def _toTileNum(self, x, y, transpose=False):
         #     raise InvalidOperationTiffError(
         #         'Tile x=%d, y=%d does not exist' % (x, y))
         if self._tiffInfo.get('istiled'):
-            tileNum = libtiff_ctypes.libtiff.TIFFComputeTile(
-                self._tiffFile, pixelX, pixelY, 0, 0).value
+            tileNum = pixelX // self._tileWidth + (pixelY // self._tileHeight) * self._tilesAcross
         else:
             # TIFFComputeStrip with sample=0 is just the row divided by the
             # strip height
-            tileNum = int(pixelY // self._stripHeight)
+            tileNum = pixelY // self._stripHeight
         return tileNum
 
     @methodcache(key=partial(strhash, '_getTileByteCountsType'))
@@ -569,34 +584,39 @@ def _getUncompressedTile(self, tileNum):
         :rtype: PIL.Image
         :raises: IOTiffError
         """
-        with self._tileLock:
-            if self._tiffInfo.get('istiled'):
-                tileSize = libtiff_ctypes.libtiff.TIFFTileSize(self._tiffFile).value
-            else:
+        if self._tiffInfo.get('istiled'):
+            if not hasattr(self, '_uncompressedTileSize'):
+                with self._tileLock:
+                    self._uncompressedTileSize = libtiff_ctypes.libtiff.TIFFTileSize(
+                        self._tiffFile).value
+            tileSize = self._uncompressedTileSize
+        else:
+            with self._tileLock:
                 stripSize = libtiff_ctypes.libtiff.TIFFStripSize(
                     self._tiffFile).value
-                stripsCount = min(self._stripsPerTile, self._stripCount - tileNum)
-                tileSize = stripSize * self._stripsPerTile
+            stripsCount = min(self._stripsPerTile, self._stripCount - tileNum)
+            tileSize = stripSize * self._stripsPerTile
         imageBuffer = ctypes.create_string_buffer(tileSize)
-        with self._tileLock:
-            if self._tiffInfo.get('istiled'):
+        if self._tiffInfo.get('istiled'):
+            with self._tileLock:
                 readSize = libtiff_ctypes.libtiff.TIFFReadEncodedTile(
                     self._tiffFile, tileNum, imageBuffer, tileSize)
-            else:
-                readSize = 0
-                for stripNum in range(stripsCount):
+        else:
+            readSize = 0
+            for stripNum in range(stripsCount):
+                with self._tileLock:
                     chunkSize = libtiff_ctypes.libtiff.TIFFReadEncodedStrip(
                         self._tiffFile,
                         tileNum + stripNum,
                         ctypes.byref(imageBuffer, stripSize * stripNum),
                         stripSize).value
-                    if chunkSize <= 0:
-                        msg = 'Read an unexpected number of bytes from an encoded strip'
-                        raise IOTiffError(msg)
-                    readSize += chunkSize
-                if readSize < tileSize:
-                    ctypes.memset(ctypes.byref(imageBuffer, readSize), 0, tileSize - readSize)
-                    readSize = tileSize
+                if chunkSize <= 0:
+                    msg = 'Read an unexpected number of bytes from an encoded strip'
+                    raise IOTiffError(msg)
+                readSize += chunkSize
+            if readSize < tileSize:
+                ctypes.memset(ctypes.byref(imageBuffer, readSize), 0, tileSize - readSize)
+                readSize = tileSize
         if readSize < tileSize:
             raise IOTiffError(
                 'Read an unexpected number of bytes from an encoded tile' if readSize >= 0 else
@@ -612,23 +632,10 @@ def _getUncompressedTile(self, tileNum):
             self._tiffInfo.get('bitspersample'),
             self._tiffInfo.get('sampleformat') if self._tiffInfo.get(
                 'sampleformat') is not None else libtiff_ctypes.SAMPLEFORMAT_UINT)
-        formattbl = {
-            (8, libtiff_ctypes.SAMPLEFORMAT_UINT): np.uint8,
-            (8, libtiff_ctypes.SAMPLEFORMAT_INT): np.int8,
-            (16, libtiff_ctypes.SAMPLEFORMAT_UINT): np.uint16,
-            (16, libtiff_ctypes.SAMPLEFORMAT_INT): np.int16,
-            (16, libtiff_ctypes.SAMPLEFORMAT_IEEEFP): np.float16,
-            (32, libtiff_ctypes.SAMPLEFORMAT_UINT): np.uint32,
-            (32, libtiff_ctypes.SAMPLEFORMAT_INT): np.int32,
-            (32, libtiff_ctypes.SAMPLEFORMAT_IEEEFP): np.float32,
-            (64, libtiff_ctypes.SAMPLEFORMAT_UINT): np.uint64,
-            (64, libtiff_ctypes.SAMPLEFORMAT_INT): np.int64,
-            (64, libtiff_ctypes.SAMPLEFORMAT_IEEEFP): np.float64,
-        }
         image = np.ctypeslib.as_array(ctypes.cast(
             imageBuffer, ctypes.POINTER(ctypes.c_uint8)), (tileSize, )).view(
-                formattbl[format]).reshape(
-                    (th, tw, self._tiffInfo.get('samplesperpixel')))
+                _ctypesFormattbl[format]).reshape(
+                    (th, tw, self._tiffInfo['samplesperpixel']))
         if (self._tiffInfo.get('samplesperpixel') == 3 and
                 self._tiffInfo.get('photometric') == libtiff_ctypes.PHOTOMETRIC_YCBCR):
             if self._tiffInfo.get('bitspersample') == 16:

test/source_geo_base.py

@@ -358,7 +358,7 @@ def testGetRegionWithProjection(self):
         ts = self.basemodule.open(imagePath, projection='EPSG:3857')
         region, _ = ts.getRegion(output=dict(maxWidth=1024, maxHeight=1024),
                                  format=constants.TILE_FORMAT_NUMPY)
-        assert region.shape == (1024, 1024, 4)
+        assert region.shape == (1024, 1024, 3)
 
     def testGCPProjection(self):
         imagePath = datastore.fetch('region_gcp.tiff')

test/test_source_base.py

@@ -369,6 +369,19 @@ def testTileOverlapWithRegionOffset():
     assert firstTile['tile_overlap']['right'] == 200
 
 
+def testGetRegionAutoOffset():
+    imagePath = datastore.fetch('sample_image.ptif')
+    source = large_image.open(imagePath)
+    region1, _ = source.getRegion(
+        region=dict(left=20, top=40, width=400, height=500),
+        format=large_image.constants.TILE_FORMAT_NUMPY)
+    region2, _ = source.getRegion(
+        region=dict(left=20, top=40, width=400, height=500),
+        tile_size=dict(width=240, height=240),
+        format=large_image.constants.TILE_FORMAT_NUMPY)
+    assert np.all(region2 == region1)
+
+
 @pytest.mark.parametrize((
     'options', 'lensrc', 'lenquads', 'frame10', 'src0', 'srclast', 'quads10',
 ), [