paprica.stitcher

Submodule containing classes and functions relative to stitching.

With this submodule the user can stitch a previously parsed dataset, typically the autofluorescence channel:

>>> import paprica
>>> tiles_autofluo = paprica.parser.tileParser(path_to_autofluo, frame_size=1024, overlap=25)
>>> stitcher = paprica.stitcher.tileStitcher(tiles_autofluo)
>>> stitcher.compute_registration_fast()

Others channel can then easily stitched using the previous one as reference:

>>> tiles_signal = paprica.parser.tileParser(path_to_data, frame_size=1024, overlap=25)
>>> stitcher_channel = paprica.stitcher.channelStitcher(stitcher, tiles_autofluo, tiles_signal)
>>> stitcher_channel.compute_rigid_registration()

Doing that each tile in the second data set will be registered to the corresponding autofluorescence tile and then their spatial position will be adjusted.

WARNING: when stitching, the expected overlap must be HIGHER than the real one. To enforce this, a margin of 20% is automatically taken (this margin can be set lower by the user for speed improvement). In order to get the best stitching quality it requires to have a good estimate of the overlap, hence why the full volume is not considered.

This submodule also contains a class for merging and reconstructing the data. It was intended to be used at lower resolution for atlasing. The generated data can quickly become out of hands, use with caution!

By using this code you agree to the terms of the software license agreement.

© Copyright 2020 Wyss Center for Bio and Neuro Engineering – All rights reserved

paprica.stitcher._compute_shift(reference_image, moving_image)

Backbone function to compute the registration and the registration error used for the global optimisation. This function can be replaced by experienced user to use their own registration and error estimation functions.

Parameters

• reference_image (array) – Reference image.

• moving_image (array) – Image to register. Must be same dimensionality as reference_image.

Returns

• d (array_like) – registration parameters found

• e (float) – error estimation for the registration (the higher the error the higher the registration uncertainty)

paprica.stitcher._get_masked_proj_shifts(proj1, proj2, threshold, upsample_factor=1)

This function computes shifts from max-projections on overlapping areas with mask on brightest area. It uses the phase cross-correlation to compute the shifts.

Parameters

• proj1 (list[ndarray]) – max-projections for tile 1

• proj2 (list[ndarray]) – max-projections for tile 2

• upsample_factor (float) – upsampling_factor for estimating the maximum phase cross-correlation position

Returns

_ – shifts in (x, y, z) and error measure (0=reliable, 1=not reliable)

Return type

array_like

paprica.stitcher._get_max_proj_apr(apr, parts, patch, patch_yx=None, plot=False)

Compute maximum projection on 3D APR data.

Parameters

• apr (pyapr.APR) – apr tree

• parts (pyapr.ParticlData) – apr particle

• patch (pyapr.ReconPatch) – patch for computing the projection only on the overlapping area.

• plot (bool) – control data plotting

Returns

_ – maximum intensity projection in each 3 dimension.

Return type

list[ndarray]

paprica.stitcher._get_proj_shifts(proj1, proj2)

This function computes shifts from max-projections on overlapping areas. It uses the phase cross-correlation to compute the shifts.

Parameters

• proj1 (list[ndarray]) – max-projections for tile 1

• proj2 (list[ndarray]) – max-projections for tile 2

• upsample_factor (float) – upsampling_factor for estimating the maximum phase cross-correlation position

Returns

_ – shifts in (x, y, z) and error measure (0=reliable, 1=not reliable)

Return type

array_like

class paprica.stitcher.baseStitcher(tiles, overlap_h: (<class 'int'>, <class 'float'>), overlap_v: (<class 'int'>, <class 'float'>))

Bases: object

Base class for stitching multi-tile data.

__init__(tiles, overlap_h: (<class 'int'>, <class 'float'>), overlap_v: (<class 'int'>, <class 'float'>))

Constructor for the baseStitcher class.

Parameters

• tiles (tileParser) – tileParser object containing the dataset to stitch.

• overlap_h (float) – expected horizontal overlap in %

• overlap_v (float) – expected vertical overlap in %

Return type

None

_load_max_projs(path)

Load the maximum intensity projection previously stored.

Parameters

path (str) – path to load the maximum intensity projection from. If None then default to max_projs folder in the acquisition folder.

Return type

None

_precompute_max_projs(progress_bar=True)

Precompute max-projections for loading the data only once during the stitching.

Return type

None

_process_GRAY_for_display(u)

Process RGB data for correctly displaying it.

Parameters

u (ndarray) – RGB data

Returns

data_to_display – RGB data displayable with correct contrast and colors.

