Source code for toasty.study

# -*- mode: python; coding: utf-8 -*-
# Copyright 2021-2022 the AAS WorldWide Telescope project
# Licensed under the MIT License.

"""Common routines for tiling images anchored to the sky in a gnomonic
(tangential) projection.

"""

__all__ = """
StudyTiling
tile_study_image
""".split()

import numpy as np

from .progress import progress_bar
from .pyramid import Pos, next_highest_power_of_2


[docs] class StudyTiling(object): """ Information about how a WWT "study" image is broken into tiles. In WWT a "study" is a large astronomical image projected onto the sky using a gnomonic (tangential or TAN) projection. The image may have many pixels, so it is broken up into tiles for efficient visualization. Note that this class doesn't know anything about how the image is projected onto the sky, or whether it is projected at all. The core tiling functionality doesn't need to care about that. """ _width = None "The width of the region in which image data are available, in pixels (int)." _height = None "The height of the region in which image data are available, in pixels (int)." _p2n = None "The size of the study when tiled, in pixels - a power of 2." _tile_size = None "The number of tiles wide and high in the tiled study - a power of 2." _tile_levels = None "The number of levels in the tiling - ``log2(tile_size)``." _img_gx0 = None """The pixel position of the image data within the global (tiled) pixelization, in pixels right from the left edge (int). Nonnegative. """ _img_gy0 = None """The pixel position of the image data within the global (tiled) pixelization, in pixels down from the top edge (int). Nonnegative. """ def __init__(self, width, height): """Set up the tiling information. Parameters ---------- width : positive integer The width of the full-resolution image, in pixels. height : positive integer The height of the full-resolution image, in pixels. """ width = int(width) height = int(height) if width <= 0: raise ValueError("bad width value: %r" % (width,)) if height <= 0: raise ValueError("bad height value: %r" % (height,)) self._width = width self._height = height p2w = next_highest_power_of_2(self._width) p2h = next_highest_power_of_2(self._height) self._p2n = max(p2w, p2h) self._tile_size = self._p2n // 256 self._tile_levels = int(np.log2(self._tile_size)) self._img_gx0 = (self._p2n - self._width) // 2 self._img_gy0 = (self._p2n - self._height) // 2
[docs] def compute_for_subimage(self, subim_ix, subim_iy, subim_width, subim_height): """ Create a new compatible tiling whose underlying image is a subset of this one. Parameters ---------- subim_ix : integer The 0-based horizontal pixel position of the left edge of the sub-image, relative to this tiling's image. subim_iy : integer The 0-based vertical pixel position of the top edge of the sub-image, relative to this tiling's image. subim_width : nonnegative integer The width of the sub-image, in pixels. subim_height : nonnegative integer The height of the sub-image, in pixels. Returns ------- A new :class:`~StudyTiling` with the same number of tile levels as this one. However, the internal information about where the available data land within that tiling will be appropriate for the specified sub-image. Methods like :meth:`count_populated_positions`, :meth:`generate_populated_positions`, and :meth:`tile_image` will behave differently. """ if subim_width < 0 or subim_width > self._width: raise ValueError("bad subimage width value {!r}".format(subim_width)) if subim_height < 0 or subim_height > self._height: raise ValueError("bad subimage height value {!r}".format(subim_height)) if subim_ix < 0 or subim_ix + subim_width > self._width: raise ValueError("bad subimage ix value {!r}".format(subim_ix)) if subim_iy < 0 or subim_iy + subim_height > self._height: raise ValueError("bad subimage iy value {!r}".format(subim_iy)) sub_tiling = StudyTiling(self._width, self._height) sub_tiling._width = subim_width sub_tiling._height = subim_height sub_tiling._img_gx0 += subim_ix sub_tiling._img_gy0 += subim_iy return sub_tiling
[docs] def n_deepest_layer_tiles(self): """Return the number of tiles in the highest-resolution layer.""" return 4**self._tile_levels
[docs] def apply_to_imageset(self, imgset): """Fill the specific ``wwt_data_formats.imageset.ImageSet`` object with parameters defined by this tiling, Parameters ---------- imgset : ``wwt_data_formats.imageset.ImageSet`` The object to modify Notes ----- The settings currently transferred are the number of tile levels and the projection type. """ from wwt_data_formats.enums import ProjectionType imgset.tile_levels = self._tile_levels if self._tile_levels == 0: imgset.projection = ProjectionType.SKY_IMAGE else: imgset.projection = ProjectionType.TAN
[docs] def image_to_tile(self, im_ix, im_iy): """Convert an image pixel position to a tiled pixel position. Parameters ---------- im_ix : integer A 0-based horizontal pixel position in the image coordinate system. im_iy : integer A 0-based vertical pixel position in the image coordinate system. Notes ----- ``(0, 0)`` is the top-left corner of the image. The input values need not lie on the image. (I.e., they may be negative.) Returns ``(tile_ix, tile_iy, subtile_ix, subtile_iy)``, where - *tile_ix* is X index of the matched tile in the tiling, between 0 and 2**tile_size - 1. Measured right from the left edge of the tiling. - *tile_iy* is Y index of the matched tile in the tiling, between 0 and 2**tile_size - 1. Measured down from the top of the tiling. - *subtile_ix* is the pixel X position within that tile, between 0 and 255. Measured right from the left edge of the tiling. - *subtile_iy* is the pixel Y position within that tile, between 0 and 255. Measured down from the top edge of the tiling. """ gx = im_ix + self._img_gx0 gy = im_iy + self._img_gy0 tile_ix = np.floor(gx // 256).astype(int) tile_iy = np.floor(gy // 256).astype(int) return (tile_ix, tile_iy, gx % 256, gy % 256)
[docs] def count_populated_positions(self): """ Count how many tiles contain image data. This is used for progress reporting. """ img_gx1 = self._img_gx0 + self._width - 1 img_gy1 = self._img_gy0 + self._height - 1 tile_start_tx = self._img_gx0 // 256 tile_start_ty = self._img_gy0 // 256 tile_end_tx = img_gx1 // 256 tile_end_ty = img_gy1 // 256 return (tile_end_ty + 1 - tile_start_ty) * (tile_end_tx + 1 - tile_start_tx)
[docs] def generate_populated_positions(self): """Generate information about tiles containing image data. Generates a sequence of tuples ``(pos, width, height, image_x, image_y, tile_x, tile_y)`` where: - *pos* is a :class:`toasty.pyramid.Pos` tuple giving parameters of a tile - *width* is the width of the rectangle of image data contained in this tile, between 1 and 256. - *height* is the height of the rectangle of image data contained in this tile, between 1 and 256. - *image_x* is the pixel X coordinate of the left edge of the image data in this tile in the image rectangle, increasing from the left edge of the tile. Between 0 and ``self._width - 1`` (inclusive). - *image_y* is the pixel Y coordinate of the *top* edge of the image data in this tile in the image rectangle, increasing from the top edge of the tile. Between 0 and ``self._height - 1`` (inclusive). - *tile_x* is the pixel X coordinate of the left edge of the image data in this tile in the tile rectangle, increasing from the left edge of the tile. Between 0 and 255 (inclusive). - *tile_y* is the pixel Y coordinate of the *top* edge of the image data in this tile in the tile rectangle, increasing from the top edge of the tile. Between 0 and 255 (inclusive). Tiles that do not overlap the image at all are not generated. Tiles that are completely filled with image data will yield tuples of the form ``(pos, 256, 256, im_x, im_y, 0, 0)``. An image that fits entirely in one tile will yield a tuple of the form ``(Pos(n=0, x=0, y=0), width, height, 0, 0, tx, ty)``. """ # Get the position of the actual image data in "global pixel # coordinates", which span the whole tiled region (a superset of the # image itself) with x=0, y=0 being the left-top corner of the tiled # region. img_gx1 = ( self._img_gx0 + self._width - 1 ) # inclusive: there are image data in this column img_gy1 = self._img_gy0 + self._height - 1 # ditto tile_start_tx = self._img_gx0 // 256 tile_start_ty = self._img_gy0 // 256 tile_end_tx = ( img_gx1 // 256 ) # inclusive; there are image data in this column of tiles tile_end_ty = img_gy1 // 256 # ditto for ity in range(tile_start_ty, tile_end_ty + 1): for itx in range(tile_start_tx, tile_end_tx + 1): # (inclusive) tile bounds in global pixel coords tile_gx0 = itx * 256 tile_gy0 = ity * 256 tile_gx1 = tile_gx0 + 255 tile_gy1 = tile_gy0 + 255 # overlap (= intersection) of the image and the tile in global pixel coords overlap_gx0 = max(tile_gx0, self._img_gx0) overlap_gy0 = max(tile_gy0, self._img_gy0) overlap_gx1 = min(tile_gx1, img_gx1) overlap_gy1 = min(tile_gy1, img_gy1) # coordinates of the overlap in image pixel coords img_overlap_x0 = overlap_gx0 - self._img_gx0 img_overlap_x1 = overlap_gx1 - self._img_gx0 img_overlap_y0 = overlap_gy0 - self._img_gy0 img_overlap_y1 = overlap_gy1 - self._img_gy0 # shape of the overlap overlap_width = img_overlap_x1 + 1 - img_overlap_x0 overlap_height = img_overlap_y1 + 1 - img_overlap_y0 # coordinates of the overlap in this tile's coordinates tile_overlap_x0 = overlap_gx0 - tile_gx0 tile_overlap_y0 = overlap_gy0 - tile_gy0 yield ( Pos(self._tile_levels, itx, ity), overlap_width, overlap_height, img_overlap_x0, img_overlap_y0, tile_overlap_x0, tile_overlap_y0, )
[docs] def tile_image(self, image, pio, cli_progress=False): """Tile an in-memory image as a study. Parameters ---------- image : :class:`toasty.image.Image` In-memory image data. The image's dimensions must match the ones for which this tiling was computed. pio : :class:`toasty.pyramid.PyramidIO` A handle for doing I/O on the tile pyramid cli_progress : optional boolean, defaults False If true, a progress bar will be printed to the terminal. Returns ------- Self. """ if image.height != self._height: raise ValueError("height of image to be sampled does not match tiling") if image.width != self._width: raise ValueError("width of image to be sampled does not match tiling") # For tiled FITS, the overall input image and tiling coordinate system # need to have negative (JPEG-like) parity, but the individal tiles need # to be saved with positive (FITS-like) parity. The upshot is that we # have to vertically flip the data layout in our tile buffers. invert_into_tiles = pio.get_default_vertical_parity_sign() == 1 # TODO: ideally make_maskable_buffer should be a method # on the Image class which could then avoid having to # manually transfer _format. buffer = image.mode.make_maskable_buffer(256, 256) buffer._default_format = image._default_format with progress_bar( total=self.count_populated_positions(), show=cli_progress ) as progress: for ( pos, width, height, image_x, image_y, tile_x, tile_y, ) in self.generate_populated_positions(): if invert_into_tiles: flip_tile_y1 = 255 - tile_y flip_tile_y0 = flip_tile_y1 - height if flip_tile_y0 == -1: flip_tile_y0 = None # with a slice, -1 does the wrong thing by_idx = slice(flip_tile_y1, flip_tile_y0, -1) else: by_idx = slice(tile_y, tile_y + height) iy_idx = slice(image_y, image_y + height) ix_idx = slice(image_x, image_x + width) bx_idx = slice(tile_x, tile_x + width) image.fill_into_maskable_buffer(buffer, iy_idx, ix_idx, by_idx, bx_idx) pio.write_image(pos, buffer) progress.update(1) return self
[docs] def tile_study_image(image, pio, cli_progress=False): """Tile an image as a study, loading the whole thing into memory. Parameters ---------- image : :class:`toasty.image.Image` The image to tile. pio : :class:`toasty.pyramid.PyramidIO` A handle for doing I/O on the tile pyramid cli_progress : optional boolean, defaults False If true, a progress bar will be printed to the terminal. Returns ------- A :class:`StudyTiling` defining the tiling of the image. """ tiling = StudyTiling(image.width, image.height) tiling.tile_image(image, pio, cli_progress=cli_progress) return tiling