Stitching Module

Stitching

Stitching module provides several functions to manage and process expansion microscopy image data, particularly for blending and stitching together image tiles, handling blob deduplication, and creating files compatible with the BigDataViewer (BDV) format.

exm.stitching.stitching.apply_new_ids(img, new_ids, tile_id)

Applies new IDs to blobs in an image based on a mapping from old IDs to new IDs.

This function creates a new image array where each blob’s ID is replaced with a new ID according to a provided mapping. The mapping is a dictionary where keys are tuples of the form (tile_id, old_blob_id), and values are the new blob IDs.

Parameters:

• img (np.ndarray) – The input image with blob IDs to be replaced.

• new_ids (Dict[Tuple[int, int], int]) – A dictionary mapping (tile_id, old_blob_id) to new_blob_id.

• tile_id (int) – The identifier of the current tile being processed.

Returns:

A new image with updated blob IDs.

Return type:

np.ndarray

exm.stitching.stitching.blend(offsets, pictures, indices=None, inverts=None)

Blends a list of volumes tiles together according to specified offsets.

This function takes a list of volumes (tiles) and their corresponding offsets, then composites them into a single volume. Optionally, some tiles can be inverted along specified axes during the blending process.

Parameters:

• offsets (np.ndarray) – An array of offsets for each image tile.

• pictures (List[np.ndarray]) – A list of image tiles to blend.

• indices (Optional[List[int]]) – Optional list of indices specifying the order in which to blend the images. If None, the images are blended in their given order.

• inverts (Optional[List[int]]) – Optional list of axes along which to invert the corresponding image tiles. If None, no inversion is performed.

Returns:

The blended image.

Return type:

np.ndarray

exm.stitching.stitching.blend2(offsets, pictures)

Blends a pair of volume tiles together into a 32-bit composite image, according to specified offsets.

Each tile is blit into separate 16 bits of the output image; one tile is placed in the upper 16 bits and the other in the lower 16 bits. This is used during the deduplication process. The blending accounts for the spatial offsets of each tile to correctly align them in the composite output.

Parameters:

• offsets (np.ndarray) – An array of offsets for each image tile, where each row corresponds to a tile’s offset.

• pictures (List[np.ndarray]) – A list containing exactly two image tiles to blend, each represented as a 3D numpy array.

Returns:

The blended image as a 32-bit numpy array, where the first and second tiles contribute to separate bit planes.

Return type:

np.ndarray

exm.stitching.stitching.deduplicate_blob_ids(tiles, offsets, progress=False)

Generates a mapping from (tile_id, blob_id) pairs to a new unique ID for each blob, merging overlapping blobs.

This function identifies blobs across tiles that should be merged based on overlaps and assigns a new unique ID to the merged set of blobs. Blobs that do not overlap with others will also receive a unique new ID.

Parameters:

• tiles (List[np.ndarray]) – A list of image tiles with blob IDs to be deduplicated.

• offsets (np.ndarray) – An array of offsets indicating the position of each tile.

• progress (bool) – Whether to display progress information.

Returns:

A dictionary mapping each (tile_id, blob_id) to a new unique ID.

Return type:

Dict[Tuple[int, int], int]

exm.stitching.stitching.find_overlapping_blobs(tiles, offsets, progress=False)

Identifies and returns a list of overlapping blobs between tiles, given their offsets.

For each pair of overlapping tiles, this function calculates the blend of the two tiles’ blob indices and determines the overlapping blobs. It returns a list where each element describes an overlap with the following format: [tile_index1, blob_index1, tile_index2, blob_index2, overlap_area].

Parameters:

• tiles (List[np.ndarray]) – A list of image volumes (tiles), each represented as a 3D numpy array.

• offsets (List[Tuple[int, int, int]]) – A list of offsets for each tile, specified as (z, x, y) tuples.

• progress (bool) – If set to True, prints progress of the operation to stdout. Default is False.

Returns:

A list where each element is a list describing an overlapping blob.