Return type

ndarray

_process_RGB_for_display(u)

Process RGB data for correctly displaying it.

Parameters

u (ndarray) – RGB data

Returns

data_to_display – RGB data displayable with correct contrast and colors.

Return type

ndarray

_reconstruct_x_slice(x=None, n_proj=0, downsample=1, color=False, debug=False, plot=True, progress_bar=True)

Reconstruct and merge the sample at a given position x.

Parameters

• x (int) – reconstruction location in x

• n_proj (int (default: 0)) – Number of planes to perform the maximum intensity projection.

• dim (int (default: 0)) – Dimension of the reconstruction, e.g. 0 will be [y, x] plane (orthogonal to z).

• downsample (int (default: 1)) – Downsample factor for the reconstruction. Must be in [1, 2, 4, 8, 16, 32].

• color (bool (default: False)) – Option to reconstruct with checkerboard color pattern. Useful to identify doubling artifacts.

• debug (bool (default: False)) – Option to add a white square for each tile, making it easy to see overlapping areas.

• plot (bool (default: True)) – Define if the function plots the results with Matplotlib or just returns an array.

Returns

merged_data – Merged frame at position x.

Return type

ndarray

_reconstruct_y_slice(y=None, n_proj=0, downsample=1, color=False, debug=False, plot=True, progress_bar=True)

Reconstruct and merge the sample at a given position y.

Parameters

• y (int) – reconstruction location in y

• n_proj (int (default: 0)) – Number of planes to perform the maximum intensity projection.

• dim (int (default: 0)) – Dimension of the reconstruction, e.g. 0 will be [y, x] plane (orthogonal to z).

• downsample (int (default: 1)) – Downsample factor for the reconstruction. Must be in [1, 2, 4, 8, 16, 32].

• color (bool (default: False)) – Option to reconstruct with checkerboard color pattern. Useful to identify doubling artifacts.

• debug (bool (default: False)) – Option to add a white square for each tile, making it easy to see overlapping areas.

• plot (bool (default: True)) – Define if the function plots the results with Matplotlib or just returns an array.

Returns

merged_data – Merged frame at position y.

Return type

ndarray

_reconstruct_z_slice(z=None, n_proj=0, downsample=1, color=False, debug=False, plot=True, seg=False, progress_bar=True)

Reconstruct and merge the sample at a given depth z.

Parameters

• z (int) – reconstruction depth (vary with downsample)

• n_proj (int (default: 0)) – Number of planes to perform the maximum intensity projection.

• dim (int (default: 0)) – Dimension of the reconstruction, e.g. 0 will be [y, x] plane (orthogonal to z).

• downsample (int (default: 1)) – Downsample factor for the reconstruction. Must be in [1, 2, 4, 8, 16, 32].

• color (bool (default: False)) – Option to reconstruct with checkerboard color pattern. Useful to identify doubling artifacts.

• debug (bool (default: False)) – Option to add a white square for each tile, making it easy to see overlapping areas.

• plot (bool (default: True)) – Define if the function plots the results with Matplotlib or just returns an array.

• seg (bool (default: False)) – Option to also reconstruct the segmentation. Only works with dim=0

Returns

merged_data – Merged frame at depth z.

Return type

ndarray

_regularize(reg, rel)

Remove too large displacement and replace them with expected one with a large uncertainty.

_save_max_projs()

Save the computed maximum intensity projection on persistent memory. This is useful to recompute the registration directly from the max. proj. but only works if the overlaps are kept the same.

Return type

None

activate_mask(threshold)

Activate the masked cross-correlation for the displacement estimation. Pixels above threshold are not taken into account.

Parameters

threshold (int) – threshold for the cross-correlation mask as a percentage of pixel to keep (e.g. 95 will create a mask removing the 5% brightest pixels).

activate_segmentation(segmenter)

Activate the segmentation. When a tile is loaded it is segmented before the stitching is done.

Parameters

segmenter (tileSegmenter) – segmenter object for segmenting each tile.

deactivate_mask()

Deactivate the masked cross-correlation and uses a classical cross correlation.

deactivate_segmentation()

Deactivate tile segmentation.

load_database(path=None, force=False)

Save database at the given path. The database must be built before calling this method.

Parameters

path (string) – path to save the database.

reconstruct_slice(loc=None, n_proj=0, dim=0, downsample=1, color=False, debug=False, plot=True, seg=False, progress_bar=True)

Reconstruct whole sample 2D section at the given location and in a given dimension. This function can also reconstruct a maximum intensity projection if n_proj>0.

