References

args = parser.parse_args() module-attribute

params = yaml.safe_load(file) module-attribute

parser = argparse.ArgumentParser(description=f'muvis-align') module-attribute

pipeline = Pipeline(params) module-attribute

muvis_align

MVSRegistration

NAPARI_PROJECT_TEMPLATE = 'ui/project_template.yaml' module-attribute

default_chunk_size = 1024 module-attribute

default_mappings_name = 'mappings.json' module-attribute

default_ome_zarr_version = '0.5' module-attribute

metrics_name = 'metrics.json' module-attribute

metrics_tabular_name = 'mappings.csv' module-attribute

original_positions_name = 'positions_original.pdf' module-attribute

prereg_mappings_name = 'prereg_mappings.csv' module-attribute

registered_positions_name = 'positions_registered.pdf' module-attribute

tiff_extension = '.ome.tiff' module-attribute

version = '0.3.0' module-attribute

zarr_extension = '.ome.zarr' module-attribute

MVSRegistration

Source code in src\muvis_align\MVSRegistration.py

class MVSRegistration:
    def __init__(self, operation='register', label='', input_path=None, output_path=None,
                 source_metadata={}, extra_metadata={},
                 global_rotation=None, global_center=None,
                 overwrite=True, clear=False, ui='', verbose=False, debug=False):
        self.state = RegState.UNINIT
        self.sims = []
        self.sources = []
        self.metrics = {}
        if input_path is not None:
            self.init(operation=operation, label=label, input_path=input_path, output_path=output_path,
                      source_metadata=source_metadata, extra_metadata=extra_metadata,
                      global_rotation=global_rotation, global_center=global_center,
                      overwrite=overwrite, clear=clear, ui=ui, verbose=verbose, debug=debug)

    def is_initialised(self):
        return self.state.value >= RegState.INIT.value

    def is_pairs_registered(self):
        return self.state.value >= RegState.PAIRS_REG.value

    def is_global_registered(self):
        return self.state.value >= RegState.GLOBAL_REG.value

    def init_params(self, params_general, params, label='', input_path=None, global_rotation=None, global_center=None):
        self.params_general = params_general
        self.params = params
        self.input_params = params.get('input')
        if isinstance(self.input_params, (str, list)):
            self.input_params = {'path': self.input_params}
        if input_path is None:
            input_path = self.input_params.get('path')
        self.output_params = params.get('output')
        if isinstance(self.output_params, str):
            self.output_params = {'path': self.output_params}
        self.preprocess_params = params.get('preprocessing', {})
        self.register_params = params.get('registration', {})
        self.fusion_params = params.get('fusion', {})

        return self.init(operation=params.get('operation'), label=label, input_path=input_path,
                         input_labels=self.input_params.get('labels'),
                         output_path=self.output_params.get('path'),
                         source_metadata=self.input_params.get('source_metadata', {}),
                         extra_metadata=self.input_params.get('extra_metadata', {}),
                         global_rotation=global_rotation, global_center=global_center,
                         overwrite=params_general.get('overwrite', False), clear=params_general.get('clear', False),
                         ui=params_general.get('ui', ''),
                         verbose=params_general.get('verbose', False), debug=params_general.get('debug', False))

    def init(self, operation='', label='', input_path=None, input_labels=None, output_path=None,
             source_metadata={}, extra_metadata={}, global_rotation=None, global_center=None,
             overwrite=True, clear=False, ui='', verbose=False, debug=False):
        self.overwrite = overwrite
        self.clear = clear
        self.ui = ui
        self.verbose = verbose
        self.debug = debug
        self.logging_dask = self.verbose
        self.logging_time = self.verbose
        self.mpl_ui = ('mpl' in self.ui or 'plot' in self.ui)
        self.operation = operation
        self.fileset_label = label
        self.global_rotation = global_rotation
        self.global_center = global_center
        self.source_transform_key = 'source_metadata'
        self.reg_transform_key = 'registered'
        self.transition_transform_key = 'transition'
        self.sims = []
        self.sources = []
        self.state = RegState.INIT

        self.input_path = input_path
        if isinstance(input_path, list):
            self.filenames = input_path
            self.input_dir = os.path.dirname(input_path[0])
        else:
            self.filenames = dir_regex(input_path)
            self.input_dir = os.path.dirname(input_path)
        if not self.filenames:
            return False

        if input_labels:
            self.file_labels = input_labels
        else:
            self.file_labels = get_unique_file_labels(self.filenames)

        self.source_metadata = source_metadata
        self.extra_metadata = extra_metadata

        output_path = output_path.format_map(split_numeric_dict(self.filenames[0]))
        self.output = os.path.join(self.input_dir, output_path)    # preserve trailing slash: do not use os.path.normpath()
        output_dir = os.path.dirname(self.output)
        if self.clear:
            shutil.rmtree(output_dir, ignore_errors=True)
        if not os.path.exists(output_dir):
            os.makedirs(output_dir)

        return True

    def run(self):
        with ProgressBar(minimum=60, dt=1) if self.logging_dask else nullcontext():
            return self._run()

    def _run(self):
        filenames = self.filenames
        file_labels = self.file_labels

        output = self.output
        operation = self.operation
        source_metadata = self.source_metadata
        extra_metadata = self.extra_metadata
        if isinstance(extra_metadata, dict):
            z_scale = extra_metadata.get('scale', {}).get('z')
            channels = extra_metadata.get('channels', [])
        else:
            z_scale = None
            channels = []
        normalise_orientation = 'norm' in source_metadata
        output_params = self.output_params
        overlap_threshold = self.register_params.get('overlap_threshold', self.params.get('overlap_threshold', 0.5))
        save_images = self.output_params.get('save_images', self.params.get('save_images', True))

        is_stack = ('stack' in operation)
        is_3d = ('3d' in operation)
        is_simple_stack = is_stack and not is_3d
        is_transition = ('transition' in operation)
        is_channel_overlay = (len(channels) > 1)

        output_format = output_params.get('format', self.params_general.get('format', zarr_extension))
        output_tile_size = output_params.get('tile_size', self.params_general.get('tile_size'))
        output_compression = output_params.get('compression', self.params_general.get('compression'))
        output_pyramid_downsample = output_params.get('pyramid_downsample', self.params_general.get('pyramid_downsample', 2))
        output_npyramid_add = output_params.get('npyramid_add', self.params_general.get('npyramid_add', 0))
        output_ome_version = output_params.get('ome_version', self.params_general.get('ome_version', default_ome_zarr_version))

        mappings_header = ['id','x_pixels', 'y_pixels', 'z_pixels', 'x', 'y', 'z', 'rotation']

        if len(filenames) == 0:
            logging.warning('Skipping (no images)')
            return False

        mappings_filename = output + output_params.get('mappings', default_mappings_name)

        output_filename = operation.split()[0] + 'ed'

        if not self.overwrite and os.path.exists(output_filename):
            logging.warning(f'Skipping existing output {output_filename}')
            return False

        with Timer('init sims', self.logging_time):
            sims = self.init_sims(target_scale=self.preprocess_params.get('scale'))
        self.sims = sims

        if not z_scale:
            z_scale = self.scales[0].get('z', 1)

        with Timer('pre-process', self.logging_time):
            register_sims, register_indices, _ = self.preprocess(sims, **self.preprocess_params)

        data = []
        for label, sim, scale in zip(file_labels, sims, self.scales):
            position, rotation = get_data_mapping(sim, transform_key=self.source_transform_key)
            position_pixels = {dim: position[dim] / float(scale.get(dim, 1)) for dim in position.keys()}
            row = [label] + dict_to_xyz(position_pixels, add_zeros=True) + dict_to_xyz(position, add_zeros=True) + [rotation]
            data.append(row)
        export_csv(output + prereg_mappings_name, data, header=mappings_header)

        if len(filenames) == 1 and save_images and not 'register' in operation and not 'stack' in operation:
            logging.warning('Skipping operation (single image)')
            self.save(output_filename, sims[0], translations0=self.positions,
                      format = output_format,
                      tile_size = output_tile_size,
                      pyramid_downsample = output_pyramid_downsample,
                      npyramid_add = output_npyramid_add,
                      ome_version = output_ome_version)
            return False

        _, has_overlaps = self.validate_overlap(sims, file_labels, is_simple_stack, is_simple_stack or is_channel_overlay)
        overall_overlap = np.mean(has_overlaps)
        if overall_overlap < overlap_threshold:
            raise ValueError(f'Not enough overlap: {overall_overlap * 100:.1f}%')

        if not self.overwrite and os.path.exists(mappings_filename):
            logging.info('Loading registration mappings...')
            # load registration mappings
            mappings = import_json(mappings_filename)
            # copy transforms to sims
            for sim, label in zip(sims, file_labels):
                mapping = param_utils.affine_to_xaffine(np.array(mappings[label]))
                if is_stack:
                    transform = param_utils.identity_transform(ndim=3)
                    transform.loc[{dim: mapping.coords[dim] for dim in mapping.dims}] = mapping
                else:
                    transform = mapping
                si_utils.set_sim_affine(sim, transform, transform_key=self.reg_transform_key)
            if is_stack:
                sims = make_sims_3d(sims, z_scale, self.positions)
        else:
            if 'register' in operation:
                with Timer('register', self.logging_time):
                    results = self.register(sims, register_sims, register_indices, self.register_params)
                reg_result = results['reg_result']
                mappings = results['mappings']
                metrics = results['metrics']

            if is_stack:
                sims = make_sims_3d(sims, z_scale, self.positions)

            if 'register' in operation:
                logging.info(metrics['summary'])
                output_mappings = {file_labels[key]: np.array(mapping.sel(t=0)).tolist() for key, mapping in mappings.items()}
                export_json(mappings_filename, output_mappings)
                output_metrics = {f'{file_labels[keys[0]]}-{file_labels[keys[1]]}': value for keys, value in metrics['pairs'].items()}
                export_json(output + metrics_name, output_metrics)
                data = []
                for label, sim, mapping, scale, position, rotation\
                        in zip(file_labels, sims, mappings.values(), self.scales, self.positions, self.rotations):
                    if not normalise_orientation:
                        # rotation already in msim affine transform
                        rotation = None
                    position, rotation = get_data_mapping(sim, transform_key=self.reg_transform_key,
                                                          transform=mapping,
                                                          translation0=position,
                                                          rotation=rotation)
                    position_pixels = {dim: position[dim] / float(scale.get(dim, 1)) for dim in position.keys()}
                    row = [label] + dict_to_xyz(position_pixels, add_zeros=True) + dict_to_xyz(position, add_zeros=True) + [rotation]
                    data.append(row)
                export_csv(output + metrics_tabular_name, data, header=mappings_header)

                for reg_label, reg_item in reg_result.items():
                    if isinstance(reg_item, dict):
                        summary_plot = reg_item.get('summary_plot')
                        if summary_plot is not None:
                            figure, axes = summary_plot
                            summary_plot_filename = output + f'{reg_label}.pdf'
                            figure.savefig(summary_plot_filename)

        self.sims = sims
        registered_positions_filename = output + registered_positions_name
        if self.reg_transform_key in sims[0].transforms:
            transform_key = self.reg_transform_key
            with Timer('plot positions', self.logging_time):
                vis_utils.plot_positions(sims, transform_key=transform_key,
                                         use_positional_colors=False, view_labels=file_labels, view_labels_size=3,
                                         show_plot=self.mpl_ui, output_filename=registered_positions_filename)
                plt_close()
        else:
            transform_key = self.source_transform_key

        logging.info('Exporting...')

        image_paths = []
        if save_images:
            if self.output_params.get('thumbnail'):
                with Timer('create thumbnail', self.logging_time):
                    self.save_thumbnail('thumb_' + output_filename,
                                        nom_sims=sims,
                                        transform_key=transform_key)

            if 'register' in operation or 'stack' in operation:
                with Timer('fuse image', self.logging_time):
                    if isinstance(self.fusion_params, dict):
                        fusion_method = self.fusion_params.get('method', '')
                        output_spacing = self.fusion_params.get('output_spacing', 'mean')
                    else:
                        fusion_method = self.fusion_params
                        output_spacing = self.params.get('output_spacing', 'mean')
                    fused_image, is_saved = self.fuse(sims, fusion_method=fusion_method, output_spacing=output_spacing,
                                                      transform_key=transform_key, output_filename=output_filename,
                                                      tile_size=output_tile_size, ome_version=output_ome_version)
            else:
                fused_image = sims
                is_saved = False

            if not is_saved or 'tif' in output_format:
                logging.info('Saving fused image...')
                with Timer('save fused image', self.logging_time):
                    self.save(output_filename, fused_image,
                              transform_key=transform_key, translations0=self.positions,
                              format = output_format,
                              tile_size = output_tile_size,
                              compression = output_compression,
                              pyramid_downsample = output_pyramid_downsample,
                              npyramid_add = output_npyramid_add,
                              ome_version = output_ome_version)

            if 'tif' in output_format:
                filename = output_filename + tiff_extension
                image_paths.append(filename)
            if 'zar' in output_format:
                filename = output_filename + zarr_extension
                image_paths.append(filename)
                create_zarr_ro_crate(data, self.output + filename)

        create_ro_crate(fused_image, self.output, image_paths)

        if is_transition:
            self.save_video(output, sims, fused_image)

        return True

    def init_sources(self):
        source_metadata0 = self.source_metadata
        source_metadata = {}
        self.sources = []
        for index, (filename, label) in enumerate(zip(self.filenames, self.file_labels)):
            if isinstance(source_metadata0, dict) and label in source_metadata0:
                source_metadata = source_metadata0[label]
                position, rotation, scale = get_properties_from_transform(param_utils.affine_to_xaffine(np.array(source_metadata)))
                source_metadata = {'position': position, 'rotation': rotation, 'scale': xyz_to_dict([scale, scale])}
            else:
                if 'position' in source_metadata0:
                    translation = source_metadata0['position']
                    if isinstance(translation, list):
                        translation = translation[index]
                    source_metadata['position'] = translation
                if 'scale' in source_metadata0:
                    scale = source_metadata0['scale']
                    if isinstance(scale, list):
                        scale = scale[index]
                    source_metadata['scale'] = scale
                if 'rotation' in source_metadata0:
                    source_metadata['rotation'] = source_metadata0['rotation']
            self.sources.append(create_dask_source(filename, source_metadata))

    def init_sims(self, source_metadata={}, extra_metadata={}, z_scale=None, chunk_size=default_chunk_size,
                  target_scale=None):
        if not source_metadata:
            source_metadata = self.source_metadata
        if not extra_metadata:
            extra_metadata = self.extra_metadata
        source_metadata = import_metadata(source_metadata, input_path=self.input_path)
        extra_metadata = import_metadata(extra_metadata, input_path=self.input_path)
        source_metadata_changed = (source_metadata != self.source_metadata)
        self.source_metadata = source_metadata
        self.extra_metadata = extra_metadata
        if isinstance(source_metadata, dict):
            z_scale = source_metadata.get('scale', {}).get('z')
        if not z_scale and isinstance(extra_metadata, dict):
            z_scale = extra_metadata.get('scale', {}).get('z')

        if len(self.filenames) == 0:
            raise ValueError('No input files')

        logging.info('Initialising sims...')
        if not self.sources or source_metadata_changed:
            self.init_sources()
        sources = self.sources
        source0 = sources[0]
        images = []
        sims = []
        scales = []
        translations = []
        rotations = []

        is_stack = ('stack' in self.operation)
        has_z_size = (source0.get_size().get('z', 0) > 0)

        output_order = 'zyx' if has_z_size else 'yx'
        ndims = len(output_order)
        if source0.get_nchannels() > 1:
            output_order += 'c'

        last_z_position = None
        different_z_positions = False
        delta_zs = []
        for filename, source in zip(self.filenames, sources):
            scale = source.get_pixel_size()
            translation = source.get_position()
            rotation = source.get_rotation()

            if 'is_center' in source_metadata:
                translation = {dim: translation[dim] - source.get_physical_size().get(dim, 0) / 2 for dim in translation}

            if 'sbem' in source_metadata:
                source_version = source.metadata.get('Creator', source.metadata.get('creator', ''))
                if '2025' in source_version:
                    path = os.path.dirname(self.filenames[0])
                    metapath = None
                    attempts = 0
                    while attempts < 3:
                        metapath = os.path.join(path, 'meta')
                        if os.path.exists(metapath):
                            break
                        path = os.path.join(path, '..')
                        attempts += 1
                    if metapath:
                        sbemimage_config = load_sbemimage_best_config(metapath, filename)
                        if sbemimage_config:
                            size = source.get_size()
                            translation, scale0 = adjust_sbemimage_properties(translation, scale, size,
                                                                              filename, sbemimage_config)
                            if scale0:
                                scale = scale0
                            elif scale.get('x') != scale.get('y'):
                                logging.warning('SBEMimage pixel size requires correction, please provide in source metadata.')
                            logging.debug(f'Adjusted SBEMimage properties for {filename}')
                        else:
                            logging.warning(f'Could not find SBEMimage config for {filename}.')

            level = 0
            rescale = 1
            if target_scale:
                # Only downscaling
                level, rescale, scale = get_level_from_scale(source.scales, scale, target_scale)
            if 'invert' in source_metadata:
                translation['x'] = -translation['x']
                translation['y'] = -translation['y']
            if 'z' in translation:
                z_position = translation['z']
            else:
                z_position = 0
            if last_z_position is not None and z_position != last_z_position:
                different_z_positions = True
                delta_zs.append(z_position - last_z_position)
            if 'rotation' in source_metadata:
                rotation = source_metadata['rotation']
            if self.global_rotation is not None:
                rotation = self.global_rotation

            dask_data = source.get_data(level=level)
            if rescale != 1:
                new_shape = [int(size / rescale) if dim in 'xy' else 1
                             for dim, size in zip(source.dimension_order, dask_data.shape)]
                dask_data = resize(dask_data, new_shape, preserve_range=True).astype(dask_data.dtype)
            image = redimension_data(dask_data, source.dimension_order, output_order)

            scales.append(scale)
            translations.append(translation)
            rotations.append(rotation)
            images.append(image)
            last_z_position = z_position

        if 'z' in output_order and z_scale is None:
            if len(delta_zs) > 0:
                z_scale = np.min(delta_zs)
            else:
                z_scale = 1

        if 'norm' in source_metadata:
            sizes = [source.get_physical_size() for source in sources]
            center = {dim: 0 for dim in output_order}
            if 'center' in source_metadata:
                if 'global' in source_metadata:
                    center = self.global_center
                else:
                    center = {dim: float(np.mean([translation[dim] for translation in translations])) for dim in translations[0]}
            translations, rotations = normalise_rotated_positions(translations, rotations, sizes, center, len(output_order))

        #translations = [np.array(translation) * 1.25 for translation in translations]

        increase_z_positions = is_stack and not different_z_positions

        z_position = 0
        final_scales = []
        final_translations = []
        for source, image, scale, translation, rotation, file_label in zip(sources, images, scales, translations, rotations, self.file_labels):
            # transform #dimensions need to match
            if 'z' in output_order:
                if len(scale) > 0 and 'z' not in scale:
                    scale['z'] = abs(z_scale)
                if (len(translation) > 0 and 'z' not in translation) or increase_z_positions:
                    translation['z'] = z_position
                if increase_z_positions:
                    z_position += z_scale
            channel_labels = [channel.get('label', '') for channel in source.get_channels()]
            if rotation is None or 'norm' in source_metadata:
                # if positions are normalised, don't use rotation
                transform = None
            else:
                transform = param_utils.invert_coordinate_order(
                    create_transform(translation, rotation, matrix_size=ndims + 1)
                )
            if file_label in extra_metadata:
                transform2 = extra_metadata[file_label]
                if transform is None:
                    transform = np.array(transform2)
                else:
                    transform = np.array(combine_transforms([transform, transform2]))

            # fix empty dictionaries args
            scale_arg = scale if scale else None
            translation_arg = translation if translation else None

            sim = si_utils.get_sim_from_array(
                image,
                dims=list(output_order),
                scale=scale_arg,
                translation=translation_arg,
                affine=transform,
                transform_key=self.source_transform_key,
                c_coords=channel_labels
            )
            if len(sim.chunksizes.get('x')) == 1 and len(sim.chunksizes.get('y')) == 1:
                if isinstance(chunk_size, int):
                    chunk_size = [chunk_size] * 2
                sim = sim.chunk(xyz_to_dict(chunk_size))
            sims.append(sim)
            final_scales.append(scale)
            final_translations.append(translation)

        self.scales = final_scales
        self.positions = final_translations
        self.rotations = rotations
        self.sims = sims
        self.state = RegState.SIMS_INIT
        return sims

    def validate_overlap(self, sims, labels, is_stack=False, expect_large_overlap=False):
        min_dists = []
        has_overlaps = []
        n = len(sims)
        positions = [get_sim_position_final(sim, get_center=True) for sim in sims]
        sizes = [float(np.linalg.norm(list(get_sim_physical_size(sim).values()))) for sim in sims]
        for i in range(n):
            norm_dists = []
            # check if only single z slices
            if is_stack:
                if i + 1 < n:
                    compare_indices = [i + 1]
                else:
                    compare_indices = []
            else:
                compare_indices = range(n)
            for j in compare_indices:
                if not j == i:
                    distance = math.dist(positions[i].values(), positions[j].values())
                    norm_dist = distance / np.mean([sizes[i], sizes[j]])
                    norm_dists.append(norm_dist)
            if len(norm_dists) > 0:
                norm_dist = min(norm_dists)
                min_dists.append(float(norm_dist))
                if norm_dist >= 1:
                    logging.warning(f'{labels[i]} has no overlap')
                    has_overlaps.append(False)
                elif expect_large_overlap and norm_dist > 0.5:
                    logging.warning(f'{labels[i]} has small overlap')
                    has_overlaps.append(False)
                else:
                    has_overlaps.append(True)
        return min_dists, has_overlaps

    def preprocess(self, sims,
                   flatfield_quantiles=None, normalisation=None, gaussian_sigma=None, filter_foreground=False):
        modified = False
        # normalise pixel size: take max pixel size
        max_scale = {dim: max(scale.get(dim, 1) for scale in self.scales) for dim in 'xy'}
        scales0 = self.scales

        if filter_foreground:
            foreground_map = calc_foreground_map(sims)
            modified = True
        else:
            foreground_map = None
        if flatfield_quantiles is not None:
            logging.info('Flat-field correction...')
            if isinstance(flatfield_quantiles, str):
                flatfield_quantiles = [float(quantile.strip()) for quantile in flatfield_quantiles.split(',')]
            new_sims = [None] * len(sims)
            for sim_indices in group_sims_by_z(sims, self.positions):
                sims_z_set = [sims[i] for i in sim_indices]
                foreground_map_z_set = [foreground_map[i] for i in sim_indices] if foreground_map is not None else None
                new_sims_z_set = flatfield_correction(sims_z_set, self.source_transform_key, flatfield_quantiles,
                                                      foreground_map=foreground_map_z_set)
                for sim_index, sim in zip(sim_indices, new_sims_z_set):
                    new_sims[sim_index] = sim
            sims = new_sims
            modified = True

        if gaussian_sigma:
            logging.info('Applying Gaussian filtering...')
            new_sims = []
            for sim, scale0 in zip(sims, scales0):
                # factor in original pixel size for gaussian sigma value
                scale = np.mean(list(scale0.values())) / np.mean(list(max_scale.values()))
                sigma = gaussian_sigma * (scale ** (1 / 3))
                new_sims.append(gaussian_filter_sim(sim, self.source_transform_key, sigma))
            sims = new_sims
            modified = True

        if normalisation is not None:
            if isinstance(normalisation, str) and normalisation.lower() in ['false', 'no', 'none', '']:
                normalisation = None
            elif isinstance(normalisation, bool) and normalisation == False:
                normalisation = None
        if normalisation:
            use_global = ('global' in str(normalisation).lower())
            if use_global:
                logging.info('Normalising (global)...')
            else:
                logging.info('Normalising (individual)...')
            sims = normalise_sims(sims, self.source_transform_key, use_global=use_global)
            modified = True

        if filter_foreground:
            logging.info('Filtering foreground images...')
            #tile_vars = np.array([np.asarray(np.std(sim)).item() for sim in sims])
            #threshold1 = np.mean(tile_vars)
            #threshold2 = np.median(tile_vars)
            #threshold3, _ = cv.threshold(np.array(tile_vars).astype(np.uint16), 0, 1, cv.THRESH_OTSU)
            #threshold = min(threshold1, threshold2, threshold3)
            #foregrounds = (tile_vars >= threshold)
            new_sims = [sim for sim, is_foreground in zip(sims, foreground_map) if is_foreground]
            logging.info(f'Foreground images: {len(new_sims)} / {len(sims)}')
            indices = np.where(foreground_map)[0]
            sims = new_sims
            modified = True
        else:
            indices = range(len(sims))
        self.register_sims = sims
        self.register_indices = indices
        return sims, indices, modified

    def create_registration_method(self, sim0, params={}, method=''):
        registration_method = None
        pairwise_reg_func_kwargs = None

        if 'registration' in params:
            params = params['registration']
        if not method:
            method = params.get('method',
                                params.get('name', ''))
        method = method.lower()

        if 'cpd' in method:
            from src.muvis_align.registration_methods.RegistrationMethodCPD import RegistrationMethodCPD
            registration_method = RegistrationMethodCPD(sim0, params, self.debug)
            pairwise_reg_func = registration_method.registration
        elif 'feature' in method or 'orb' in method or 'sift' in method:
            if 'cv' in method:
                from src.muvis_align.registration_methods.RegistrationMethodCvFeatures import RegistrationMethodCvFeatures
                registration_method = RegistrationMethodCvFeatures(sim0, params, self.debug)
            else:
                from src.muvis_align.registration_methods.RegistrationMethodSkFeatures import RegistrationMethodSkFeatures
                registration_method = RegistrationMethodSkFeatures(sim0, params, self.debug)
            pairwise_reg_func = registration_method.registration
        elif 'ant' in method:
            pairwise_reg_func = registration.registration_ANTsPy
            # args for ANTsPy registration: used internally by ANYsPy algorithm
            pairwise_reg_func_kwargs = {
                'transform_types': ['Rigid'],
                "aff_random_sampling_rate": 0.5,
                "aff_iterations": (2000, 2000, 1000, 1000),
                "aff_smoothing_sigmas": (4, 2, 1, 0),
                "aff_shrink_factors": (16, 8, 2, 1),
            }
        else:
            pairwise_reg_func = registration.phase_correlation_registration

        self.registration_method = registration_method

        return method, pairwise_reg_func, pairwise_reg_func_kwargs

    def create_fusion_method(self, fusion_method, sim0):
        if fusion_method is None:
            fusion_method = ''
        if 'compos' in fusion_method:
            fusion_method = None
            fuse_func = None
        elif 'exclus' in fusion_method:
            from src.muvis_align.fusion_methods.FusionMethodExclusive import FusionMethodExclusive
            fusion_method = FusionMethodExclusive(sim0, self.debug)
            fuse_func = fusion_method.fusion
        elif 'add' in fusion_method:
            from src.muvis_align.fusion_methods.FusionMethodAdditive import FusionMethodAdditive
            fusion_method = FusionMethodAdditive(sim0, self.debug)
            fuse_func = fusion_method.fusion
        else:
            fuse_func = fusion.simple_average_fusion

        self.fusion_method = fusion_method

        return fuse_func

    def register(self, sims, register_sims=None, register_indices=None, params=None):
        self.register_pairs(sims, register_sims=register_sims, params=params)
        results = self.register_global(sims, self.msims, register_indices=register_indices, params=params)
        return results

    def register_pairs(self, sims, register_sims=None, params=None):
        operation = self.operation
        pairing = params.get('pairing',
                             params.get('registration', {}).get('pairing', '')).lower()
        n_parallel_pairwise_regs = params.get('n_parallel_pairwise_regs',
                                              params.get('registration', {}).get('n_parallel_pairwise_regs'))
        if n_parallel_pairwise_regs is not None and n_parallel_pairwise_regs == '0':
            n_parallel_pairwise_regs = None

        is_stack = ('stack' in operation)
        is_3d = ('3d' in operation)

        reg_channel = params.get('channel', 0)
        if isinstance(reg_channel, int):
            reg_channel_index = reg_channel
            reg_channel = None
        else:
            reg_channel_index = None

        if register_sims is None:
            register_sims = sims
        if is_stack and not is_3d:
            # register in 2d; pairwise consecutive views
            register_sims = [si_utils.max_project_sim(sim, dim='z') if 'z' in sim.dims else sim
                             for sim in register_sims]
            pairs = [(index, index + 1) for index in range(len(register_sims) - 1)]
        elif 'ortho' in pairing or 'overla' in pairing:
            origins = np.array([get_sim_position_final(sim, position, get_center=True)
                                for sim, position in zip(sims, self.positions)])
            sizes = [get_sim_physical_size(sim) for sim in sims]
            pairs, _ = get_pairs(origins, sizes, pairing)
            logging.info(f'#pairs: {len(pairs)}')
            #for pair in pairs:
            #    print(f'{self.file_labels[pair[0]]} - {self.file_labels[pair[1]]}')
        else:
            pairs = None

        reg_method, pairwise_reg_func, pairwise_reg_func_kwargs = self.create_registration_method(register_sims[0],
                                                                                                  params=params)
        logging.info(f'Registration method: {reg_method}')
        logging.info('Registering...')
        register_msims = [msi_utils.get_msim_from_sim(sim) for sim in register_sims]

        overlap_tolerance = 0

        # ******* start MVS registration functions

        if "c" in msi_utils.get_dims(register_msims[0]):
            if reg_channel is None:
                if reg_channel_index is None:
                    for msim in register_msims:
                        if "c" in msi_utils.get_dims(msim):
                            raise (
                                Exception("Please choose a registration channel.")
                            )
                else:
                    reg_channel = sims[0].coords["c"][reg_channel_index]

            msims_reg = [
                msi_utils.multiscale_sel_coords(msim, {"c": reg_channel})
                if "c" in msi_utils.get_dims(msim)
                else msim
                for imsim, msim in enumerate(register_msims)
            ]
        else:
            msims_reg = register_msims

        try:
            with dask.config.set(scheduler='threads'):
                g_reg = mv_graph.build_view_adjacency_graph_from_msims(
                    msims_reg,
                    transform_key=self.source_transform_key,
                    pairs=pairs,
                    overlap_tolerance=overlap_tolerance,
                )

                g_reg_computed = compute_pairwise_registrations(
                    msims_reg,
                    g_reg,
                    transform_key=self.source_transform_key,
                    overlap_tolerance=overlap_tolerance,
                    pairwise_reg_func=pairwise_reg_func,
                    pairwise_reg_func_kwargs=pairwise_reg_func_kwargs,
                    n_parallel_pairwise_regs=n_parallel_pairwise_regs,
                )

                # ******* end MVS registration functions

                # reg_result = registration.register(
                #     register_msims,
                #     reg_channel=reg_channel,
                #     reg_channel_index=reg_channel_index,
                #     transform_key=self.source_transform_key,
                #     new_transform_key=self.reg_transform_key,
                #
                #     pairs=pairs,
                #     pre_registration_pruning_method=None,
                #
                #     pairwise_reg_func=pairwise_reg_func,
                #     pairwise_reg_func_kwargs=pairwise_reg_func_kwargs,
                #
                #     groupwise_resolution_method=groupwise_resolution_method,
                #     groupwise_resolution_kwargs=groupwise_resolution_kwargs,
                #
                #     post_registration_do_quality_filter=(post_registration_quality_threshold is not None),
                #     post_registration_quality_threshold=post_registration_quality_threshold,
                #
                #     n_parallel_pairwise_regs=n_parallel_pairwise_regs,
                #
                #     plot_summary=self.mpl_ui,
                #     return_dict=return_dict,
                # )

        except NotEnoughOverlapError:
            g_reg_computed = g_reg

        metrics = calc_pair_metrics(msims_reg, g_reg_computed, params.get('metrics', []), self.source_transform_key,
                                    reg_channel=reg_channel_index, n_parallel_pairs=n_parallel_pairwise_regs)

        self.pairs_graph = g_reg_computed
        self.msims = msims_reg
        self.pairs = pairs
        self.metrics = metrics
        self.state = RegState.PAIRS_REG
        return {
            'pairs_graph': self.pairs_graph,
            'msims': msims_reg,
            'pairs': pairs,
            'metrics': metrics
        }

    def register_global(self, sims, msims, register_indices=None, params=None,
                        pairs_graph=None):
        if pairs_graph is not None:
            g_reg_computed = pairs_graph
        else:
            g_reg_computed = self.pairs_graph

        sim0 = sims[0]
        ndims = si_utils.get_ndim_from_sim(sim0)

        groupwise_resolution_method = params.get('groupwise_resolution_method',
                                                 params.get('registration', {}).get('groupwise_resolution_method', 'global_optimization'))
        groupwise_resolution_kwargs = {}
        if groupwise_resolution_method == 'global_optimization':
           groupwise_resolution_kwargs['transform'] = params.get('transform_type',
                                                                 params.get('registration', {}).get('transform_type'))
           # transform_type options include 'translation', 'rigid', 'affine', 'similarity'

        post_registration_quality_threshold = params.get('post_registration_quality_threshold',
                                                         params.get('registration', {}).get('post_registration_quality_threshold'))
        post_registration_do_quality_filter = (post_registration_quality_threshold is not None)

        n_parallel_pairwise_regs = params.get('n_parallel_pairwise_regs',
                                              params.get('registration', {}).get('n_parallel_pairwise_regs'))
        if n_parallel_pairwise_regs is not None and n_parallel_pairwise_regs == '0':
            n_parallel_pairwise_regs = None

        plot_summary = self.mpl_ui

        # ******* start MVS registration functions

        if post_registration_do_quality_filter:
            # filter edges by quality
            g_reg_computed = mv_graph.filter_edges(
                g_reg_computed,
                threshold=post_registration_quality_threshold,
                weight_key="quality",
            )

        with dask.config.set(scheduler='threads'):
            transforms_dict, groupwise_resolution_info_dict = groupwise_resolution(
                g_reg_computed,
                method=groupwise_resolution_method,
                **groupwise_resolution_kwargs,
            )

        transforms = [
            transforms_dict[iview] for iview in sorted(g_reg_computed.nodes())
        ]

        for imsim, msim in enumerate(msims):
            msi_utils.set_affine_transform(
                msim,
                transforms[imsim],
                transform_key=self.reg_transform_key,
                base_transform_key=self.source_transform_key,
            )

        if plot_summary:
            plot_info = _plot_registration_summaries(
                msims,
                self.source_transform_key,
                self.reg_transform_key,
                g_reg_computed,
                groupwise_resolution_info_dict,
                show_plot=plot_summary,
            )
        else:
            plot_info = {}

        reg_result = {
            "params": transforms,
            "pairwise_registration": {
                "graph": g_reg_computed,
                "metrics": {
                    "qualities": nx.get_edge_attributes(
                        g_reg_computed, "quality"
                    )
                },
                "summary_plot": None if plot_summary is False
                else (
                    plot_info['fig_pair_reg'],
                    plot_info['ax_pair_reg']
                )
            },
            "groupwise_resolution": {
                "metrics": groupwise_resolution_info_dict,
                "summary_plot": None if plot_summary is False
                else (
                    plot_info['fig_group_res'],
                    plot_info['ax_group_res']
                )
            },
        }

        # ******* end MVS registration functions

        if register_indices is None:
            register_indices = range(len(msims))

        # copy transforms from register sims to unmodified sims
        for reg_msim, index in zip(msims, register_indices):
            si_utils.set_sim_affine(
                sims[index],
                msi_utils.get_transform_from_msim(reg_msim, transform_key=self.reg_transform_key),
                transform_key=self.reg_transform_key)

        # set missing transforms
        for sim in sims:
            if self.reg_transform_key not in si_utils.get_tranform_keys_from_sim(sim):
                si_utils.set_sim_affine(
                    sim,
                    param_utils.identity_transform(ndim=ndims, t_coords=[0]),
                    transform_key=self.reg_transform_key)

        mappings = reg_result['params']
        # re-index from subset of sims
        residual_error_dict = reg_result.get('groupwise_resolution', {}).get('metrics', {}).get('residuals', {})
        residual_error_dict = {(register_indices[key[0]], register_indices[key[1]]): value.item()
                               for key, value in residual_error_dict.items()}
        registration_qualities_dict = reg_result.get('pairwise_registration', {}).get('metrics', {}).get('qualities', {})
        registration_qualities_dict = {(register_indices[key[0]], register_indices[key[1]]): value
                                       for key, value in registration_qualities_dict.items()}

        # re-index from subset of sims
        mappings_dict = {index: mapping for index, mapping in zip(register_indices, mappings)}

        reg_channel = params.get('channel', 0)
        metrics = calc_global_metrics(msims, self.source_transform_key, self.reg_transform_key,
                                      params.get('metrics', []), reg_channel=reg_channel, reg_results=reg_result,
                                      n_parallel_pairs=n_parallel_pairwise_regs)

        self.metrics = metrics
        self.state = RegState.GLOBAL_REG
        return {'reg_result': reg_result,
                'mappings': mappings_dict,
                'residual_errors': residual_error_dict,
                'registration_qualities': registration_qualities_dict,
                'metrics': metrics}

    def fuse(self, sims, fusion_method=None, output_spacing='mean', transform_key=None,
             output_filename=None, tile_size=None, ome_version=default_ome_zarr_version):
        if output_filename is not None:
            output_filename = self.output + output_filename
        sim0 = sims[0]
        if transform_key is None:
            transform_key = self.reg_transform_key
        if isinstance(self.extra_metadata, dict):
            z_scale = self.extra_metadata.get('scale', {}).get('z')
            channels = self.extra_metadata.get('channels', [])
        else:
            z_scale = None
            channels = []
        is_channel_overlay = (len(channels) > 1)

        if z_scale is None and self.scales is not None:
            z_scale0 = np.mean([scale.get('z', 0) for scale in self.scales])
            if z_scale0 > 0:
                z_scale = z_scale0
        if z_scale is None:
            if 'z' in sim0.dims:
                diffs = np.diff(sorted(set([si_utils.get_origin_from_sim(sim).get('z', 0) for sim in sims])))
                if len(diffs) > 0:
                    z_scale = min(diffs)

        output_stack_properties = calc_output_properties(sims, transform_key,
                                                         output_spacing=output_spacing, z_scale=z_scale)

        if self.verbose:
            logging.info(f'Output stack: {numpy_to_native(output_stack_properties)}')
        data_size = np.prod(list(output_stack_properties['shape'].values())) * sim0.dtype.itemsize
        logging.info(f'Fusing {print_hbytes(data_size)}')

        saving_zarr = False
        if is_channel_overlay:
            # convert to multichannel images
            channel_sims = [fusion.fuse(
                [sim],
                transform_key=transform_key,
                output_stack_properties=output_stack_properties
            ) for sim in sims]
            channel_sims = [sim.assign_coords({'c': [channels[simi]['label']]}) for simi, sim in enumerate(channel_sims)]
            fused_image = xr.combine_nested([sim.rename() for sim in channel_sims], concat_dim='c', combine_attrs='override')
        else:
            fuse_func = self.create_fusion_method(fusion_method, sim0)
            if fuse_func:
                saving_zarr = output_filename is not None
                output_chunksize = None
                if saving_zarr and not output_filename.lower().endswith('.zarr'):
                    output_filename += '.ome.zarr'
                    zarr_options = {'ome_zarr': saving_zarr, 'ngff_version': ome_version}
                    if tile_size is not None:
                        if not isinstance(tile_size, (list, tuple)):
                            tile_size = [tile_size] * 2
                        output_chunksize = xyz_to_dict(tile_size)
                        if 'z' in output_stack_properties['shape'] and 'z' not in output_chunksize:
                            output_chunksize['z'] = 1
                else:
                    zarr_options = None
                with dask.config.set(scheduler='threads'):
                    fused_image = fusion.fuse(
                        sims,
                        fusion_func=fuse_func,
                        transform_key=transform_key,
                        output_stack_properties=output_stack_properties,
                        output_zarr_url=output_filename,
                        zarr_options=zarr_options,
                        output_chunksize=output_chunksize
                    )
                if saving_zarr:
                    open(output_filename.rstrip('.zarr').rstrip('.ome'), 'w')
            else:
                fused_image = sims
        return fused_image, saving_zarr

    def save_thumbnail(self, output_filename, nom_sims=None, transform_key=None):
        output_params = self.params_general['output']
        thumbnail_scale = output_params.get('thumbnail_scale', 16)
        is_stack = ('stack' in self.operation)
        if isinstance(self.extra_metadata, dict):
            z_scale = self.extra_metadata.get('scale', {}).get('z')
        else:
            z_scale = None

        sims = self.init_sims(target_scale=thumbnail_scale)
        if is_stack:
            sims = make_sims_3d(sims, z_scale, self.positions)

        if nom_sims is not None:
            if sims[0].sizes['x'] >= nom_sims[0].sizes['x']:
                logging.warning('Unable to generate scaled down thumbnail due to lack of source pyramid sizes')
                return

            if transform_key is not None and transform_key != self.source_transform_key:
                for nom_sim, sim in zip(nom_sims, sims):
                    si_utils.set_sim_affine(sim,
                                            si_utils.get_affine_from_sim(nom_sim, transform_key=transform_key),
                                            transform_key=transform_key)
        fused_image, is_saved = self.fuse(sims, transform_key=transform_key, output_spacing='max',
                                          output_filename=output_filename)
        if not is_saved or 'tif' in output_params.get('thumbnail'):
            self.save(output_filename, fused_image.squeeze(), transform_key=transform_key,
                      format=output_params.get('thumbnail'), ome_version=output_params.get('ome_version'))
        self.state = RegState.FUSED

    def save(self, output_filename, data, format=zarr_extension, transform_key=None, translations0=None,
             tile_size=None, compression=None, pyramid_downsample=2, npyramid_add=0, ome_version=default_ome_zarr_version):
        if output_filename is not None:
            output_filename = self.output + output_filename
        if isinstance(self.extra_metadata, dict):
            channels = self.extra_metadata.get('channels', [])
        else:
            channels = []
        save_image(output_filename, data, format,
                   transform_key=transform_key, channels=channels, translations0=translations0,
                   tile_size=tile_size, compression=compression,
                   pyramid_downsample=pyramid_downsample, npyramid_add=npyramid_add,
                   ome_version=ome_version,
                   verbose=self.verbose)

    def save_video(self, output, sims, fused_image):
        logging.info('Creating transition video...')
        pixel_size = [si_utils.get_spacing_from_sim(sims[0]).get(dim, 1) for dim in 'xy']
        params = self.params
        nframes = params.get('frames', 1)
        spacing = params.get('spacing', [1.1, 1])
        scale = params.get('scale', 1)
        transition_filename = output + 'transition'
        video = Video(transition_filename + '.mp4', fps=params.get('fps', 1))
        positions0 = np.array([si_utils.get_origin_from_sim(sim, asarray=True) for sim in sims])
        center = np.mean(positions0, 0)
        window = get_image_window(fused_image)

        max_size = None
        acum = 0
        for framei in range(nframes):
            c = (1 - np.cos(framei / (nframes - 1) * 2 * math.pi)) / 2
            acum += c / (nframes / 2)
            spacing1 = spacing[0] + (spacing[1] - spacing[0]) * acum
            for sim, position0 in zip(sims, positions0):
                transform = param_utils.identity_transform(ndim=2, t_coords=[0])
                transform[0][:2, 2] += (position0 - center) * spacing1
                si_utils.set_sim_affine(sim, transform, transform_key=self.transition_transform_key)
            frame = fusion.fuse(sims, transform_key=self.transition_transform_key).squeeze()
            frame = float2int_image(normalise_values(frame, window[0], window[1]))
            frame = cv.resize(np.asarray(frame), None, fx=scale, fy=scale)
            if max_size is None:
                max_size = frame.shape[1], frame.shape[0]
                video.size = max_size
            frame = image_reshape(frame, max_size)
            save_tiff(transition_filename + f'{framei:04d}.tiff', frame, None, pixel_size)
            video.write(frame)

        video.close()

    def get_metrics(self, metric=None, pair=None):
        if pair is not None:
            if isinstance(pair, np.ndarray):
                pair = pair.tolist()
                pair = tuple(pair)
            metrics = self.metrics.get(pair, {})
        else:
            metrics = self.metrics
        if metric is not None:
            return metrics.get(metric, 0)
        else:
            return metrics

