diff --git a/CHANGELOG.md b/CHANGELOG.md index aadb92f77..1577f04c2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ - Add a convert endpoint to the Girder plugin (#578) - Added support for creating Aperio svs files (#580) - Added support for geospatial files with GCP (#588) +- Regions can be output as tiled tiffs with scale or geospatial metadata (#594) ### Improvements - More untiled tiff files are handles by the bioformats reader (#569) diff --git a/girder/girder_large_image/rest/tiles.py b/girder/girder_large_image/rest/tiles.py index e7bd2c123..54592fe46 100644 --- a/girder/girder_large_image/rest/tiles.py +++ b/girder/girder_large_image/rest/tiles.py @@ -18,6 +18,7 @@ import hashlib import math import os +import pathlib import re import urllib @@ -29,6 +30,7 @@ from girder.exceptions import RestException from girder.models.file import File from girder.models.item import Item +from girder.utility.progress import setResponseTimeLimit from large_image.constants import TileInputUnits from large_image.exceptions import TileGeneralException @@ -627,7 +629,6 @@ def getDZITile(self, item, level, xandy, params): @loadmodel(model='item', map={'itemId': 'item'}, level=AccessType.WRITE) def deleteTiles(self, item, params): deleted = self.imageItemModel.delete(item) - # TODO: a better response return { 'deleted': deleted } @@ -654,8 +655,11 @@ def deleteTiles(self, item, params): .param('frame', 'For multiframe images, the 0-based frame number. ' 'This is ignored on non-multiframe images.', required=False, dataType='int') - .param('encoding', 'Thumbnail output encoding', required=False, - enum=['JPEG', 'PNG', 'TIFF'], default='JPEG') + .param('encoding', 'Output image encoding. TILED generates a tiled ' + 'tiff without the upper limit on image size the other options ' + 'have. For geospatial sources, TILED will also have ' + 'appropriate tagging.', required=False, + enum=['JPEG', 'PNG', 'TIFF', 'TILED'], default='JPEG') .param('contentDisposition', 'Specify the Content-Disposition response ' 'header disposition-type value.', required=False, enum=['inline', 'attachment']) @@ -752,8 +756,11 @@ def getTilesThumbnail(self, item, params): .param('frame', 'For multiframe images, the 0-based frame number. ' 'This is ignored on non-multiframe images.', required=False, dataType='int') - .param('encoding', 'Output image encoding', required=False, - enum=['JPEG', 'PNG', 'TIFF'], default='JPEG') + .param('encoding', 'Output image encoding. TILED generates a tiled ' + 'tiff without the upper limit on image size the other options ' + 'have. For geospatial sources, TILED will also have ' + 'appropriate tagging.', required=False, + enum=['JPEG', 'PNG', 'TIFF', 'TILED'], default='JPEG') .param('jpegQuality', 'Quality used for generating JPEG images', required=False, dataType='int', default=95) .param('jpegSubsampling', 'Chroma subsampling used for generating ' @@ -809,6 +816,7 @@ def getTilesRegion(self, item, params): ('contentDisposition', str), ]) _handleETag('getTilesRegion', item, params) + setResponseTimeLimit(86400) try: regionData, regionMime = self.imageItemModel.getRegion( item, **params) @@ -819,6 +827,20 @@ def getTilesRegion(self, item, params): self._setContentDisposition( item, params.get('contentDisposition'), regionMime, 'region') setResponseHeader('Content-Type', regionMime) + if isinstance(regionData, pathlib.Path): + BUF_SIZE = 65536 + + def stream(): + try: + with regionData.open('rb') as f: + while True: + data = f.read(BUF_SIZE) + if not data: + break + yield data + finally: + regionData.unlink() + return stream setRawResponse() return regionData diff --git a/girder/test_girder/girder_utilities.py b/girder/test_girder/girder_utilities.py index 9870514bc..4f5740ed1 100644 --- a/girder/test_girder/girder_utilities.py +++ b/girder/test_girder/girder_utilities.py @@ -4,7 +4,7 @@ from girder.models.upload import Upload from test.datastore import datastore -from test.utilities import JFIFHeader, JPEGHeader, PNGHeader # noqa +from test.utilities import JFIFHeader, JPEGHeader, PNGHeader, TIFFHeader, BigTIFFHeader # noqa def namedFolder(user, folderName='Public'): diff --git a/girder/test_girder/test_tiles_rest.py b/girder/test_girder/test_tiles_rest.py index 1f8dcdd9f..43fc40fb7 100644 --- a/girder/test_girder/test_tiles_rest.py +++ b/girder/test_girder/test_tiles_rest.py @@ -783,6 +783,15 @@ def testRegions(server, admin, fsAssetstore): assert width == 500 assert height == 375 + # Get a tiled image + params = {'regionWidth': 1000, 'regionHeight': 1000, + 'left': 48000, 'top': 3000, 'encoding': 'TILED'} + resp = server.request(path='/item/%s/tiles/region' % itemId, + user=admin, isJson=False, params=params) + assert utilities.respStatus(resp) == 200 + image = origImage = utilities.getBody(resp, text=False) + assert image[:len(utilities.BigTIFFHeader)] == utilities.BigTIFFHeader + @pytest.mark.usefixtures('unbindLargeImage') @pytest.mark.plugin('large_image') diff --git a/large_image/constants.py b/large_image/constants.py index fa3277669..a0cf99472 100644 --- a/large_image/constants.py +++ b/large_image/constants.py @@ -42,12 +42,13 @@ class SourcePriority: 'JPEG': 'image/jpeg', 'PNG': 'image/png', 'TIFF': 'image/tiff', + # TILED indicates the region output should be generated as a tiled TIFF + 'TILED': 'image/tiff', } TileOutputPILFormat = { 'JFIF': 'JPEG' } - TileInputUnits = { None: 'base_pixels', 'base': 'base_pixels', @@ -64,3 +65,17 @@ class SourcePriority: 'millimeters': 'mm', 'fraction': 'fraction', } + +# numpy dtype to pyvips GValue +dtypeToGValue = { + 'b': 'char', + 'B': 'uchar', + 'd': 'double', + 'D': 'dpcomplex', + 'f': 'float', + 'F': 'complex', + 'h': 'short', + 'H': 'ushort', + 'i': 'int', + 'I': 'uint', +} diff --git a/large_image/tilesource/base.py b/large_image/tilesource/base.py index 678f4b963..2f278a034 100644 --- a/large_image/tilesource/base.py +++ b/large_image/tilesource/base.py @@ -2,18 +2,21 @@ import json import math import numpy +import os +import pathlib import PIL import PIL.Image import PIL.ImageColor import PIL.ImageDraw +import tempfile import threading import xml.etree.ElementTree from collections import defaultdict from ..cache_util import getTileCache, strhash, methodcache -from ..constants import SourcePriority, \ - TILE_FORMAT_IMAGE, TILE_FORMAT_NUMPY, TILE_FORMAT_PIL, \ - TileOutputMimeTypes, TileOutputPILFormat, TileInputUnits +from ..constants import ( + SourcePriority, TILE_FORMAT_IMAGE, TILE_FORMAT_NUMPY, TILE_FORMAT_PIL, + TileOutputMimeTypes, TileOutputPILFormat, TileInputUnits, dtypeToGValue) from .. import config from .. import exceptions @@ -68,7 +71,7 @@ def _encodeImage(image, encoding='JPEG', jpegQuality=95, jpegSubsampling=0, if encoding == 'JPEG': params['quality'] = jpegQuality params['subsampling'] = jpegSubsampling - elif encoding == 'TIFF': + elif encoding in {'TIFF', 'TILED'}: params['compression'] = { 'none': 'raw', 'lzw': 'tiff_lzw', @@ -163,6 +166,152 @@ def _letterboxImage(image, width, height, fill): return result +def _vipsCast(image, mustBe8Bit=False, originalScale=None): + """ + Cast a vips image to a format we want. + + :param image: a vips image + :param mustBe9Bit: if True, then always cast to unsigned 8-bit. + :param originalScale: + # ##DWM:: + :returns: a vips image + """ + import pyvips + + formats = { + pyvips.BandFormat.CHAR: (pyvips.BandFormat.UCHAR, 2**7, 1), + pyvips.BandFormat.COMPLEX: (pyvips.BandFormat.USHORT, 0, 65535), + pyvips.BandFormat.DOUBLE: (pyvips.BandFormat.USHORT, 0, 65535), + pyvips.BandFormat.DPCOMPLEX: (pyvips.BandFormat.USHORT, 0, 65535), + pyvips.BandFormat.FLOAT: (pyvips.BandFormat.USHORT, 0, 65535), + pyvips.BandFormat.INT: (pyvips.BandFormat.USHORT, 2**31, 2**-16), + pyvips.BandFormat.USHORT: (pyvips.BandFormat.UCHAR, 0, 2**-8), + pyvips.BandFormat.SHORT: (pyvips.BandFormat.USHORT, 2**15, 1), + pyvips.BandFormat.UINT: (pyvips.BandFormat.USHORT, 0, 2**-16), + } + if image.format not in formats or (image.format == pyvips.BandFormat.USHORT and not mustBe8Bit): + return image + target, offset, multiplier = formats[image.format] + if image.format == pyvips.BandFormat.DOUBLE: + maxVal = image.max() + # These thresholds are higher than 256 and 65536 because bicubic and + # other interpolations can cause value spikes + if maxVal >= 2 and maxVal < 2**9: + multiplier = 256 + elif maxVal >= 256 and maxVal < 2**17: + multiplier = 1 + if mustBe8Bit and target != pyvips.BandFormat.UCHAR: + target = pyvips.BandFormat.UCHAR + multiplier /= 256 + # logger.debug('Casting image from %r to %r', image.format, target) + image = ((image.cast(pyvips.BandFormat.DOUBLE) + offset) * multiplier).cast(target) + return image + + +def _gdalParameters(defaultCompression=None, **kwargs): + """ + Return an array of gdal translation parameters. + + :param defaultCompression: if not specified, use this value. + + Optional parameters that can be specified in kwargs: + + :param tileSize: the horizontal and vertical tile size. + :param compression: one of 'jpeg', 'deflate' (zip), 'lzw', 'packbits', + 'zstd', or 'none'. + :param quality: a jpeg quality passed to gdal. 0 is small, 100 is high + quality. 90 or above is recommended. + :param level: compression level for zstd, 1-22 (default is 10). + :param predictor: one of 'none', 'horizontal', or 'float' used for lzw and + deflate. + :returns: a dictionary of parameters. + """ + options = { + 'tileSize': 256, + 'compression': 'lzw', + 'quality': 90, + 'predictor': 'yes', + } + predictor = { + 'none': 'NO', + 'horizontal': 'STANDARD', + 'float': 'FLOATING_POINT', + 'yes': 'YES', + } + options.update({k: v for k, v in kwargs.items() if v not in (None, '')}) + cmdopt = ['-of', 'COG', '-co', 'BIGTIFF=IF_SAFER'] + cmdopt += ['-co', 'BLOCKSIZE=%d' % options['tileSize']] + cmdopt += ['-co', 'COMPRESS=%s' % options['compression'].upper()] + cmdopt += ['-co', 'QUALITY=%s' % options['quality']] + cmdopt += ['-co', 'PREDICTOR=%s' % predictor[options['predictor']]] + if 'level' in options: + cmdopt += ['-co', 'LEVEL=%s' % options['level']] + return cmdopt + + +def _vipsParameters(forTiled=True, defaultCompression=None, **kwargs): + """ + Return a dictionary of vips conversion parameters. + + :param forTiled: True if this is for a tiled image. False for an + associated image. + :param defaultCompression: if not specified, use this value. + + Optional parameters that can be specified in kwargs: + + :param tileSize: the horizontal and vertical tile size. + :param compression: one of 'jpeg', 'deflate' (zip), 'lzw', 'packbits', + 'zstd', or 'none'. + :param quality: a jpeg quality passed to vips. 0 is small, 100 is high + quality. 90 or above is recommended. + :param level: compression level for zstd, 1-22 (default is 10). + :param predictor: one of 'none', 'horizontal', or 'float' used for lzw and + deflate. + :returns: a dictionary of parameters. + """ + if not forTiled: + convertParams = { + 'compression': defaultCompression or 'jpeg', + 'Q': 90, + 'predictor': 'horizontal', + 'tile': False, + } + if 'mime' in kwargs and kwargs.get('mime') != 'image/jpeg': + convertParams['compression'] = 'lzw' + return convertParams + convertParams = { + 'tile': True, + 'tile_width': 256, + 'tile_height': 256, + 'pyramid': True, + 'bigtiff': True, + 'compression': defaultCompression or 'jpeg', + 'Q': 90, + 'predictor': 'horizontal', + } + for vkey, kwkeys in { + 'tile_width': {'tileSize'}, + 'tile_height': {'tileSize'}, + 'compression': {'compression', 'tiffCompression'}, + 'Q': {'quality', 'jpegQuality'}, + 'level': {'level'}, + 'predictor': {'predictor'}, + }.items(): + for kwkey in kwkeys: + if kwkey in kwargs and kwargs[kwkey] not in {None, ''}: + convertParams[vkey] = kwargs[kwkey] + if convertParams['compression'] == 'jp2k': + convertParams['compression'] = 'none' + if convertParams['compression'] == 'webp' and kwargs.get('quality') == 0: + convertParams['lossless'] = True + convertParams.pop('Q', None) + if convertParams['predictor'] == 'yes': + convertParams['predictor'] = 'horizontal' + if convertParams['compression'] == 'jpeg': + convertParams['rgbjpeg'] = True + return convertParams + + def etreeToDict(t): """ Convert an xml etree to a nested dictionary without schema names in the @@ -485,7 +634,7 @@ def __init__(self, encoding='JPEG', jpegQuality=95, jpegSubsampling=0, :param jpegQuality: when serving jpegs, use this quality. :param jpegSubsampling: when serving jpegs, use this subsampling (0 is full chroma, 1 is half, 2 is quarter). - :param encoding: 'JPEG', 'PNG', or 'TIFF'. + :param encoding: 'JPEG', 'PNG', 'TIFF', or 'TILED'. :param edge: False to leave edge tiles whole, True or 'crop' to crop edge tiles, otherwise, an #rrggbb color to fill edges. :param tiffCompression: the compression format to use when encoding a @@ -969,7 +1118,7 @@ def _tileIteratorInfo(self, **kwargs): 'ymin': ymin, 'ymax': ymax}) # Use RGB for JPEG, RGBA for PNG - mode = 'RGBA' if kwargs.get('encoding') in ('PNG', 'TIFF') else 'RGB' + mode = 'RGBA' if kwargs.get('encoding') in {'PNG', 'TIFF', 'TILED'} else 'RGB' info = { 'region': { @@ -1542,7 +1691,7 @@ def _outputTile(self, tile, tileEncoding, x, y, z, pilImageAllowed=False, if encoding == 'JPEG': params['quality'] = self.jpegQuality params['subsampling'] = self.jpegSubsampling - elif encoding == 'TIFF': + elif encoding in {'TIFF', 'TILED'}: params['compression'] = self.tiffCompression tile.save(output, encoding, **params) return output.getvalue() @@ -1855,19 +2004,11 @@ 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']) - if image is None: - try: - image = numpy.zeros( - (regionHeight, regionWidth, subimage.shape[2]), - dtype=subimage.dtype) - except MemoryError: - raise exceptions.TileSourceException( - 'Insufficient memory to get region of %d x %d pixels.' % ( - regionWidth, regionHeight)) x0, y0 = tile['x'] - left, tile['y'] - top if x0 < 0: subimage = subimage[:, -x0:] @@ -1877,15 +2018,13 @@ def getRegion(self, format=(TILE_FORMAT_IMAGE, ), **kwargs): y0 = 0 subimage = subimage[:min(subimage.shape[0], regionHeight - y0), :min(subimage.shape[1], regionWidth - x0)] - if subimage.shape[2] > image.shape[2]: - newimage = numpy.ones((image.shape[0], image.shape[1], subimage.shape[2])) - newimage[:, :, :image.shape[2]] = image - image = newimage - image[y0:y0 + subimage.shape[0], x0:x0 + subimage.shape[1], - :subimage.shape[2]] = subimage + image = self._addRegionTileToImage( + image, subimage, x0, y0, regionWidth, regionHeight, tiled, tile, **kwargs) # Scale if we need to outWidth = int(math.floor(outWidth)) outHeight = int(math.floor(outHeight)) + if tiled: + return self._encodeTiledImage(image, outWidth, outHeight, iterInfo, **kwargs) if outWidth != regionWidth or outHeight != regionHeight: image = _imageToPIL(image, mode).resize( (outWidth, outHeight), @@ -1897,6 +2036,200 @@ 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, tiled=False, tile=None, **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. + :param tiled: true to generate a tiled output image. + :param tile: the original tile record with the current scale, etc. + :returns: the output image record. + """ + if tiled: + return self._addRegionTileToTiled(image, subimage, x, y, width, height, tile, **kwargs) + if image is None: + try: + image = numpy.zeros( + (height, width, subimage.shape[2]), + dtype=subimage.dtype) + except MemoryError: + raise exceptions.TileSourceException( + 'Insufficient memory to get region of %d x %d pixels.' % ( + width, height)) + if subimage.shape[2] > image.shape[2]: + newimage = numpy.ones((image.shape[0], image.shape[1], subimage.shape[2])) + newimage[:, :, :image.shape[2]] = image + image = newimage + image[y:y + subimage.shape[0], x:x + subimage.shape[1], + :subimage.shape[2]] = subimage + return image + + def _vipsAddAlphaBand(self, vimg, *otherImages): + """ + Add an alpha band to a vips image. The alpha value is either 1, 255, + or 65535 depending on the max value in the image and any other images + passed for reference. + + :param vimg: the image to modify. + :param otherImages: a list of other images to use for determining the + alpha value. + :returns: the original image with an alpha band. + """ + maxValue = vimg.max() + for img in otherImages: + maxValue = max(maxValue, img.max()) + alpha = 1 + if maxValue >= 2 and maxValue < 2**9: + alpha = 255 + elif maxValue >= 2**8 and maxValue < 2**17: + alpha = 65535 + return vimg.bandjoin(alpha) + + def _addRegionTileToTiled(self, image, subimage, x, y, width, height, tile=None, **kwargs): + """ + Add a subtile to a vips image. + + :param image: an object with information on the output. + :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. + :param tile: the original tile record with the current scale, etc. + :returns: the output object. + """ + import pyvips + + if subimage.dtype.char not in dtypeToGValue: + subimage = subimage.astype('d') + vimgMem = pyvips.Image.new_from_memory( + numpy.ascontiguousarray(subimage).data, + subimage.shape[1], subimage.shape[0], subimage.shape[2], + dtypeToGValue[subimage.dtype.char]) + vimg = pyvips.Image.new_temp_file('%s.v') + vimgMem.write(vimg) + if image is None: + image = { + 'width': width, + 'height': height, + 'mm_x': tile.get('mm_x') if tile else None, + 'mm_y': tile.get('mm_y') if tile else None, + 'magnification': tile.get('magnification') if tile else None, + 'channels': subimage.shape[2], + 'strips': {}, + } + if y not in image['strips']: + image['strips'][y] = vimg + if not x: + return image + if image['strips'][y].bands + 1 == vimg.bands: + image['strips'][y] = self._vipsAddAlphaBand(image['strips'][y], vimg) + elif vimg.bands + 1 == image['strips'][y].bands: + vimg = self._vipsAddAlphaBand(vimg, image['strips'][y]) + image['strips'][y] = image['strips'][y].insert(vimg, x, 0, expand=True) + return image + + def _encodeTiledImage(self, image, outWidth, outHeight, iterInfo, **kwargs): + """ + Given an image record of a set of vips image strips, generate a tiled + tiff file at the specified output size. + + :param image: a record with partial vips images and the current output + size. + :param outWidth: the output size after scaling and before any + letterboxing. + :param outHeight: the output size after scaling and before any + letterboxing. + :param iterInfo: information about the region based on the tile + iterator. + + Additional parameters are available. + + :param fill: a color to use in letterboxing. + :param maxWidth: the output size if letterboxing is applied. + :param maxHeight: the output size if letterboxing is applied. + :param compression: the internal compression format. This can handle + a variety of options similar to the converter utility. + :returns: a pathlib.Path of the output file and the output mime type. + """ + vimg = image['strips'][0] + for y in sorted(image['strips'].keys())[1:]: + if image['strips'][y].bands + 1 == vimg.bands: + image['strips'][y] = self._vipsAddAlphaBand(image['strips'][y], vimg) + elif vimg.bands + 1 == image['strips'][y].bands: + vimg = self._vipsAddAlphaBand(vimg, image['strips'][y]) + vimg = vimg.insert(image['strips'][y], 0, y, expand=True) + + if outWidth != image['width'] or outHeight != image['height']: + scale = outWidth / image['width'] + vimg = vimg.resize(outWidth / image['width'], vscale=outHeight / image['height']) + image['width'] = outWidth + image['height'] = outHeight + image['mm_x'] = image['mm_x'] / scale if image['mm_x'] else image['mm_x'] + image['mm_y'] = image['mm_y'] / scale if image['mm_y'] else image['mm_y'] + image['magnification'] = ( + image['magnification'] * scale + if image['magnification'] else image['magnification']) + return self._encodeTiledImageFromVips(vimg, iterInfo, image, **kwargs) + + def _encodeTiledImageFromVips(self, vimg, iterInfo, image, **kwargs): + """ + Save a vips image as a tiled tiff. + + :param vimg: a vips image. + :param iterInfo: information about the region based on the tile + iterator. + :param image: a record with partial vips images and the current output + size. + + Additional parameters are available. + + :param compression: the internal compression format. This can handle + a variety of options similar to the converter utility. + :returns: a pathlib.Path of the output file and the output mime type. + """ + import pyvips + + convertParams = _vipsParameters(defaultCompression='lzw', **kwargs) + vimg = _vipsCast(vimg, convertParams['compression'] in {'webp', 'jpeg'}) + maxWidth = kwargs.get('output', {}).get('maxWidth') + maxHeight = kwargs.get('output', {}).get('maxHeight') + if (kwargs.get('fill') and str(kwargs.get('fill')).lower() != 'none' and + maxWidth and maxHeight and + (maxWidth > image['width'] or maxHeight > image['height'])): + color = PIL.ImageColor.getcolor( + kwargs.get('fill'), ['L', 'LA', 'RGB', 'RGBA'][vimg.bands - 1]) + lbimage = pyvips.Image.black(maxWidth, maxHeight, bands=vimg.bands) + lbimage = lbimage.cast(vimg.format) + lbimage = lbimage.draw_rect( + [c * (257 if vimg.format == pyvips.BandFormat.USHORT else 1) for c in color], + 0, 0, maxWidth, maxHeight, fill=True) + vimg = lbimage.insert( + vimg, (maxWidth - image['width']) // 2, (maxHeight - image['height']) // 2) + if image['mm_x'] and image['mm_y']: + vimg = vimg.copy(xres=1 / image['mm_x'], yres=1 / image['mm_y']) + fd, outputPath = tempfile.mkstemp('.tiff', 'tiledRegion_') + os.close(fd) + try: + vimg.write_to_file(outputPath, **convertParams) + return pathlib.Path(outputPath), TileOutputMimeTypes['TILED'] + except Exception as exc: + try: + pathlib.Path(outputPath).unlink() + except Exception: + pass + raise exc + def getRegionAtAnotherScale(self, sourceRegion, sourceScale=None, targetScale=None, targetUnits=None, **kwargs): """ @@ -2166,8 +2499,8 @@ def tileIterator(self, format=(TILE_FORMAT_NUMPY, ), resample=True, with the non-overlapped area of each as 0, 1, 2, 3, 4, 5, 6, 7. :param encoding: if format includes TILE_FORMAT_IMAGE, a valid PIL - encoding (typically 'PNG', 'JPEG', or 'TIFF'). Must also be in the - TileOutputMimeTypes map. + encoding (typically 'PNG', 'JPEG', or 'TIFF') or 'TILED' (identical + to TIFF). Must also be in the TileOutputMimeTypes map. :param jpegQuality: the quality to use when encoding a JPEG. :param jpegSubsampling: the subsampling level to use when encoding a JPEG. diff --git a/sources/gdal/large_image_source_gdal/__init__.py b/sources/gdal/large_image_source_gdal/__init__.py index 53e70fec7..09a18c0d0 100644 --- a/sources/gdal/large_image_source_gdal/__init__.py +++ b/sources/gdal/large_image_source_gdal/__init__.py @@ -16,10 +16,13 @@ import math import numpy +import os import palettable +import pathlib import PIL.Image import pyproj import struct +import tempfile import threading from operator import attrgetter from osgeo import gdal @@ -28,9 +31,11 @@ from osgeo import osr from pkg_resources import DistributionNotFound, get_distribution +import large_image from large_image import config from large_image.cache_util import LruCacheMetaclass, methodcache, CacheProperties -from large_image.constants import SourcePriority, TileInputUnits, TILE_FORMAT_NUMPY, TILE_FORMAT_PIL +from large_image.constants import ( + SourcePriority, TileInputUnits, TileOutputMimeTypes, TILE_FORMAT_NUMPY, TILE_FORMAT_PIL) from large_image.exceptions import TileSourceException from large_image.tilesource import FileTileSource @@ -878,6 +883,31 @@ def _getRegionBounds(self, metadata, left=None, top=None, right=None, return super()._getRegionBounds( metadata, left, top, right, bottom, width, height, units, **kwargs) + def pixelToProjection(self, x, y, level=None): + """ + Convert from pixels back to projection coordinates. + + :param x, y: base pixel coordinates. + :param level: the level of the pixel. None for maximum level. + :returns: x, y in projection coordinates. + """ + if level is None: + level = self.levels - 1 + if not self.projection: + x *= 2 ** (self.levels - 1 - level) + y *= 2 ** (self.levels - 1 - level) + gt = self._getGeoTransform() + px = gt[0] + gt[1] * x + gt[2] * y + py = gt[3] + gt[4] * x + gt[5] * y + return px, py + xScale = 2 ** level * self.tileWidth + yScale = 2 ** level * self.tileHeight + x = x / xScale - 0.5 + y = 0.5 - y / yScale + x = x * self.unitsAcrossLevel0 + self.projectionOrigin[0] + y = y * self.unitsAcrossLevel0 + self.projectionOrigin[1] + return x, y + @methodcache() def getThumbnail(self, width=None, height=None, **kwargs): """ @@ -974,6 +1004,66 @@ def getPixel(self, **kwargs): pass return pixel + def _encodeTiledImageFromVips(self, vimg, iterInfo, image, **kwargs): + """ + Save a vips image as a tiled tiff. + + :param vimg: a vips image. + :param iterInfo: information about the region based on the tile + iterator. + :param image: a record with partial vips images and the current output + size. + + Additional parameters are available. + + :param compression: the internal compression format. This can handle + a variety of options similar to the converter utility. + :returns: a pathlib.Path of the output file and the output mime type. + """ + convertParams = large_image.tilesource.base._vipsParameters( + defaultCompression='lzw', **kwargs) + convertParams.pop('pyramid', None) + vimg = large_image.tilesource.base._vipsCast( + vimg, convertParams['compression'] in {'webp', 'jpeg'}) + gdalParams = large_image.tilesource.base._gdalParameters( + defaultCompression='lzw', **kwargs) + for ch in range(image['channels']): + gdalParams += [ + '-b' if ch not in (1, 3) or ch + 1 != image['channels'] else '-mask', str(ch + 1)] + tl = self.pixelToProjection( + iterInfo['region']['left'], iterInfo['region']['top'], iterInfo['level']) + br = self.pixelToProjection( + iterInfo['region']['right'], iterInfo['region']['bottom'], iterInfo['level']) + gdalParams += [ + '-a_srs', + iterInfo['metadata']['bounds']['srs'], + '-a_ullr', + str(tl[0]), + str(tl[1]), + str(br[0]), + str(br[1]), + ] + fd, tempPath = tempfile.mkstemp('.tiff', 'tiledRegion_') + os.close(fd) + fd, outputPath = tempfile.mkstemp('.tiff', 'tiledGeoRegion_') + os.close(fd) + try: + vimg.write_to_file(tempPath, **convertParams) + ds = gdal.Open(tempPath, gdalconst.GA_ReadOnly) + gdal.Translate(outputPath, ds, options=gdalParams) + os.unlink(tempPath) + except Exception as exc: + try: + os.unlink(tempPath) + except Exception: + pass + try: + os.unlink(outputPath) + except Exception: + pass + raise exc + return pathlib.Path(outputPath), TileOutputMimeTypes['TILED'] + def open(*args, **kwargs): """ diff --git a/test/test_source_gdal.py b/test/test_source_gdal.py index dade5c244..64281a6e5 100644 --- a/test/test_source_gdal.py +++ b/test/test_source_gdal.py @@ -374,3 +374,64 @@ def testGCPProjection(): assert tileMetadata['bounds']['ymin'] == pytest.approx(3899358, 1) assert tileMetadata['bounds']['srs'] == '+init=epsg:3857' assert tileMetadata['geospatial'] + + +def testGetTiledRegion(): + imagePath = datastore.fetch('landcover_sample_1000.tif') + ts = large_image_source_gdal.open(imagePath) + region, _ = ts.getRegion(output=dict(maxWidth=1024, maxHeight=1024), + encoding='TILED') + result = large_image_source_gdal.open(str(region)) + tileMetadata = result.getMetadata() + assert tileMetadata['bounds']['xmax'] == pytest.approx(2006547, 1) + assert tileMetadata['bounds']['xmin'] == pytest.approx(1319547, 1) + assert tileMetadata['bounds']['ymax'] == pytest.approx(2658548, 1) + assert tileMetadata['bounds']['ymin'] == pytest.approx(2149548, 1) + assert '+proj=aea' in tileMetadata['bounds']['srs'] + region.unlink() + + +def testGetTiledRegionWithProjection(): + imagePath = datastore.fetch('landcover_sample_1000.tif') + ts = large_image_source_gdal.open(imagePath, projection='EPSG:3857') + # This gets the whole world + region, _ = ts.getRegion(output=dict(maxWidth=1024, maxHeight=1024), + encoding='TILED') + result = large_image_source_gdal.open(str(region)) + tileMetadata = result.getMetadata() + assert tileMetadata['bounds']['xmax'] == pytest.approx(20037508, 1) + assert tileMetadata['bounds']['xmin'] == pytest.approx(-20037508, 1) + assert tileMetadata['bounds']['ymax'] == pytest.approx(20037508, 1) + assert tileMetadata['bounds']['ymin'] == pytest.approx(-20037508, 1) + assert '+proj=merc' in tileMetadata['bounds']['srs'] + region.unlink() + + # Ask for a smaller part + region, _ = ts.getRegion( + output=dict(maxWidth=1024, maxHeight=1024), + region=dict(left=-8622811, right=-8192317, bottom=5294998, + top=5477835, units='projection'), + encoding='TILED') + result = large_image_source_gdal.open(str(region)) + tileMetadata = result.getMetadata() + assert tileMetadata['bounds']['xmax'] == pytest.approx(-8192215, 1) + assert tileMetadata['bounds']['xmin'] == pytest.approx(-8622708, 1) + assert tileMetadata['bounds']['ymax'] == pytest.approx(5477783, 1) + assert tileMetadata['bounds']['ymin'] == pytest.approx(5294946, 1) + assert '+proj=merc' in tileMetadata['bounds']['srs'] + region.unlink() + + +def testGetTiledRegion16Bit(): + imagePath = datastore.fetch('region_gcp.tiff') + ts = large_image_source_gdal.open(imagePath) + region, _ = ts.getRegion(output=dict(maxWidth=1024, maxHeight=1024), + encoding='TILED') + result = large_image_source_gdal.open(str(region)) + tileMetadata = result.getMetadata() + assert tileMetadata['bounds']['xmax'] == pytest.approx(-10753925, 1) + assert tileMetadata['bounds']['xmin'] == pytest.approx(-10871650, 1) + assert tileMetadata['bounds']['ymax'] == pytest.approx(3949393, 1) + assert tileMetadata['bounds']['ymin'] == pytest.approx(3899358, 1) + assert '+proj=merc' in tileMetadata['bounds']['srs'] + region.unlink() diff --git a/test/test_source_tiff.py b/test/test_source_tiff.py index bd9eeacc1..e4e9e2bc3 100644 --- a/test/test_source_tiff.py +++ b/test/test_source_tiff.py @@ -1,8 +1,10 @@ +import io import json import numpy import os import pytest import struct +import tifftools from large_image import constants import large_image_source_tiff @@ -337,35 +339,37 @@ def testThumbnails(): source.getThumbnail(**entry[0]) +@pytest.mark.parametrize('badParams,errMessage', [ + ({'encoding': 'invalid', 'width': 10}, 'Invalid encoding'), + ({'output': {'maxWidth': 'invalid'}}, 'ValueError'), + ({'output': {'maxWidth': -5}}, 'Invalid output width or height'), + ({'output': {'maxWidth': None, 'maxHeight': 'invalid'}}, 'ValueError'), + ({'output': {'maxWidth': None, 'maxHeight': -5}}, 'Invalid output width or height'), + ({'jpegQuality': 'invalid'}, 'ValueError'), + ({'jpegSubsampling': 'invalid'}, 'TypeError'), + ({'region': {'left': 'invalid'}}, 'TypeError'), + ({'region': {'right': 'invalid'}}, 'TypeError'), + ({'region': {'top': 'invalid'}}, 'TypeError'), + ({'region': {'bottom': 'invalid'}}, 'TypeError'), + ({'region': {'width': 'invalid'}}, 'TypeError'), + ({'region': {'height': 'invalid'}}, 'TypeError'), + ({'region': {'units': 'invalid'}}, 'Invalid units'), + ({'region': {'unitsWH': 'invalid'}}, 'Invalid units'), +]) +def testRegionBadParameters(badParams, errMessage): + imagePath = datastore.fetch('sample_image.ptif') + source = large_image_source_tiff.open(imagePath) + with pytest.raises(Exception): + params = {'output': {'maxWidth': 400}} + nestedUpdate(params, badParams) + source.getRegion(**params) + + def testRegions(): imagePath = datastore.fetch('sample_image.ptif') source = large_image_source_tiff.open(imagePath) tileMetadata = source.getMetadata() - # Test bad parameters - badParams = [ - ({'encoding': 'invalid', 'width': 10}, 'Invalid encoding'), - ({'output': {'maxWidth': 'invalid'}}, 'ValueError'), - ({'output': {'maxWidth': -5}}, 'Invalid output width or height'), - ({'output': {'maxWidth': None, 'maxHeight': 'invalid'}}, 'ValueError'), - ({'output': {'maxWidth': None, 'maxHeight': -5}}, 'Invalid output width or height'), - ({'jpegQuality': 'invalid'}, 'ValueError'), - ({'jpegSubsampling': 'invalid'}, 'TypeError'), - ({'region': {'left': 'invalid'}}, 'TypeError'), - ({'region': {'right': 'invalid'}}, 'TypeError'), - ({'region': {'top': 'invalid'}}, 'TypeError'), - ({'region': {'bottom': 'invalid'}}, 'TypeError'), - ({'region': {'width': 'invalid'}}, 'TypeError'), - ({'region': {'height': 'invalid'}}, 'TypeError'), - ({'region': {'units': 'invalid'}}, 'Invalid units'), - ({'region': {'unitsWH': 'invalid'}}, 'Invalid units'), - ] - for entry in badParams: - with pytest.raises(Exception): - params = {'output': {'maxWidth': 400}} - nestedUpdate(params, entry[0]) - source.getRegion(**params) - # Get a small region for testing. Our test file is sparse, so # initially get a region where there is full information. params = {'region': {'width': 1000, 'height': 1000, 'left': 48000, 'top': 3000}} @@ -453,6 +457,53 @@ def testRegions(): assert image != nextimage +def testRegionTiledOutputIsTiled(): + imagePath = datastore.fetch('sample_image.ptif') + source = large_image_source_tiff.open(imagePath) + + # TIFF isn't tiled and has only one layer + params = {'output': {'maxWidth': 500, 'maxHeight': 500}, + 'encoding': 'TIFF'} + image, mimeType = source.getRegion(**params) + info = tifftools.read_tiff(io.BytesIO(image)) + assert len(info['ifds']) == 1 + assert tifftools.Tag.StripOffsets.value in info['ifds'][0]['tags'] + assert tifftools.Tag.TileOffsets.value not in info['ifds'][0]['tags'] + + # TILED is tiled and has multiple layers + params = {'output': {'maxWidth': 500, 'maxHeight': 500}, + 'encoding': 'TILED'} + image, mimeType = source.getRegion(**params) + info = tifftools.read_tiff(image) + assert len(info['ifds']) == 2 + assert tifftools.Tag.StripOffsets.value not in info['ifds'][0]['tags'] + assert tifftools.Tag.TileOffsets.value in info['ifds'][0]['tags'] + os.unlink(image) + + # Bigger outputs should have more layers + params = {'output': {'maxWidth': 3000, 'maxHeight': 3000}, + 'encoding': 'TILED'} + image, mimeType = source.getRegion(**params) + info = tifftools.read_tiff(image) + assert len(info['ifds']) == 5 + assert tifftools.Tag.StripOffsets.value not in info['ifds'][0]['tags'] + assert tifftools.Tag.TileOffsets.value in info['ifds'][0]['tags'] + os.unlink(image) + + +def testRegionTiledOutputLetterbox(): + imagePath = datastore.fetch('sample_image.ptif') + source = large_image_source_tiff.open(imagePath) + params = {'output': {'maxWidth': 500, 'maxHeight': 500}, + 'fill': 'pink', + 'encoding': 'TILED'} + image, mimeType = source.getRegion(**params) + result = large_image_source_tiff.open(str(image)) + assert result.sizeX == 500 + assert result.sizeY == 500 + os.unlink(image) + + def testPixel(): imagePath = datastore.fetch('sample_image.ptif') source = large_image_source_tiff.open(imagePath) diff --git a/test/utilities.py b/test/utilities.py index f5085167a..d59fec1e1 100644 --- a/test/utilities.py +++ b/test/utilities.py @@ -6,6 +6,7 @@ JPEGHeader = b'\xff\xd8\xff' PNGHeader = b'\x89PNG' TIFFHeader = b'II\x2a\x00' +BigTIFFHeader = b'II\x2b\x00' def checkTilesZXY(source, metadata, tileParams=None, imgHeader=JPEGHeader): diff --git a/tox.ini b/tox.ini index 1599333be..8b35283e8 100644 --- a/tox.ini +++ b/tox.ini @@ -43,6 +43,12 @@ setenv = GDAL_PAM_ENABLED=no PIP_FIND_LINKS=https://girder.github.io/large_image_wheels +# To run just some non-web-client tests, do tox -e server -- -k +# Don't use this for CI or full tests. +[testenv:server] +commands = + pytest {posargs} + [testenv:flake8] skipsdist = true skip_install = true diff --git a/utilities/converter/large_image_converter/__init__.py b/utilities/converter/large_image_converter/__init__.py index 1d2d61620..54183f37e 100644 --- a/utilities/converter/large_image_converter/__init__.py +++ b/utilities/converter/large_image_converter/__init__.py @@ -98,26 +98,7 @@ def _generate_geotiff(inputPath, outputPath, **kwargs): from osgeo import gdal from osgeo import gdalconst - options = { - 'tileSize': 256, - 'compression': 'lzw', - 'quality': 90, - 'predictor': 'yes', - } - predictor = { - 'none': 'NO', - 'horizontal': 'STANDARD', - 'float': 'FLOATING_POINT', - 'yes': 'YES', - } - options.update({k: v for k, v in kwargs.items() if v not in (None, '')}) - cmdopt = ['-of', 'COG', '-co', 'BIGTIFF=YES'] - cmdopt += ['-co', 'BLOCKSIZE=%d' % options['tileSize']] - cmdopt += ['-co', 'COMPRESS=%s' % options['compression'].upper()] - cmdopt += ['-co', 'QUALITY=%s' % options['quality']] - cmdopt += ['-co', 'PREDICTOR=%s' % predictor[options['predictor']]] - if 'level' in options: - cmdopt += ['-co', 'LEVEL=%s' % options['level']] + cmdopt = large_image.tilesource.base._gdalParameters(**kwargs) cmd = ['gdal_translate', inputPath, outputPath] + cmdopt logger.info('Convert to geotiff: %r', cmd) try: @@ -245,10 +226,10 @@ def _convert_via_vips(inputPathOrBuffer, outputPath, tempPath, forTiled=True, also stores files in TMPDIR :param forTiled: True if the output should be tiled, false if not. :param status: an optional additional string to add to log messages. - :param kwargs: addition arguments that get passed to _vips_parameters + :param kwargs: addition arguments that get passed to _vipsParameters and _convert_to_jp2k. """ - convertParams = _vips_parameters(forTiled, **kwargs) + convertParams = large_image.tilesource.base._vipsParameters(forTiled, **kwargs) status = (', ' + status) if status else '' if type(inputPathOrBuffer) == pyvips.vimage.Image: source = 'vips image' @@ -271,7 +252,7 @@ def _convert_via_vips(inputPathOrBuffer, outputPath, tempPath, forTiled=True, image.interpretation != pyvips.Interpretation.SCRGB): # jp2k compression supports more than 8-bits per sample, but the # decompressor claims this is unsupported. - image = _vips_cast( + image = large_image.tilesource.base._vipsCast( image, convertParams['compression'] in {'webp', 'jpeg'} or kwargs.get('compression') in {'jp2k'}) @@ -462,24 +443,12 @@ def _convert_large_image_tile(tilelock, strips, tile): :param tile: a tileIterator tile. """ data = tile['tile'] - dtypeToGValue = { - 'b': 'char', - 'B': 'uchar', - 'd': 'double', - 'D': 'dpcomplex', - 'f': 'float', - 'F': 'complex', - 'h': 'short', - 'H': 'ushort', - 'i': 'int', - 'I': 'uint', - } - if data.dtype.char not in dtypeToGValue: + if data.dtype.char not in large_image.constants.dtypeToGValue: data = data.astype('d') vimg = pyvips.Image.new_from_memory( numpy.ascontiguousarray(data).data, data.shape[1], data.shape[0], data.shape[2], - dtypeToGValue[data.dtype.char]) + large_image.constants.dtypeToGValue[data.dtype.char]) vimgTemp = pyvips.Image.new_temp_file('%s.v') vimg.write(vimgTemp) vimg = vimgTemp @@ -810,97 +779,6 @@ def _make_li_description( return json.dumps(results, separators=(',', ':'), sort_keys=True, default=json_serial) -def _vips_cast(image, mustBe8Bit=False): - """ - Cast a vips image to a format we want. - - :param image: a vips image - :param mustBe9Bit: if True, then always cast to unsigned 8-bit. - :returns: a vips image - """ - formats = { - pyvips.BandFormat.CHAR: (pyvips.BandFormat.UCHAR, 2**7, 1), - pyvips.BandFormat.COMPLEX: (pyvips.BandFormat.USHORT, 0, 65535), - pyvips.BandFormat.DOUBLE: (pyvips.BandFormat.USHORT, 0, 65535), - pyvips.BandFormat.DPCOMPLEX: (pyvips.BandFormat.USHORT, 0, 65535), - pyvips.BandFormat.FLOAT: (pyvips.BandFormat.USHORT, 0, 65535), - pyvips.BandFormat.INT: (pyvips.BandFormat.USHORT, 2**31, 2**-16), - pyvips.BandFormat.USHORT: (pyvips.BandFormat.UCHAR, 0, 2**-8), - pyvips.BandFormat.SHORT: (pyvips.BandFormat.USHORT, 2**15, 1), - pyvips.BandFormat.UINT: (pyvips.BandFormat.USHORT, 0, 2**-16), - } - if image.format not in formats or (image.format == pyvips.BandFormat.USHORT and not mustBe8Bit): - return image - target, offset, multiplier = formats[image.format] - if mustBe8Bit and target != pyvips.BandFormat.UCHAR: - target = pyvips.BandFormat.UCHAR - multiplier /= 256 - logger.debug('Casting image from %r to %r', image.format, target) - image = ((image.cast(pyvips.BandFormat.DOUBLE) + offset) * multiplier).cast(target) - return image - - -def _vips_parameters(forTiled=True, **kwargs): - """ - Return a dictionary of vips conversion parameters. - - :param forTiled: True if this is for a tiled image. False for an - associated image. - - Optional parameters that can be specified in kwargs: - - :param tileSize: the horizontal and vertical tile size. - :param compression: one of 'jpeg', 'deflate' (zip), 'lzw', 'packbits', - 'zstd', or 'none'. - :param quality: a jpeg quality passed to vips. 0 is small, 100 is high - quality. 90 or above is recommended. - :param level: compression level for zstd, 1-22 (default is 10). - :param predictor: one of 'none', 'horizontal', or 'float' used for lzw and - deflate. - :returns: a dictionary of parameters. - """ - if not forTiled: - convertParams = { - 'compression': 'jpeg', - 'Q': 90, - 'predictor': 'horizontal', - 'tile': False, - } - if 'mime' in kwargs and kwargs.get('mime') != 'image/jpeg': - convertParams['compression'] = 'lzw' - return convertParams - convertParams = { - 'tile': True, - 'tile_width': 256, - 'tile_height': 256, - 'pyramid': True, - 'bigtiff': True, - 'compression': 'jpeg', - 'Q': 90, - 'predictor': 'horizontal', - } - for vkey, kwkey in { - 'tile_width': 'tileSize', - 'tile_height': 'tileSize', - 'compression': 'compression', - 'Q': 'quality', - 'level': 'level', - 'predictor': 'predictor', - }.items(): - if kwkey in kwargs and kwargs[kwkey] not in {None, ''}: - convertParams[vkey] = kwargs[kwkey] - if convertParams['compression'] == 'jp2k': - convertParams['compression'] = 'none' - if convertParams['compression'] == 'webp' and kwargs.get('quality') == 0: - convertParams['lossless'] = True - convertParams.pop('Q', None) - if convertParams['predictor'] == 'yes': - convertParams['predictor'] = 'horizontal' - if convertParams['compression'] == 'jpeg': - convertParams['rgbjpeg'] = True - return convertParams - - def format_hook(funcname, *args, **kwargs): """ Call a function specific to a file format.