Parameters

• loc (int (default: middle of the sample)) – Position of the plane where the reconstruction should be done. The location varies depending on the downsample parameter and should be adapted.

• n_proj (int (default: 0)) – Number of planes to perform the maximum intensity projection.

• dim (int (default: 0)) – Dimension of the reconstruction, e.g. 0 will be [y, x] plane (orthogonal to z).

• downsample (int (default: 1)) – Downsample factor for the reconstruction. Must be in [1, 2, 4, 8, 16, 32].

• color (bool (default: False)) – Option to reconstruct with checkerboard color pattern. Useful to identify doubling artifacts.

• debug (bool (default: False)) – Option to add a white square for each tile, making it easy to see overlapping areas.

• plot (bool (default: True)) – Define if the function plots the results with Matplotlib or just returns an array.

• seg (bool (default: False)) – Option to also reconstruct the segmentation. Only works with dim=0

Returns

_ – Array containing the reconstructed data.

Return type

ndarray

reconstruct_z_color(z=None, n_proj=10, downsample=1, debug=False, plot=True, progress_bar=True)

Reconstruct and merge the sample at a given depth z.

Parameters

• z (int) – reconstruction depth

• downsample (int) – downsample for reconstruction (must be a power of 2)

• debug (bool) – if true the border of each tile will be highlighted

Returns

merged_data – Merged frame at depth z.

Return type

ndarray

save_database(path=None)

Save database at the given path. The database must be built before calling this method.

Parameters

path (string) – path to save the database.

set_regularization(reg_x, reg_y, reg_z)

Set the regularization for the stitching to prevent aberrant displacements.

Parameters

• reg_x (int) – if the horizontal displacement computed in the pairwise registration for any tile is greater than reg_x (in pixel unit) then the expected displacement (from motor position) is taken.

• reg_y (int) – if the horizontal displacement computed in the pairwise registration for any tile is greater than reg_z (in pixel unit) then the expected displacement (from motor position) is taken.

• reg_z (int) – if the horizontal displacement computed in the pairwise registration for any tile is greater than reg_z (in pixel unit) then the expected displacement (from motor position) is taken.

Return type

None

class paprica.stitcher.channelStitcher(stitcher, ref, moving)

Bases: paprica.stitcher.baseStitcher

Class used to perform the stitching between different channels. The registration must be performed first a single channel (typically auto-fluorescence)

The stitching is performed between each corresponding tile and the relative displacement is added to the previously determined stitching parameters.

The number and position of tile should matched for bot dataset.

__init__(stitcher, ref, moving)

Constructor for the channelStitcher class.

Parameters

• stitcher (tileStitcher) – tileStitcher object with the multitile registration parameters

• tiles_stitched (tileParser) – tiles corresponding to the stitcher

• tiles_channel (tileParser) – tiles to be registered to tiles_stitched

_update_database(row, col, d)

Update database after the registration.

Parameters

• row (int) – row number

• col (int) – col number

• d (int) – computed displacement

compute_rigid_registration(progress_bar=True)

Compute the rigid registration between each pair of tiles across different channels.

Return type

None

set_lim(x_begin=None, x_end=None, y_begin=None, y_end=None, z_begin=None, z_end=None)

Define spatial limits to compute the maximum intensity projection.

Parameters

• x_begin (int) –

• x_end (int) –

• y_begin (int) –

• y_end (int) –

• z_begin (int) –

• z_end (int) –

Return type

None

paprica.stitcher.max_sum_over_single_max(reference_image, moving_image, d)

This function is a reliability metric which works well for sparse data. It computes the 99 percentile of the sum of reference and shifted image divided by twice the 99 percentile of the reference image.

Parameters

• reference_image (ndarray) – 2D array of the reference image

• moving_image (ndarray) – 2D array of the moving image image

• d (array_like) – registration parameters

Returns

e – error estimation of the registration. The lower e, the more reliable the registration.

Return type

float

paprica.stitcher.mse(reference_image, moving_image, d)

Normalized root mean square error.

Parameters

• reference_image (ndarray) – 2D array of the reference image

• moving_image (ndarray) – 2D array of the moving image image

• d (array_like) – registration parameters

Returns

_ – error estimation of the registration. The lower e, the more reliable the registration.

Return type

float

paprica.stitcher.phase_cross_correlation(reference_image, moving_image, upsample_factor=1, return_error=True)

Phase cross correlation. Because skimage function compute the NORMAL cross correlation to estimate the shift I modified it to compute the TRUE phase cross correlation, as per the standard definition.

Parameters