metrics = {} instance-attribute

sims = [] instance-attribute

sources = [] instance-attribute

state = RegState.UNINIT instance-attribute

__init__(operation='register', label='', input_path=None, output_path=None, source_metadata={}, extra_metadata={}, global_rotation=None, global_center=None, overwrite=True, clear=False, ui='', verbose=False, debug=False)

Source code in src\muvis_align\MVSRegistration.py

def __init__(self, operation='register', label='', input_path=None, output_path=None,
             source_metadata={}, extra_metadata={},
             global_rotation=None, global_center=None,
             overwrite=True, clear=False, ui='', verbose=False, debug=False):
    self.state = RegState.UNINIT
    self.sims = []
    self.sources = []
    self.metrics = {}
    if input_path is not None:
        self.init(operation=operation, label=label, input_path=input_path, output_path=output_path,
                  source_metadata=source_metadata, extra_metadata=extra_metadata,
                  global_rotation=global_rotation, global_center=global_center,
                  overwrite=overwrite, clear=clear, ui=ui, verbose=verbose, debug=debug)

create_fusion_method(fusion_method, sim0)

Source code in src\muvis_align\MVSRegistration.py

def create_fusion_method(self, fusion_method, sim0):
    if fusion_method is None:
        fusion_method = ''
    if 'compos' in fusion_method:
        fusion_method = None
        fuse_func = None
    elif 'exclus' in fusion_method:
        from src.muvis_align.fusion_methods.FusionMethodExclusive import FusionMethodExclusive
        fusion_method = FusionMethodExclusive(sim0, self.debug)
        fuse_func = fusion_method.fusion
    elif 'add' in fusion_method:
        from src.muvis_align.fusion_methods.FusionMethodAdditive import FusionMethodAdditive
        fusion_method = FusionMethodAdditive(sim0, self.debug)
        fuse_func = fusion_method.fusion
    else:
        fuse_func = fusion.simple_average_fusion

    self.fusion_method = fusion_method

    return fuse_func