Return type:

List[List[int]]

exm.stitching.stitching.get_all_blobs(tiles)

Extracts all unique blob identifiers from a list of tiles.

This function iterates through a list of images (tiles) and extracts all unique (tile index, blob identifier) pairs present across the tiles.

Parameters:

tiles (List[np.ndarray]) – A list of image tiles containing blob identifiers.

Returns:

A set containing all unique (tile index, blob identifier) pairs.

Return type:

Set[Tuple[int, int]]

exm.stitching.stitching.get_bottom_surface(img, thresh=None)

Identifies the bottom surface of the visible cells in an image stack.

This function is a wrapper for get_contact_surface that calculates the surface starting from the bottom of the stack.

Parameters:

• img (np.ndarray) – 3D image stack.

• thresh (Optional[np.ndarray]) – Optional threshold value(s); if not provided, it’s calculated from the image.

Returns:

Tuple containing the image surface, the indices of the surface, and the threshold used.

Return type:

Tuple[np.ndarray, np.ndarray, np.ndarray]

exm.stitching.stitching.get_contact_surface(img, direction=1, thresh=None)

Calculates the surface where cells start being visible in an image stack.

The function moves through the stack in the specified direction, identifying the first z-slice where the intensity is above a calculated threshold, which can also be provided.

Parameters:

• img (np.ndarray) – 3D image stack.

• direction (int) – Direction to process the stack; 1 from bottom, -1 from top.

• thresh (Optional[np.ndarray]) – Optional threshold value(s); if not provided, it’s calculated from the image.

Returns:

Tuple containing the image surface, the indices of the surface, and the threshold used.

Return type:

Tuple[np.ndarray, np.ndarray, np.ndarray]

exm.stitching.stitching.get_offsets(filename)

Extracts the transformation offsets from a BigDataViewer (BDV) or Hierarchical Data Format 5 (H5) XML file. The offsets correspond to the affine transformations applied to the views in the dataset.

Each view’s transformation is represented as a 4x4 matrix, which is accumulated (multiplied) if multiple transformations are present. The function returns these as an array of matrices, one for each view.

Parameters:

filename (str) – Path to the XML file containing view registration and transformation or raw postioning data.

Returns:

An array of corresponding to the offsets of each view as (N,3) array in the ((X,Y,Z),…) order.

Return type:

np.ndarray

exm.stitching.stitching.get_offsets_nd2(filename)

Extracts the transformation offsets from a Nikon Digital Sight 2 (ND2) file. The offsets correspond to the XYZ positions of each acquisition FOV.

The function reads metadata from the ND2 file to extract the positions and returns these as an array of (X, Y, Z) coordinates.

Parameters:

filename (str) – Path to the ND2 file containing acquisition point data.

Returns:

An array of (X, Y, Z) coordinates for each acquisition point.

Return type:

np.ndarray

exm.stitching.stitching.get_tiles_overlaps(transforms, tile_size)

Identifies and returns a set of pairs of indices of transforms that cause their corresponding tiles to overlap.

Given a list of transforms describing the spatial positioning of tiles and the size of the tiles, this function determines all pairs of tiles that overlap each other.

Parameters:

• transforms (List[Tuple[int, int, int]]) – A list of 3-tuples where each tuple represents the transform of a tile (z, x, y offsets).

• tile_size (Tuple[int, int, int]) – A 3-tuple representing the size of the tiles (z, x, y dimensions).

Returns:

A set of unique pairs (as tuples) of indices that have overlapping tiles.

Return type:

Set[Tuple[int, int]]

exm.stitching.stitching.get_top_surface(img, thresh=None)

Identifies the top surface of the visible cells in an image stack.

This function is a wrapper for get_contact_surface that calculates the surface starting from the top of the stack.

Parameters:

• img (np.ndarray) – 3D image stack.

• thresh (Optional[np.ndarray]) – Optional threshold value(s); if not provided, it’s calculated from the image.