• reference_image (array) – Reference image.

• moving_image (array) – Image to register. Must be same dimensionality as reference_image.

• upsample_factor (int, optional) – Upsampling factor. Images will be registered to within 1 / upsample_factor of a pixel. For example upsample_factor == 20 means the images will be registered within 1/20th of a pixel. Default is 1 (no upsampling). Not used if any of reference_mask or moving_mask is not None.

• return_error (bool, optional) – Returns error and phase difference if on, otherwise only shifts are returned. Has noeffect if any of reference_mask or moving_mask is not None. In this case only shifts is returned.

Returns

• shifts (ndarray) – Shift vector (in pixels) required to register moving_image with reference_image. Axis ordering is consistent with numpy (e.g. Z, Y, X)

• error (float) – Translation invariant normalized RMS error between reference_image and moving_image.

• phasediff (float) – Global phase difference between the two images (should be zero if images are non-negative).

paprica.stitcher.phase_cross_correlation_cv(reference_image, moving_image)

Compute openCV to compute the phase cross correlation. It is around 16 times faster than the implementation using numpy FFT (same as skimage).

Parameters

• reference_image (array) – Reference image.

• moving_image (array) – Image to register. Must be same dimensionality as reference_image.

Returns

• shifts (ndarray) – Shift vector (in pixels) required to register moving_image with reference_image. Axis ordering is consistent with numpy (e.g. Z, Y, X)