create_registration_method(sim0, params={}, method='')

Source code in src\muvis_align\MVSRegistration.py

def create_registration_method(self, sim0, params={}, method=''):
    registration_method = None
    pairwise_reg_func_kwargs = None

    if 'registration' in params:
        params = params['registration']
    if not method:
        method = params.get('method',
                            params.get('name', ''))
    method = method.lower()

    if 'cpd' in method:
        from src.muvis_align.registration_methods.RegistrationMethodCPD import RegistrationMethodCPD
        registration_method = RegistrationMethodCPD(sim0, params, self.debug)
        pairwise_reg_func = registration_method.registration
    elif 'feature' in method or 'orb' in method or 'sift' in method:
        if 'cv' in method:
            from src.muvis_align.registration_methods.RegistrationMethodCvFeatures import RegistrationMethodCvFeatures
            registration_method = RegistrationMethodCvFeatures(sim0, params, self.debug)
        else:
            from src.muvis_align.registration_methods.RegistrationMethodSkFeatures import RegistrationMethodSkFeatures
            registration_method = RegistrationMethodSkFeatures(sim0, params, self.debug)
        pairwise_reg_func = registration_method.registration
    elif 'ant' in method:
        pairwise_reg_func = registration.registration_ANTsPy
        # args for ANTsPy registration: used internally by ANYsPy algorithm
        pairwise_reg_func_kwargs = {
            'transform_types': ['Rigid'],
            "aff_random_sampling_rate": 0.5,
            "aff_iterations": (2000, 2000, 1000, 1000),
            "aff_smoothing_sigmas": (4, 2, 1, 0),
            "aff_shrink_factors": (16, 8, 2, 1),
        }
    else:
        pairwise_reg_func = registration.phase_correlation_registration

    self.registration_method = registration_method

    return method, pairwise_reg_func, pairwise_reg_func_kwargs

fuse(sims, fusion_method=None, output_spacing='mean', transform_key=None, output_filename=None, tile_size=None, ome_version=default_ome_zarr_version)

Source code in src\muvis_align\MVSRegistration.py

def fuse(self, sims, fusion_method=None, output_spacing='mean', transform_key=None,
         output_filename=None, tile_size=None, ome_version=default_ome_zarr_version):
    if output_filename is not None:
        output_filename = self.output + output_filename
    sim0 = sims[0]
    if transform_key is None:
        transform_key = self.reg_transform_key
    if isinstance(self.extra_metadata, dict):
        z_scale = self.extra_metadata.get('scale', {}).get('z')
        channels = self.extra_metadata.get('channels', [])
    else:
        z_scale = None
        channels = []
    is_channel_overlay = (len(channels) > 1)

    if z_scale is None and self.scales is not None:
        z_scale0 = np.mean([scale.get('z', 0) for scale in self.scales])
        if z_scale0 > 0:
            z_scale = z_scale0
    if z_scale is None:
        if 'z' in sim0.dims:
            diffs = np.diff(sorted(set([si_utils.get_origin_from_sim(sim).get('z', 0) for sim in sims])))
            if len(diffs) > 0:
                z_scale = min(diffs)

    output_stack_properties = calc_output_properties(sims, transform_key,
                                                     output_spacing=output_spacing, z_scale=z_scale)

    if self.verbose:
        logging.info(f'Output stack: {numpy_to_native(output_stack_properties)}')
    data_size = np.prod(list(output_stack_properties['shape'].values())) * sim0.dtype.itemsize
    logging.info(f'Fusing {print_hbytes(data_size)}')

    saving_zarr = False
    if is_channel_overlay:
        # convert to multichannel images
        channel_sims = [fusion.fuse(
            [sim],
            transform_key=transform_key,
            output_stack_properties=output_stack_properties
        ) for sim in sims]
        channel_sims = [sim.assign_coords({'c': [channels[simi]['label']]}) for simi, sim in enumerate(channel_sims)]
        fused_image = xr.combine_nested([sim.rename() for sim in channel_sims], concat_dim='c', combine_attrs='override')
    else:
        fuse_func = self.create_fusion_method(fusion_method, sim0)
        if fuse_func:
            saving_zarr = output_filename is not None
            output_chunksize = None
            if saving_zarr and not output_filename.lower().endswith('.zarr'):
                output_filename += '.ome.zarr'
                zarr_options = {'ome_zarr': saving_zarr, 'ngff_version': ome_version}
                if tile_size is not None:
                    if not isinstance(tile_size, (list, tuple)):
                        tile_size = [tile_size] * 2
                    output_chunksize = xyz_to_dict(tile_size)
                    if 'z' in output_stack_properties['shape'] and 'z' not in output_chunksize:
                        output_chunksize['z'] = 1
            else:
                zarr_options = None
            with dask.config.set(scheduler='threads'):
                fused_image = fusion.fuse(
                    sims,
                    fusion_func=fuse_func,
                    transform_key=transform_key,
                    output_stack_properties=output_stack_properties,
                    output_zarr_url=output_filename,
                    zarr_options=zarr_options,
                    output_chunksize=output_chunksize
                )
            if saving_zarr:
                open(output_filename.rstrip('.zarr').rstrip('.ome'), 'w')
        else:
            fused_image = sims
    return fused_image, saving_zarr

get_metrics(metric=None, pair=None)

Source code in src\muvis_align\MVSRegistration.py

def get_metrics(self, metric=None, pair=None):
    if pair is not None:
        if isinstance(pair, np.ndarray):
            pair = pair.tolist()
            pair = tuple(pair)
        metrics = self.metrics.get(pair, {})
    else:
        metrics = self.metrics
    if metric is not None:
        return metrics.get(metric, 0)
    else:
        return metrics

init(operation='', label='', input_path=None, input_labels=None, output_path=None, source_metadata={}, extra_metadata={}, global_rotation=None, global_center=None, overwrite=True, clear=False, ui='', verbose=False, debug=False)

Source code in src\muvis_align\MVSRegistration.py

def init(self, operation='', label='', input_path=None, input_labels=None, output_path=None,
         source_metadata={}, extra_metadata={}, global_rotation=None, global_center=None,
         overwrite=True, clear=False, ui='', verbose=False, debug=False):
    self.overwrite = overwrite
    self.clear = clear
    self.ui = ui
    self.verbose = verbose
    self.debug = debug
    self.logging_dask = self.verbose
    self.logging_time = self.verbose
    self.mpl_ui = ('mpl' in self.ui or 'plot' in self.ui)
    self.operation = operation
    self.fileset_label = label
    self.global_rotation = global_rotation
    self.global_center = global_center
    self.source_transform_key = 'source_metadata'
    self.reg_transform_key = 'registered'
    self.transition_transform_key = 'transition'
    self.sims = []
    self.sources = []
    self.state = RegState.INIT

    self.input_path = input_path
    if isinstance(input_path, list):
        self.filenames = input_path
        self.input_dir = os.path.dirname(input_path[0])
    else:
        self.filenames = dir_regex(input_path)
        self.input_dir = os.path.dirname(input_path)
    if not self.filenames:
        return False

    if input_labels:
        self.file_labels = input_labels
    else:
        self.file_labels = get_unique_file_labels(self.filenames)

    self.source_metadata = source_metadata
    self.extra_metadata = extra_metadata

    output_path = output_path.format_map(split_numeric_dict(self.filenames[0]))
    self.output = os.path.join(self.input_dir, output_path)    # preserve trailing slash: do not use os.path.normpath()
    output_dir = os.path.dirname(self.output)
    if self.clear:
        shutil.rmtree(output_dir, ignore_errors=True)
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    return True

init_params(params_general, params, label='', input_path=None, global_rotation=None, global_center=None)

Source code in src\muvis_align\MVSRegistration.py

def init_params(self, params_general, params, label='', input_path=None, global_rotation=None, global_center=None):
    self.params_general = params_general
    self.params = params
    self.input_params = params.get('input')
    if isinstance(self.input_params, (str, list)):
        self.input_params = {'path': self.input_params}
    if input_path is None:
        input_path = self.input_params.get('path')
    self.output_params = params.get('output')
    if isinstance(self.output_params, str):
        self.output_params = {'path': self.output_params}
    self.preprocess_params = params.get('preprocessing', {})
    self.register_params = params.get('registration', {})
    self.fusion_params = params.get('fusion', {})

    return self.init(operation=params.get('operation'), label=label, input_path=input_path,
                     input_labels=self.input_params.get('labels'),
                     output_path=self.output_params.get('path'),
                     source_metadata=self.input_params.get('source_metadata', {}),
                     extra_metadata=self.input_params.get('extra_metadata', {}),
                     global_rotation=global_rotation, global_center=global_center,
                     overwrite=params_general.get('overwrite', False), clear=params_general.get('clear', False),
                     ui=params_general.get('ui', ''),
                     verbose=params_general.get('verbose', False), debug=params_general.get('debug', False))

init_sims(source_metadata={}, extra_metadata={}, z_scale=None, chunk_size=default_chunk_size, target_scale=None)

Source code in src\muvis_align\MVSRegistration.py

def init_sims(self, source_metadata={}, extra_metadata={}, z_scale=None, chunk_size=default_chunk_size,
              target_scale=None):
    if not source_metadata:
        source_metadata = self.source_metadata
    if not extra_metadata:
        extra_metadata = self.extra_metadata
    source_metadata = import_metadata(source_metadata, input_path=self.input_path)
    extra_metadata = import_metadata(extra_metadata, input_path=self.input_path)
    source_metadata_changed = (source_metadata != self.source_metadata)
    self.source_metadata = source_metadata
    self.extra_metadata = extra_metadata
    if isinstance(source_metadata, dict):
        z_scale = source_metadata.get('scale', {}).get('z')
    if not z_scale and isinstance(extra_metadata, dict):
        z_scale = extra_metadata.get('scale', {}).get('z')

    if len(self.filenames) == 0:
        raise ValueError('No input files')

    logging.info('Initialising sims...')
    if not self.sources or source_metadata_changed:
        self.init_sources()
    sources = self.sources
    source0 = sources[0]
    images = []
    sims = []
    scales = []
    translations = []
    rotations = []

    is_stack = ('stack' in self.operation)
    has_z_size = (source0.get_size().get('z', 0) > 0)

    output_order = 'zyx' if has_z_size else 'yx'
    ndims = len(output_order)
    if source0.get_nchannels() > 1:
        output_order += 'c'

    last_z_position = None
    different_z_positions = False
    delta_zs = []
    for filename, source in zip(self.filenames, sources):
        scale = source.get_pixel_size()
        translation = source.get_position()
        rotation = source.get_rotation()

        if 'is_center' in source_metadata:
            translation = {dim: translation[dim] - source.get_physical_size().get(dim, 0) / 2 for dim in translation}

        if 'sbem' in source_metadata:
            source_version = source.metadata.get('Creator', source.metadata.get('creator', ''))
            if '2025' in source_version:
                path = os.path.dirname(self.filenames[0])
                metapath = None
                attempts = 0
                while attempts < 3:
                    metapath = os.path.join(path, 'meta')
                    if os.path.exists(metapath):
                        break
                    path = os.path.join(path, '..')
                    attempts += 1
                if metapath:
                    sbemimage_config = load_sbemimage_best_config(metapath, filename)
                    if sbemimage_config:
                        size = source.get_size()
                        translation, scale0 = adjust_sbemimage_properties(translation, scale, size,
                                                                          filename, sbemimage_config)
                        if scale0:
                            scale = scale0
                        elif scale.get('x') != scale.get('y'):
                            logging.warning('SBEMimage pixel size requires correction, please provide in source metadata.')
                        logging.debug(f'Adjusted SBEMimage properties for {filename}')
                    else:
                        logging.warning(f'Could not find SBEMimage config for {filename}.')

        level = 0
        rescale = 1
        if target_scale:
            # Only downscaling
            level, rescale, scale = get_level_from_scale(source.scales, scale, target_scale)
        if 'invert' in source_metadata:
            translation['x'] = -translation['x']
            translation['y'] = -translation['y']
        if 'z' in translation:
            z_position = translation['z']
        else:
            z_position = 0
        if last_z_position is not None and z_position != last_z_position:
            different_z_positions = True
            delta_zs.append(z_position - last_z_position)
        if 'rotation' in source_metadata:
            rotation = source_metadata['rotation']
        if self.global_rotation is not None:
            rotation = self.global_rotation

        dask_data = source.get_data(level=level)
        if rescale != 1:
            new_shape = [int(size / rescale) if dim in 'xy' else 1
                         for dim, size in zip(source.dimension_order, dask_data.shape)]
            dask_data = resize(dask_data, new_shape, preserve_range=True).astype(dask_data.dtype)
        image = redimension_data(dask_data, source.dimension_order, output_order)

        scales.append(scale)
        translations.append(translation)
        rotations.append(rotation)
        images.append(image)
        last_z_position = z_position

    if 'z' in output_order and z_scale is None:
        if len(delta_zs) > 0:
            z_scale = np.min(delta_zs)
        else:
            z_scale = 1

    if 'norm' in source_metadata:
        sizes = [source.get_physical_size() for source in sources]
        center = {dim: 0 for dim in output_order}
        if 'center' in source_metadata:
            if 'global' in source_metadata:
                center = self.global_center
            else:
                center = {dim: float(np.mean([translation[dim] for translation in translations])) for dim in translations[0]}
        translations, rotations = normalise_rotated_positions(translations, rotations, sizes, center, len(output_order))

    #translations = [np.array(translation) * 1.25 for translation in translations]

    increase_z_positions = is_stack and not different_z_positions

    z_position = 0
    final_scales = []
    final_translations = []
    for source, image, scale, translation, rotation, file_label in zip(sources, images, scales, translations, rotations, self.file_labels):
        # transform #dimensions need to match
        if 'z' in output_order:
            if len(scale) > 0 and 'z' not in scale:
                scale['z'] = abs(z_scale)
            if (len(translation) > 0 and 'z' not in translation) or increase_z_positions:
                translation['z'] = z_position
            if increase_z_positions:
                z_position += z_scale
        channel_labels = [channel.get('label', '') for channel in source.get_channels()]
        if rotation is None or 'norm' in source_metadata:
            # if positions are normalised, don't use rotation
            transform = None
        else:
            transform = param_utils.invert_coordinate_order(
                create_transform(translation, rotation, matrix_size=ndims + 1)
            )
        if file_label in extra_metadata:
            transform2 = extra_metadata[file_label]
            if transform is None:
                transform = np.array(transform2)
            else:
                transform = np.array(combine_transforms([transform, transform2]))

        # fix empty dictionaries args
        scale_arg = scale if scale else None
        translation_arg = translation if translation else None

        sim = si_utils.get_sim_from_array(
            image,
            dims=list(output_order),
            scale=scale_arg,
            translation=translation_arg,
            affine=transform,
            transform_key=self.source_transform_key,
            c_coords=channel_labels
        )
        if len(sim.chunksizes.get('x')) == 1 and len(sim.chunksizes.get('y')) == 1:
            if isinstance(chunk_size, int):
                chunk_size = [chunk_size] * 2
            sim = sim.chunk(xyz_to_dict(chunk_size))
        sims.append(sim)
        final_scales.append(scale)
        final_translations.append(translation)

    self.scales = final_scales
    self.positions = final_translations
    self.rotations = rotations
    self.sims = sims
    self.state = RegState.SIMS_INIT
    return sims

init_sources()

Source code in src\muvis_align\MVSRegistration.py