Returns:

Tuple containing the image surface, the indices of the surface, and the threshold used.

Return type:

Tuple[np.ndarray, np.ndarray, np.ndarray]

exm.stitching.stitching.make_stitched_dedup_h5(target_file, tiles, trans, new_ids)

Creates a BigDataViewer (BDV) H5 file with deduplicated and stitched tiles.

This function takes a list of image tiles, their transformations (translation offsets), and a dictionary of new IDs for deduplication. It stitches these tiles into a single volume, remaps their blob IDs according to the new_ids, and saves the volume as an H5 file readable by the BigDataViewer.

Parameters:

• target_file (str) – The path to the output H5 file.

• tiles (List[np.ndarray]) – A list of 3D image tiles (numpy arrays).

• trans (List[Tuple[float, float, float]]) – A list of translation offsets for each tile.

• new_ids (dict) – A dictionary mapping old blob IDs to new unique blob IDs.

Returns:

A list of new image tiles with deduplicated IDs.

Return type:

List[np.ndarray]

exm.stitching.stitching.make_stitched_h5(target_file, tiles, transforms)

Creates a BigDataViewer (BDV) H5 file with stitched tiles.

This function stitches a list of image tiles together using the provided transformations and saves the volume as an H5 file readable by BigDataViewer.

Parameters:

• target_file (str) – The path to the output H5 file.

• tiles (List[np.ndarray]) – A list of 3D image tiles (numpy arrays).

• transforms (List[Tuple[float, float, float]]) – A list of translation offsets for each tile.

Return type:

None

exm.stitching.stitching.overlapping(t1, t2, tile_size)

Determines whether the tiles positioned by two transforms will overlap, given the size of the tiles.

This function checks if the spatial positioning of two tiles described by their transform vectors (t1, t2) will result in any overlap between them, considering the dimensions of the tiles.

Parameters:

• t1 (Tuple[int, int, int]) – A 3-tuple representing the transform of the first tile (z, x, y offsets).

• t2 (Tuple[int, int, int]) – A 3-tuple representing the transform of the second tile (z, x, y offsets).

• tile_size (Tuple[int, int, int]) – A 3-tuple representing the size of the tiles (z, x, y dimensions).

Returns:

True if the tiles overlap, False otherwise.

Return type:

bool

Tileset

class exm.stitching.tileset.Tileset(voxel_size, orig_size=None)

A class for managing a set of image tiles from various microscopy formats.

This class supports initializing tilesets from JP2, ND2, BDV H5, and TIFF files, with methods to handle tiles, preview data, and manipulate the intensity and offsets of the tiles. It is designed to work with image data in a voxel grid format, with methods provided for scaling, previewing, and exporting the data.

Attributes:

• tiles (list): A list of Tile objects that make up the tileset.

• voxel_size (tuple): The size of a voxel in the dataset.

• original_xyz_size (tuple): Original size of the dataset in XYZ dimensions.

• nd2 (ND2Reader) : An ND2Reader object for ND2 file operations, if applicable.

Methods:

• init_from_jp2(files_list, offsets_file, downscale): Initialize from a list of JP2 files.

• init_from_nd2(nd2): Initialize from an ND2 file.

• init_from_bdv(filename): Initialize from a BDV H5 file.

• init_from_h5(filename, downscale, progress): Initialize from a non-BDV H5 file.

• init_from_tiff_files(filelist, downscale, progress): Initialize from a list of TIFF files.

• scale_offset(scale): Scale all tile offsets by a given factor.

• __getitem__(i): Access individual tiles using the [] operator.

• get_slice(fovnum, z): Get a specific Z-slice from a specific FOV.

• preview_nd2(z_layer, down_sample): Preview a slice of the ND2 file.

• get_centroids(): Calculate and return centroids of tiles.

• create_blank_h5bdv(target_file, tile_size, dtype, noise_scale): Create a blank H5/BDV file with random noise.