• error (float) – Peak response (see opencv description here: https://docs.opencv.org/4.5.3/d7/df3/group__imgproc__motion.html#ga552420a2ace9ef3fb053cd630fdb4952)

class paprica.stitcher.tileMerger(tiles, database)

Bases: object

Class to merge tiles and create a stitched volume. Typically used at a lower resolution for registering the sample to an Atlas.

__init__(tiles, database)

Constructor for the tileMerger class.

Parameters

• tiles (tileParser) – tileParser object containing the dataset to merge.

• database ((pd.DataFrame, string)) – database or path to the database containing the registered tile position

• n_planes (int) – number of planes per files.

Return type

None

_find_if_lazy()

Function to test if all tile can be lazy loaded.

Returns

_ – True if the tile can be lazy loaded, false if not.

Return type

bool

_get_n_planes()

Load a tile and check the number of planes per tile.

Returns

_ – Number of planes per tile;

Return type

int

_get_nx()

Compute the merged array size for x dimension.

Returns

_ – x size for merged array

Return type

int

_get_ny()

Compute the merged array size for y dimension.

Returns

_ – y size for merged array

Return type

int

_get_nz()

Compute the merged array size for y dimension.

Returns

_ – y size for merged array

Return type

int

_initialize_merged_array()

Initialize the merged array in accordance with the asked downsampling.

Return type

None

_initialize_merged_segmentation()

Initialize the merged array in accordance with the asked downsampling.

Return type

None

crop(background=0, xlim=None, ylim=None, zlim=None)

Add a black mask around the brain (rather than really cropping which makes the overlays complicated in a later stage).

Parameters

• background (int) – constant value to replace the cropped area with.

• xlim (array_like) – x limits for cropping

• ylim (array_like) – y limits for cropping

• zlim (array_like) – z limits for cropping

Return type

None

equalize_hist(method='opencv')

Perform histogram equalization to improve the contrast on merged data. Both OpenCV (only 2D) and Skimage (3D but 10 times slower) are available.

Parameters

method (string) – method for performing histogram equalization among ‘skimage’ and ‘opencv’.

Return type

None

merge_additive(reconstruction_mode='constant', tree_mode='mean', progress_bar=True)

Perform merging with a mean algorithm for overlapping areas. Maximum merging should be preferred to avoid integer overflowing and higher signals on the overlapping areas.

Parameters

mode (string) – APR reconstruction type among (‘constant’, ‘smooth’, ‘level’)

Return type

None

merge_max(reconstruction_mode='constant', tree_mode='mean', debug=False, progress_bar=True)

Perform merging with a maximum algorithm for overlapping areas.

Parameters

• mode (string) – APR reconstruction type among (‘constant’, ‘smooth’, ‘level’)

• debug (bool) – add white border on the edge of each tile to see where it was overlapping.

Return type

None

merge_segmentation(reconstruction_mode='constant', tree_mode='max', debug=False, progress_bar=True)

Perform merging with a maximum algorithm for overlapping areas.

Parameters

• mode (string) – APR reconstruction type among (‘constant’, ‘smooth’, ‘level’)

• debug (bool) – add white border on the edge of each tile to see where it was overlapping.

Return type

None

set_downsample(downsample)

Set the downsampling value for the merging reconstruction.

Parameters

downsample (int) – downsample factor

Return type

None

class paprica.stitcher.tileStitcher(tiles, overlap_h: (<class 'int'>, <class 'float'>), overlap_v: (<class 'int'>, <class 'float'>))

Bases: paprica.stitcher.baseStitcher

Class used to perform the stitching. The stitching is performed in 4 steps:

1. The pairwise registration parameters of each neighboring tile is computed on the max-projection

2. A sparse graph (edges = tiles and vertex = registration between neighboring tiles) is constructed to store the registration parameters (displacements and reliability)

3. The sparse graph is optimized to satisfy the constraints (every loop in the graph should sum to 0) using the maximum spanning tree on the reliability estimation.

4. The maximum spanning tree is parsed to extract optimal tile positions solution.

The beauty of this method is that it scales well with increasing dataset sizes and because the final optimization is very fast and does not require to reload the data.

__init__(tiles, overlap_h: (<class 'int'>, <class 'float'>), overlap_v: (<class 'int'>, <class 'float'>))

Constructor for the tileStitcher class.

Parameters

• tiles (tileParser) – tileParser object containing the dataset to stitch.

• overlap_h (float) – expected horizontal overlap in %

• overlap_v (float) – expected vertical overlap in %

_build_database()

Build the database for storing the registration parameters. This method needs to be called after the registration map has been produced.

Return type

None

_build_sparse_graphs()

Build the sparse graph from the reliability and (row, col). This method needs to be called after the pair-wise registration has been performed for all neighbors pair.

Return type

None

_compute_east_registration(apr_1, parts_1, apr_2, parts_2)

Compute the registration between the current tile and its eastern neighbor.

Parameters

• u (list) – current tile

• v (list) – neighboring tile

Return type

None

_compute_registration_old()

Compute the pair-wise registration for all tiles. This implementation loads the data twice and is therefore not efficient.

_compute_south_registration(apr_1, parts_1, apr_2, parts_2)

Compute the registration between the current tile and its southern neighbor.

Parameters

• u (list) – current tile

• v (list) – neighboring tile

Return type

None

_get_ind(ind_from, ind_to)

Returns the ind in the original graph which corresponds to (ind_from, ind_to) in the minimum spanning tree.

Parameters

• ind_from (int) – starting node in the directed graph

• ind_to (int) – ending node in the directed graph

Returns

ind – corresponding ind in the original graph

Return type

int

_optimize_sparse_graphs()

Optimize the sparse graph by computing the minimum spanning tree for each direction (H, D, V). This method needs to be called after the sparse graphs have been built.

Return type

None

_print_info()

Display stitching result information.

_produce_registration_map()

Produce the registration map where reg_rel_map[d, row, col] (d = H,V,D) is the relative tile position in pixel from the expected one. This method needs to be called after the optimization has been done.

Return type

None

compute_expected_registration()

Compute the expected registration if the expected overlap are correct.

compute_registration(on_disk=False, progress_bar=True)

Compute the pair-wise registration for all tiles. This implementation loads the data once by precomputing the max-proj and is therefore efficient.

compute_registration_from_max_projs(path=None)

Compute the registration directly from the max-projections. Max-projections must have been computed before.

dump_stitcher(path)

Use dill to store a tgraph object.

Parameters

path (string) – path to save the database.

plot_graph(annotate=False)

Plot the graph for each direction (H, D, V). This method needs to be called after the graph optimization.

Parameters

annotate (bool) – control if annotation are drawn on the graph

plot_min_trees(annotate=False)

Plot the minimum spanning tree for each direction (H, D, V). This method needs to be called after the graph optimization.

Parameters

annotate (bool) – control if annotation are drawn on the graph

plot_registration_map()

Display the registration map using matplotlib.

plot_stitching_info()

Plot pair-wise registration error for each axis [H, V, D].

Returns

• rel_map (array) – error matrix

• d_map (array) – shift matrix

set_overlap_margin(margin)

Modify the overlaping area size. If the overlaping area is smaller than the true one, the stitching can’t be performed properly. If the overlaping area area is more than twice the size of the true one it will also fail (due to the circular FFT in the phase cross correlation).

Parameters

margin (float) – safety margin in % to take the overlaping area.

Return type

None

set_z_range(z_begin, z_end)

Set a range of depth fo computing the stitching.

Parameters

• z_begin (int) – first depth to be included in the max-proj

• z_end (int) – last depth to be included in the max-proj

Return type

None