def init_sources(self):
    source_metadata0 = self.source_metadata
    source_metadata = {}
    self.sources = []
    for index, (filename, label) in enumerate(zip(self.filenames, self.file_labels)):
        if isinstance(source_metadata0, dict) and label in source_metadata0:
            source_metadata = source_metadata0[label]
            position, rotation, scale = get_properties_from_transform(param_utils.affine_to_xaffine(np.array(source_metadata)))
            source_metadata = {'position': position, 'rotation': rotation, 'scale': xyz_to_dict([scale, scale])}
        else:
            if 'position' in source_metadata0:
                translation = source_metadata0['position']
                if isinstance(translation, list):
                    translation = translation[index]
                source_metadata['position'] = translation
            if 'scale' in source_metadata0:
                scale = source_metadata0['scale']
                if isinstance(scale, list):
                    scale = scale[index]
                source_metadata['scale'] = scale
            if 'rotation' in source_metadata0:
                source_metadata['rotation'] = source_metadata0['rotation']
        self.sources.append(create_dask_source(filename, source_metadata))

is_global_registered()

Source code in src\muvis_align\MVSRegistration.py

def is_global_registered(self):
    return self.state.value >= RegState.GLOBAL_REG.value

is_initialised()

Source code in src\muvis_align\MVSRegistration.py

def is_initialised(self):
    return self.state.value >= RegState.INIT.value

is_pairs_registered()

Source code in src\muvis_align\MVSRegistration.py

def is_pairs_registered(self):
    return self.state.value >= RegState.PAIRS_REG.value

preprocess(sims, flatfield_quantiles=None, normalisation=None, gaussian_sigma=None, filter_foreground=False)

Source code in src\muvis_align\MVSRegistration.py

def preprocess(self, sims,
               flatfield_quantiles=None, normalisation=None, gaussian_sigma=None, filter_foreground=False):
    modified = False
    # normalise pixel size: take max pixel size
    max_scale = {dim: max(scale.get(dim, 1) for scale in self.scales) for dim in 'xy'}
    scales0 = self.scales

    if filter_foreground:
        foreground_map = calc_foreground_map(sims)
        modified = True
    else:
        foreground_map = None
    if flatfield_quantiles is not None:
        logging.info('Flat-field correction...')
        if isinstance(flatfield_quantiles, str):
            flatfield_quantiles = [float(quantile.strip()) for quantile in flatfield_quantiles.split(',')]
        new_sims = [None] * len(sims)
        for sim_indices in group_sims_by_z(sims, self.positions):
            sims_z_set = [sims[i] for i in sim_indices]
            foreground_map_z_set = [foreground_map[i] for i in sim_indices] if foreground_map is not None else None
            new_sims_z_set = flatfield_correction(sims_z_set, self.source_transform_key, flatfield_quantiles,
                                                  foreground_map=foreground_map_z_set)
            for sim_index, sim in zip(sim_indices, new_sims_z_set):
                new_sims[sim_index] = sim
        sims = new_sims
        modified = True

    if gaussian_sigma:
        logging.info('Applying Gaussian filtering...')
        new_sims = []
        for sim, scale0 in zip(sims, scales0):
            # factor in original pixel size for gaussian sigma value
            scale = np.mean(list(scale0.values())) / np.mean(list(max_scale.values()))
            sigma = gaussian_sigma * (scale ** (1 / 3))
            new_sims.append(gaussian_filter_sim(sim, self.source_transform_key, sigma))
        sims = new_sims
        modified = True

    if normalisation is not None:
        if isinstance(normalisation, str) and normalisation.lower() in ['false', 'no', 'none', '']:
            normalisation = None
        elif isinstance(normalisation, bool) and normalisation == False:
            normalisation = None
    if normalisation:
        use_global = ('global' in str(normalisation).lower())
        if use_global:
            logging.info('Normalising (global)...')
        else:
            logging.info('Normalising (individual)...')
        sims = normalise_sims(sims, self.source_transform_key, use_global=use_global)
        modified = True

    if filter_foreground:
        logging.info('Filtering foreground images...')
        #tile_vars = np.array([np.asarray(np.std(sim)).item() for sim in sims])
        #threshold1 = np.mean(tile_vars)
        #threshold2 = np.median(tile_vars)
        #threshold3, _ = cv.threshold(np.array(tile_vars).astype(np.uint16), 0, 1, cv.THRESH_OTSU)
        #threshold = min(threshold1, threshold2, threshold3)
        #foregrounds = (tile_vars >= threshold)
        new_sims = [sim for sim, is_foreground in zip(sims, foreground_map) if is_foreground]
        logging.info(f'Foreground images: {len(new_sims)} / {len(sims)}')
        indices = np.where(foreground_map)[0]
        sims = new_sims
        modified = True
    else:
        indices = range(len(sims))
    self.register_sims = sims
    self.register_indices = indices
    return sims, indices, modified

register(sims, register_sims=None, register_indices=None, params=None)

Source code in src\muvis_align\MVSRegistration.py

def register(self, sims, register_sims=None, register_indices=None, params=None):
    self.register_pairs(sims, register_sims=register_sims, params=params)
    results = self.register_global(sims, self.msims, register_indices=register_indices, params=params)
    return results

register_global(sims, msims, register_indices=None, params=None, pairs_graph=None)

Source code in src\muvis_align\MVSRegistration.py

def register_global(self, sims, msims, register_indices=None, params=None,
                    pairs_graph=None):
    if pairs_graph is not None:
        g_reg_computed = pairs_graph
    else:
        g_reg_computed = self.pairs_graph

    sim0 = sims[0]
    ndims = si_utils.get_ndim_from_sim(sim0)

    groupwise_resolution_method = params.get('groupwise_resolution_method',
                                             params.get('registration', {}).get('groupwise_resolution_method', 'global_optimization'))
    groupwise_resolution_kwargs = {}
    if groupwise_resolution_method == 'global_optimization':
       groupwise_resolution_kwargs['transform'] = params.get('transform_type',
                                                             params.get('registration', {}).get('transform_type'))
       # transform_type options include 'translation', 'rigid', 'affine', 'similarity'

    post_registration_quality_threshold = params.get('post_registration_quality_threshold',
                                                     params.get('registration', {}).get('post_registration_quality_threshold'))
    post_registration_do_quality_filter = (post_registration_quality_threshold is not None)

    n_parallel_pairwise_regs = params.get('n_parallel_pairwise_regs',
                                          params.get('registration', {}).get('n_parallel_pairwise_regs'))
    if n_parallel_pairwise_regs is not None and n_parallel_pairwise_regs == '0':
        n_parallel_pairwise_regs = None

    plot_summary = self.mpl_ui

    # ******* start MVS registration functions

    if post_registration_do_quality_filter:
        # filter edges by quality
        g_reg_computed = mv_graph.filter_edges(
            g_reg_computed,
            threshold=post_registration_quality_threshold,
            weight_key="quality",
        )

    with dask.config.set(scheduler='threads'):
        transforms_dict, groupwise_resolution_info_dict = groupwise_resolution(
            g_reg_computed,
            method=groupwise_resolution_method,
            **groupwise_resolution_kwargs,
        )

    transforms = [
        transforms_dict[iview] for iview in sorted(g_reg_computed.nodes())
    ]

    for imsim, msim in enumerate(msims):
        msi_utils.set_affine_transform(
            msim,
            transforms[imsim],
            transform_key=self.reg_transform_key,
            base_transform_key=self.source_transform_key,
        )

    if plot_summary:
        plot_info = _plot_registration_summaries(
            msims,
            self.source_transform_key,
            self.reg_transform_key,
            g_reg_computed,
            groupwise_resolution_info_dict,
            show_plot=plot_summary,
        )
    else:
        plot_info = {}

    reg_result = {
        "params": transforms,
        "pairwise_registration": {
            "graph": g_reg_computed,
            "metrics": {
                "qualities": nx.get_edge_attributes(
                    g_reg_computed, "quality"
                )
            },
            "summary_plot": None if plot_summary is False
            else (
                plot_info['fig_pair_reg'],
                plot_info['ax_pair_reg']
            )
        },
        "groupwise_resolution": {
            "metrics": groupwise_resolution_info_dict,
            "summary_plot": None if plot_summary is False
            else (
                plot_info['fig_group_res'],
                plot_info['ax_group_res']
            )
        },
    }

    # ******* end MVS registration functions

    if register_indices is None:
        register_indices = range(len(msims))

    # copy transforms from register sims to unmodified sims
    for reg_msim, index in zip(msims, register_indices):
        si_utils.set_sim_affine(
            sims[index],
            msi_utils.get_transform_from_msim(reg_msim, transform_key=self.reg_transform_key),
            transform_key=self.reg_transform_key)

    # set missing transforms
    for sim in sims:
        if self.reg_transform_key not in si_utils.get_tranform_keys_from_sim(sim):
            si_utils.set_sim_affine(
                sim,
                param_utils.identity_transform(ndim=ndims, t_coords=[0]),
                transform_key=self.reg_transform_key)

    mappings = reg_result['params']
    # re-index from subset of sims
    residual_error_dict = reg_result.get('groupwise_resolution', {}).get('metrics', {}).get('residuals', {})
    residual_error_dict = {(register_indices[key[0]], register_indices[key[1]]): value.item()
                           for key, value in residual_error_dict.items()}
    registration_qualities_dict = reg_result.get('pairwise_registration', {}).get('metrics', {}).get('qualities', {})
    registration_qualities_dict = {(register_indices[key[0]], register_indices[key[1]]): value
                                   for key, value in registration_qualities_dict.items()}

    # re-index from subset of sims
    mappings_dict = {index: mapping for index, mapping in zip(register_indices, mappings)}

    reg_channel = params.get('channel', 0)
    metrics = calc_global_metrics(msims, self.source_transform_key, self.reg_transform_key,
                                  params.get('metrics', []), reg_channel=reg_channel, reg_results=reg_result,
                                  n_parallel_pairs=n_parallel_pairwise_regs)

    self.metrics = metrics
    self.state = RegState.GLOBAL_REG
    return {'reg_result': reg_result,
            'mappings': mappings_dict,
            'residual_errors': residual_error_dict,
            'registration_qualities': registration_qualities_dict,
            'metrics': metrics}

register_pairs(sims, register_sims=None, params=None)

Source code in src\muvis_align\MVSRegistration.py

def register_pairs(self, sims, register_sims=None, params=None):
    operation = self.operation
    pairing = params.get('pairing',
                         params.get('registration', {}).get('pairing', '')).lower()
    n_parallel_pairwise_regs = params.get('n_parallel_pairwise_regs',
                                          params.get('registration', {}).get('n_parallel_pairwise_regs'))
    if n_parallel_pairwise_regs is not None and n_parallel_pairwise_regs == '0':
        n_parallel_pairwise_regs = None

    is_stack = ('stack' in operation)
    is_3d = ('3d' in operation)

    reg_channel = params.get('channel', 0)
    if isinstance(reg_channel, int):
        reg_channel_index = reg_channel
        reg_channel = None
    else:
        reg_channel_index = None

    if register_sims is None:
        register_sims = sims
    if is_stack and not is_3d:
        # register in 2d; pairwise consecutive views
        register_sims = [si_utils.max_project_sim(sim, dim='z') if 'z' in sim.dims else sim
                         for sim in register_sims]
        pairs = [(index, index + 1) for index in range(len(register_sims) - 1)]
    elif 'ortho' in pairing or 'overla' in pairing:
        origins = np.array([get_sim_position_final(sim, position, get_center=True)
                            for sim, position in zip(sims, self.positions)])
        sizes = [get_sim_physical_size(sim) for sim in sims]
        pairs, _ = get_pairs(origins, sizes, pairing)
        logging.info(f'#pairs: {len(pairs)}')
        #for pair in pairs:
        #    print(f'{self.file_labels[pair[0]]} - {self.file_labels[pair[1]]}')
    else:
        pairs = None

    reg_method, pairwise_reg_func, pairwise_reg_func_kwargs = self.create_registration_method(register_sims[0],
                                                                                              params=params)
    logging.info(f'Registration method: {reg_method}')
    logging.info('Registering...')
    register_msims = [msi_utils.get_msim_from_sim(sim) for sim in register_sims]

    overlap_tolerance = 0

    # ******* start MVS registration functions

    if "c" in msi_utils.get_dims(register_msims[0]):
        if reg_channel is None:
            if reg_channel_index is None:
                for msim in register_msims:
                    if "c" in msi_utils.get_dims(msim):
                        raise (
                            Exception("Please choose a registration channel.")
                        )
            else:
                reg_channel = sims[0].coords["c"][reg_channel_index]

        msims_reg = [
            msi_utils.multiscale_sel_coords(msim, {"c": reg_channel})
            if "c" in msi_utils.get_dims(msim)
            else msim
            for imsim, msim in enumerate(register_msims)
        ]
    else:
        msims_reg = register_msims

    try:
        with dask.config.set(scheduler='threads'):
            g_reg = mv_graph.build_view_adjacency_graph_from_msims(
                msims_reg,
                transform_key=self.source_transform_key,
                pairs=pairs,
                overlap_tolerance=overlap_tolerance,
            )

            g_reg_computed = compute_pairwise_registrations(
                msims_reg,
                g_reg,
                transform_key=self.source_transform_key,
                overlap_tolerance=overlap_tolerance,
                pairwise_reg_func=pairwise_reg_func,
                pairwise_reg_func_kwargs=pairwise_reg_func_kwargs,
                n_parallel_pairwise_regs=n_parallel_pairwise_regs,
            )

            # ******* end MVS registration functions

            # reg_result = registration.register(
            #     register_msims,
            #     reg_channel=reg_channel,
            #     reg_channel_index=reg_channel_index,
            #     transform_key=self.source_transform_key,
            #     new_transform_key=self.reg_transform_key,
            #
            #     pairs=pairs,
            #     pre_registration_pruning_method=None,
            #
            #     pairwise_reg_func=pairwise_reg_func,
            #     pairwise_reg_func_kwargs=pairwise_reg_func_kwargs,
            #
            #     groupwise_resolution_method=groupwise_resolution_method,
            #     groupwise_resolution_kwargs=groupwise_resolution_kwargs,
            #
            #     post_registration_do_quality_filter=(post_registration_quality_threshold is not None),
            #     post_registration_quality_threshold=post_registration_quality_threshold,
            #
            #     n_parallel_pairwise_regs=n_parallel_pairwise_regs,
            #
            #     plot_summary=self.mpl_ui,
            #     return_dict=return_dict,
            # )

    except NotEnoughOverlapError:
        g_reg_computed = g_reg

    metrics = calc_pair_metrics(msims_reg, g_reg_computed, params.get('metrics', []), self.source_transform_key,
                                reg_channel=reg_channel_index, n_parallel_pairs=n_parallel_pairwise_regs)

    self.pairs_graph = g_reg_computed
    self.msims = msims_reg
    self.pairs = pairs
    self.metrics = metrics
    self.state = RegState.PAIRS_REG
    return {
        'pairs_graph': self.pairs_graph,
        'msims': msims_reg,
        'pairs': pairs,
        'metrics': metrics
    }

run()

Source code in src\muvis_align\MVSRegistration.py

def run(self):
    with ProgressBar(minimum=60, dt=1) if self.logging_dask else nullcontext():
        return self._run()

save(output_filename, data, format=zarr_extension, transform_key=None, translations0=None, tile_size=None, compression=None, pyramid_downsample=2, npyramid_add=0, ome_version=default_ome_zarr_version)

Source code in src\muvis_align\MVSRegistration.py

def save(self, output_filename, data, format=zarr_extension, transform_key=None, translations0=None,
         tile_size=None, compression=None, pyramid_downsample=2, npyramid_add=0, ome_version=default_ome_zarr_version):
    if output_filename is not None:
        output_filename = self.output + output_filename
    if isinstance(self.extra_metadata, dict):
        channels = self.extra_metadata.get('channels', [])
    else:
        channels = []
    save_image(output_filename, data, format,
               transform_key=transform_key, channels=channels, translations0=translations0,
               tile_size=tile_size, compression=compression,
               pyramid_downsample=pyramid_downsample, npyramid_add=npyramid_add,
               ome_version=ome_version,
               verbose=self.verbose)

save_thumbnail(output_filename, nom_sims=None, transform_key=None)

Source code in src\muvis_align\MVSRegistration.py

def save_thumbnail(self, output_filename, nom_sims=None, transform_key=None):
    output_params = self.params_general['output']
    thumbnail_scale = output_params.get('thumbnail_scale', 16)
    is_stack = ('stack' in self.operation)
    if isinstance(self.extra_metadata, dict):
        z_scale = self.extra_metadata.get('scale', {}).get('z')
    else:
        z_scale = None

    sims = self.init_sims(target_scale=thumbnail_scale)
    if is_stack:
        sims = make_sims_3d(sims, z_scale, self.positions)

    if nom_sims is not None:
        if sims[0].sizes['x'] >= nom_sims[0].sizes['x']:
            logging.warning('Unable to generate scaled down thumbnail due to lack of source pyramid sizes')
            return

        if transform_key is not None and transform_key != self.source_transform_key:
            for nom_sim, sim in zip(nom_sims, sims):
                si_utils.set_sim_affine(sim,
                                        si_utils.get_affine_from_sim(nom_sim, transform_key=transform_key),
                                        transform_key=transform_key)
    fused_image, is_saved = self.fuse(sims, transform_key=transform_key, output_spacing='max',
                                      output_filename=output_filename)
    if not is_saved or 'tif' in output_params.get('thumbnail'):
        self.save(output_filename, fused_image.squeeze(), transform_key=transform_key,
                  format=output_params.get('thumbnail'), ome_version=output_params.get('ome_version'))
    self.state = RegState.FUSED

save_video(output, sims, fused_image)

Source code in src\muvis_align\MVSRegistration.py

def save_video(self, output, sims, fused_image):
    logging.info('Creating transition video...')
    pixel_size = [si_utils.get_spacing_from_sim(sims[0]).get(dim, 1) for dim in 'xy']
    params = self.params
    nframes = params.get('frames', 1)
    spacing = params.get('spacing', [1.1, 1])
    scale = params.get('scale', 1)
    transition_filename = output + 'transition'
    video = Video(transition_filename + '.mp4', fps=params.get('fps', 1))
    positions0 = np.array([si_utils.get_origin_from_sim(sim, asarray=True) for sim in sims])
    center = np.mean(positions0, 0)
    window = get_image_window(fused_image)

    max_size = None
    acum = 0
    for framei in range(nframes):
        c = (1 - np.cos(framei / (nframes - 1) * 2 * math.pi)) / 2
        acum += c / (nframes / 2)
        spacing1 = spacing[0] + (spacing[1] - spacing[0]) * acum
        for sim, position0 in zip(sims, positions0):
            transform = param_utils.identity_transform(ndim=2, t_coords=[0])
            transform[0][:2, 2] += (position0 - center) * spacing1
            si_utils.set_sim_affine(sim, transform, transform_key=self.transition_transform_key)
        frame = fusion.fuse(sims, transform_key=self.transition_transform_key).squeeze()
        frame = float2int_image(normalise_values(frame, window[0], window[1]))
        frame = cv.resize(np.asarray(frame), None, fx=scale, fy=scale)
        if max_size is None:
            max_size = frame.shape[1], frame.shape[0]
            video.size = max_size
        frame = image_reshape(frame, max_size)
        save_tiff(transition_filename + f'{framei:04d}.tiff', frame, None, pixel_size)
        video.write(frame)

    video.close()

validate_overlap(sims, labels, is_stack=False, expect_large_overlap=False)

Source code in src\muvis_align\MVSRegistration.py

def validate_overlap(self, sims, labels, is_stack=False, expect_large_overlap=False):
    min_dists = []
    has_overlaps = []
    n = len(sims)
    positions = [get_sim_position_final(sim, get_center=True) for sim in sims]
    sizes = [float(np.linalg.norm(list(get_sim_physical_size(sim).values()))) for sim in sims]
    for i in range(n):
        norm_dists = []
        # check if only single z slices
        if is_stack:
            if i + 1 < n:
                compare_indices = [i + 1]
            else:
                compare_indices = []
        else:
            compare_indices = range(n)
        for j in compare_indices:
            if not j == i:
                distance = math.dist(positions[i].values(), positions[j].values())
                norm_dist = distance / np.mean([sizes[i], sizes[j]])
                norm_dists.append(norm_dist)
        if len(norm_dists) > 0:
            norm_dist = min(norm_dists)
            min_dists.append(float(norm_dist))
            if norm_dist >= 1:
                logging.warning(f'{labels[i]} has no overlap')
                has_overlaps.append(False)
            elif expect_large_overlap and norm_dist > 0.5:
                logging.warning(f'{labels[i]} has small overlap')
                has_overlaps.append(False)
            else:
                has_overlaps.append(True)
    return min_dists, has_overlaps

RegState

Bases: Enum

Source code in src\muvis_align\MVSRegistration.py

