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)[source]
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)[source]
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)[source]
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)[source]
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'>))[source]
Bases:
object
Base class for stitching multi-tile data.
- __init__(tiles, overlap_h: (<class 'int'>, <class 'float'>), overlap_v: (<class 'int'>, <class 'float'>))[source]
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)[source]
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)[source]
Precompute max-projections for loading the data only once during the stitching.
- Return type
None
- _process_GRAY_for_display(u)[source]
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)[source]
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)[source]
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)[source]
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)[source]
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)[source]
Remove too large displacement and replace them with expected one with a large uncertainty.
- _save_max_projs()[source]
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)[source]
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)[source]
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()[source]
Deactivate the masked cross-correlation and uses a classical cross correlation.
- load_database(path=None, force=False)[source]
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)[source]
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)[source]
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)[source]
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)[source]
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)[source]
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)[source]
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)[source]
Update database after the registration.
- Parameters
row (int) – row number
col (int) – col number
d (int) – computed displacement
- paprica.stitcher.max_sum_over_single_max(reference_image, moving_image, d)[source]
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)[source]
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)[source]
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 exampleupsample_factor == 20
means the images will be registered within 1/20th of a pixel. Default is 1 (no upsampling). Not used if any ofreference_mask
ormoving_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
ormoving_mask
is not None. In this case only shifts is returned.
- Returns
shifts (ndarray) – Shift vector (in pixels) required to register
moving_image
withreference_image
. Axis ordering is consistent with numpy (e.g. Z, Y, X)error (float) – Translation invariant normalized RMS error between
reference_image
andmoving_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)[source]
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
withreference_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)[source]
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)[source]
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()[source]
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()[source]
Load a tile and check the number of planes per tile.
- Returns
_ – Number of planes per tile;
- Return type
int
- _get_nx()[source]
Compute the merged array size for x dimension.
- Returns
_ – x size for merged array
- Return type
int
- _get_ny()[source]
Compute the merged array size for y dimension.
- Returns
_ – y size for merged array
- Return type
int
- _get_nz()[source]
Compute the merged array size for y dimension.
- Returns
_ – y size for merged array
- Return type
int
- _initialize_merged_array()[source]
Initialize the merged array in accordance with the asked downsampling.
- Return type
None
- _initialize_merged_segmentation()[source]
Initialize the merged array in accordance with the asked downsampling.
- Return type
None
- crop(background=0, xlim=None, ylim=None, zlim=None)[source]
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')[source]
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)[source]
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)[source]
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)[source]
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
- class paprica.stitcher.tileStitcher(tiles, overlap_h: (<class 'int'>, <class 'float'>), overlap_v: (<class 'int'>, <class 'float'>))[source]
Bases:
paprica.stitcher.baseStitcher
Class used to perform the stitching. The stitching is performed in 4 steps:
The pairwise registration parameters of each neighboring tile is computed on the max-projection
A sparse graph (edges = tiles and vertex = registration between neighboring tiles) is constructed to store the registration parameters (displacements and reliability)
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.
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'>))[source]
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()[source]
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()[source]
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)[source]
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()[source]
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)[source]
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)[source]
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()[source]
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
- _produce_registration_map()[source]
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()[source]
Compute the expected registration if the expected overlap are correct.
- compute_registration(on_disk=False, progress_bar=True)[source]
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)[source]
Compute the registration directly from the max-projections. Max-projections must have been computed before.
- dump_stitcher(path)[source]
Use dill to store a tgraph object.
- Parameters
path (string) – path to save the database.
- plot_graph(annotate=False)[source]
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)[source]
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_stitching_info()[source]
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)[source]
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