• load_tile(fovnums, downscale, method): Load image data for a specified tile from the ND2 file.

• load_all(downscale, method): Load image data for all tiles from the ND2 file.

• find_intensity_scale(progress): Find the 1st and 99th percentile intensity values.

• scale_intensity(p1p99, progress): Scale the intensity of the images in the tileset.

• write_into_h5bdv(filename, dtype): Write the tileset into an H5BDV file format.

• show_slice(zslice, down_sample): Show a single Z-slice of the stitched tileset.

• update_offsets(xml_file, scale_factor): Update the offsets for the tiles based on an XML file.

• produce_output_volume(): Create and return a stitched volume from the tiles.

• local_to_global(coords): Transform local tile coordinates to global coordinates.

• dedup_segmentation_ids(progress): De-duplicate segmentation IDs across the tileset.

Note:

Some methods include TODOs for future improvements and considerations for handling specific data formats and conditions.

create_blank_h5bdv(target_file, tile_size=None, dtype=<class 'numpy.int16'>, noise_scale=200.0)

Creates a H5/BDV file pair. It is hard to add new tiles to an existing file so instead the H5BDV files are created with all their tiles filled with random noise. Why noise instead of zeros? Well Fiji’s BigStitcher does not like it when some tiles are totally empty and saturates the other tiles. That’s why we fill it with noise. To deactivate that put noise_scale at \(0.0\)

Parameters:

• target_file (str) – Must be the path to the desired H5 file. An XML file will be created at the same place with the same name and the extension changed.

• tile_size (tuple(int)) – The size we want the tiles to be.

• dtype (numpy.dtype) – The type to initialize the arrays with.

• noise_scale (float) – The intensity of the white noise the tiles are initialized with.

dedup_segmentation_ids(progress=False)

If the tileset is a set of segmented FOVs, this function replaces the tiles IDs by identifiers that are unique accross the dataset

find_intensity_scale(progress=False)

Finds the \(1^{st}\) and \(99^{th}\) percentile of intensity values in the dataset. Can be long.

get_centroids()

Returns a list of centroids as an array of XYZ coordinates. It can take a while. As this function internally calls produce_output_volume() it is recommended to scale down the tiles first

get_slice(fovnum, z)

Returns a specific Z-slice from a specific FOV. That’s a fast function that does not load the whole tile

init_from_bdv(filename)

Loads a tileset from a BDV H5 (BigDataViewer format). The argument should be either the .h5 or the .xml file

init_from_h5(filename, downscale=[1, 1, 1], progress=False)

Loads a tileset from a non-BDV H5, will expect transforms from another source

init_from_jp2(files_list, offsets_file, downscale=[1, 1, 1])

Initializes the tileset from a list of JP2 files, each representing a single tile, and an offset files, which is .npy file storing a 2D array. The first row of this array gives the original pixel size of the tile and the other rows the offset of the tiles.

init_from_nd2(nd2)

Initializes the tileset from a ND2 file. As these files can be slow to load it does not actually read the content of the tiles, but only reads metadata and will lazily load the necessary information on calls to load_fov or preview_nd2

init_from_tiff_files(filelist, downscale=[1, 1, 1], progress=False)

Init tiles from a list of TIFF files but will not initiialize offsets

load_all(downscale=[1, 1, 1], method='fast')

Loads all tiles from the ND2 file. Before that call, tiles exist but are placeholders without image data Method can be “fast” or “safe”. The fast method makes the loading faster but makes some assumptions about the data format Note that the downscale parameter is in the axis order of the tiles, typically ZYX

load_tile(fovnums, downscale=[1, 1, 1], method='fast')

Loads a tile from the ND2 file. Before that call, tiles exist but are placeholders without image data Method can be “fast” or “safe”. The fast method makes the loading faster but makes some assumptions about the data format

local_to_global(coords)

(See the notebook Tileset 02) Transforms an array of coordinates that are local to tilesets into an array of coordinates into the global coorddinates system, that is, into the reconstructed volume.