class RegState(Enum):
    UNINIT = auto()
    INIT = auto()
    SIMS_INIT = auto()
    PAIRS_REG = auto()
    GLOBAL_REG = auto()
    FUSED = auto()

FUSED = auto() class-attribute instance-attribute

GLOBAL_REG = auto() class-attribute instance-attribute

INIT = auto() class-attribute instance-attribute

PAIRS_REG = auto() class-attribute instance-attribute

SIMS_INIT = auto() class-attribute instance-attribute

UNINIT = auto() class-attribute instance-attribute

adjust_sbemimage_properties(translation, scale, size, filename, sbemimage_config)

Source code in src\muvis_align\util.py

def adjust_sbemimage_properties(translation, scale, size, filename, sbemimage_config):
    cfg = ConfigParser()
    cfg.read_string(sbemimage_config)
    if bool(cfg['sys'].get('use_microtome')):
        props = cfg['microtome']
    else:
        props = cfg['sem']
    scale_x = float(props.get('stage_scale_factor_x'))
    scale_y = float(props.get('stage_scale_factor_y'))
    rotation_x = float(props.get('stage_rotation_angle_x'))
    rotation_y = float(props.get('stage_rotation_angle_y'))
    rotation_diff = rotation_x - rotation_y
    rot_mat_a = math.cos(rotation_y) / math.cos(rotation_diff)
    rot_mat_b = -math.sin(rotation_y) / math.cos(rotation_diff)
    rot_mat_c = math.sin(rotation_x) / math.cos(rotation_diff)
    rot_mat_d = math.cos(rotation_x) / math.cos(rotation_diff)
    rot_mat_determinant = rot_mat_a * rot_mat_d - rot_mat_b * rot_mat_c

    pixel_size = None
    parts = split_numeric_dict(filename)
    if 't' in parts:
        grids = cfg['grids']
        if 'r' in parts:
            grid_index = json.loads(grids.get('roi_index')).index(int(parts['r']))
        else:
            grid_index = int(parts['g'])
        pixel_size = json.loads(grids.get('pixel_size'))[grid_index] * 1e-3
    else:
        ov = cfg['overviews']
        if 'ov' in parts:
            ov_index = int(parts['ov'])
            pixel_size = json.loads(ov.get('ov_pixel_size'))[ov_index] * 1e-3
        elif '_stubov_' in filename.lower():
            pixel_size = float(ov.get('stub_ov_pixel_size')) * 1e-3

    if pixel_size:
        scale = {'x': pixel_size, 'y': pixel_size}
    else:
        scale = None

    stage_x, stage_y = translation['x'], translation['y']
    stage_x /= scale_x
    stage_y /= scale_y
    dx = ((rot_mat_d * stage_x - rot_mat_b * stage_y) / rot_mat_determinant)
    dy = ((-rot_mat_c * stage_x + rot_mat_a * stage_y) / rot_mat_determinant)
    # convert center to top/left
    physical_size = {dim: size[dim] * scale[dim] for dim in size}
    translation['x'] = dx - physical_size['x'] / 2
    translation['y'] = dy - physical_size['y'] / 2

    return translation, scale

apply_transform(points, transform)

Source code in src\muvis_align\util.py

def apply_transform(points, transform):
    new_points = []
    for point in points:
        point_len = len(point)
        while len(point) < len(transform):
            point = list(point) + [1]
        new_point = np.dot(point, np.transpose(np.array(transform)))
        new_points.append(new_point[:point_len])
    return new_points

apply_transform_dict(points, transform, transform_dims='xyz')

Source code in src\muvis_align\util.py

def apply_transform_dict(points, transform, transform_dims='xyz'):
    new_points = []
    for point in points:
        point = dict_to_xyz(point, dims=transform_dims)
        while len(point) < max(len(transform), 3):
            point = list(point) + [1]
        new_point = np.dot(point, np.transpose(np.array(transform)))
        new_point = xyz_to_dict(new_point, dims=transform_dims)
        new_point.pop('1', None)
        new_points.append(new_point)
    return new_points

check_round_significants(a, significant_digits)

Source code in src\muvis_align\util.py

def check_round_significants(a: float, significant_digits: int) -> float:
    rounded = round_significants(a, significant_digits)
    if a != 0:
        dif = 1 - rounded / a
    else:
        dif = rounded - a
    if abs(dif) < 10 ** -significant_digits:
        return rounded
    return a

convert_rational_value(value)

Source code in src\muvis_align\util.py

def convert_rational_value(value) -> float:
    if value is not None and isinstance(value, tuple):
        if value[0] == value[1]:
            value = value[0]
        else:
            value = value[0] / value[1]
    return value

convert_to_um(value, unit)

Source code in src\muvis_align\util.py

def convert_to_um(value, unit):
    conversions = {
        'nm': 1e-3,
        'µm': 1, 'um': 1, 'micrometer': 1, 'micron': 1,
        'mm': 1e3, 'millimeter': 1e3,
        'cm': 1e4, 'centimeter': 1e4,
        'm': 1e6, 'meter': 1e6
    }
    return value * conversions.get(unit, 1)

create_transform(center, angle, matrix_size=3)

Source code in src\muvis_align\util.py

def create_transform(center, angle, matrix_size=3):
    if isinstance(center, dict):
        center = dict_to_xyz(center)
    if len(center) == 2:
        center = np.array(list(center) + [0])
    if angle is None:
        angle = 0
    r = Rotation.from_euler('z', angle, degrees=True)
    t = center - r.apply(center, inverse=True)
    transform = np.eye(matrix_size)
    transform[:3, :3] = np.transpose(r.as_matrix())
    transform[:3, -1] += t
    return transform

create_transform0(center=(0, 0), angle=0, scale=1, translate=(0, 0))

Source code in src\muvis_align\util.py

def create_transform0(center=(0, 0), angle=0, scale=1, translate=(0, 0)):
    transform = cv.getRotationMatrix2D(center[:2], angle, scale)
    transform[:, 2] += translate
    if len(transform) == 2:
        transform = np.vstack([transform, [0, 0, 1]])   # create 3x3 matrix
    return transform

desc_to_dict(desc)

Source code in src\muvis_align\util.py

def desc_to_dict(desc: str) -> dict:
    desc_dict = {}
    if desc.startswith('{'):
        try:
            metadata = ast.literal_eval(desc)
            return metadata
        except:
            pass
    for item in re.split(r'[\r\n\t|]', desc):
        item_sep = '='
        if ':' in item:
            item_sep = ':'
        if item_sep in item:
            items = item.split(item_sep)
            key = items[0].strip()
            value = items[1].strip()
            for dtype in (int, float, bool):
                try:
                    value = dtype(value)
                    break
                except:
                    pass
            desc_dict[key] = value
    return desc_dict

dict_to_xyz(dct, dims='xyz', add_zeros=False)

Source code in src\muvis_align\util.py

def dict_to_xyz(dct, dims='xyz', add_zeros=False):
    array = [dct[dim] for dim in dims if dim in dct]
    if len(array) < len(dims) and add_zeros:
        array = array + [0] * (len(dims) - len(array))
    return array

dir_regex(pattern)

Source code in src\muvis_align\util.py

def dir_regex(pattern):
    files = []
    for pattern_item in ensure_list(pattern):
        files.extend(glob.glob(pattern_item, recursive=True))
    files_sorted = sorted(files, key=lambda file: find_all_numbers(get_filetitle(file)))
    return files_sorted

draw_edge_filter(bounds)

Source code in src\muvis_align\util.py

def draw_edge_filter(bounds):
    out_image = np.zeros(np.flip(bounds))
    y, x = np.where(out_image == 0)
    points = np.transpose([x, y])

    center = np.array(bounds) / 2
    dist_center = np.abs(points / center - 1)
    position_weights = np.clip((1 - np.max(dist_center, axis=-1)) * 10, 0, 1)
    return position_weights.reshape(np.flip(bounds))

ensure_list(x)

Source code in src\muvis_align\util.py

def ensure_list(x) -> list:
    if x is None:
        return []
    elif isinstance(x, list):
        return x
    else:
        return [x]

eval_context(data, key, default_value, context)

Source code in src\muvis_align\util.py

def eval_context(data, key, default_value, context):
    value = data.get(key, default_value)
    if isinstance(value, str):
        try:
            value = value.format_map(context)
        except:
            pass
        try:
            value = eval(value, context)
        except:
            pass
    if not isinstance(value, (float, int)):
        value = default_value
    return value

export_csv(filename, data, header=None)

Source code in src\muvis_align\util.py

def export_csv(filename, data, header=None):
    with open(filename, 'w', encoding='utf8', newline='') as file:
        csvwriter = csv.writer(file)
        if header is not None:
            csvwriter.writerow(header)
        for row in data:
            csvwriter.writerow(row)

export_json(filename, data)

Source code in src\muvis_align\util.py

def export_json(filename, data):
    with open(filename, 'w', encoding='utf8') as file:
        json.dump(data, file, indent=4)

filter_dict(dict0)

Source code in src\muvis_align\util.py

def filter_dict(dict0: dict) -> dict:
    new_dict = {}
    for key, value0 in dict0.items():
        if value0 is not None:
            values = []
            for value in ensure_list(value0):
                if isinstance(value, dict):
                    value = filter_dict(value)
                values.append(value)
            if len(values) == 1:
                values = values[0]
            new_dict[key] = values
    return new_dict

filter_edge_points(points, bounds, filter_factor=0.1, threshold=0.5)

Source code in src\muvis_align\util.py

def filter_edge_points(points, bounds, filter_factor=0.1, threshold=0.5):
    center = np.array(bounds) / 2
    dist_center = np.abs(points / center - 1)
    position_weights = np.clip((1 - np.max(dist_center, axis=-1)) / filter_factor, 0, 1)
    order_weights = 1 - np.array(range(len(points))) / len(points) / 2
    weights = position_weights * order_weights
    return weights > threshold

find_all_numbers(text)

Source code in src\muvis_align\util.py

def find_all_numbers(text: str) -> list:
    return list(map(int, re.findall(r'\d+', text)))

find_target_numeric(text, target)

Source code in src\muvis_align\util.py

def find_target_numeric(text: str, target: str) -> int|None:
    parts = split_path_parts(text)
    for part in parts:
        if part.startswith(target):
            part = part.lstrip(target)
            if part.isdecimal():
                return int(part)
    return None

get_center(data, offset=(0, 0))

Source code in src\muvis_align\util.py

def get_center(data, offset=(0, 0)):
    moments = get_moments(data, offset=offset)
    if moments['m00'] != 0:
        center = get_moments_center(moments)
    else:
        center = np.mean(data, 0).flatten()  # close approximation
    return center.astype(np.float32)

get_center_from_transform(transform)

Source code in src\muvis_align\util.py

def get_center_from_transform(transform):
    # from opencv:
    # t0 = (1-alpha) * cx - beta * cy
    # t1 = beta * cx + (1-alpha) * cy
    # where
    # alpha = cos(angle) * scale
    # beta = sin(angle) * scale
    # isolate cx and cy:
    t0, t1 = transform[:2, 2]
    scale = 1
    angle = np.arctan2(transform[0][1], transform[0][0])
    alpha = np.cos(angle) * scale
    beta = np.sin(angle) * scale
    cx = (t1 + t0 * (1 - alpha) / beta) / (beta + (1 - alpha) ** 2 / beta)
    cy = ((1 - alpha) * cx - t0) / beta
    return cx, cy

get_default(x, default)

Source code in src\muvis_align\util.py

def get_default(x, default):
    return default if x is None else x

get_filetitle(filename)

Source code in src\muvis_align\util.py

def get_filetitle(filename: str) -> str:
    filebase = os.path.basename(filename)
    title = os.path.splitext(filebase)[0].rstrip('.ome')
    return title

get_label_element(elements, label)

Source code in src\muvis_align\util.py

def get_label_element(elements, label):
    for element in elements:
        if element.get('label') == label:
            return element
    return None

get_mean_nn_distance(points1, points2)

Source code in src\muvis_align\util.py

def get_mean_nn_distance(points1, points2):
    return np.mean([get_nn_distance(points1), get_nn_distance(points2)])

get_moments(data, offset=(0, 0))

Source code in src\muvis_align\util.py

def get_moments(data, offset=(0, 0)):
    moments = cv.moments((np.array(data) + offset).astype(np.float32))    # doesn't work for float64!
    return moments

get_moments_center(moments, offset=(0, 0))

Source code in src\muvis_align\util.py

def get_moments_center(moments, offset=(0, 0)):
    return np.array([moments['m10'], moments['m01']]) / moments['m00'] + np.array(offset)

get_nn_distance(points0)

Source code in src\muvis_align\util.py

def get_nn_distance(points0):
    points = list(set(map(tuple, points0)))     # get unique points
    if len(points) >= 2:
        tree = KDTree(points, leaf_size=2)
        dist, ind = tree.query(points, k=2)
        nn_distance = np.median(dist[:, 1])
    else:
        nn_distance = 1
    return nn_distance

get_pairs(positions, sizes, pairing=None)

Get pairs of orthogonal neighbors from a list of tiles. Tiles don't have to be placed on a regular grid.

Source code in src\muvis_align\util.py

def get_pairs(positions, sizes, pairing=None):
    """
    Get pairs of orthogonal neighbors from a list of tiles.
    Tiles don't have to be placed on a regular grid.
    """
    pairs = []
    angles = []
    z_positions = [position['z'] for position in positions if 'z' in position]
    ordered_z = sorted(set(z_positions))
    is_mixed_3dstack = len(ordered_z) < len(z_positions)
    for i, j in np.transpose(np.triu_indices(len(positions), 1)):
        posi, posj = positions[i], positions[j]
        sizei, sizej = sizes[i], sizes[j]
        if is_mixed_3dstack:
            # ignore z value for distance
            distance = math.dist([posi[dim] for dim in 'xy'], [posj[dim] for dim in 'xy'])
            min_distance = max([size[dim] for size in [sizei, sizej] for dim in 'xy'])
            is_same_z = (posi['z'] == posj['z'])
            is_close_z = abs(ordered_z.index(posi['z']) - ordered_z.index(posj['z'])) <= 1
            if not is_close_z:
                # if not close, discard as pair
                min_distance = 0
            elif not is_same_z:
                # for tiles in different z stack, require greater overlap
                min_distance *= 0.8
        else:
            distance = math.dist(posi.values(), posj.values())
            min_distance = max(list(sizei.values()) + list(sizej.values()))

        if pairing and 'overla' in pairing:
            ok = (distance / min_distance < 0.5)
        else:
            ok = (distance < min_distance)
        if ok:
            pairs.append((int(i), int(j)))
            vector = np.array(list(posi.values())) - np.array(list(posj.values()))
            angle = math.degrees(math.atan2(vector[1], vector[0]))
            if distance < min(list(sizei.values()) + list(sizej.values())):
                angle += 90
            while angle < -90:
                angle += 180
            while angle > 90:
                angle -= 180
            angles.append(angle)
    return pairs, angles

get_rotation_from_transform(transform, dims='xyz')

Source code in src\muvis_align\util.py

def get_rotation_from_transform(transform, dims='xyz'):
    # TODO: assume 2D rotation, expand to 3D
    # Rotation.from_matrix(transform).as_euler() only works for simple rotation matrices
    if isinstance(transform, DataArray):
        dims = transform['x_in'].data.tolist()
    x_index, y_index = dims.index('x'), dims.index('y')
    transform = np.array(transform)
    if y_index > x_index:
        rotation = np.arctan2(transform[0][1], transform[0][0])
    else:
        rotation = np.arctan2(transform[1][0], transform[1][1])
    return float(np.rad2deg(rotation))

get_scale_from_transform(transform)

Source code in src\muvis_align\util.py

def get_scale_from_transform(transform):
    scale = np.mean(np.linalg.norm(transform, axis=0)[:-1])
    return float(scale)

get_translation_from_transform(transform)

Source code in src\muvis_align\util.py

def get_translation_from_transform(transform):
    ndim = len(transform) - 1
    #translation = transform[:ndim, ndim]
    translation = apply_transform([[0] * ndim], transform)[0]
    return translation

get_unique_file_labels(filenames)

Source code in src\muvis_align\util.py

def get_unique_file_labels(filenames: list) -> list:
    file_labels = []

    ntot = len(filenames)
    parts_dic = get_unique_nums([split_numeric_dict(get_filetitle(filename)) for filename in filenames])
    parts_dic_full = get_unique_nums([split_numeric_dict(filename) for filename in filenames])
    parts_num_full = get_unique_nums([{index: value for index, value in enumerate(split_numeric(filename))} for filename in filenames])

    dic_ok = (len(set(['_'.join(parts.values()) for parts in parts_dic])) == ntot)
    full_dic_ok = (len(set(['_'.join(parts.values()) for parts in parts_dic_full])) == ntot)
    full_num_ok = (len(set(['_'.join(parts.values()) for parts in parts_num_full])) == ntot)

    if dic_ok:
        all_parts = parts_dic
    elif full_dic_ok:
        all_parts = parts_dic_full
    elif full_num_ok:
        all_parts = parts_num_full
    else:
        all_parts = [{0: filename} for filename in filenames]

    for parts in all_parts:
        file_label = '_'.join([key + part if isinstance(key, str) else part for key, part in parts.items()])
        file_labels.append(file_label)

    if len(set(file_labels)) < len(file_labels):
        # fallback for duplicate labels
        file_labels = [get_filetitle(filename) for filename in filenames]

    return file_labels

get_unique_nums(all_parts)

Source code in src\muvis_align\util.py

def get_unique_nums(all_parts: list) -> list:
    keys = []
    for parts in all_parts:
        for key in parts:
            if key not in keys:
                keys.append(key)

    changing_keys = []
    for key in keys:
        values = [parts.get(key) for parts in all_parts]
        if len(set(values)) > 1:
            changing_keys.append(key)

    final_parts = [{key: parts[key] for key in changing_keys if key in parts} for parts in all_parts]
    return final_parts

get_value_units_micrometer(value_units0)

Source code in src\muvis_align\util.py

def get_value_units_micrometer(value_units0: list|dict) -> list|dict|None:
    conversions = {
        'nm': 1e-3,
        'µm': 1, 'um': 1, 'micrometer': 1, 'micron': 1,
        'mm': 1e3, 'millimeter': 1e3,
        'cm': 1e4, 'centimeter': 1e4,
        'm': 1e6, 'meter': 1e6
    }
    if value_units0 is None:
        return None

    if isinstance(value_units0, dict):
        values_um = {}
        for dim, value_unit in value_units0.items():
            if isinstance(value_unit, (list, tuple)):
                value_um = value_unit[0] * conversions.get(value_unit[1], 1)
            else:
                value_um = value_unit
            values_um[dim] = value_um
    else:
        values_um = []
        for value_unit in value_units0:
            if isinstance(value_unit, (list, tuple)):
                value_um = value_unit[0] * conversions.get(value_unit[1], 1)
            else:
                value_um = value_unit
            values_um.append(value_um)
    return values_um

import_csv(filename)

Source code in src\muvis_align\util.py

def import_csv(filename):
    with open(filename, encoding='utf8') as file:
        data = csv.reader(file)
    return data

import_json(filename)

Source code in src\muvis_align\util.py

def import_json(filename):
    with open(filename, encoding='utf8') as file:
        data = json.load(file)
    return data

import_metadata(content, fields=None, input_path=None)

Source code in src\muvis_align\util.py

def import_metadata(content, fields=None, input_path=None):
    # return dict[id] = {values}
    if isinstance(content, str):
        ext = os.path.splitext(content)[1].lower()
        if input_path:
            if isinstance(input_path, list):
                input_path = input_path[0]
            content = os.path.normpath(os.path.join(os.path.dirname(input_path), content))
        if ext == '.csv':
            content = import_csv(content)
        elif ext in ['.json', '.ome.json']:
            content = import_json(content)
    if fields is not None:
        content = [[data[field] for field in fields] for data in content]
    return content

is_valid_value(value)

Source code in src\muvis_align\util.py

def is_valid_value(value):
    return value is not None and value != ''

load_sbemimage_best_config(metapath, filename)

Source code in src\muvis_align\util.py

def load_sbemimage_best_config(metapath, filename):
    target_datetime = datetime.fromtimestamp(os.path.getmtime(filename))

    for config_filename in sorted(glob.glob(os.path.join(metapath, 'logs/config_*.txt')), reverse=True):
        match = re.split(r'config_(\d+-\d+-\d+).txt', config_filename)
        if len(match) >= 2:
            file_date = datetime.strptime(match[1], '%Y-%m-%d%H%M%S%f')
            if file_date <= target_datetime:
                with open(config_filename, 'r') as file:
                    sbemimage_config = file.read()
                return sbemimage_config
    return None

metric_to_color(value)

Source code in src\muvis_align\util.py

def metric_to_color(value):
    # metric range 0...1 to traffic-light color
    if value > 0.5:
        color = 'green'
    elif value > 0.25:
        color = 'gold'
    elif value > 0.1:
        color = 'orange'
    else:
        color = 'red'
    return color

metric_to_rgb(value, min_light=0, max_light=1, range=1.0)

Source code in src\muvis_align\util.py

def metric_to_rgb(value, min_light=0, max_light=1, range=1.0):
    # metric range 0...1 to red-yellow-green ranged rgb
    colormap = colormaps.get('RdYlGn')
    index = int(value * colormap.N)
    r, g, b, a = [float(value) for value in colormap(index)]
    light = 0.2125 * r + 0.7154 * g + 0.0721 * b
    if light < min_light:
        factor = light / min_light
        r = 1 - (1 - r) * factor
        g = 1 - (1 - g) * factor
        b = 1 - (1 - b) * factor
    elif light > max_light:
        factor = max_light / light
        r *= factor
        g *= factor
        b *= factor
    r *= range
    g *= range
    b *= range
    if isinstance(range, int):
        r, g, b = int(r), int(g), int(b)
    return r, g, b

normalise_rotated_positions(centers0, rotations0, sizes, center, ndims)

Source code in src\muvis_align\util.py

def normalise_rotated_positions(centers0, rotations0, sizes, center, ndims):
    # in [xy(z)]
    centers = []
    rotations = []
    _, angles = get_pairs(centers0, sizes)
    for center0, rotation in zip(centers0, rotations0):
        if rotation is None and len(angles) > 0:
            rotation = -float(np.mean(angles))
        angle = -rotation if rotation is not None else None
        transform = create_transform(center=center, angle=angle, matrix_size=ndims + 1)
        center = apply_transform_dict([center0], transform)[0]
        centers.append(center)
        rotations.append(rotation)
    return centers, rotations

normalise_rotation(rotation)

Normalise rotation to be in the range [-180, 180].

Source code in src\muvis_align\util.py

def normalise_rotation(rotation):
    """
    Normalise rotation to be in the range [-180, 180].
    """
    while rotation < -180:
        rotation += 360
    while rotation > 180:
        rotation -= 360
    return rotation

numpy_to_native(value)

Source code in src\muvis_align\util.py

def numpy_to_native(value):
    if isinstance(value, np.ndarray):
        return value.tolist()
    elif isinstance(value, list):
        return [numpy_to_native(v) for v in value]
    elif isinstance(value, tuple):
        return (numpy_to_native(v) for v in value)
    elif isinstance(value, dict):
        return {k: numpy_to_native(v) for k, v in value.items()}
    elif hasattr(value, 'dtype'):
        return value.item()
    else:
        return value

points_to_3d(points)

Source code in src\muvis_align\util.py

def points_to_3d(points):
    return [list(point) + [0] for point in points]

print_dict(dct, indent=0)

Source code in src\muvis_align\util.py

def print_dict(dct: dict, indent: int = 0) -> str:
    s = ''
    if isinstance(dct, dict):
        for key, value in dct.items():
            s += '\n'
            if not isinstance(value, list):
                s += '\t' * indent + str(key) + ': '
            if isinstance(value, dict):
                s += print_dict(value, indent=indent + 1)
            elif isinstance(value, list):
                for v in value:
                    s += print_dict(v)
            else:
                s += str(value)
    else:
        s += str(dct)
    return s

print_dict_simple(dct)

Source code in src\muvis_align\util.py

def print_dict_simple(dct: dict) -> str:
    items = []
    for key, value in dct.items():
        if isinstance(value, float):
            value = f'{value:.3f}'
        items.append(f'{key}: {value}')
    return ' '.join(items)

print_dict_xyz(dct, dims='xyz', decimals=3, as_tuple=False)

Source code in src\muvis_align\util.py

def print_dict_xyz(dct: dict, dims='xyz', decimals=3, as_tuple=False) -> str:
    s = ''
    for dim in dims:
        if dim in dct:
            if s:
                s += ' '
            if as_tuple:
                s += f'{dct[dim]:.{decimals}f}'
            else:
                s += f'{dim}:{dct[dim]:.{decimals}f}'
    return s

print_hbytes(nbytes)

Source code in src\muvis_align\util.py

def print_hbytes(nbytes: int) -> str:
    exps = ['', 'K', 'M', 'G', 'T', 'P', 'E']
    div = 1024
    exp = 0

    while nbytes > div:
        nbytes /= div
        exp += 1
    if exp < len(exps):
        e = exps[exp]
    else:
        e = f'e{exp * 3}'
    return f'{nbytes:.1f}{e}B'

reorder(items, old_order, new_order, default_value=0)

Source code in src\muvis_align\util.py

def reorder(items: list, old_order: str, new_order: str, default_value: int = 0) -> list:
    new_items = []
    for label in new_order:
        if label in old_order:
            item = items[old_order.index(label)]
        else:
            item = default_value
        new_items.append(item)
    return new_items

retuple(chunks, shape)

Expand chunks to match shape.

E.g. if chunks is (64, 64) and shape is (3, 4, 5, 1028, 1028) return (3, 4, 5, 64, 64)

If chunks is an integer, it is applied to all dimensions, to match the behaviour of zarr-python.

Source code in src\muvis_align\util.py

def retuple(chunks, shape):
    # from ome-zarr-py
    """
    Expand chunks to match shape.

    E.g. if chunks is (64, 64) and shape is (3, 4, 5, 1028, 1028)
    return (3, 4, 5, 64, 64)

    If chunks is an integer, it is applied to all dimensions, to match
    the behaviour of zarr-python.
    """

    if isinstance(chunks, int):
        return tuple([chunks] * len(shape))

    dims_to_add = len(shape) - len(chunks)
    return *shape[:dims_to_add], *chunks

round_significants(a, significant_digits)

Source code in src\muvis_align\util.py

def round_significants(a: float, significant_digits: int) -> float:
    if a != 0:
        round_decimals = significant_digits - int(np.floor(np.log10(abs(a)))) - 1
        return round(a, round_decimals)
    return a

set_dict_value(dct, keys, value)

Source code in src\muvis_align\util.py

def set_dict_value(dct, keys, value):
    try:
        value = float(value)
    except:
        pass
    for index, key in enumerate(keys):
        if index == len(keys) - 1:
            dct[key] = value
        else:
            if key not in dct:
                dct[key] = {}
            dct = dct[key]

split_num_text(text)

Source code in src\muvis_align\util.py

def split_num_text(text: str) -> list:
    num_texts = []
    block = ''
    is_num0 = None
    if text is None:
        return []

    for c in text:
        is_num = (c.isnumeric() or c == '.')
        if is_num0 is not None and is_num != is_num0:
            num_texts.append(block)
            block = ''
        block += c
        is_num0 = is_num
    if block != '':
        num_texts.append(block)

    num_texts2 = []
    for block in num_texts:
        block = block.strip()
        try:
            block = float(block)
        except:
            pass
        if block not in [' ', ',', '|']:
            num_texts2.append(block)
    return num_texts2

split_numeric(text)

Source code in src\muvis_align\util.py

def split_numeric(text: str) -> list:
    num_parts = []
    parts = split_path_parts(text)
    for part in parts:
        num_span = re.search(r'\d+', part)
        if num_span:
            num_parts.append(part)
    return num_parts

split_numeric_dict(text)

Source code in src\muvis_align\util.py

def split_numeric_dict(text: str) -> dict:
    num_parts = {}
    parts = split_path_parts(text)
    parti = 0
    for part in parts:
        num_span = re.search(r'\d+', part)
        if num_span:
            index = num_span.start()
            label = part[:index]
            if label == '':
                label = parti
            num_parts[label] = num_span.group()
            parti += 1
    return num_parts

split_path(path)

Source code in src\muvis_align\util.py

def split_path(path: str) -> list:
    return os.path.normpath(path).split(os.path.sep)

split_path_parts(text)

Source code in src\muvis_align\util.py

def split_path_parts(text: str) -> list:
    return text.replace('/', '_').replace('\\', '_').replace('.', '_').split('_')

split_value_unit_list(text)

Source code in src\muvis_align\util.py

def split_value_unit_list(text: str) -> list:
    value_units = []
    if text is None:
        return None

    items = split_num_text(text)
    if isinstance(items[-1], str):
        def_unit = items[-1]
    else:
        def_unit = ''

    i = 0
    while i < len(items):
        value = items[i]
        if i + 1 < len(items):
            unit = items[i + 1]
        else:
            unit = ''
        if not isinstance(value, str):
            if isinstance(unit, str):
                i += 1
            else:
                unit = def_unit
            value_units.append((value, unit))
        i += 1
    return value_units

validate_transform(transform, max_scale=1.25, max_rotation=None)

Source code in src\muvis_align\util.py

def validate_transform(transform, max_scale = 1.25, max_rotation=None):
    if transform is None:
        return False
    transform = np.array(transform)
    if np.any(np.isnan(transform)):
        return False
    if np.any(np.isinf(transform)):
        return False
    if np.linalg.det(transform) == 0:
        return False
    scale = get_scale_from_transform(transform)
    if scale < 1 / max_scale or scale > max_scale:
        return False
    if  max_rotation is not None and abs(normalise_rotation(get_rotation_from_transform(transform))) > max_rotation:
        return False
    return True

xyz_to_dict(xyz, dims='xyz')

Source code in src\muvis_align\util.py

def xyz_to_dict(xyz, dims='xyz'):
    dct = {dim: float(value) for dim, value in zip(dims, xyz)}
    return dct

Pipeline

Pipeline

Bases: Thread

Source code in src\muvis_align\Pipeline.py

class Pipeline(Thread):
    def __init__(self, params):
        super().__init__()
        self.params = params

        self.params_general = params['general']
        params_logging = self.params_general.get('logging', {})
        log_filename = params_logging.get('filename', 'log/muvis-align.log')
        log_format = params_logging.get('format', '%(asctime)s %(levelname)s: %(message)s')
        self.verbose = params_logging.get('verbose', False)
        init_logging(log_filename, log_format, self.verbose)

    def run(self):
        break_on_error = self.params_general.get('break_on_error', False)
        for operation_params in tqdm(self.params['operations']):
            error = False
            input_path = operation_params['input']
            if isinstance(input_path, dict):
                input_path = input_path.get('path')
            logging.info(f'Input: {input_path}')
            try:
                self.run_operation(operation_params)
            except Exception as e:
                logging.exception(f'Error processing: {input_path}')
                print(f'Error processing: {input_path}: {e}')
                error = True

            if error and break_on_error:
                break

        logging.info('Done!')

    def run_operation(self, params):
        operation = params['operation']
        use_global_metadata = 'global' in params.get('source_metadata', '')
        metadata_summary = self.params_general.get('metadata_summary', False)

        if isinstance(params['input'], dict):
            path = params['input'].get('path')
        else:
            path = params['input']
        filenames = sorted(dir_regex(path), key=lambda file: list(find_all_numbers(file)))    # sort first key first
        if len(filenames) == 0:
            logging.warning(f'Skipping operation {operation} (no files)')
            return False
        elif self.verbose:
            logging.info(f'# total files: {len(filenames)}')

        operation_parts = operation.split()
        if 'match' in operation_parts:
            # check if match label provided
            index = operation_parts.index('match') + 1
            if index < len(operation_parts):
                match_label = operation_parts[index]
            else:
                match_label = None
            matches = {}
            for filename in filenames:
                match_value = find_target_numeric(filename, match_label)
                if match_value is not None:
                    if match_value not in matches:
                        matches[match_value] = []
                    matches[match_value].append(filename)
            if len(matches) == 0:
                matches[0] = filenames
            filesets = []
            fileset_labels = []
            for label in sorted(matches):
                filesets.append(matches[label])
                fileset_labels.append(f'{match_label}:{label}')
            logging.info(f'# matched file sets: {len(filesets)}')
        else:
            filesets = [filenames]
            fileset_labels = [get_filetitle(filename) for filename in filenames]

        metadatas = []
        rotations = []
        global_center = None
        if metadata_summary or use_global_metadata:
            for fileset, fileset_label in zip(filesets, fileset_labels):
                metadata = get_images_metadata(fileset, params.get('source_metadata'))
                if metadata_summary:
                    logging.info(f'File set: {fileset_label} metadata:\n' + metadata['summary'])
                metadatas.append(metadata)
            if use_global_metadata:
                global_center = {dim: np.mean([metadata['center'][dim] for metadata in metadatas]) for dim in metadatas[0]['center']}
                rotations = [metadata['rotation'] for metadata in metadatas]
                # fix missing rotation values
                rotations = pd.Series(rotations).interpolate(limit_direction='both').to_numpy()

        n_set_workers = params.get('n_set_workers', 1)
        if n_set_workers > 1:
            value_sets = []
            for index, (fileset, fileset_label) in enumerate(zip(filesets, fileset_labels)):
                center = global_center if use_global_metadata else None
                rotation = rotations[index] if use_global_metadata else None
                value_sets.append({'fileset_label': fileset_label,
                                   'fileset': fileset,
                                   'params': params,
                                   'center': center,
                                   'rotation': rotation})

            with ThreadPoolExecutor(max_workers=n_set_workers) as executor:
                oks = executor.map(lambda kwargs: self.run_operation_thread(**kwargs), value_sets)
            return np.all([ok for ok in oks])
        else:
            ok = False
            for index, (fileset, fileset_label) in enumerate(zip(filesets, fileset_labels)):
                center = global_center if use_global_metadata else None
                rotation = rotations[index] if use_global_metadata else None
                ok |= self.run_operation_thread(fileset_label=fileset_label, fileset=fileset, params=params,
                                                 center=center, rotation=rotation)
            return ok

    def run_operation_thread(self, fileset_label, fileset, params, center, rotation):
        if fileset_label:
            logging.info(f'File set: {fileset_label}')

        mvs_registration = MVSRegistration()
        mvs_registration.init_params(params_general=self.params_general, params=params,
                                     label=fileset_label, input_path=fileset,
                                     global_center=center, global_rotation=rotation)
        return mvs_registration.run()

params = params instance-attribute

params_general = params['general'] instance-attribute

verbose = params_logging.get('verbose', False) instance-attribute

__init__(params)

Source code in src\muvis_align\Pipeline.py

def __init__(self, params):
    super().__init__()
    self.params = params

    self.params_general = params['general']
    params_logging = self.params_general.get('logging', {})
    log_filename = params_logging.get('filename', 'log/muvis-align.log')
    log_format = params_logging.get('format', '%(asctime)s %(levelname)s: %(message)s')
    self.verbose = params_logging.get('verbose', False)
    init_logging(log_filename, log_format, self.verbose)

run()

Source code in src\muvis_align\Pipeline.py

def run(self):
    break_on_error = self.params_general.get('break_on_error', False)
    for operation_params in tqdm(self.params['operations']):
        error = False
        input_path = operation_params['input']
        if isinstance(input_path, dict):
            input_path = input_path.get('path')
        logging.info(f'Input: {input_path}')
        try:
            self.run_operation(operation_params)
        except Exception as e:
            logging.exception(f'Error processing: {input_path}')
            print(f'Error processing: {input_path}: {e}')
            error = True

        if error and break_on_error:
            break

    logging.info('Done!')

run_operation(params)

Source code in src\muvis_align\Pipeline.py

def run_operation(self, params):
    operation = params['operation']
    use_global_metadata = 'global' in params.get('source_metadata', '')
    metadata_summary = self.params_general.get('metadata_summary', False)

    if isinstance(params['input'], dict):
        path = params['input'].get('path')
    else:
        path = params['input']
    filenames = sorted(dir_regex(path), key=lambda file: list(find_all_numbers(file)))    # sort first key first
    if len(filenames) == 0:
        logging.warning(f'Skipping operation {operation} (no files)')
        return False
    elif self.verbose:
        logging.info(f'# total files: {len(filenames)}')

    operation_parts = operation.split()
    if 'match' in operation_parts:
        # check if match label provided
        index = operation_parts.index('match') + 1
        if index < len(operation_parts):
            match_label = operation_parts[index]
        else:
            match_label = None
        matches = {}
        for filename in filenames:
            match_value = find_target_numeric(filename, match_label)
            if match_value is not None:
                if match_value not in matches:
                    matches[match_value] = []
                matches[match_value].append(filename)
        if len(matches) == 0:
            matches[0] = filenames
        filesets = []
        fileset_labels = []
        for label in sorted(matches):
            filesets.append(matches[label])
            fileset_labels.append(f'{match_label}:{label}')
        logging.info(f'# matched file sets: {len(filesets)}')
    else:
        filesets = [filenames]
        fileset_labels = [get_filetitle(filename) for filename in filenames]

    metadatas = []
    rotations = []
    global_center = None
    if metadata_summary or use_global_metadata:
        for fileset, fileset_label in zip(filesets, fileset_labels):
            metadata = get_images_metadata(fileset, params.get('source_metadata'))
            if metadata_summary:
                logging.info(f'File set: {fileset_label} metadata:\n' + metadata['summary'])
            metadatas.append(metadata)
        if use_global_metadata:
            global_center = {dim: np.mean([metadata['center'][dim] for metadata in metadatas]) for dim in metadatas[0]['center']}
            rotations = [metadata['rotation'] for metadata in metadatas]
            # fix missing rotation values
            rotations = pd.Series(rotations).interpolate(limit_direction='both').to_numpy()

    n_set_workers = params.get('n_set_workers', 1)
    if n_set_workers > 1:
        value_sets = []
        for index, (fileset, fileset_label) in enumerate(zip(filesets, fileset_labels)):
            center = global_center if use_global_metadata else None
            rotation = rotations[index] if use_global_metadata else None
            value_sets.append({'fileset_label': fileset_label,
                               'fileset': fileset,
                               'params': params,
                               'center': center,
                               'rotation': rotation})

        with ThreadPoolExecutor(max_workers=n_set_workers) as executor:
            oks = executor.map(lambda kwargs: self.run_operation_thread(**kwargs), value_sets)
        return np.all([ok for ok in oks])
    else:
        ok = False
        for index, (fileset, fileset_label) in enumerate(zip(filesets, fileset_labels)):
            center = global_center if use_global_metadata else None
            rotation = rotations[index] if use_global_metadata else None
            ok |= self.run_operation_thread(fileset_label=fileset_label, fileset=fileset, params=params,
                                             center=center, rotation=rotation)
        return ok

run_operation_thread(fileset_label, fileset, params, center, rotation)

Source code in src\muvis_align\Pipeline.py

def run_operation_thread(self, fileset_label, fileset, params, center, rotation):
    if fileset_label:
        logging.info(f'File set: {fileset_label}')

    mvs_registration = MVSRegistration()
    mvs_registration.init_params(params_general=self.params_general, params=params,
                                 label=fileset_label, input_path=fileset,
                                 global_center=center, global_rotation=rotation)
    return mvs_registration.run()

Timer

Timer

Bases: object

Source code in src\muvis_align\Timer.py

class Timer(object):
    def __init__(self, title, auto_unit=True, verbose=True):
        self.title = title
        self.auto_unit = auto_unit
        self.verbose = verbose

    def __enter__(self):
        self.ptime_start = time.process_time()
        self.time_start = time.time()

    def __exit__(self, type, value, traceback):
        if self.verbose:
            ptime_end = time.process_time()
            time_end = time.time()
            pelapsed = ptime_end - self.ptime_start
            elapsed = time_end - self.time_start
            unit = 'seconds'
            if self.auto_unit and elapsed >= 60:
                pelapsed /= 60
                elapsed /= 60
                unit = 'minutes'
                if elapsed >= 60:
                    pelapsed /= 60
                    elapsed /= 60
                    unit = 'hours'
            logging.info(f'Time {self.title}: {elapsed:.1f} ({pelapsed:.1f}) {unit}')

auto_unit = auto_unit instance-attribute

title = title instance-attribute

verbose = verbose instance-attribute

__enter__()

Source code in src\muvis_align\Timer.py

def __enter__(self):
    self.ptime_start = time.process_time()
    self.time_start = time.time()

__exit__(type, value, traceback)

Source code in src\muvis_align\Timer.py

def __exit__(self, type, value, traceback):
    if self.verbose:
        ptime_end = time.process_time()
        time_end = time.time()
        pelapsed = ptime_end - self.ptime_start
        elapsed = time_end - self.time_start
        unit = 'seconds'
        if self.auto_unit and elapsed >= 60:
            pelapsed /= 60
            elapsed /= 60
            unit = 'minutes'
            if elapsed >= 60:
                pelapsed /= 60
                elapsed /= 60
                unit = 'hours'
        logging.info(f'Time {self.title}: {elapsed:.1f} ({pelapsed:.1f}) {unit}')

__init__(title, auto_unit=True, verbose=True)

Source code in src\muvis_align\Timer.py

def __init__(self, title, auto_unit=True, verbose=True):
    self.title = title
    self.auto_unit = auto_unit
    self.verbose = verbose