Parameters:

coords (numpy.ndarray) – Array of shape \((n,4)\) describing \(n\) points. The rows are ordered in the \((i, x, y, z)\) order, \(i\) being the index of the tile the point belongs to and \((x, y, z)\) its coordinates in \(\mu m\) within the tile.

Returns:

Array of shape \((n,3)\) containing the transformed coordinates in \(\mu m\)

Return type:

numpy.ndarray

preview_nd2(z_layer=0, down_sample='auto')

Previews a slice of a stitched version of the ND2 without loading all the tiles

produce_output_volume(data_format=None)

Creates and returns a stitched volume. Careful, if you never downsampled the tiles it easily fills up all memory.

scale_intensity(p1p99=None, progress=False)

Scales the intensities over the whole dataset and convert tiles to 8bpp.

Parameters:

• p1p99 (tuple(int, int)) – lowest and highest value between which to scale. If None, the function will calculate the \(1^{st}\) and \(99^{th}\) percentiles of intensities and use these.

• progress (bool) – If True, displays the progress of the percentile calculation (which is long)

scale_offset(scale=(1, 1, 1))

Scales all offsets by a factor. Expects a X, Y, Z order in the argument

show_slice(zslice, down_sample='auto')

Shows a single \(Z\)-slice of the stitched tileset.

Parameters:

• zslice (int) – \(z\)-index to display

• down_sample (list[int, int], optional) – integer factor by which to downsample each dimension of the result. If set to "auto", will make the biggest dimension as close as possible to \(1500\)

update_offsets(xml_file, scale_factor=None)

Replaces the current offsets by those in a specific XML file. Does not touch the pixels information.

Parameters:

• xml_file (str) – the file to read new offsets from.

• scale_factor (numpy.ndarray) – the vector by which to multiply offsets if necessary, defaults to None. Usually unnecessary

write_into_h5bdv(filename, dtype=<class 'numpy.int16'>)

Exports the tileset as a H5BDV file format. It will create a .h5 file containing the image information and a XML file containing the offsets.

Tiles

class exm.stitching.tileset.Tile(fovnum)

A Tile represents a FOV in the dataset. It contains the FOV number in the dataset, methods to retrieve the original image and to process it, with the ability to cache various intermediate results.

Initializes a Tile object representing a Field Of View (FOV) with an identifier.

apply_mask()

Applies a previously computed mask to the tile image, setting unmasked pixels to zero.

fast_nd2_read_downscale(nd2, downscale)

Reads the FOV specified by fovnum from the specified nd2 file. This is a fast function that needs a several assumptions to be true to work. Notably dtype needs to be int16.

Parameters:

downscale (list(int)) – steps to apply while reading the image, in \((z, y, x)\) order

find_intensity_scale()

Returns the value of the first and \(99^{th}\) percentile of intensities in the tile

from_array(array)

Initializes the tile’s image from a given array.

from_jp2(filename, downscale=[1, 1, 1])

Loads the tile image from a JPEG 2000 file with optional downscaling.

from_nd2(nd2, downsample=[1, 1, 1])

Loads a tile from a specific FOV

Parameters:

• nd2 (str) – can be a filename or an existing ND2Reader

• fovnum (int) – the FOV (along the \(v\) axis) to extract

• downsample (tuple(int)) – a tuple of \(3\) items describing a stride to apply in the \((x, y, z)\) axes, resulting in a downscale along these axes

Returns:

An initialized tile with image and offset information

Return type:

tileset.Tile

from_tiff(filename, cached=True)

Loads the tile image from a TIFF file.

save_as_jp2(filename)

Saves the tile image as a JPEG 2000 file at the specified path. (Note: method currently lacks implementation details).

save_as_tiff(filename)

Saves the tile image as a TIFF file at the specified path.

scale_intensity(p1, p99, convert_to_8bpp=True)

Scales the image intensities between the provided percentiles, optionally converting it to 8 bits per pixel.