constants

NAPARI_PROJECT_TEMPLATE = 'ui/project_template.yaml' module-attribute

default_chunk_size = 1024 module-attribute

default_mappings_name = 'mappings.json' module-attribute

default_ome_zarr_version = '0.5' module-attribute

metrics_name = 'metrics.json' module-attribute

metrics_tabular_name = 'mappings.csv' module-attribute

original_positions_name = 'positions_original.pdf' module-attribute

prereg_mappings_name = 'prereg_mappings.csv' module-attribute

registered_positions_name = 'positions_registered.pdf' module-attribute

tiff_extension = '.ome.tiff' module-attribute

version = '0.3.0' module-attribute

zarr_extension = '.ome.zarr' module-attribute

logging

init_logging(log_filename='log/muvis-align.log', log_format='%(asctime)s %(levelname)s: %(message)s', verbose=False)

Source code in src\muvis_align\logging.py

def init_logging(log_filename='log/muvis-align.log', log_format='%(asctime)s %(levelname)s: %(message)s', verbose=False):
    basepath = os.path.dirname(log_filename)
    if basepath and not os.path.exists(basepath):
        os.makedirs(basepath)

    handlers = [logging.FileHandler(log_filename, encoding='utf-8')]
    if verbose:
        handlers += [logging.StreamHandler()]
    logging.basicConfig(level=logging.INFO, format=log_format, handlers=handlers, encoding='utf-8')

    # verbose external modules
    if verbose:
        # expose multiview_stitcher.registration logger and make more verbose
        mvsr_logger = logging.getLogger('multiview_stitcher.registration')
        mvsr_logger.setLevel(logging.DEBUG)
        if len(mvsr_logger.handlers) == 0:
            mvsr_logger.addHandler(logging.StreamHandler())
    else:
        # reduce verbose level
        for module in ['multiview_stitcher', 'multiview_stitcher.registration', 'multiview_stitcher.fusion']:
            logging.getLogger(module).setLevel(logging.WARNING)

    for module in ['ome_zarr']:
        logging.getLogger(module).setLevel(logging.WARNING)

    logging.info(f'muvis-align version {version}')
    logging.info(f'Multiview-stitcher version: {mvs_version}')

metrics

calc_frc(image1, image2)

Source code in src\muvis_align\metrics.py

def calc_frc(image1, image2):
    pixel_size1 = si_utils.get_spacing_from_sim(image1)
    pixel_size2 = si_utils.get_spacing_from_sim(image2)
    pixel_size = np.mean([pixel_size1['x'], pixel_size1['y'], pixel_size2['x'], pixel_size2['y']])
    max_size = np.flip(np.max([image1.shape, image2.shape], 0))
    image1 = frc.util.square_image(image_reshape(image1, max_size), add_padding=True)
    image2 = frc.util.square_image(image_reshape(image2, max_size), add_padding=True)

    frc_curve = frc.two_frc(image1, image2)
    xs_pix = np.arange(len(frc_curve)) / max(max_size)
    # scale has units [pixels <length unit>^-1] corresponding to original image
    xs_nm_freq = xs_pix / pixel_size
    frc_res, res_y, thres = frc.frc_res(xs_nm_freq, frc_curve, max_size)
    #plt.plot(xs_nm_freq, thres(xs_nm_freq))
    #plt.plot(xs_nm_freq, frc_curve)
    #plt.show()
    return frc_res

calc_global_metrics(msims, base_transform_key, reg_transform_key, metric_methods, reg_channel=None, reg_results=None, n_parallel_pairs=None)

Source code in src\muvis_align\metrics.py

def calc_global_metrics(msims, base_transform_key, reg_transform_key, metric_methods, reg_channel=None,
                        reg_results=None, n_parallel_pairs=None):
    metric_funcs = create_metric_methods(metric_methods, msims[0], reg_channel=reg_channel)
    with dask.config.set(scheduler='single-threaded'):
        metric_results = multiview_stitcher.metrics.tile_pair_image_metrics(
            msims,
            base_transform_key=base_transform_key,  # defines overlap region
            query_transform_keys=[
                base_transform_key,
                reg_transform_key
            ],
            metric_funcs=metric_funcs,
            n_parallel_pairs=n_parallel_pairs
        )

    if reg_results is not None:
        qualities = reg_results['pairwise_registration']['metrics']['qualities']

        quality_values = []
        for pair_key, value in qualities.items():
            if isinstance(value, DataArray):
                if 't' in value.dims:
                    value = value.sel(t=0)
                value = value.item()
            if value:
                metric_results['pairs'][pair_key][reg_transform_key]['quality'] = value
                quality_values.append(value)

        metric_results['summary'][reg_transform_key]['quality'] = float(np.nanmean(quality_values))

    return metric_results

calc_match_metrics(points1, points2, transform, threshold, lowe_ratio=None)

Source code in src\muvis_align\metrics.py

def calc_match_metrics(points1, points2, transform, threshold, lowe_ratio=None):
    metrics = {}
    transformed_points1 = apply_transform(points1, transform)
    npoints1, npoints2 = len(points1), len(points2)
    npoints = min(npoints1, npoints2)
    if npoints1 == 0 or npoints2 == 0:
        return metrics

    swapped = (npoints1 > npoints2)
    if swapped:
        points1, points2 = points2, points1

    distance_matrix = euclidean_distances(transformed_points1, points2)
    matching_distances = np.diag(distance_matrix)
    if npoints1 == npoints2 and np.mean(matching_distances < threshold) > 0.5:
        # already matching points lists
        nmatches = np.sum(matching_distances < threshold)
    else:
        matches = []
        distances0 = []
        for rowi, row in enumerate(distance_matrix):
            sorted_indices = np.argsort(row)
            index0 = sorted_indices[0]
            distance0 = row[index0]
            matches.append((rowi, sorted_indices))
            distances0.append(distance0)
        sorted_matches = np.argsort(distances0)

        done = []
        nmatches = 0
        matching_distances = []
        for sorted_match in sorted_matches:
            i, match = matches[sorted_match]
            for ji, j in enumerate(match):
                if j not in done:
                    # found best, available match
                    distance0 = distance_matrix[i, j]
                    distance1 = distance_matrix[i, match[ji + 1]] if ji + 1 < len(match) else np.inf
                    matching_distances.append(distance0)    # use all distances to also weigh in the non-matches
                    if distance0 < threshold and (lowe_ratio is None or distance0 < lowe_ratio * distance1):
                        done.append(j)
                        nmatches += 1
                    break

    metrics['nmatches'] = nmatches
    metrics['match_rate'] = nmatches / npoints if npoints > 0 else 0
    distance = np.mean(matching_distances) if nmatches > 0 else np.inf
    metrics['distance'] = float(distance)
    metrics['norm_distance'] = float(distance / threshold)
    return metrics

calc_ncc(image1, image2)

Source code in src\muvis_align\metrics.py

def calc_ncc(image1, image2):
    max_size = np.flip(np.max([image1.shape, image2.shape], 0))
    image1 = image_reshape(image1, max_size)
    image2 = image_reshape(image2, max_size)

    normimage1 = np.array(image1 - np.mean(image1))
    normimage2 = np.array(image2 - np.mean(image2))
    ncc = np.sum(normimage1 * normimage2) / (np.linalg.norm(normimage1) * np.linalg.norm(normimage2))
    return float(ncc)

calc_ncc2(image1, image2)

Source code in src\muvis_align\metrics.py

def calc_ncc2(image1, image2):
    max_size = np.flip(np.max([image1.shape, image2.shape], 0))
    image1 = image_reshape(image1, max_size)
    image2 = image_reshape(image2, max_size)

    normimage1 = (image1 - np.mean(image1)) / np.std(image1)
    normimage2 = (image2 - np.mean(image2)) / np.std(image2)
    array1 = np.array(normimage1).reshape(-1)
    array2 = np.array(normimage2).reshape(-1)
    ncc = (np.correlate(array1, array2) / max(len(array1), len(array2)))[0]
    return float(ncc)

calc_pair_metrics(msims, pairs_graph, metric_methods, base_transform_key, reg_channel=None, n_parallel_pairs=None)

Source code in src\muvis_align\metrics.py

def calc_pair_metrics(msims, pairs_graph, metric_methods, base_transform_key, reg_channel=None,
                      n_parallel_pairs=None):
    metric_funcs = create_metric_methods(metric_methods, msims[0], reg_channel=reg_channel)
    with dask.config.set(scheduler='single-threaded'):
        metric_results = multiview_stitcher.metrics.tile_pair_image_metrics(
            msims,
            base_transform_key=base_transform_key,  # defines overlap region
            pairs_graph=pairs_graph,
            metric_funcs=metric_funcs,
            n_parallel_pairs=n_parallel_pairs
        )

    qualities = nx.get_edge_attributes(pairs_graph, 'quality')

    quality_values = []
    for pair_key, value in qualities.items():
        if isinstance(value, DataArray):
            if 't' in value.dims:
                value = value.sel(t=0)
            value = value.item()
        if value:
            metric_results['pairs'][pair_key]['transform']['quality'] = value
            quality_values.append(value)

    value = float(np.nanmean(quality_values)) if quality_values else None
    metric_results['summary']['transform']['quality'] = value

    return metric_results

calc_sims_metrics(sims, pair_transforms, qualities, base_transform_key=None, metric_methods='all', reg_channel=None, n_parallel_pairs=None)

Source code in src\muvis_align\metrics.py

def calc_sims_metrics(sims, pair_transforms, qualities, base_transform_key=None, metric_methods='all',
                      reg_channel=None, n_parallel_pairs=None):
    if base_transform_key is None:
        reg_keys = si_utils.get_tranform_keys_from_sim(sims[0])
        base_transform_key = reg_keys[0]
    msims = [msi_utils.get_msim_from_sim(sim) for sim in sims]
    with dask.config.set(scheduler='single-threaded'):
        pairs_graph = mv_graph.build_view_adjacency_graph_from_msims(
            msims,
            transform_key=base_transform_key,
            pairs=list(pair_transforms.keys())
        )
    nx.set_edge_attributes(pairs_graph, pair_transforms, 'transform')
    nx.set_edge_attributes(pairs_graph, qualities, 'quality')
    return calc_pair_metrics(msims=msims, pairs_graph=pairs_graph, base_transform_key=base_transform_key,
                             metric_methods=metric_methods, reg_channel=reg_channel, n_parallel_pairs=n_parallel_pairs)

calc_ssim(image1, image2)

Source code in src\muvis_align\metrics.py

def calc_ssim(image1, image2):
    dtype = image1.dtype
    maxval = 2 ** (8 * dtype.itemsize) - 1
    max_size = np.flip(np.max([image1.shape, image2.shape], 0))
    image1 = image_reshape(image1, max_size)
    image2 = image_reshape(image2, max_size)
    try:
        ssim = structural_similarity(np.array(image1), np.array(image2), data_range=maxval)
    except ValueError:
        ssim = np.nan
    return float(ssim)

create_metric_methods(metric_methods, msim, reg_channel=None)

Source code in src\muvis_align\metrics.py

def create_metric_methods(metric_methods, msim, reg_channel=None):
    data_range = np.iinfo(msim["scale0/image"].dtype).max
    all_metric_funcs = {
        'ncc': multiview_stitcher.metrics.normalized_cross_correlation,
        'ssim': lambda im1, im2: structural_similarity(np.nan_to_num(im1), np.nan_to_num(im2),
                                                       data_range=data_range, channel_axis=reg_channel),
        'onmi': lambda im1, im2: normalized_mutual_information(np.nan_to_num(im1), np.nan_to_num(im2)) - 1,
        'mse': lambda im1, im2: 1 / mean_squared_error(im1, im2),
    }
    if metric_methods == 'all':
        metric_funcs = all_metric_funcs
    else:
        metric_funcs = {metric_method: all_metric_funcs[metric_method] for metric_method in metric_methods}
    return metric_funcs

util

adjust_sbemimage_properties(translation, scale, size, filename, sbemimage_config)

Source code in src\muvis_align\util.py

def adjust_sbemimage_properties(translation, scale, size, filename, sbemimage_config):
    cfg = ConfigParser()
    cfg.read_string(sbemimage_config)
    if bool(cfg['sys'].get('use_microtome')):
        props = cfg['microtome']
    else:
        props = cfg['sem']
    scale_x = float(props.get('stage_scale_factor_x'))
    scale_y = float(props.get('stage_scale_factor_y'))
    rotation_x = float(props.get('stage_rotation_angle_x'))
    rotation_y = float(props.get('stage_rotation_angle_y'))
    rotation_diff = rotation_x - rotation_y
    rot_mat_a = math.cos(rotation_y) / math.cos(rotation_diff)
    rot_mat_b = -math.sin(rotation_y) / math.cos(rotation_diff)
    rot_mat_c = math.sin(rotation_x) / math.cos(rotation_diff)
    rot_mat_d = math.cos(rotation_x) / math.cos(rotation_diff)
    rot_mat_determinant = rot_mat_a * rot_mat_d - rot_mat_b * rot_mat_c

    pixel_size = None
    parts = split_numeric_dict(filename)
    if 't' in parts:
        grids = cfg['grids']
        if 'r' in parts:
            grid_index = json.loads(grids.get('roi_index')).index(int(parts['r']))
        else:
            grid_index = int(parts['g'])
        pixel_size = json.loads(grids.get('pixel_size'))[grid_index] * 1e-3
    else:
        ov = cfg['overviews']
        if 'ov' in parts:
            ov_index = int(parts['ov'])
            pixel_size = json.loads(ov.get('ov_pixel_size'))[ov_index] * 1e-3
        elif '_stubov_' in filename.lower():
            pixel_size = float(ov.get('stub_ov_pixel_size')) * 1e-3

    if pixel_size:
        scale = {'x': pixel_size, 'y': pixel_size}
    else:
        scale = None

    stage_x, stage_y = translation['x'], translation['y']
    stage_x /= scale_x
    stage_y /= scale_y
    dx = ((rot_mat_d * stage_x - rot_mat_b * stage_y) / rot_mat_determinant)
    dy = ((-rot_mat_c * stage_x + rot_mat_a * stage_y) / rot_mat_determinant)
    # convert center to top/left
    physical_size = {dim: size[dim] * scale[dim] for dim in size}
    translation['x'] = dx - physical_size['x'] / 2
    translation['y'] = dy - physical_size['y'] / 2

    return translation, scale

apply_transform(points, transform)

Source code in src\muvis_align\util.py

def apply_transform(points, transform):
    new_points = []
    for point in points:
        point_len = len(point)
        while len(point) < len(transform):
            point = list(point) + [1]
        new_point = np.dot(point, np.transpose(np.array(transform)))
        new_points.append(new_point[:point_len])
    return new_points

apply_transform_dict(points, transform, transform_dims='xyz')

Source code in src\muvis_align\util.py

def apply_transform_dict(points, transform, transform_dims='xyz'):
    new_points = []
    for point in points:
        point = dict_to_xyz(point, dims=transform_dims)
        while len(point) < max(len(transform), 3):
            point = list(point) + [1]
        new_point = np.dot(point, np.transpose(np.array(transform)))
        new_point = xyz_to_dict(new_point, dims=transform_dims)
        new_point.pop('1', None)
        new_points.append(new_point)
    return new_points

check_round_significants(a, significant_digits)

Source code in src\muvis_align\util.py

def check_round_significants(a: float, significant_digits: int) -> float:
    rounded = round_significants(a, significant_digits)
    if a != 0:
        dif = 1 - rounded / a
    else:
        dif = rounded - a
    if abs(dif) < 10 ** -significant_digits:
        return rounded
    return a

convert_rational_value(value)

Source code in src\muvis_align\util.py

def convert_rational_value(value) -> float:
    if value is not None and isinstance(value, tuple):
        if value[0] == value[1]:
            value = value[0]
        else:
            value = value[0] / value[1]
    return value

convert_to_um(value, unit)

Source code in src\muvis_align\util.py

def convert_to_um(value, unit):
    conversions = {
        'nm': 1e-3,
        'µm': 1, 'um': 1, 'micrometer': 1, 'micron': 1,
        'mm': 1e3, 'millimeter': 1e3,
        'cm': 1e4, 'centimeter': 1e4,
        'm': 1e6, 'meter': 1e6
    }
    return value * conversions.get(unit, 1)

create_transform(center, angle, matrix_size=3)

Source code in src\muvis_align\util.py

def create_transform(center, angle, matrix_size=3):
    if isinstance(center, dict):
        center = dict_to_xyz(center)
    if len(center) == 2:
        center = np.array(list(center) + [0])
    if angle is None:
        angle = 0
    r = Rotation.from_euler('z', angle, degrees=True)
    t = center - r.apply(center, inverse=True)
    transform = np.eye(matrix_size)
    transform[:3, :3] = np.transpose(r.as_matrix())
    transform[:3, -1] += t
    return transform

create_transform0(center=(0, 0), angle=0, scale=1, translate=(0, 0))

Source code in src\muvis_align\util.py

def create_transform0(center=(0, 0), angle=0, scale=1, translate=(0, 0)):
    transform = cv.getRotationMatrix2D(center[:2], angle, scale)
    transform[:, 2] += translate
    if len(transform) == 2:
        transform = np.vstack([transform, [0, 0, 1]])   # create 3x3 matrix
    return transform

desc_to_dict(desc)

Source code in src\muvis_align\util.py

def desc_to_dict(desc: str) -> dict:
    desc_dict = {}
    if desc.startswith('{'):
        try:
            metadata = ast.literal_eval(desc)
            return metadata
        except:
            pass
    for item in re.split(r'[\r\n\t|]', desc):
        item_sep = '='
        if ':' in item:
            item_sep = ':'
        if item_sep in item:
            items = item.split(item_sep)
            key = items[0].strip()
            value = items[1].strip()
            for dtype in (int, float, bool):
                try:
                    value = dtype(value)
                    break
                except:
                    pass
            desc_dict[key] = value
    return desc_dict

dict_to_xyz(dct, dims='xyz', add_zeros=False)

Source code in src\muvis_align\util.py

def dict_to_xyz(dct, dims='xyz', add_zeros=False):
    array = [dct[dim] for dim in dims if dim in dct]
    if len(array) < len(dims) and add_zeros:
        array = array + [0] * (len(dims) - len(array))
    return array

dir_regex(pattern)

Source code in src\muvis_align\util.py

def dir_regex(pattern):
    files = []
    for pattern_item in ensure_list(pattern):
        files.extend(glob.glob(pattern_item, recursive=True))
    files_sorted = sorted(files, key=lambda file: find_all_numbers(get_filetitle(file)))
    return files_sorted

draw_edge_filter(bounds)

Source code in src\muvis_align\util.py

def draw_edge_filter(bounds):
    out_image = np.zeros(np.flip(bounds))
    y, x = np.where(out_image == 0)
    points = np.transpose([x, y])

    center = np.array(bounds) / 2
    dist_center = np.abs(points / center - 1)
    position_weights = np.clip((1 - np.max(dist_center, axis=-1)) * 10, 0, 1)
    return position_weights.reshape(np.flip(bounds))

ensure_list(x)

Source code in src\muvis_align\util.py

def ensure_list(x) -> list:
    if x is None:
        return []
    elif isinstance(x, list):
        return x
    else:
        return [x]

eval_context(data, key, default_value, context)

Source code in src\muvis_align\util.py

def eval_context(data, key, default_value, context):
    value = data.get(key, default_value)
    if isinstance(value, str):
        try:
            value = value.format_map(context)
        except:
            pass
        try:
            value = eval(value, context)
        except:
            pass
    if not isinstance(value, (float, int)):
        value = default_value
    return value

export_csv(filename, data, header=None)

Source code in src\muvis_align\util.py

def export_csv(filename, data, header=None):
    with open(filename, 'w', encoding='utf8', newline='') as file:
        csvwriter = csv.writer(file)
        if header is not None:
            csvwriter.writerow(header)
        for row in data:
            csvwriter.writerow(row)

export_json(filename, data)

Source code in src\muvis_align\util.py

def export_json(filename, data):
    with open(filename, 'w', encoding='utf8') as file:
        json.dump(data, file, indent=4)

filter_dict(dict0)

Source code in src\muvis_align\util.py

def filter_dict(dict0: dict) -> dict:
    new_dict = {}
    for key, value0 in dict0.items():
        if value0 is not None:
            values = []
            for value in ensure_list(value0):
                if isinstance(value, dict):
                    value = filter_dict(value)
                values.append(value)
            if len(values) == 1:
                values = values[0]
            new_dict[key] = values
    return new_dict

filter_edge_points(points, bounds, filter_factor=0.1, threshold=0.5)

Source code in src\muvis_align\util.py

def filter_edge_points(points, bounds, filter_factor=0.1, threshold=0.5):
    center = np.array(bounds) / 2
    dist_center = np.abs(points / center - 1)
    position_weights = np.clip((1 - np.max(dist_center, axis=-1)) / filter_factor, 0, 1)
    order_weights = 1 - np.array(range(len(points))) / len(points) / 2
    weights = position_weights * order_weights
    return weights > threshold

find_all_numbers(text)

Source code in src\muvis_align\util.py

def find_all_numbers(text: str) -> list:
    return list(map(int, re.findall(r'\d+', text)))

find_target_numeric(text, target)

Source code in src\muvis_align\util.py

def find_target_numeric(text: str, target: str) -> int|None:
    parts = split_path_parts(text)
    for part in parts:
        if part.startswith(target):
            part = part.lstrip(target)
            if part.isdecimal():
                return int(part)
    return None

get_center(data, offset=(0, 0))

Source code in src\muvis_align\util.py

def get_center(data, offset=(0, 0)):
    moments = get_moments(data, offset=offset)
    if moments['m00'] != 0:
        center = get_moments_center(moments)
    else:
        center = np.mean(data, 0).flatten()  # close approximation
    return center.astype(np.float32)

get_center_from_transform(transform)

Source code in src\muvis_align\util.py

def get_center_from_transform(transform):
    # from opencv:
    # t0 = (1-alpha) * cx - beta * cy
    # t1 = beta * cx + (1-alpha) * cy
    # where
    # alpha = cos(angle) * scale
    # beta = sin(angle) * scale
    # isolate cx and cy:
    t0, t1 = transform[:2, 2]
    scale = 1
    angle = np.arctan2(transform[0][1], transform[0][0])
    alpha = np.cos(angle) * scale
    beta = np.sin(angle) * scale
    cx = (t1 + t0 * (1 - alpha) / beta) / (beta + (1 - alpha) ** 2 / beta)
    cy = ((1 - alpha) * cx - t0) / beta
    return cx, cy

get_default(x, default)

Source code in src\muvis_align\util.py

def get_default(x, default):
    return default if x is None else x

get_filetitle(filename)

Source code in src\muvis_align\util.py

def get_filetitle(filename: str) -> str:
    filebase = os.path.basename(filename)
    title = os.path.splitext(filebase)[0].rstrip('.ome')
    return title

get_label_element(elements, label)

Source code in src\muvis_align\util.py

def get_label_element(elements, label):
    for element in elements:
        if element.get('label') == label:
            return element
    return None

get_mean_nn_distance(points1, points2)

Source code in src\muvis_align\util.py

def get_mean_nn_distance(points1, points2):
    return np.mean([get_nn_distance(points1), get_nn_distance(points2)])

get_moments(data, offset=(0, 0))

Source code in src\muvis_align\util.py

def get_moments(data, offset=(0, 0)):
    moments = cv.moments((np.array(data) + offset).astype(np.float32))    # doesn't work for float64!
    return moments

get_moments_center(moments, offset=(0, 0))

Source code in src\muvis_align\util.py

def get_moments_center(moments, offset=(0, 0)):
    return np.array([moments['m10'], moments['m01']]) / moments['m00'] + np.array(offset)

get_nn_distance(points0)

Source code in src\muvis_align\util.py

def get_nn_distance(points0):
    points = list(set(map(tuple, points0)))     # get unique points
    if len(points) >= 2:
        tree = KDTree(points, leaf_size=2)
        dist, ind = tree.query(points, k=2)
        nn_distance = np.median(dist[:, 1])
    else:
        nn_distance = 1
    return nn_distance

get_pairs(positions, sizes, pairing=None)

Get pairs of orthogonal neighbors from a list of tiles. Tiles don't have to be placed on a regular grid.

Source code in src\muvis_align\util.py

def get_pairs(positions, sizes, pairing=None):
    """
    Get pairs of orthogonal neighbors from a list of tiles.
    Tiles don't have to be placed on a regular grid.
    """
    pairs = []
    angles = []
    z_positions = [position['z'] for position in positions if 'z' in position]
    ordered_z = sorted(set(z_positions))
    is_mixed_3dstack = len(ordered_z) < len(z_positions)
    for i, j in np.transpose(np.triu_indices(len(positions), 1)):
        posi, posj = positions[i], positions[j]
        sizei, sizej = sizes[i], sizes[j]
        if is_mixed_3dstack:
            # ignore z value for distance
            distance = math.dist([posi[dim] for dim in 'xy'], [posj[dim] for dim in 'xy'])
            min_distance = max([size[dim] for size in [sizei, sizej] for dim in 'xy'])
            is_same_z = (posi['z'] == posj['z'])
            is_close_z = abs(ordered_z.index(posi['z']) - ordered_z.index(posj['z'])) <= 1
            if not is_close_z:
                # if not close, discard as pair
                min_distance = 0
            elif not is_same_z:
                # for tiles in different z stack, require greater overlap
                min_distance *= 0.8
        else:
            distance = math.dist(posi.values(), posj.values())
            min_distance = max(list(sizei.values()) + list(sizej.values()))

        if pairing and 'overla' in pairing:
            ok = (distance / min_distance < 0.5)
        else:
            ok = (distance < min_distance)
        if ok:
            pairs.append((int(i), int(j)))
            vector = np.array(list(posi.values())) - np.array(list(posj.values()))
            angle = math.degrees(math.atan2(vector[1], vector[0]))
            if distance < min(list(sizei.values()) + list(sizej.values())):
                angle += 90
            while angle < -90:
                angle += 180
            while angle > 90:
                angle -= 180
            angles.append(angle)
    return pairs, angles

get_rotation_from_transform(transform, dims='xyz')

Source code in src\muvis_align\util.py

def get_rotation_from_transform(transform, dims='xyz'):
    # TODO: assume 2D rotation, expand to 3D
    # Rotation.from_matrix(transform).as_euler() only works for simple rotation matrices
    if isinstance(transform, DataArray):
        dims = transform['x_in'].data.tolist()
    x_index, y_index = dims.index('x'), dims.index('y')
    transform = np.array(transform)
    if y_index > x_index:
        rotation = np.arctan2(transform[0][1], transform[0][0])
    else:
        rotation = np.arctan2(transform[1][0], transform[1][1])
    return float(np.rad2deg(rotation))

get_scale_from_transform(transform)

Source code in src\muvis_align\util.py

def get_scale_from_transform(transform):
    scale = np.mean(np.linalg.norm(transform, axis=0)[:-1])
    return float(scale)

get_translation_from_transform(transform)

Source code in src\muvis_align\util.py

def get_translation_from_transform(transform):
    ndim = len(transform) - 1
    #translation = transform[:ndim, ndim]
    translation = apply_transform([[0] * ndim], transform)[0]
    return translation

get_unique_file_labels(filenames)

Source code in src\muvis_align\util.py

def get_unique_file_labels(filenames: list) -> list:
    file_labels = []

    ntot = len(filenames)
    parts_dic = get_unique_nums([split_numeric_dict(get_filetitle(filename)) for filename in filenames])
    parts_dic_full = get_unique_nums([split_numeric_dict(filename) for filename in filenames])
    parts_num_full = get_unique_nums([{index: value for index, value in enumerate(split_numeric(filename))} for filename in filenames])

    dic_ok = (len(set(['_'.join(parts.values()) for parts in parts_dic])) == ntot)
    full_dic_ok = (len(set(['_'.join(parts.values()) for parts in parts_dic_full])) == ntot)
    full_num_ok = (len(set(['_'.join(parts.values()) for parts in parts_num_full])) == ntot)

    if dic_ok:
        all_parts = parts_dic
    elif full_dic_ok:
        all_parts = parts_dic_full
    elif full_num_ok:
        all_parts = parts_num_full
    else:
        all_parts = [{0: filename} for filename in filenames]

    for parts in all_parts:
        file_label = '_'.join([key + part if isinstance(key, str) else part for key, part in parts.items()])
        file_labels.append(file_label)

    if len(set(file_labels)) < len(file_labels):
        # fallback for duplicate labels
        file_labels = [get_filetitle(filename) for filename in filenames]

    return file_labels

get_unique_nums(all_parts)

Source code in src\muvis_align\util.py

def get_unique_nums(all_parts: list) -> list:
    keys = []
    for parts in all_parts:
        for key in parts:
            if key not in keys:
                keys.append(key)

    changing_keys = []
    for key in keys:
        values = [parts.get(key) for parts in all_parts]
        if len(set(values)) > 1:
            changing_keys.append(key)

    final_parts = [{key: parts[key] for key in changing_keys if key in parts} for parts in all_parts]
    return final_parts

get_value_units_micrometer(value_units0)

Source code in src\muvis_align\util.py

def get_value_units_micrometer(value_units0: list|dict) -> list|dict|None:
    conversions = {
        'nm': 1e-3,
        'µm': 1, 'um': 1, 'micrometer': 1, 'micron': 1,
        'mm': 1e3, 'millimeter': 1e3,
        'cm': 1e4, 'centimeter': 1e4,
        'm': 1e6, 'meter': 1e6
    }
    if value_units0 is None:
        return None

    if isinstance(value_units0, dict):
        values_um = {}
        for dim, value_unit in value_units0.items():
            if isinstance(value_unit, (list, tuple)):
                value_um = value_unit[0] * conversions.get(value_unit[1], 1)
            else:
                value_um = value_unit
            values_um[dim] = value_um
    else:
        values_um = []
        for value_unit in value_units0:
            if isinstance(value_unit, (list, tuple)):
                value_um = value_unit[0] * conversions.get(value_unit[1], 1)
            else:
                value_um = value_unit
            values_um.append(value_um)
    return values_um

import_csv(filename)

Source code in src\muvis_align\util.py

def import_csv(filename):
    with open(filename, encoding='utf8') as file:
        data = csv.reader(file)
    return data

import_json(filename)

Source code in src\muvis_align\util.py

def import_json(filename):
    with open(filename, encoding='utf8') as file:
        data = json.load(file)
    return data

import_metadata(content, fields=None, input_path=None)

Source code in src\muvis_align\util.py

def import_metadata(content, fields=None, input_path=None):
    # return dict[id] = {values}
    if isinstance(content, str):
        ext = os.path.splitext(content)[1].lower()
        if input_path:
            if isinstance(input_path, list):
                input_path = input_path[0]
            content = os.path.normpath(os.path.join(os.path.dirname(input_path), content))
        if ext == '.csv':
            content = import_csv(content)
        elif ext in ['.json', '.ome.json']:
            content = import_json(content)
    if fields is not None:
        content = [[data[field] for field in fields] for data in content]
    return content

is_valid_value(value)

Source code in src\muvis_align\util.py

def is_valid_value(value):
    return value is not None and value != ''

load_sbemimage_best_config(metapath, filename)

Source code in src\muvis_align\util.py

def load_sbemimage_best_config(metapath, filename):
    target_datetime = datetime.fromtimestamp(os.path.getmtime(filename))

    for config_filename in sorted(glob.glob(os.path.join(metapath, 'logs/config_*.txt')), reverse=True):
        match = re.split(r'config_(\d+-\d+-\d+).txt', config_filename)
        if len(match) >= 2:
            file_date = datetime.strptime(match[1], '%Y-%m-%d%H%M%S%f')
            if file_date <= target_datetime:
                with open(config_filename, 'r') as file:
                    sbemimage_config = file.read()
                return sbemimage_config
    return None

metric_to_color(value)

Source code in src\muvis_align\util.py

def metric_to_color(value):
    # metric range 0...1 to traffic-light color
    if value > 0.5:
        color = 'green'
    elif value > 0.25:
        color = 'gold'
    elif value > 0.1:
        color = 'orange'
    else:
        color = 'red'
    return color

metric_to_rgb(value, min_light=0, max_light=1, range=1.0)

Source code in src\muvis_align\util.py

def metric_to_rgb(value, min_light=0, max_light=1, range=1.0):
    # metric range 0...1 to red-yellow-green ranged rgb
    colormap = colormaps.get('RdYlGn')
    index = int(value * colormap.N)
    r, g, b, a = [float(value) for value in colormap(index)]
    light = 0.2125 * r + 0.7154 * g + 0.0721 * b
    if light < min_light:
        factor = light / min_light
        r = 1 - (1 - r) * factor
        g = 1 - (1 - g) * factor
        b = 1 - (1 - b) * factor
    elif light > max_light:
        factor = max_light / light
        r *= factor
        g *= factor
        b *= factor
    r *= range
    g *= range
    b *= range
    if isinstance(range, int):
        r, g, b = int(r), int(g), int(b)
    return r, g, b

normalise_rotated_positions(centers0, rotations0, sizes, center, ndims)

Source code in src\muvis_align\util.py

def normalise_rotated_positions(centers0, rotations0, sizes, center, ndims):
    # in [xy(z)]
    centers = []
    rotations = []
    _, angles = get_pairs(centers0, sizes)
    for center0, rotation in zip(centers0, rotations0):
        if rotation is None and len(angles) > 0:
            rotation = -float(np.mean(angles))
        angle = -rotation if rotation is not None else None
        transform = create_transform(center=center, angle=angle, matrix_size=ndims + 1)
        center = apply_transform_dict([center0], transform)[0]
        centers.append(center)
        rotations.append(rotation)
    return centers, rotations

normalise_rotation(rotation)

Normalise rotation to be in the range [-180, 180].

Source code in src\muvis_align\util.py

def normalise_rotation(rotation):
    """
    Normalise rotation to be in the range [-180, 180].
    """
    while rotation < -180:
        rotation += 360
    while rotation > 180:
        rotation -= 360
    return rotation

numpy_to_native(value)

Source code in src\muvis_align\util.py

def numpy_to_native(value):
    if isinstance(value, np.ndarray):
        return value.tolist()
    elif isinstance(value, list):
        return [numpy_to_native(v) for v in value]
    elif isinstance(value, tuple):
        return (numpy_to_native(v) for v in value)
    elif isinstance(value, dict):
        return {k: numpy_to_native(v) for k, v in value.items()}
    elif hasattr(value, 'dtype'):
        return value.item()
    else:
        return value

points_to_3d(points)

Source code in src\muvis_align\util.py

def points_to_3d(points):
    return [list(point) + [0] for point in points]

print_dict(dct, indent=0)

Source code in src\muvis_align\util.py

def print_dict(dct: dict, indent: int = 0) -> str:
    s = ''
    if isinstance(dct, dict):
        for key, value in dct.items():
            s += '\n'
            if not isinstance(value, list):
                s += '\t' * indent + str(key) + ': '
            if isinstance(value, dict):
                s += print_dict(value, indent=indent + 1)
            elif isinstance(value, list):
                for v in value:
                    s += print_dict(v)
            else:
                s += str(value)
    else:
        s += str(dct)
    return s

print_dict_simple(dct)

Source code in src\muvis_align\util.py

def print_dict_simple(dct: dict) -> str:
    items = []
    for key, value in dct.items():
        if isinstance(value, float):
            value = f'{value:.3f}'
        items.append(f'{key}: {value}')
    return ' '.join(items)

print_dict_xyz(dct, dims='xyz', decimals=3, as_tuple=False)

Source code in src\muvis_align\util.py

def print_dict_xyz(dct: dict, dims='xyz', decimals=3, as_tuple=False) -> str:
    s = ''
    for dim in dims:
        if dim in dct:
            if s:
                s += ' '
            if as_tuple:
                s += f'{dct[dim]:.{decimals}f}'
            else:
                s += f'{dim}:{dct[dim]:.{decimals}f}'
    return s

print_hbytes(nbytes)

Source code in src\muvis_align\util.py

def print_hbytes(nbytes: int) -> str:
    exps = ['', 'K', 'M', 'G', 'T', 'P', 'E']
    div = 1024
    exp = 0

    while nbytes > div:
        nbytes /= div
        exp += 1
    if exp < len(exps):
        e = exps[exp]
    else:
        e = f'e{exp * 3}'
    return f'{nbytes:.1f}{e}B'

reorder(items, old_order, new_order, default_value=0)

Source code in src\muvis_align\util.py

def reorder(items: list, old_order: str, new_order: str, default_value: int = 0) -> list:
    new_items = []
    for label in new_order:
        if label in old_order:
            item = items[old_order.index(label)]
        else:
            item = default_value
        new_items.append(item)
    return new_items

retuple(chunks, shape)

Expand chunks to match shape.

E.g. if chunks is (64, 64) and shape is (3, 4, 5, 1028, 1028) return (3, 4, 5, 64, 64)

If chunks is an integer, it is applied to all dimensions, to match the behaviour of zarr-python.

Source code in src\muvis_align\util.py

def retuple(chunks, shape):
    # from ome-zarr-py
    """
    Expand chunks to match shape.

    E.g. if chunks is (64, 64) and shape is (3, 4, 5, 1028, 1028)
    return (3, 4, 5, 64, 64)

    If chunks is an integer, it is applied to all dimensions, to match
    the behaviour of zarr-python.
    """

    if isinstance(chunks, int):
        return tuple([chunks] * len(shape))

    dims_to_add = len(shape) - len(chunks)
    return *shape[:dims_to_add], *chunks

round_significants(a, significant_digits)

Source code in src\muvis_align\util.py

def round_significants(a: float, significant_digits: int) -> float:
    if a != 0:
        round_decimals = significant_digits - int(np.floor(np.log10(abs(a)))) - 1
        return round(a, round_decimals)
    return a

set_dict_value(dct, keys, value)

Source code in src\muvis_align\util.py

def set_dict_value(dct, keys, value):
    try:
        value = float(value)
    except:
        pass
    for index, key in enumerate(keys):
        if index == len(keys) - 1:
            dct[key] = value
        else:
            if key not in dct:
                dct[key] = {}
            dct = dct[key]

split_num_text(text)

Source code in src\muvis_align\util.py

def split_num_text(text: str) -> list:
    num_texts = []
    block = ''
    is_num0 = None
    if text is None:
        return []

    for c in text:
        is_num = (c.isnumeric() or c == '.')
        if is_num0 is not None and is_num != is_num0:
            num_texts.append(block)
            block = ''
        block += c
        is_num0 = is_num
    if block != '':
        num_texts.append(block)

    num_texts2 = []
    for block in num_texts:
        block = block.strip()
        try:
            block = float(block)
        except:
            pass
        if block not in [' ', ',', '|']:
            num_texts2.append(block)
    return num_texts2

split_numeric(text)

Source code in src\muvis_align\util.py

def split_numeric(text: str) -> list:
    num_parts = []
    parts = split_path_parts(text)
    for part in parts:
        num_span = re.search(r'\d+', part)
        if num_span:
            num_parts.append(part)
    return num_parts

split_numeric_dict(text)

Source code in src\muvis_align\util.py

def split_numeric_dict(text: str) -> dict:
    num_parts = {}
    parts = split_path_parts(text)
    parti = 0
    for part in parts:
        num_span = re.search(r'\d+', part)
        if num_span:
            index = num_span.start()
            label = part[:index]
            if label == '':
                label = parti
            num_parts[label] = num_span.group()
            parti += 1
    return num_parts

split_path(path)

Source code in src\muvis_align\util.py

def split_path(path: str) -> list:
    return os.path.normpath(path).split(os.path.sep)

split_path_parts(text)

Source code in src\muvis_align\util.py

def split_path_parts(text: str) -> list:
    return text.replace('/', '_').replace('\\', '_').replace('.', '_').split('_')

split_value_unit_list(text)

Source code in src\muvis_align\util.py

def split_value_unit_list(text: str) -> list:
    value_units = []
    if text is None:
        return None

    items = split_num_text(text)
    if isinstance(items[-1], str):
        def_unit = items[-1]
    else:
        def_unit = ''

    i = 0
    while i < len(items):
        value = items[i]
        if i + 1 < len(items):
            unit = items[i + 1]
        else:
            unit = ''
        if not isinstance(value, str):
            if isinstance(unit, str):
                i += 1
            else:
                unit = def_unit
            value_units.append((value, unit))
        i += 1
    return value_units

validate_transform(transform, max_scale=1.25, max_rotation=None)

Source code in src\muvis_align\util.py

def validate_transform(transform, max_scale = 1.25, max_rotation=None):
    if transform is None:
        return False
    transform = np.array(transform)
    if np.any(np.isnan(transform)):
        return False
    if np.any(np.isinf(transform)):
        return False
    if np.linalg.det(transform) == 0:
        return False
    scale = get_scale_from_transform(transform)
    if scale < 1 / max_scale or scale > max_scale:
        return False
    if  max_rotation is not None and abs(normalise_rotation(get_rotation_from_transform(transform))) > max_rotation:
        return False
    return True

xyz_to_dict(xyz, dims='xyz')

Source code in src\muvis_align\util.py

def xyz_to_dict(xyz, dims='xyz'):
    dct = {dim: float(value) for dim, value in zip(dims, xyz)}
    return dct