cylindra.widgets.sta.SubtomogramAveraging

Methods are available in the namespace ui.sta.

Widget for subtomogram averaging.

Source code in cylindra/widgets/sta.py

@magicclass(widget_type="scrollable", use_native_menubar=False)
@_shared_doc.update_cls
class SubtomogramAveraging(ChildWidget):
    """Widget for subtomogram averaging."""

    AveragingMenu = field(Averaging, name="Averaging")
    STAnalysisMenu = field(STAnalysis, name="Analysis")
    AlignmentMenu = field(Alignment, name="Alignment")
    LandscapeMenu = field(LandscapeMenu, name="Landscape")
    params = field(StaParameters)

    @property
    def sub_viewer(self) -> "napari.Viewer | None":
        """The napari viewer for subtomogram averaging."""
        return StaParameters._viewer

    def _template_param(self, *_) -> Path | None:  # for bind
        return self.params._get_template_input(allow_multiple=False)

    def _template_params(self, *_) -> Path | list[Path] | None:  # for bind
        return self.params._get_template_input(allow_multiple=True)

    def _get_mask_params(self, *_):  # for bind
        return self.params._get_mask_params()

    def _get_template_image(self) -> ip.ImgArray:
        scale = self._get_dummy_loader().scale
        template = self.params._norm_template_param(
            self._template_params(),
            allow_none=False,
            allow_multiple=True,
        ).provide(scale)
        if isinstance(template, list):
            template = ip.asarray(np.stack(template, axis=0), axes="zyx")
        else:
            template = ip.asarray(template, axes="zyx")
        return template.set_scale(zyx=scale, unit="nm")

    def _get_mask_image(self, template_params) -> ip.ImgArray:
        loader = self._get_dummy_loader()
        _, mask = loader.normalize_input(
            self.params._norm_template_param(
                template_params, allow_none=True, allow_multiple=True
            ),
            self.params._get_mask(),
        )
        if mask is None:
            raise ValueError("Mask is None.")
        return ip.asarray(mask, axes="zyx").set_scale(zyx=loader.scale, unit="nm")

    @magictoolbar
    class STATools(MagicTemplate):
        show_template = abstractapi()
        show_template_original = abstractapi()
        show_mask = abstractapi()

    @set_design(icon="ic:baseline-view-in-ar", location=STATools)
    @do_not_record
    def show_template(self):
        """Load and show template image in the scale of the tomogram."""
        template = self._get_template_image()
        self._show_rec(template, name="Template image", store=False)

    @set_design(icon="material-symbols:view-in-ar", location=STATools)
    @do_not_record
    def show_template_original(self):
        """Load and show template image in the original scale."""
        _input = self._template_params()
        if _input is None:
            raise ValueError("No template path provided.")
        elif isinstance(_input, Path):
            self._show_rec(ip.imread(_input), name="Template image", store=False)
        else:
            for i, fp in enumerate(_input):
                img = ip.imread(fp)
                self._show_rec(img, name=f"Template image [{i}]", store=False)

    @set_design(icon="fluent:shape-organic-20-filled", location=STATools)
    @do_not_record
    def show_mask(self):
        """Load and show mask image in the scale of the tomogram."""
        mask = self._get_mask_image(self._template_params())
        self._show_rec(mask, name="Mask image", store=False, threshold=0.5)

    @property
    def last_average(self) -> "ip.ImgArray | None":
        """Last averaged image if exists."""
        return StaParameters._last_average

    def _get_shape_in_nm(self, default: int | None = None) -> tuple[nm, nm, nm]:
        if default is None:
            tmp = self._get_template_image()
            return tuple(np.array(tmp.sizesof("zyx")) * tmp.scale.x)
        else:
            return (default,) * 3

    @thread_worker.callback
    def _show_rec(self, img: ip.ImgArray, name: str, store=True, threshold=None):
        return self.params._show_reconstruction(img, name, store, threshold=threshold)

    def _get_loader(
        self,
        binsize: int,
        molecules: Molecules,
        shape: tuple[nm, nm, nm] = None,
        order: int = 1,
    ) -> SubtomogramLoader:
        """
        Returns proper subtomogram loader, template image and mask image that matche the
        bin size.
        """
        return self._get_main().tomogram.get_subtomogram_loader(
            molecules,
            binsize=binsize,
            order=order,
            output_shape=shape,
        )

    def _get_dummy_loader(self):
        return self._get_loader(binsize=1, molecules=Molecules.empty())

    def _get_available_binsize(self, _=None) -> list[int]:
        parent = self._get_main()
        out = [x[0] for x in parent.tomogram.multiscaled]
        if 1 not in out:
            out = [1] + out
        return out

    @set_design(text="Average all molecules", location=Averaging)
    @dask_worker.with_progress(desc=_pdesc.fmt_layers("Subtomogram averaging of {!r}"))
    def average_all(
        self,
        layers: MoleculesLayersType,
        size: _SubVolumeSize = None,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 1,
        bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
    ):
        """
        Subtomogram averaging using all the molecules in the selected layer(s).

        If multiple layers are selected, subtomograms around all the molecules will
        be averaged.

        Parameters
        ----------
        {layers}{size}{interpolation}{bin_size}
        """
        t0 = timer()
        layers = assert_list_of_layers(layers, self.parent_viewer)
        parent = self._get_main()
        tomo = parent.tomogram
        shape = self._get_shape_in_nm(size)
        loader = tomo.get_subtomogram_loader(
            _concat_molecules(layers), shape, binsize=bin_size, order=interpolation
        )
        img = ip.asarray(loader.average(), axes="zyx")
        img.set_scale(zyx=loader.scale, unit="nm")
        t0.toc()
        return self._show_rec.with_args(img, f"[AVG]{_avg_name(layers)}")

    @set_design(text="Average subset of molecules", location=Averaging)
    @dask_worker.with_progress(desc=_pdesc.fmt_layers("Subtomogram averaging (subset) of {!r}"))  # fmt: skip
    def average_subset(
        self,
        layers: MoleculesLayersType,
        size: _SubVolumeSize = None,
        method: Literal["steps", "first", "last", "random"] = "steps",
        number: int = 64,
        bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
    ):
        """
        Subtomogram averaging using a subset of subvolumes.

        If multiple layers are selected, subtomograms around all the molecules will
        be concatenated before choosing a subset.

        Parameters
        ----------
        {layers}{size}
        method : str, optional
            How to choose subtomogram subset.
            (1) steps: Each 'steps' subtomograms from the tip of spline.
            (2) first: First subtomograms.
            (3) last: Last subtomograms.
            (4) random: choose randomly.
        number : int, default
            Number of subtomograms to use.
        {bin_size}
        """
        t0 = timer()
        layers = assert_list_of_layers(layers, self.parent_viewer)
        parent = self._get_main()
        molecules = _concat_molecules(layers)
        nmole = len(molecules)
        shape = self._get_shape_in_nm(size)
        sl = _get_slice_for_average_subset(method, nmole, number)
        mole = molecules.subset(sl)
        loader = parent.tomogram.get_subtomogram_loader(
            mole, shape, binsize=bin_size, order=1
        )
        img = ip.asarray(loader.average(), axes="zyx").set_scale(zyx=loader.scale)
        t0.toc()
        return self._show_rec.with_args(img, f"[AVG(n={number})]{_avg_name(layers)}")

    @set_design(text="Average group-wise", location=Averaging)
    @dask_worker.with_progress(desc=_pdesc.fmt_layers("Grouped subtomogram averaging of {!r}"))  # fmt: skip
    def average_groups(
        self,
        layers: MoleculesLayersType,
        size: _SubVolumeSize = None,
        by: PolarsExprStr = "col('pf-id')",
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 1,
        bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
    ):
        """
        Group-wise subtomogram averaging using molecules grouped by the given expression.

        This method first group molecules by its features, and then average each group.
        This method is useful for such as get average of each protofilament and segmented
        subtomogram averaging.

        Parameters
        ----------
        {layers}{size}
        by : str or polars expression
            Expression to group molecules.
        {interpolation}{bin_size}
        """
        t0 = timer()
        layers = assert_list_of_layers(layers, self.parent_viewer)
        parent = self._get_main()
        tomo = parent.tomogram
        shape = self._get_shape_in_nm(size)
        loader = tomo.get_subtomogram_loader(
            _concat_molecules(layers), shape, binsize=bin_size, order=interpolation
        )
        expr = widget_utils.norm_expr(by)
        avg_dict = loader.groupby(expr).average()
        avgs = np.stack([avg_dict[k] for k in sorted(avg_dict.keys())], axis=0)
        img = ip.asarray(avgs, axes="pzyx")
        img.set_scale(zyx=loader.scale, unit="nm")
        t0.toc()
        return self._show_rec.with_args(img, f"[AVG]{_avg_name(layers)}", store=False)

    @set_design(text="Average filtered", location=Averaging)
    @dask_worker.with_progress(desc=_pdesc.fmt_layers("Filtered subtomogram averaging of {!r}"))  # fmt: skip
    def average_filtered(
        self,
        layers: MoleculesLayersType,
        size: _SubVolumeSize = None,
        predicate: PolarsExprStr = "col('pf-id') == 0",
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 1,
        bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
    ):
        """
        Subtomogram averaging using molecules filtered by the given expression.

        This method first concatenate molecules in the selected layers, and then filter them
        by the predicate.

        Parameters
        ----------
        {layers}{size}
        predicate : str or polars expression
            Filter expression to select molecules.
        {interpolation}{bin_size}
        """
        t0 = timer()
        layers = assert_list_of_layers(layers, self.parent_viewer)
        parent = self._get_main()
        tomo = parent.tomogram
        shape = self._get_shape_in_nm(size)
        loader = tomo.get_subtomogram_loader(
            _concat_molecules(layers), shape, binsize=bin_size, order=interpolation
        )
        avg = loader.filter(widget_utils.norm_expr(predicate)).average()
        img = ip.asarray(avg, axes="zyx")
        img.set_scale(zyx=loader.scale, unit="nm")
        t0.toc()
        return self._show_rec.with_args(img, f"[AVG]{_avg_name(layers)}")

    @set_design(text="Split and average molecules", location=Averaging)
    @dask_worker.with_progress(desc=_pdesc.fmt_layers("Split-and-averaging of {!r}"))  # fmt: skip
    def split_and_average(
        self,
        layers: MoleculesLayersType,
        n_pairs: Annotated[int, {"min": 1, "label": "number of image pairs"}] = 1,
        size: _SubVolumeSize = None,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 1,
        bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
    ):
        """
        Split molecules into two groups and average separately.

        Parameters
        ----------
        {layers}
        n_pairs : int, default 1
            How many pairs of average will be calculated.
        {size}{interpolation}{bin_size}
        """
        t0 = timer()
        layers = assert_list_of_layers(layers, self.parent_viewer)
        parent = self._get_main()
        molecules = _concat_molecules(layers)
        shape = self._get_shape_in_nm(size)
        loader = parent.tomogram.get_subtomogram_loader(
            molecules, shape, binsize=bin_size, order=interpolation
        )
        axes = "ipzyx" if n_pairs > 1 else "pzyx"
        img = ip.asarray(loader.average_split(n_set=n_pairs), axes=axes)
        img.set_scale(zyx=loader.scale)
        t0.toc()
        return self._show_rec.with_args(img, f"[Split]{_avg_name(layers)}", store=False)

    @set_design(text="Align average to template", location=Alignment)
    @dask_worker.with_progress(descs=_pdesc.align_averaged_fmt)
    def align_averaged(
        self,
        layers: MoleculesLayersType,
        template_path: Annotated[_PathOrNone, {"bind": _template_param}],
        mask_params: Annotated[Any, {"bind": _get_mask_params}],
        max_shifts: Optional[_MaxShifts] = None,
        rotations: _Rotations = ((0.0, 0.0), (15.0, 1.0), (3.0, 1.0)),
        bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
        method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
    ):  # fmt: skip
        """
        Align the averaged image at current monomers to the template image.

        This function creates a new layer with transformed monomers, which should
        align well with template image.

        Parameters
        ----------
        {layers}{template_path}{mask_params}{max_shifts}{rotations}{bin_size}{method}
        """
        t0 = timer()
        layers = assert_list_of_layers(layers, self.parent_viewer)
        parent = self._get_main()

        new_layers = list[MoleculesLayer]()

        @thread_worker.callback
        def _on_yield(mole_trans: Molecules, layer: MoleculesLayer):
            points = parent.add_molecules(
                mole_trans,
                name=_coerce_aligned_name(layer.name, self.parent_viewer),
                source=layer.source_component,
            )
            new_layers.append(points)
            layer.visible = False
            _Logger.print_html(f"{layer.name!r} → {points.name!r}")

        mole = layers[0].molecules
        loader = self._get_loader(bin_size, mole, order=1)
        template, mask = loader.normalize_input(
            template=self.params._norm_template_param(
                template_path, allow_multiple=False
            ),
            mask=self.params._get_mask(params=mask_params),
        )
        temp_norm = utils.normalize_image(template)

        _scale = parent.tomogram.scale * bin_size

        if max_shifts is None:
            max_shifts = _default_align_averaged_shifts(mole)

        model = _get_alignment(method)(
            template,
            mask,
            rotations=rotations,
            tilt=None,  # NOTE: because input is an average
        )
        _spl_globs = list[
            tuple[weakref.ReferenceType["CylSpline"], pl.DataFrame, pl.DataFrame]
        ]()
        for layer in layers:
            mole = layer.molecules
            loader = self._get_loader(bin_size, mole, order=1)
            _img_trans, result = model.fit(
                loader.average(template.shape),
                max_shifts=[_s / _scale for _s in max_shifts],
            )

            rotator = Rotation.from_quat(result.quat)
            svec = result.shift * _scale
            _mole_trans = mole.linear_transform(result.shift * _scale, rotator)

            # write offsets to spline globalprops if available
            if spl := layer.source_spline:
                _mole_trans = _update_mole_pos(_mole_trans, mole, spl)
                if spl.radius is None:
                    _radius: nm = cylmeasure.calc_radius(mole, spl).mean()
                else:
                    _radius = spl.radius
                _glob_old = spl.props.glob.clone()
                _glob_new = _update_offset(spl, rotator.apply(svec), _radius)
                spl.props.glob = _glob_new
                _spl_globs.append((weakref.ref(spl), _glob_old, _glob_new))

            yield _on_yield.with_args(_mole_trans, layer)

            # create images for visualization in the logger. Image is magenta, template is green
            img_norm = utils.normalize_image(_img_trans)
            merge = np.stack([img_norm, temp_norm, img_norm], axis=-1)
            with _Logger.set_plt():
                widget_utils.plot_projections(merge)

            # logging
            rvec = rotator.as_rotvec()
            _fmt = "  {:.2f}  ".format
            _Logger.print_table(
                [
                    ["", "X", "Y", "Z"],
                    ["Shift (nm)", _fmt(svec[2]), _fmt(svec[1]), _fmt(svec[0])],
                    ["Rot vector", _fmt(rvec[2]), _fmt(rvec[1]), _fmt(rvec[0])],
                ],
                header=False,
                index=False,
            )

        t0.toc()

        @thread_worker.callback
        def _align_averaged_on_return():
            @undo_callback
            def _out():
                parent._try_removing_layers(new_layers)
                for spl_ref, old, _ in _spl_globs:
                    if spl := spl_ref():
                        spl.props.glob = old

            @_out.with_redo
            def _out():
                parent._add_layers_future(new_layers)()
                for spl_ref, _, new in _spl_globs:
                    if spl := spl_ref():
                        spl.props.glob = new

            return _out

        return _align_averaged_on_return

    sep0 = Separator

    @set_design(text="Align all molecules", location=Alignment)
    @dask_worker.with_progress(descs=_pdesc.align_all_fmt)
    def align_all(
        self,
        layers: MoleculesLayersType,
        template_path: Annotated[_PathOrPathsOrNone, {"bind": _template_params}],
        mask_params: Annotated[Any, {"bind": _get_mask_params}],
        max_shifts: _MaxShifts = (1.0, 1.0, 1.0),
        rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
        cutoff: _CutoffFreq = 0.5,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
        method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
        bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
    ):  # fmt: skip
        """
        Align the input template image to all the molecules.

        Parameters
        ----------
        {layers}{template_path}{mask_params}{max_shifts}{rotations}{cutoff}
        {interpolation}{method}{bin_size}
        """
        t0 = timer()
        layers = assert_list_of_layers(layers, self.parent_viewer)
        main = self._get_main()

        combiner = MoleculesCombiner()

        loader = self._get_loader(
            binsize=bin_size,
            molecules=combiner.concat(layer.molecules for layer in layers),
            order=interpolation,
        )
        aligned_loader = loader.align(
            template=self.params._norm_template_param(
                template_path, allow_multiple=True
            ),
            mask=self.params._get_mask(params=mask_params),
            max_shifts=max_shifts,
            rotations=rotations,
            cutoff=cutoff,
            alignment_model=_get_alignment(method),
            tilt=main.tomogram.tilt_model,
        )
        molecules = combiner.split(aligned_loader.molecules, layers)
        t0.toc()
        return self._align_all_on_return.with_args(molecules, layers)

    @set_design(text="Align all (template-free)", location=Alignment)
    @dask_worker.with_progress(descs=_pdesc.align_template_free_fmt)
    def align_all_template_free(
        self,
        layers: MoleculesLayersType,
        mask_params: Annotated[Any, {"bind": _get_mask_params}],
        size: _SubVolumeSize = 12.0,
        max_shifts: _MaxShifts = (1.0, 1.0, 1.0),
        rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
        cutoff: _CutoffFreq = 0.5,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
        method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
        bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
    ):  # fmt: skip
        """
        Run template-free alignment for the given layers (EXPERIMENTAL).

        Parameters
        ----------
        {layers}{mask_params}{size}{max_shifts}{rotations}{cutoff}{interpolation}
        {method}{bin_size}
        """
        t0 = timer()
        layers = assert_list_of_layers(layers, self.parent_viewer)
        main = self._get_main()
        combiner = MoleculesCombiner()
        molecules = combiner.concat(layer.molecules for layer in layers)
        mask = self.params._get_mask(params=mask_params)
        if size is None:
            raise NotImplementedError("'size' must be given.")
        else:
            shape = tuple(
                main.tomogram.nm2pixel(self._get_shape_in_nm(size), binsize=bin_size)
            )

        aligned_loader = (
            self._get_loader(binsize=bin_size, molecules=molecules, order=interpolation)
            .reshape(shape=shape)
            .align_no_template(
                mask=mask,
                max_shifts=max_shifts,
                rotations=rotations,
                cutoff=cutoff,
                alignment_model=_get_alignment(method),
                tilt=main.tomogram.tilt_model,
            )
        )
        molecules = combiner.split(aligned_loader.molecules, layers)
        t0.toc()
        return self._align_all_on_return.with_args(molecules, layers)

    sep1 = Separator

    @set_design(text="Viterbi Alignment", location=Alignment)
    @dask_worker.with_progress(descs=_pdesc.align_viterbi_fmt)
    def align_all_viterbi(
        self,
        layer: MoleculesLayerType,
        template_path: Annotated[_PathOrPathsOrNone, {"bind": _template_params}],
        mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
        max_shifts: _MaxShifts = (0.8, 0.8, 0.8),
        rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
        cutoff: _CutoffFreq = 0.5,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
        range_long: _DistRangeLon = (4.0, 4.28),
        angle_max: _AngleMaxLon = 5.0,
        bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
        upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 5,
    ):  # fmt: skip
        """
        Subtomogram alignment using 1D Viterbi alignment.

        1D Viterbi alignment is an alignment algorithm that considers the distance and
        the skew angle between every longitudinally adjacent monomers. The classical
        Viterbi algorithm is used to find the global optimal solution of the alignment.
        Note that Viterbi alignment is data size dependent, i.e. the alignment result
        of a molecule may vary depending on the total number of molecules in the dataset.

        Parameters
        ----------
        {layer}{template_path}{mask_params}{max_shifts}{rotations}{cutoff}
        {interpolation}{range_long}{angle_max}{bin_size}{upsample_factor}
        """
        t0 = timer()
        layer = assert_layer(layer, self.parent_viewer)
        landscape = self._construct_landscape(
            molecules=layer.molecules,
            template_path=template_path,
            mask_params=mask_params,
            max_shifts=max_shifts,
            rotations=rotations,
            cutoff=cutoff,
            order=interpolation,
            upsample_factor=upsample_factor,
            bin_size=bin_size,
        )

        yield
        mole = landscape.run_viterbi_along_spline(
            spl=layer.source_spline,
            range_long=range_long,
            angle_max=angle_max,
        )
        t0.toc()
        return self._align_all_on_return.with_args([mole], [layer])

    @set_design(text="Simulated annealing", location=Alignment)
    @dask_worker.with_progress(descs=_pdesc.align_annealing_fmt)
    def align_all_annealing(
        self,
        layer: MoleculesLayerType,
        template_path: Annotated[_PathOrPathsOrNone, {"bind": _template_params}],
        mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
        max_shifts: _MaxShifts = (0.8, 0.8, 0.8),
        rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
        cutoff: _CutoffFreq = 0.5,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
        range_long: _DistRangeLon = (4.0, 4.28),
        range_lat: _DistRangeLat = (5.1, 5.3),
        angle_max: _AngleMaxLon = 5.0,
        bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
        temperature_time_const: Annotated[float, {"min": 0.01, "max": 10.0}] = 1.0,
        upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 5,
        random_seeds: _RandomSeeds = (0, 1, 2, 3, 4),
    ):  # fmt: skip
        """
        2D-constrained subtomogram alignment using simulated annealing.

        This alignment method considers the distance between every adjacent monomers.
        Two-dimensionally connected optimization can be approximated by the simulated
        annealing algorithm.

        Parameters
        ----------
        {layer}{template_path}{mask_params}{max_shifts}{rotations}{cutoff}
        {interpolation}{range_long}{range_lat}{angle_max}{bin_size}
        {temperature_time_const}{upsample_factor}{random_seeds}
        """
        t0 = timer()
        layer = assert_layer(layer, self.parent_viewer)
        if layer.source_spline is None:
            raise ValueError("RMA requires a spline.")
        main = self._get_main()
        landscape = self._construct_landscape(
            molecules=layer.molecules,
            template_path=template_path,
            mask_params=mask_params,
            max_shifts=max_shifts,
            rotations=rotations,
            cutoff=cutoff,
            order=interpolation,
            bin_size=bin_size,
            upsample_factor=upsample_factor,
        )
        yield
        mole, results = landscape.run_annealing_along_spline(
            layer.source_spline,
            range_long=range_long,
            range_lat=range_lat,
            angle_max=angle_max,
            temperature_time_const=temperature_time_const,
            random_seeds=random_seeds,
        )
        t0.toc()

        @thread_worker.callback
        def _on_return():
            points = main.add_molecules(
                mole,
                name=_coerce_aligned_name(layer.name, self.parent_viewer),
                source=layer.source_component,
                metadata={ANNEALING_RESULT: results[0]},
            )
            layer.visible = False
            with _Logger.set_plt():
                _annealing.plot_annealing_result(results)
            return self._undo_for_new_layer([layer.name], [points])

        return _on_return

    @set_design(text=capitalize, location=LandscapeMenu)
    @dask_worker.with_progress(descs=_pdesc.construct_landscape_fmt)
    def construct_landscape(
        self,
        layer: MoleculesLayerType,
        template_path: Annotated[_PathOrPathsOrNone, {"bind": _template_params}],
        mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
        max_shifts: _MaxShifts = (0.8, 0.8, 0.8),
        rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
        cutoff: _CutoffFreq = 0.5,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
        bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
        upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 5,
        method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
        norm: bool = True,
    ):
        """
        Construct a cross-correlation landscape for subtomogram alignment.

        Parameters
        ----------
        {layer}{template_path}{mask_params}{max_shifts}{rotations}{cutoff}
        {interpolation}{bin_size}{upsample_factor}{method}
        norm: bool, default True
            If true, each landscape will be normalized by its mean and standard
            deviation.
        """
        layer = assert_layer(layer, self.parent_viewer)
        lnd = self._construct_landscape(
            molecules=layer.molecules,
            template_path=template_path,
            mask_params=mask_params,
            max_shifts=max_shifts,
            rotations=rotations,
            cutoff=cutoff,
            order=interpolation,
            bin_size=bin_size,
            upsample_factor=upsample_factor,
            norm=norm,
            method=method,
        )
        surf = LandscapeSurface(lnd, name=f"{LANDSCAPE_PREFIX}{layer.name}")
        surf.source_component = layer.source_component

        @thread_worker.callback
        def _on_return():
            self.parent_viewer.add_layer(surf)
            self._get_main()._reserved_layers.to_be_removed.add(surf)
            layer.visible = False

        return _on_return

    @set_design(text="Run alignment on landscape", location=LandscapeMenu)
    @dask_worker.with_progress(desc="Peak detection on landscape")
    def run_align_on_landscape(self, landscape_layer: _LandscapeLayer):
        """Find the optimal displacement for each molecule on the landscape."""
        landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
        landscape = landscape_layer.landscape
        spl = landscape_layer.source_spline
        mole_opt, _ = landscape.run_min_energy(spl)
        return self._align_on_landscape_on_return.with_args(
            mole_opt, landscape_layer.name, spl
        )

    @set_design(text="Run Viterbi alignment on landscape", location=LandscapeMenu)
    @dask_worker.with_progress(desc="Running Viterbi alignment")
    def run_viterbi_on_landscape(
        self,
        landscape_layer: _LandscapeLayer,
        range_long: _DistRangeLon = (4.0, 4.28),
        angle_max: _AngleMaxLon = 5.0,
    ):
        """
        Run Viterbi alignment on the landscape.

        Parameters
        ----------
        {landscape_layer}{range_long}{angle_max}
        """
        t0 = timer()
        landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
        spl = landscape_layer.source_spline
        mole = landscape_layer.landscape.run_viterbi_along_spline(
            spl=spl,
            range_long=range_long,
            angle_max=angle_max,
        )
        t0.toc()
        return self._align_on_landscape_on_return.with_args(
            mole, landscape_layer.name, spl
        )

    @set_design(text="Run annealing on landscape", location=LandscapeMenu)
    @dask_worker.with_progress(desc="Running simulated annealing")
    def run_annealing_on_landscape(
        self,
        landscape_layer: _LandscapeLayer,
        range_long: _DistRangeLon = (4.0, 4.28),
        range_lat: _DistRangeLat = (5.1, 5.3),
        angle_max: _AngleMaxLon = 5.0,
        temperature_time_const: Annotated[float, {"min": 0.01, "max": 10.0}] = 1.0,
        random_seeds: _RandomSeeds = (0, 1, 2, 3, 4),
    ):
        """
        Run simulated annealing on the landscape.

        Parameters
        ----------
        {landscape_layer}{range_long}{range_lat}{angle_max}{temperature_time_const}
        {random_seeds}
        """
        t0 = timer()
        landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
        spl = landscape_layer.source_spline
        if spl is None:
            raise ValueError("RMA requires a spline.")
        mole, results = landscape_layer.landscape.run_annealing_along_spline(
            spl=spl,
            range_long=range_long,
            range_lat=range_lat,
            angle_max=angle_max,
            temperature_time_const=temperature_time_const,
            random_seeds=random_seeds,
        )
        t0.toc()

        @thread_worker.callback
        def _plot_result():
            with _Logger.set_plt():
                _annealing.plot_annealing_result(results)

        yield _plot_result
        return self._align_on_landscape_on_return.with_args(
            mole,
            landscape_layer.name,
            source=spl,
            metadata={ANNEALING_RESULT: results[0]},
        )

    @set_design(text=capitalize, location=LandscapeMenu)
    def remove_landscape_outliers(
        self,
        landscape_layer: _LandscapeLayer,
        lower: Annotated[Optional[float], {"text": "Do not process lower outliers"}] = None,
        upper: Annotated[Optional[float], {"text": "Do not process upper outliers"}] = None,
    ):  # fmt: skip
        """
        Remove outliers from the landscape.

        This method will replace energy (inverse score) outliers with the thresholds.
        This method is useful for lattice with such as defects or strong artifacts.

        Parameters
        ----------
        {landscape_layer}
        lower : float, optional
            Lower limit of the energy.
        upper : float, optional
            Upper limit of the energy.
        """
        landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
        new = landscape_layer.landscape.clip_energies(lower, upper)
        surf = LandscapeSurface(new, name=f"{landscape_layer}-Clip")
        return self._add_new_landscape_layer(landscape_layer, surf)

    @set_design(text=capitalize, location=LandscapeMenu)
    def normalize_landscape(
        self,
        landscape_layer: _LandscapeLayer,
        norm_sd: bool = True,
    ):
        """
        Normalize the landscape.

        Parameters
        ----------
        {landscape_layer}
        norm_sd : bool, default True
            If true, each landscape will also be normalized by its standard deviation.
        """
        landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
        new = landscape_layer.landscape.normed(sd=norm_sd)
        surf = LandscapeSurface(new, name=f"{landscape_layer}-Norm")
        return self._add_new_landscape_layer(landscape_layer, surf)

    def _add_new_landscape_layer(self, old: LandscapeSurface, new: LandscapeSurface):
        new.source_component = old.source_component

        self.parent_viewer.add_layer(new)
        self._get_main()._reserved_layers.to_be_removed.add(new)
        old.visible = False
        return None

    def _get_layers_with_annealing_result(self, *_) -> list[MoleculesLayer]:
        if self.parent_viewer is None:
            return []
        return [
            (layer.name, layer)
            for layer in self.parent_viewer.layers
            if ANNEALING_RESULT in layer.metadata
        ]

    @set_design(text=capitalize, location=Alignment)
    @do_not_record
    def save_annealing_scores(
        self,
        layer: Annotated[MoleculesLayer, {"choices": _get_layers_with_annealing_result}],
        path: Path.Save[FileFilter.CSV],
    ):  # fmt: skip
        """Save RMA scores to a CSV file."""
        layer = assert_layer(layer, self.parent_viewer)
        try:
            result: AnnealingResult = layer.metadata[ANNEALING_RESULT]
        except KeyError:
            raise ValueError(
                f"Layer {layer!r} does not have annealing result."
            ) from None
        x = result.batch_size * np.arange(result.energies.size)
        df = pl.DataFrame({"iteration": x, "score": -result.energies})
        return df.write_csv(path, include_header=False)

    @set_design(text=capitalize, location=STAnalysis)
    @dask_worker.with_progress(desc=_pdesc.fmt_layers("Calculating correlations of {!r}"))  # fmt: skip
    def calculate_correlation(
        self,
        layers: MoleculesLayersType,
        template_path: Annotated[_PathOrPathsOrNone, {"bind": _template_params}],
        mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
        cutoff: _CutoffFreq = 0.5,
        bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
        method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
        column_prefix: str = "score",
    ):
        """
        Calculate correlation between template images and the subtomograms.

        This method will load every subtomograms, calculate the correlation between
        the template images and each subtomogram, and save the correlation values
        as new columns in the molecules features.

        Parameters
        ----------
        {layers}{template_path}{mask_params}{interpolation}{cutoff}{bin_size}{method}
        column_prefix : str, default "score"
            Prefix of the column names of the calculated correlations.
        """
        layers = assert_list_of_layers(layers, self.parent_viewer)
        main = self._get_main()
        combiner = MoleculesCombiner()

        if isinstance(template_path, (Path, str)):
            template_path = [template_path]
        mask = self.params._get_mask(mask_params)
        all_mole = combiner.concat(layer.molecules for layer in layers)

        out = main.tomogram.get_subtomogram_loader(
            all_mole,
            order=interpolation,
            binsize=bin_size,
        ).score(
            templates=[pipe.from_file(t) for t in template_path],
            mask=mask,
            alignment_model=_get_alignment(method),
            cutoff=cutoff,
            tilt=main.tomogram.tilt_model,
        )
        all_mole = all_mole.with_features(
            pl.Series(f"{column_prefix}_{i}", col) for i, col in enumerate(out)
        )

        @thread_worker.callback
        def _on_return():
            moles = combiner.split(all_mole, layers)
            for layer, each_mole in zip(layers, moles, strict=True):
                features = each_mole.features.select(
                    [f"{column_prefix}_{i}" for i in range(len(out))]
                )
                layer.set_molecules_with_new_features(
                    layer.molecules.with_features(features)
                )

        return _on_return

    @set_design(text="Calculate FSC", location=STAnalysis)
    @dask_worker.with_progress(desc=_pdesc.fmt_layers("Calculating FSC of {!r}"))
    def calculate_fsc(
        self,
        layers: MoleculesLayersType,
        template_path: Annotated[_PathOrNone, {"bind": _template_param}] = None,
        mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
        size: _SubVolumeSize = None,
        seed: Annotated[Optional[int], {"text": "Do not use random seed."}] = 0,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 1,
        n_pairs: Annotated[int, {"min": 1, "label": "number of image pairs"}] = 1,
        show_average: bool = True,
        dfreq: FSCFreq = None,
    ):
        """
        Calculate Fourier Shell Correlation using the selected monomer layer.

        Parameters
        ----------
        {layers}
        template_path : template input type
            Used only when soft-Otsu mask parameters are given.
        {mask_params}{size}
        seed : int, optional
            Random seed used for subtomogram sampling.
        {interpolation}
        n_pairs : int, default 1
            How many sets of image pairs will be generated to average FSC.
        show_average : bool, default True
            If true, subtomogram average will be shown after FSC calculation.
        dfreq : float, default 0.02
            Precision of frequency to calculate FSC. "0.02" means that FSC will be
            calculated at frequency 0.01, 0.03, 0.05, ..., 0.45.
        """
        t0 = timer()
        layers = assert_list_of_layers(layers, self.parent_viewer)
        main = self._get_main()
        mole = _concat_molecules(layers)

        loader = main.tomogram.get_subtomogram_loader(mole, order=interpolation)
        template, mask = loader.normalize_input(
            template=self.params._norm_template_param(template_path, allow_none=True),
            mask=self.params._get_mask(params=mask_params),
        )
        fsc, (img_0, img_1), img_mask = loader.reshape(
            template=template if size is None else None,
            mask=mask,
            shape=None if size is None else (main.tomogram.nm2pixel(size),) * 3,
        ).fsc_with_halfmaps(mask, seed=seed, n_set=n_pairs, dfreq=dfreq, squeeze=False)

        def _as_imgarray(im, axes: str = "zyx") -> ip.ImgArray | None:
            if np.isscalar(im):
                return None
            return ip.asarray(im, axes=axes).set_scale(zyx=loader.scale, unit="nm")

        if show_average:
            avg = (img_0[0] + img_1[0]) / 2
            img_avg = _as_imgarray(avg)
        else:
            img_avg = None

        result = widget_utils.FscResult.from_dataframe(fsc, loader.scale)
        criteria = [0.5, 0.143]
        _name = _avg_name(layers)
        t0.toc()

        @thread_worker.callback
        def _calculate_fsc_on_return():
            _Logger.print_html(f"<b>Fourier Shell Correlation of {_name!r}</b>")
            with _Logger.set_plt():
                result.plot(criteria)
            for _c in criteria:
                _r = result.get_resolution(_c)
                _Logger.print_html(f"Resolution at FSC={_c:.3f} ... <b>{_r:.3f} nm</b>")

            if img_avg is not None:
                _imlayer: "Image" = self._show_rec(img_avg, name=f"[AVG]{_name}")
                _imlayer.metadata["fsc"] = result
                _imlayer.metadata["fsc_halfmaps"] = (
                    _as_imgarray(img_0, axes="izyx"),
                    _as_imgarray(img_1, axes="izyx"),
                )
                _imlayer.metadata["fsc_mask"] = _as_imgarray(img_mask)

        return _calculate_fsc_on_return

    @set_design(text="PCA/K-means classification", location=STAnalysis)
    @dask_worker.with_progress(descs=_pdesc.classify_pca_fmt)
    def classify_pca(
        self,
        layer: MoleculesLayerType,
        template_path: Annotated[_PathOrNone, {"bind": _template_param}] = None,
        mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
        size: _SubVolumeSize = None,
        cutoff: _CutoffFreq = 0.5,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
        bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
        n_components: Annotated[int, {"min": 2, "max": 20}] = 2,
        n_clusters: Annotated[int, {"min": 2, "max": 100}] = 2,
        seed: Annotated[Optional[int], {"text": "Do not use random seed."}] = 0,
    ):  # fmt: skip
        """
        Classify molecules in a layer using PCA and K-means clustering.

        Parameters
        ----------
        {layer}
        template_path : template input type
            Used only when soft-Otsu mask parameters are given.
        {mask_params}{size}{cutoff}{interpolation}{bin_size}
        n_components : int, default 2
            The number of PCA dimensions.
        n_clusters : int, default
            The number of clusters.
        seed : int, default 0
            Random seed.
        """
        from cylindra.widgets.subwidgets import PcaViewer

        t0 = timer()
        layer = assert_layer(layer, self.parent_viewer)
        tomo = self._get_main().tomogram
        loader = self._get_loader(
            binsize=bin_size, molecules=layer.molecules, order=interpolation
        )
        template, mask = loader.normalize_input(
            template=self.params._norm_template_param(template_path, allow_none=True),
            mask=self.params._get_mask(params=mask_params),
        )
        shape = None
        if size is not None and mask is None:
            shape = (tomo.nm2pixel(size, binsize=bin_size),) * 3
        out, pca = loader.reshape(
            template=template if mask is None and shape is None else None,
            mask=mask,
            shape=shape,
        ).classify(
            mask=mask,
            seed=seed,
            cutoff=cutoff,
            n_components=n_components,
            n_clusters=n_clusters,
            label_name="cluster",
        )

        avgs_dict = out.groupby("cluster").average()
        avgs = ip.asarray(
            np.stack(list(avgs_dict.values()), axis=0), axes=["cluster", "z", "y", "x"]
        ).set_scale(zyx=loader.scale, unit="nm")
        t0.toc()

        @thread_worker.callback
        def _on_return():
            layer.molecules = out.molecules  # update features
            pca_viewer = PcaViewer(pca)
            pca_viewer.native.setParent(self.native, pca_viewer.native.windowFlags())
            pca_viewer.show()
            self._show_rec(avgs, name=f"[PCA]{layer.name}", store=False)
            ACTIVE_WIDGETS.add(pca_viewer)

        return _on_return

    @set_design(text="Seam search by correlation", location=STAnalysis.SeamSearch)
    @dask_worker.with_progress(desc=_pdesc.fmt_layer("Seam search of {!r}"))
    def seam_search(
        self,
        layer: MoleculesLayerType,
        template_path: Annotated[_PathOrNone, {"bind": _template_param}],
        mask_params: Annotated[Any, {"bind": _get_mask_params}],
        anti_template_path: Annotated[Optional[Path.Read[FileFilter.IMAGE]], {"text": "Do not use anti-template", "label": "anti-template path"}] = None,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
        npf: Annotated[Optional[int], {"text": "use global properties"}] = None,
        show_average: Annotated[str, {"label": "show averages as", "choices": [None, "Raw", "Filtered"]}] = "Filtered",
        cutoff: _CutoffFreq = 0.25,
    ):  # fmt: skip
        """
        Search for the best seam position.

        Try all patterns of seam positions and compare cross correlation values. If
        molecule assembly has 13 protofilaments, this method will try 26 patterns.

        Parameters
        ----------
        {layer}{template_path}{mask_params}
        anti_template_path : Path, optional
            The anti-template used for seam search. For instance, if the template is
            beta-tubulin, the anti-template is alpha-tubulin.
        {interpolation}
        npf : int, optional
            Number of protofilaments. By default the global properties stored in the
            corresponding spline will be used.
        show_average : bool, default True
            If true, all the subtomogram averages will be shown.
        {cutoff}
        """
        t0 = timer()
        layer = assert_layer(layer, self.parent_viewer)
        loader, npf = self._seam_search_input(layer, npf, interpolation)
        template, mask = loader.normalize_input(
            template=self.params._norm_template_param(template_path),
            mask=self.params._get_mask(params=mask_params),
        )
        if anti_template_path is not None:
            anti_template = ip.asarray(
                pipe.from_file(anti_template_path).provide(loader.scale),
                axes="zyx",
            )
        else:
            anti_template = None

        seam_searcher = CorrelationSeamSearcher(npf)
        result = seam_searcher.search(
            loader=loader,
            template=ip.asarray(template, axes="zyx"),
            anti_template=anti_template,
            mask=mask,
            cutoff=cutoff,
        )

        t0.toc()

        @thread_worker.callback
        def _seam_search_on_return():
            new_feat = result.as_series(loader.molecules.count())
            layer.features = layer.molecules.features.with_columns(new_feat)
            layer.metadata[SEAM_SEARCH_RESULT] = result
            if show_average is not None:
                if show_average == "Filtered":
                    sigma = 0.25 / loader.scale
                    result.averages.gaussian_filter(sigma=sigma, update=True)
                _imlayer: "Image" = self._show_rec(
                    result.averages, layer.name, store=False
                )
                _imlayer.metadata[SEAM_SEARCH_RESULT] = result

            # plot all the correlation
            _Logger.print_html("<code>seam_search</code>")
            with _Logger.set_plt():
                _Logger.print(f"layer = {layer.name!r}")
                _Logger.print(f"template = {Path(template_path).as_posix()!r}")
                if anti_template_path is not None:
                    _Logger.print(
                        f"anti_template = {Path(anti_template_path).as_posix()!r}"
                    )
                widget_utils.plot_seam_search_result(result.scores, npf)

        return _seam_search_on_return

    @set_design(text=capitalize, location=STAnalysis.SeamSearch)
    def seam_search_by_feature(
        self,
        layer: MoleculesLayerType,
        by: Annotated[str, {"choices": _choice_getter("seam_search_by_feature")}],
    ):
        """
        Search for seams by a feature.

        Parameters
        ----------
        {layer}
        by : str
            Name of the feature that will be used for seam search.
        """
        layer = assert_layer(layer, self.parent_viewer)
        feat = layer.molecules.features
        if by not in feat.columns:
            raise ValueError(f"Column {by} does not exist.")
        npf = utils.roundint(layer.molecules.features[Mole.pf].max() + 1)
        seam_searcher = BooleanSeamSearcher(npf)
        result = seam_searcher.search(feat[by])
        new_feat = result.as_series(feat.shape[0])
        layer.features = layer.molecules.features.with_columns(new_feat)
        return undo_callback(layer.feature_setter(feat, layer.colormap_info))

    @set_design(text=capitalize, location=STAnalysis.SeamSearch)
    def seam_search_manually(
        self,
        layer: MoleculesLayerType,
        location: int = 0,
    ):
        """
        Search for seams manually.

        Seam location is represented by a number in the range [0, 2 * npf - 1].

        Parameters
        ----------
        {layer}
        location : int
            Seam location.
        """
        layer = assert_layer(layer, self.parent_viewer)
        feat = layer.molecules.features
        npf = utils.roundint(layer.molecules.features[Mole.pf].max() + 1)
        seam_searcher = ManualSeamSearcher(npf)
        result = seam_searcher.search(location)
        new_feat = result.as_series(feat.shape[0])
        layer.features = layer.molecules.features.with_columns(new_feat)
        return undo_callback(layer.feature_setter(feat, layer.colormap_info))

    def _seam_search_input(
        self, layer: MoleculesLayer, npf: int | None, order: int
    ) -> tuple[SubtomogramLoader, int]:
        parent = self._get_main()
        mole = layer.molecules
        loader = parent.tomogram.get_subtomogram_loader(mole, order=order)
        if npf is None:
            npf = mole.features[Mole.pf].unique().len()
        return loader, npf

    @set_design(text="Save last average", location=STAnalysis)
    def save_last_average(self, path: Path.Save[FileFilter.IMAGE]):
        """Save the lastly generated average image."""
        path = Path(path)
        img = self.last_average
        if img is None:
            raise ValueError(
                "No average image is available. You have to average subtomograms first."
            )
        return img.imsave(path)

    @average_all.started.connect
    @align_averaged.started.connect
    @align_all.started.connect
    @calculate_fsc.started.connect
    def _show_subtomogram_averaging(self):
        return self.show()

    @thread_worker.callback
    def _align_all_on_return(
        self, molecules: list[Molecules], old_layers: list[MoleculesLayer]
    ):
        """The return callback function for alignment methods."""
        main = self._get_main()
        new_layers = list[MoleculesLayer]()
        for mole, layer in zip(molecules, old_layers, strict=True):
            points = main.add_molecules(
                mole,
                name=_coerce_aligned_name(layer.name, self.parent_viewer),
                source=layer.source_component,
            )
            new_layers.append(points)
            layer.visible = False
            _Logger.print_html(f"{layer.name!r} &#8594; {points.name!r}")
        return self._undo_for_new_layer([l.name for l in old_layers], new_layers)

    @thread_worker.callback
    def _align_on_landscape_on_return(
        self,
        mole: Molecules,
        name: str,
        source=None,
        metadata: dict[str, Any] = {},
    ):
        main = self._get_main()
        if name.startswith(LANDSCAPE_PREFIX):
            nchars = len(LANDSCAPE_PREFIX)
            mole_name = name[nchars:].strip()
        else:
            mole_name = name
        points = main.add_molecules(
            mole,
            name=_coerce_aligned_name(mole_name, self.parent_viewer),
            source=source,
            metadata=metadata,
        )
        _Logger.print_html(f"{name!r} &#8594; {points.name!r}")
        return mole

    def _undo_for_new_layer(
        self,
        old_names: list[str],
        new_layers: list[MoleculesLayer],
    ):
        @undo_callback
        def out():
            main = self._get_main()
            main._try_removing_layers(new_layers)
            for name in old_names:
                if name not in main.parent_viewer.layers:
                    continue
                main.parent_viewer.layers[name].visible = True

        @out.with_redo
        def out():
            main = self._get_main()
            for points in new_layers:
                main.parent_viewer.add_layer(points)

        return out

    @nogui
    @do_not_record
    def get_template(
        self, template_path: str | Path, scale: float | None = None
    ) -> ip.ImgArray:
        """A non-GUI method to get the template"""
        if scale is None:
            scale = self._get_main().tomogram.scale
        img = self.params._norm_template_param(template_path).provide(scale)
        return ip.asarray(img, axes="zyx").set_scale(zyx=scale)

    @nogui
    @do_not_record
    def get_mask(
        self,
        mask_params: Any,
        scale: float | None = None,
        template_path: str | Path | None = None,
    ) -> ip.ImgArray:
        """A non-GUI method to get the mask."""
        if scale is None:
            scale = self._get_main().tomogram.scale
        if isinstance(mask_params, tuple):
            if template_path is None:
                raise ValueError("Template path is required when using soft-Otsu mask.")
            template = self.params._norm_template_param(template_path).provide(scale)
            radius, sigma = mask_params
            mask = pipe.soft_otsu(radius=radius, sigma=sigma).convert(template, scale)
        elif isinstance(mask_params, (str, Path)):
            mask = pipe.from_file(mask_params).provide(scale)
        else:
            raise TypeError(
                f"Cannot create mask image using parameter: {mask_params!r}"
            )
        return ip.asarray(mask, axes="zyx").set_scale(zyx=scale)

    @nogui
    @do_not_record
    def get_subtomograms(
        self,
        layers: str | MoleculesLayer | list[str | MoleculesLayer],
        shape: tuple[nm, nm, nm],
        bin_size: int = 1,
        order: int = 3,
    ) -> "Array":
        """
        A non-GUI method to get all the subtomograms as a dask array.

        Parameters
        ----------
        layers : str, MoleculesLayer or list of them
            All the layers that will be used to construct the subtomogram array.
        shape : (nm, nm, nm)
            Shape of output subtomograms.
        bin_size : int, default
            Bin size of the subtomograms.
        order : int, default 3
            Interpolation order.

        Returns
        -------
        Array
            4D Dask array.
        """
        layers = assert_list_of_layers(layers, self.parent_viewer)
        parent = self._get_main()
        tomo = parent.tomogram
        loader = tomo.get_subtomogram_loader(
            _concat_molecules(layers), shape, binsize=bin_size, order=order
        )
        return loader.construct_dask()

    def _get_simple_annealing_model(self, layer: MoleculesLayer):
        # TODO: This method should finally be moved to some utils module since
        # this analysis is independent of annealing. Currently annealing and
        # graph construction cannot be separated.
        parent = self._get_main()
        scale = parent.tomogram.scale
        return _annealing.get_annealing_model(
            layer.molecules,
            layer.source_spline,
            (0, 0, 0),
            scale,
        )

    def _construct_landscape(
        self,
        molecules: Molecules,
        template_path: Any,
        mask_params=None,
        max_shifts: tuple[nm, nm, nm] = (0.8, 0.8, 0.8),
        rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
        cutoff: float = 0.5,
        order: int = 3,
        upsample_factor: int = 5,
        bin_size: int = 1,
        method: str = "zncc",
        norm: bool = True,
    ):  # fmt: skip
        parent = self._get_main()
        loader = parent.tomogram.get_subtomogram_loader(
            molecules, binsize=bin_size, order=order
        )
        model = _get_alignment(method)
        landscape = Landscape.from_loader(
            loader=loader,
            template=template_path,
            mask=self.params._get_mask(params=mask_params),
            max_shifts=max_shifts,
            upsample_factor=upsample_factor,
            alignment_model=model.with_params(
                rotations=rotations,
                cutoff=cutoff,
                tilt=parent.tomogram.tilt_model,
            ),
        )
        return landscape.normed() if norm else landscape

last_average: ip.ImgArray | None property

Last averaged image if exists.

sub_viewer: napari.Viewer | None property

The napari viewer for subtomogram averaging.

align_all(layers, template_path, mask_params, max_shifts=(1.0, 1.0, 1.0), rotations=((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)), cutoff=0.5, interpolation=3, method='zncc', bin_size=1)

Align the input template image to all the molecules.

Parameters:

Name	Type	Description	Default
layers	list of MoleculesLayer	All the points layers of molecules to be used.	required
template_path	Path or str	Path to template image.	required
mask_params	str or (float, float)	Mask image path or dilation/Gaussian blur parameters. If a path is given, image must in the same shape as the template.	required
max_shifts	int or tuple of int	Maximum shift between subtomograms and template in nm. ZYX order.	(1.0, 1.0, 1.0)
rotations	((float, float), (float, float), (float, float))	Rotation in external degree around each axis.	((0.0, 0.0), (0.0, 0.0), (0.0, 0.0))
cutoff	float	Cutoff frequency of low-pass filter applied in each subtomogram.	0.5
interpolation	int	Interpolation order.	3
method	str	Correlation metrics for alignment.	'zncc'
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1

Source code in cylindra/widgets/sta.py

@set_design(text="Align all molecules", location=Alignment)
@dask_worker.with_progress(descs=_pdesc.align_all_fmt)
def align_all(
    self,
    layers: MoleculesLayersType,
    template_path: Annotated[_PathOrPathsOrNone, {"bind": _template_params}],
    mask_params: Annotated[Any, {"bind": _get_mask_params}],
    max_shifts: _MaxShifts = (1.0, 1.0, 1.0),
    rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
    cutoff: _CutoffFreq = 0.5,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
    method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
    bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
):  # fmt: skip
    """
    Align the input template image to all the molecules.

    Parameters
    ----------
    {layers}{template_path}{mask_params}{max_shifts}{rotations}{cutoff}
    {interpolation}{method}{bin_size}
    """
    t0 = timer()
    layers = assert_list_of_layers(layers, self.parent_viewer)
    main = self._get_main()

    combiner = MoleculesCombiner()

    loader = self._get_loader(
        binsize=bin_size,
        molecules=combiner.concat(layer.molecules for layer in layers),
        order=interpolation,
    )
    aligned_loader = loader.align(
        template=self.params._norm_template_param(
            template_path, allow_multiple=True
        ),
        mask=self.params._get_mask(params=mask_params),
        max_shifts=max_shifts,
        rotations=rotations,
        cutoff=cutoff,
        alignment_model=_get_alignment(method),
        tilt=main.tomogram.tilt_model,
    )
    molecules = combiner.split(aligned_loader.molecules, layers)
    t0.toc()
    return self._align_all_on_return.with_args(molecules, layers)

align_all_annealing(layer, template_path, mask_params=None, max_shifts=(0.8, 0.8, 0.8), rotations=((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)), cutoff=0.5, interpolation=3, range_long=(4.0, 4.28), range_lat=(5.1, 5.3), angle_max=5.0, bin_size=1, temperature_time_const=1.0, upsample_factor=5, random_seeds=(0, 1, 2, 3, 4))

2D-constrained subtomogram alignment using simulated annealing.

This alignment method considers the distance between every adjacent monomers. Two-dimensionally connected optimization can be approximated by the simulated annealing algorithm.

Parameters:

Name	Type	Description	Default
layer	MoleculesLayer	Points layer of molecules to be used.	required
template_path	Path or str	Path to template image.	required
mask_params	str or (float, float)	Mask image path or dilation/Gaussian blur parameters. If a path is given, image must in the same shape as the template.	None
max_shifts	int or tuple of int	Maximum shift between subtomograms and template in nm. ZYX order.	(0.8, 0.8, 0.8)
rotations	((float, float), (float, float), (float, float))	Rotation in external degree around each axis.	((0.0, 0.0), (0.0, 0.0), (0.0, 0.0))
cutoff	float	Cutoff frequency of low-pass filter applied in each subtomogram.	0.5
interpolation	int	Interpolation order.	3
range_long	(float, float)	Minimum and maximum allowed distances between longitudinally consecutive monomers	(4.0, 4.28)
range_lat	(float, float)	Minimum and maximum allowed distances between laterally consecutive monomers	(5.1, 5.3)
angle_max	float	Maximum allowed angle between longitudinally consecutive monomers and the Y axis.	5.0
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1
temperature_time_const	float	Time constant of the temperature decay during annealing. Larger value results in slower annealing. 1.0 is a moderate value.	1.0
upsample_factor	int	Upsampling factor of ZNCC landscape. Be careful not to set this parameter too large. Calculation will take much longer for larger upsample_factor.	5
random_seeds	iterable of int	Random seed integers. Number of integers will be the number of trials.	(0, 1, 2, 3, 4)

Source code in cylindra/widgets/sta.py

@set_design(text="Simulated annealing", location=Alignment)
@dask_worker.with_progress(descs=_pdesc.align_annealing_fmt)
def align_all_annealing(
    self,
    layer: MoleculesLayerType,
    template_path: Annotated[_PathOrPathsOrNone, {"bind": _template_params}],
    mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
    max_shifts: _MaxShifts = (0.8, 0.8, 0.8),
    rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
    cutoff: _CutoffFreq = 0.5,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
    range_long: _DistRangeLon = (4.0, 4.28),
    range_lat: _DistRangeLat = (5.1, 5.3),
    angle_max: _AngleMaxLon = 5.0,
    bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
    temperature_time_const: Annotated[float, {"min": 0.01, "max": 10.0}] = 1.0,
    upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 5,
    random_seeds: _RandomSeeds = (0, 1, 2, 3, 4),
):  # fmt: skip
    """
    2D-constrained subtomogram alignment using simulated annealing.

    This alignment method considers the distance between every adjacent monomers.
    Two-dimensionally connected optimization can be approximated by the simulated
    annealing algorithm.

    Parameters
    ----------
    {layer}{template_path}{mask_params}{max_shifts}{rotations}{cutoff}
    {interpolation}{range_long}{range_lat}{angle_max}{bin_size}
    {temperature_time_const}{upsample_factor}{random_seeds}
    """
    t0 = timer()
    layer = assert_layer(layer, self.parent_viewer)
    if layer.source_spline is None:
        raise ValueError("RMA requires a spline.")
    main = self._get_main()
    landscape = self._construct_landscape(
        molecules=layer.molecules,
        template_path=template_path,
        mask_params=mask_params,
        max_shifts=max_shifts,
        rotations=rotations,
        cutoff=cutoff,
        order=interpolation,
        bin_size=bin_size,
        upsample_factor=upsample_factor,
    )
    yield
    mole, results = landscape.run_annealing_along_spline(
        layer.source_spline,
        range_long=range_long,
        range_lat=range_lat,
        angle_max=angle_max,
        temperature_time_const=temperature_time_const,
        random_seeds=random_seeds,
    )
    t0.toc()

    @thread_worker.callback
    def _on_return():
        points = main.add_molecules(
            mole,
            name=_coerce_aligned_name(layer.name, self.parent_viewer),
            source=layer.source_component,
            metadata={ANNEALING_RESULT: results[0]},
        )
        layer.visible = False
        with _Logger.set_plt():
            _annealing.plot_annealing_result(results)
        return self._undo_for_new_layer([layer.name], [points])

    return _on_return

align_all_template_free(layers, mask_params, size=12.0, max_shifts=(1.0, 1.0, 1.0), rotations=((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)), cutoff=0.5, interpolation=3, method='zncc', bin_size=1)

Run template-free alignment for the given layers (EXPERIMENTAL).

Parameters:

Name	Type	Description	Default
layers	list of MoleculesLayer	All the points layers of molecules to be used.	required
mask_params	str or (float, float)	Mask image path or dilation/Gaussian blur parameters. If a path is given, image must in the same shape as the template.	required
size	nm	Size of the template in nm. Use the size of template image by default.	12.0
max_shifts	int or tuple of int	Maximum shift between subtomograms and template in nm. ZYX order.	(1.0, 1.0, 1.0)
rotations	((float, float), (float, float), (float, float))	Rotation in external degree around each axis.	((0.0, 0.0), (0.0, 0.0), (0.0, 0.0))
cutoff	float	Cutoff frequency of low-pass filter applied in each subtomogram.	0.5
interpolation	int	Interpolation order.	3
method	str	Correlation metrics for alignment.	'zncc'
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1

Source code in cylindra/widgets/sta.py

@set_design(text="Align all (template-free)", location=Alignment)
@dask_worker.with_progress(descs=_pdesc.align_template_free_fmt)
def align_all_template_free(
    self,
    layers: MoleculesLayersType,
    mask_params: Annotated[Any, {"bind": _get_mask_params}],
    size: _SubVolumeSize = 12.0,
    max_shifts: _MaxShifts = (1.0, 1.0, 1.0),
    rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
    cutoff: _CutoffFreq = 0.5,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
    method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
    bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
):  # fmt: skip
    """
    Run template-free alignment for the given layers (EXPERIMENTAL).

    Parameters
    ----------
    {layers}{mask_params}{size}{max_shifts}{rotations}{cutoff}{interpolation}
    {method}{bin_size}
    """
    t0 = timer()
    layers = assert_list_of_layers(layers, self.parent_viewer)
    main = self._get_main()
    combiner = MoleculesCombiner()
    molecules = combiner.concat(layer.molecules for layer in layers)
    mask = self.params._get_mask(params=mask_params)
    if size is None:
        raise NotImplementedError("'size' must be given.")
    else:
        shape = tuple(
            main.tomogram.nm2pixel(self._get_shape_in_nm(size), binsize=bin_size)
        )

    aligned_loader = (
        self._get_loader(binsize=bin_size, molecules=molecules, order=interpolation)
        .reshape(shape=shape)
        .align_no_template(
            mask=mask,
            max_shifts=max_shifts,
            rotations=rotations,
            cutoff=cutoff,
            alignment_model=_get_alignment(method),
            tilt=main.tomogram.tilt_model,
        )
    )
    molecules = combiner.split(aligned_loader.molecules, layers)
    t0.toc()
    return self._align_all_on_return.with_args(molecules, layers)

align_all_viterbi(layer, template_path, mask_params=None, max_shifts=(0.8, 0.8, 0.8), rotations=((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)), cutoff=0.5, interpolation=3, range_long=(4.0, 4.28), angle_max=5.0, bin_size=1, upsample_factor=5)

Subtomogram alignment using 1D Viterbi alignment.

1D Viterbi alignment is an alignment algorithm that considers the distance and the skew angle between every longitudinally adjacent monomers. The classical Viterbi algorithm is used to find the global optimal solution of the alignment. Note that Viterbi alignment is data size dependent, i.e. the alignment result of a molecule may vary depending on the total number of molecules in the dataset.

Parameters:

Name	Type	Description	Default
layer	MoleculesLayer	Points layer of molecules to be used.	required
template_path	Path or str	Path to template image.	required
mask_params	str or (float, float)	Mask image path or dilation/Gaussian blur parameters. If a path is given, image must in the same shape as the template.	None
max_shifts	int or tuple of int	Maximum shift between subtomograms and template in nm. ZYX order.	(0.8, 0.8, 0.8)
rotations	((float, float), (float, float), (float, float))	Rotation in external degree around each axis.	((0.0, 0.0), (0.0, 0.0), (0.0, 0.0))
cutoff	float	Cutoff frequency of low-pass filter applied in each subtomogram.	0.5
interpolation	int	Interpolation order.	3
range_long	(float, float)	Minimum and maximum allowed distances between longitudinally consecutive monomers	(4.0, 4.28)
angle_max	float	Maximum allowed angle between longitudinally consecutive monomers and the Y axis.	5.0
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1
upsample_factor	int	Upsampling factor of ZNCC landscape. Be careful not to set this parameter too large. Calculation will take much longer for larger upsample_factor.	5

Source code in cylindra/widgets/sta.py

@set_design(text="Viterbi Alignment", location=Alignment)
@dask_worker.with_progress(descs=_pdesc.align_viterbi_fmt)
def align_all_viterbi(
    self,
    layer: MoleculesLayerType,
    template_path: Annotated[_PathOrPathsOrNone, {"bind": _template_params}],
    mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
    max_shifts: _MaxShifts = (0.8, 0.8, 0.8),
    rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
    cutoff: _CutoffFreq = 0.5,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
    range_long: _DistRangeLon = (4.0, 4.28),
    angle_max: _AngleMaxLon = 5.0,
    bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
    upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 5,
):  # fmt: skip
    """
    Subtomogram alignment using 1D Viterbi alignment.

    1D Viterbi alignment is an alignment algorithm that considers the distance and
    the skew angle between every longitudinally adjacent monomers. The classical
    Viterbi algorithm is used to find the global optimal solution of the alignment.
    Note that Viterbi alignment is data size dependent, i.e. the alignment result
    of a molecule may vary depending on the total number of molecules in the dataset.

    Parameters
    ----------
    {layer}{template_path}{mask_params}{max_shifts}{rotations}{cutoff}
    {interpolation}{range_long}{angle_max}{bin_size}{upsample_factor}
    """
    t0 = timer()
    layer = assert_layer(layer, self.parent_viewer)
    landscape = self._construct_landscape(
        molecules=layer.molecules,
        template_path=template_path,
        mask_params=mask_params,
        max_shifts=max_shifts,
        rotations=rotations,
        cutoff=cutoff,
        order=interpolation,
        upsample_factor=upsample_factor,
        bin_size=bin_size,
    )

    yield
    mole = landscape.run_viterbi_along_spline(
        spl=layer.source_spline,
        range_long=range_long,
        angle_max=angle_max,
    )
    t0.toc()
    return self._align_all_on_return.with_args([mole], [layer])

align_averaged(layers, template_path, mask_params, max_shifts=None, rotations=((0.0, 0.0), (15.0, 1.0), (3.0, 1.0)), bin_size=1, method='zncc')

Align the averaged image at current monomers to the template image.

This function creates a new layer with transformed monomers, which should align well with template image.

Parameters:

Name	Type	Description	Default
layers	list of MoleculesLayer	All the points layers of molecules to be used.	required
template_path	Path or str	Path to template image.	required
mask_params	str or (float, float)	Mask image path or dilation/Gaussian blur parameters. If a path is given, image must in the same shape as the template.	required
max_shifts	int or tuple of int	Maximum shift between subtomograms and template in nm. ZYX order.	None
rotations	((float, float), (float, float), (float, float))	Rotation in external degree around each axis.	((0.0, 0.0), (15.0, 1.0), (3.0, 1.0))
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1
method	str	Correlation metrics for alignment.	'zncc'

Source code in cylindra/widgets/sta.py

@set_design(text="Align average to template", location=Alignment)
@dask_worker.with_progress(descs=_pdesc.align_averaged_fmt)
def align_averaged(
    self,
    layers: MoleculesLayersType,
    template_path: Annotated[_PathOrNone, {"bind": _template_param}],
    mask_params: Annotated[Any, {"bind": _get_mask_params}],
    max_shifts: Optional[_MaxShifts] = None,
    rotations: _Rotations = ((0.0, 0.0), (15.0, 1.0), (3.0, 1.0)),
    bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
    method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
):  # fmt: skip
    """
    Align the averaged image at current monomers to the template image.

    This function creates a new layer with transformed monomers, which should
    align well with template image.

    Parameters
    ----------
    {layers}{template_path}{mask_params}{max_shifts}{rotations}{bin_size}{method}
    """
    t0 = timer()
    layers = assert_list_of_layers(layers, self.parent_viewer)
    parent = self._get_main()

    new_layers = list[MoleculesLayer]()

    @thread_worker.callback
    def _on_yield(mole_trans: Molecules, layer: MoleculesLayer):
        points = parent.add_molecules(
            mole_trans,
            name=_coerce_aligned_name(layer.name, self.parent_viewer),
            source=layer.source_component,
        )
        new_layers.append(points)
        layer.visible = False
        _Logger.print_html(f"{layer.name!r} → {points.name!r}")

    mole = layers[0].molecules
    loader = self._get_loader(bin_size, mole, order=1)
    template, mask = loader.normalize_input(
        template=self.params._norm_template_param(
            template_path, allow_multiple=False
        ),
        mask=self.params._get_mask(params=mask_params),
    )
    temp_norm = utils.normalize_image(template)

    _scale = parent.tomogram.scale * bin_size

    if max_shifts is None:
        max_shifts = _default_align_averaged_shifts(mole)

    model = _get_alignment(method)(
        template,
        mask,
        rotations=rotations,
        tilt=None,  # NOTE: because input is an average
    )
    _spl_globs = list[
        tuple[weakref.ReferenceType["CylSpline"], pl.DataFrame, pl.DataFrame]
    ]()
    for layer in layers:
        mole = layer.molecules
        loader = self._get_loader(bin_size, mole, order=1)
        _img_trans, result = model.fit(
            loader.average(template.shape),
            max_shifts=[_s / _scale for _s in max_shifts],
        )

        rotator = Rotation.from_quat(result.quat)
        svec = result.shift * _scale
        _mole_trans = mole.linear_transform(result.shift * _scale, rotator)

        # write offsets to spline globalprops if available
        if spl := layer.source_spline:
            _mole_trans = _update_mole_pos(_mole_trans, mole, spl)
            if spl.radius is None:
                _radius: nm = cylmeasure.calc_radius(mole, spl).mean()
            else:
                _radius = spl.radius
            _glob_old = spl.props.glob.clone()
            _glob_new = _update_offset(spl, rotator.apply(svec), _radius)
            spl.props.glob = _glob_new
            _spl_globs.append((weakref.ref(spl), _glob_old, _glob_new))

        yield _on_yield.with_args(_mole_trans, layer)

        # create images for visualization in the logger. Image is magenta, template is green
        img_norm = utils.normalize_image(_img_trans)
        merge = np.stack([img_norm, temp_norm, img_norm], axis=-1)
        with _Logger.set_plt():
            widget_utils.plot_projections(merge)

        # logging
        rvec = rotator.as_rotvec()
        _fmt = "  {:.2f}  ".format
        _Logger.print_table(
            [
                ["", "X", "Y", "Z"],
                ["Shift (nm)", _fmt(svec[2]), _fmt(svec[1]), _fmt(svec[0])],
                ["Rot vector", _fmt(rvec[2]), _fmt(rvec[1]), _fmt(rvec[0])],
            ],
            header=False,
            index=False,
        )

    t0.toc()

    @thread_worker.callback
    def _align_averaged_on_return():
        @undo_callback
        def _out():
            parent._try_removing_layers(new_layers)
            for spl_ref, old, _ in _spl_globs:
                if spl := spl_ref():
                    spl.props.glob = old

        @_out.with_redo
        def _out():
            parent._add_layers_future(new_layers)()
            for spl_ref, _, new in _spl_globs:
                if spl := spl_ref():
                    spl.props.glob = new

        return _out

    return _align_averaged_on_return

average_all(layers, size=None, interpolation=1, bin_size=1)

Subtomogram averaging using all the molecules in the selected layer(s).

If multiple layers are selected, subtomograms around all the molecules will be averaged.

Parameters:

Name	Type	Description	Default
layers	list of MoleculesLayer	All the points layers of molecules to be used.	required
size	nm	Size of the template in nm. Use the size of template image by default.	None
interpolation	int	Interpolation order.	1
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1

Source code in cylindra/widgets/sta.py

@set_design(text="Average all molecules", location=Averaging)
@dask_worker.with_progress(desc=_pdesc.fmt_layers("Subtomogram averaging of {!r}"))
def average_all(
    self,
    layers: MoleculesLayersType,
    size: _SubVolumeSize = None,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 1,
    bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
):
    """
    Subtomogram averaging using all the molecules in the selected layer(s).

    If multiple layers are selected, subtomograms around all the molecules will
    be averaged.

    Parameters
    ----------
    {layers}{size}{interpolation}{bin_size}
    """
    t0 = timer()
    layers = assert_list_of_layers(layers, self.parent_viewer)
    parent = self._get_main()
    tomo = parent.tomogram
    shape = self._get_shape_in_nm(size)
    loader = tomo.get_subtomogram_loader(
        _concat_molecules(layers), shape, binsize=bin_size, order=interpolation
    )
    img = ip.asarray(loader.average(), axes="zyx")
    img.set_scale(zyx=loader.scale, unit="nm")
    t0.toc()
    return self._show_rec.with_args(img, f"[AVG]{_avg_name(layers)}")

average_filtered(layers, size=None, predicate="col('pf-id') == 0", interpolation=1, bin_size=1)

Subtomogram averaging using molecules filtered by the given expression.

This method first concatenate molecules in the selected layers, and then filter them by the predicate.

Parameters:

Name	Type	Description	Default
layers	list of MoleculesLayer	All the points layers of molecules to be used.	required
size	nm	Size of the template in nm. Use the size of template image by default.	None
predicate	str or polars expression	Filter expression to select molecules.	"col('pf-id') == 0"
interpolation	int	Interpolation order.	1
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1

Source code in cylindra/widgets/sta.py

@set_design(text="Average filtered", location=Averaging)
@dask_worker.with_progress(desc=_pdesc.fmt_layers("Filtered subtomogram averaging of {!r}"))  # fmt: skip
def average_filtered(
    self,
    layers: MoleculesLayersType,
    size: _SubVolumeSize = None,
    predicate: PolarsExprStr = "col('pf-id') == 0",
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 1,
    bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
):
    """
    Subtomogram averaging using molecules filtered by the given expression.

    This method first concatenate molecules in the selected layers, and then filter them
    by the predicate.

    Parameters
    ----------
    {layers}{size}
    predicate : str or polars expression
        Filter expression to select molecules.
    {interpolation}{bin_size}
    """
    t0 = timer()
    layers = assert_list_of_layers(layers, self.parent_viewer)
    parent = self._get_main()
    tomo = parent.tomogram
    shape = self._get_shape_in_nm(size)
    loader = tomo.get_subtomogram_loader(
        _concat_molecules(layers), shape, binsize=bin_size, order=interpolation
    )
    avg = loader.filter(widget_utils.norm_expr(predicate)).average()
    img = ip.asarray(avg, axes="zyx")
    img.set_scale(zyx=loader.scale, unit="nm")
    t0.toc()
    return self._show_rec.with_args(img, f"[AVG]{_avg_name(layers)}")

average_groups(layers, size=None, by="col('pf-id')", interpolation=1, bin_size=1)

Group-wise subtomogram averaging using molecules grouped by the given expression.

This method first group molecules by its features, and then average each group. This method is useful for such as get average of each protofilament and segmented subtomogram averaging.

Parameters:

Name	Type	Description	Default
layers	list of MoleculesLayer	All the points layers of molecules to be used.	required
size	nm	Size of the template in nm. Use the size of template image by default.	None
by	str or polars expression	Expression to group molecules.	"col('pf-id')"
interpolation	int	Interpolation order.	1
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1

Source code in cylindra/widgets/sta.py

@set_design(text="Average group-wise", location=Averaging)
@dask_worker.with_progress(desc=_pdesc.fmt_layers("Grouped subtomogram averaging of {!r}"))  # fmt: skip
def average_groups(
    self,
    layers: MoleculesLayersType,
    size: _SubVolumeSize = None,
    by: PolarsExprStr = "col('pf-id')",
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 1,
    bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
):
    """
    Group-wise subtomogram averaging using molecules grouped by the given expression.

    This method first group molecules by its features, and then average each group.
    This method is useful for such as get average of each protofilament and segmented
    subtomogram averaging.

    Parameters
    ----------
    {layers}{size}
    by : str or polars expression
        Expression to group molecules.
    {interpolation}{bin_size}
    """
    t0 = timer()
    layers = assert_list_of_layers(layers, self.parent_viewer)
    parent = self._get_main()
    tomo = parent.tomogram
    shape = self._get_shape_in_nm(size)
    loader = tomo.get_subtomogram_loader(
        _concat_molecules(layers), shape, binsize=bin_size, order=interpolation
    )
    expr = widget_utils.norm_expr(by)
    avg_dict = loader.groupby(expr).average()
    avgs = np.stack([avg_dict[k] for k in sorted(avg_dict.keys())], axis=0)
    img = ip.asarray(avgs, axes="pzyx")
    img.set_scale(zyx=loader.scale, unit="nm")
    t0.toc()
    return self._show_rec.with_args(img, f"[AVG]{_avg_name(layers)}", store=False)

average_subset(layers, size=None, method='steps', number=64, bin_size=1)

Subtomogram averaging using a subset of subvolumes.

If multiple layers are selected, subtomograms around all the molecules will be concatenated before choosing a subset.

Parameters:

Name	Type	Description	Default
layers	list of MoleculesLayer	All the points layers of molecules to be used.	required
size	nm	Size of the template in nm. Use the size of template image by default.	None
method	str	How to choose subtomogram subset. (1) steps: Each 'steps' subtomograms from the tip of spline. (2) first: First subtomograms. (3) last: Last subtomograms. (4) random: choose randomly.	'steps'
number	(int, default)	Number of subtomograms to use.	64
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1

Source code in cylindra/widgets/sta.py

@set_design(text="Average subset of molecules", location=Averaging)
@dask_worker.with_progress(desc=_pdesc.fmt_layers("Subtomogram averaging (subset) of {!r}"))  # fmt: skip
def average_subset(
    self,
    layers: MoleculesLayersType,
    size: _SubVolumeSize = None,
    method: Literal["steps", "first", "last", "random"] = "steps",
    number: int = 64,
    bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
):
    """
    Subtomogram averaging using a subset of subvolumes.

    If multiple layers are selected, subtomograms around all the molecules will
    be concatenated before choosing a subset.

    Parameters
    ----------
    {layers}{size}
    method : str, optional
        How to choose subtomogram subset.
        (1) steps: Each 'steps' subtomograms from the tip of spline.
        (2) first: First subtomograms.
        (3) last: Last subtomograms.
        (4) random: choose randomly.
    number : int, default
        Number of subtomograms to use.
    {bin_size}
    """
    t0 = timer()
    layers = assert_list_of_layers(layers, self.parent_viewer)
    parent = self._get_main()
    molecules = _concat_molecules(layers)
    nmole = len(molecules)
    shape = self._get_shape_in_nm(size)
    sl = _get_slice_for_average_subset(method, nmole, number)
    mole = molecules.subset(sl)
    loader = parent.tomogram.get_subtomogram_loader(
        mole, shape, binsize=bin_size, order=1
    )
    img = ip.asarray(loader.average(), axes="zyx").set_scale(zyx=loader.scale)
    t0.toc()
    return self._show_rec.with_args(img, f"[AVG(n={number})]{_avg_name(layers)}")

calculate_correlation(layers, template_path, mask_params=None, interpolation=3, cutoff=0.5, bin_size=1, method='zncc', column_prefix='score')

Calculate correlation between template images and the subtomograms.

This method will load every subtomograms, calculate the correlation between the template images and each subtomogram, and save the correlation values as new columns in the molecules features.

Parameters:

Name	Type	Description	Default
layers	list of MoleculesLayer	All the points layers of molecules to be used.	required
template_path	Path or str	Path to template image.	required
mask_params	str or (float, float)	Mask image path or dilation/Gaussian blur parameters. If a path is given, image must in the same shape as the template.	None
interpolation	int	Interpolation order.	3
cutoff	float	Cutoff frequency of low-pass filter applied in each subtomogram.	0.5
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1
method	str	Correlation metrics for alignment.	'zncc'
column_prefix	str	Prefix of the column names of the calculated correlations.	"score"

Source code in cylindra/widgets/sta.py

@set_design(text=capitalize, location=STAnalysis)
@dask_worker.with_progress(desc=_pdesc.fmt_layers("Calculating correlations of {!r}"))  # fmt: skip
def calculate_correlation(
    self,
    layers: MoleculesLayersType,
    template_path: Annotated[_PathOrPathsOrNone, {"bind": _template_params}],
    mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
    cutoff: _CutoffFreq = 0.5,
    bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
    method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
    column_prefix: str = "score",
):
    """
    Calculate correlation between template images and the subtomograms.

    This method will load every subtomograms, calculate the correlation between
    the template images and each subtomogram, and save the correlation values
    as new columns in the molecules features.

    Parameters
    ----------
    {layers}{template_path}{mask_params}{interpolation}{cutoff}{bin_size}{method}
    column_prefix : str, default "score"
        Prefix of the column names of the calculated correlations.
    """
    layers = assert_list_of_layers(layers, self.parent_viewer)
    main = self._get_main()
    combiner = MoleculesCombiner()

    if isinstance(template_path, (Path, str)):
        template_path = [template_path]
    mask = self.params._get_mask(mask_params)
    all_mole = combiner.concat(layer.molecules for layer in layers)

    out = main.tomogram.get_subtomogram_loader(
        all_mole,
        order=interpolation,
        binsize=bin_size,
    ).score(
        templates=[pipe.from_file(t) for t in template_path],
        mask=mask,
        alignment_model=_get_alignment(method),
        cutoff=cutoff,
        tilt=main.tomogram.tilt_model,
    )
    all_mole = all_mole.with_features(
        pl.Series(f"{column_prefix}_{i}", col) for i, col in enumerate(out)
    )

    @thread_worker.callback
    def _on_return():
        moles = combiner.split(all_mole, layers)
        for layer, each_mole in zip(layers, moles, strict=True):
            features = each_mole.features.select(
                [f"{column_prefix}_{i}" for i in range(len(out))]
            )
            layer.set_molecules_with_new_features(
                layer.molecules.with_features(features)
            )

    return _on_return

calculate_fsc(layers, template_path=None, mask_params=None, size=None, seed=0, interpolation=1, n_pairs=1, show_average=True, dfreq=None)

Calculate Fourier Shell Correlation using the selected monomer layer.

Parameters:

Name	Type	Description	Default
layers	list of MoleculesLayer	All the points layers of molecules to be used.	required
template_path	template input type	Used only when soft-Otsu mask parameters are given.	None
mask_params	str or (float, float)	Mask image path or dilation/Gaussian blur parameters. If a path is given, image must in the same shape as the template.	None
size	nm	Size of the template in nm. Use the size of template image by default.	None
seed	int	Random seed used for subtomogram sampling.	0
interpolation	int	Interpolation order.	1
n_pairs	int	How many sets of image pairs will be generated to average FSC.	1
show_average	bool	If true, subtomogram average will be shown after FSC calculation.	True
dfreq	float	Precision of frequency to calculate FSC. "0.02" means that FSC will be calculated at frequency 0.01, 0.03, 0.05, ..., 0.45.	0.02

Source code in cylindra/widgets/sta.py

@set_design(text="Calculate FSC", location=STAnalysis)
@dask_worker.with_progress(desc=_pdesc.fmt_layers("Calculating FSC of {!r}"))
def calculate_fsc(
    self,
    layers: MoleculesLayersType,
    template_path: Annotated[_PathOrNone, {"bind": _template_param}] = None,
    mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
    size: _SubVolumeSize = None,
    seed: Annotated[Optional[int], {"text": "Do not use random seed."}] = 0,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 1,
    n_pairs: Annotated[int, {"min": 1, "label": "number of image pairs"}] = 1,
    show_average: bool = True,
    dfreq: FSCFreq = None,
):
    """
    Calculate Fourier Shell Correlation using the selected monomer layer.

    Parameters
    ----------
    {layers}
    template_path : template input type
        Used only when soft-Otsu mask parameters are given.
    {mask_params}{size}
    seed : int, optional
        Random seed used for subtomogram sampling.
    {interpolation}
    n_pairs : int, default 1
        How many sets of image pairs will be generated to average FSC.
    show_average : bool, default True
        If true, subtomogram average will be shown after FSC calculation.
    dfreq : float, default 0.02
        Precision of frequency to calculate FSC. "0.02" means that FSC will be
        calculated at frequency 0.01, 0.03, 0.05, ..., 0.45.
    """
    t0 = timer()
    layers = assert_list_of_layers(layers, self.parent_viewer)
    main = self._get_main()
    mole = _concat_molecules(layers)

    loader = main.tomogram.get_subtomogram_loader(mole, order=interpolation)
    template, mask = loader.normalize_input(
        template=self.params._norm_template_param(template_path, allow_none=True),
        mask=self.params._get_mask(params=mask_params),
    )
    fsc, (img_0, img_1), img_mask = loader.reshape(
        template=template if size is None else None,
        mask=mask,
        shape=None if size is None else (main.tomogram.nm2pixel(size),) * 3,
    ).fsc_with_halfmaps(mask, seed=seed, n_set=n_pairs, dfreq=dfreq, squeeze=False)

    def _as_imgarray(im, axes: str = "zyx") -> ip.ImgArray | None:
        if np.isscalar(im):
            return None
        return ip.asarray(im, axes=axes).set_scale(zyx=loader.scale, unit="nm")

    if show_average:
        avg = (img_0[0] + img_1[0]) / 2
        img_avg = _as_imgarray(avg)
    else:
        img_avg = None

    result = widget_utils.FscResult.from_dataframe(fsc, loader.scale)
    criteria = [0.5, 0.143]
    _name = _avg_name(layers)
    t0.toc()

    @thread_worker.callback
    def _calculate_fsc_on_return():
        _Logger.print_html(f"<b>Fourier Shell Correlation of {_name!r}</b>")
        with _Logger.set_plt():
            result.plot(criteria)
        for _c in criteria:
            _r = result.get_resolution(_c)
            _Logger.print_html(f"Resolution at FSC={_c:.3f} ... <b>{_r:.3f} nm</b>")

        if img_avg is not None:
            _imlayer: "Image" = self._show_rec(img_avg, name=f"[AVG]{_name}")
            _imlayer.metadata["fsc"] = result
            _imlayer.metadata["fsc_halfmaps"] = (
                _as_imgarray(img_0, axes="izyx"),
                _as_imgarray(img_1, axes="izyx"),
            )
            _imlayer.metadata["fsc_mask"] = _as_imgarray(img_mask)

    return _calculate_fsc_on_return

classify_pca(layer, template_path=None, mask_params=None, size=None, cutoff=0.5, interpolation=3, bin_size=1, n_components=2, n_clusters=2, seed=0)

Classify molecules in a layer using PCA and K-means clustering.

Parameters:

Name	Type	Description	Default
layer	MoleculesLayer	Points layer of molecules to be used.	required
template_path	template input type	Used only when soft-Otsu mask parameters are given.	None
mask_params	str or (float, float)	Mask image path or dilation/Gaussian blur parameters. If a path is given, image must in the same shape as the template.	None
size	nm	Size of the template in nm. Use the size of template image by default.	None
cutoff	float	Cutoff frequency of low-pass filter applied in each subtomogram.	0.5
interpolation	int	Interpolation order.	3
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1
n_components	int	The number of PCA dimensions.	2
n_clusters	(int, default)	The number of clusters.	2
seed	int	Random seed.	0

Source code in cylindra/widgets/sta.py

@set_design(text="PCA/K-means classification", location=STAnalysis)
@dask_worker.with_progress(descs=_pdesc.classify_pca_fmt)
def classify_pca(
    self,
    layer: MoleculesLayerType,
    template_path: Annotated[_PathOrNone, {"bind": _template_param}] = None,
    mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
    size: _SubVolumeSize = None,
    cutoff: _CutoffFreq = 0.5,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
    bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
    n_components: Annotated[int, {"min": 2, "max": 20}] = 2,
    n_clusters: Annotated[int, {"min": 2, "max": 100}] = 2,
    seed: Annotated[Optional[int], {"text": "Do not use random seed."}] = 0,
):  # fmt: skip
    """
    Classify molecules in a layer using PCA and K-means clustering.

    Parameters
    ----------
    {layer}
    template_path : template input type
        Used only when soft-Otsu mask parameters are given.
    {mask_params}{size}{cutoff}{interpolation}{bin_size}
    n_components : int, default 2
        The number of PCA dimensions.
    n_clusters : int, default
        The number of clusters.
    seed : int, default 0
        Random seed.
    """
    from cylindra.widgets.subwidgets import PcaViewer

    t0 = timer()
    layer = assert_layer(layer, self.parent_viewer)
    tomo = self._get_main().tomogram
    loader = self._get_loader(
        binsize=bin_size, molecules=layer.molecules, order=interpolation
    )
    template, mask = loader.normalize_input(
        template=self.params._norm_template_param(template_path, allow_none=True),
        mask=self.params._get_mask(params=mask_params),
    )
    shape = None
    if size is not None and mask is None:
        shape = (tomo.nm2pixel(size, binsize=bin_size),) * 3
    out, pca = loader.reshape(
        template=template if mask is None and shape is None else None,
        mask=mask,
        shape=shape,
    ).classify(
        mask=mask,
        seed=seed,
        cutoff=cutoff,
        n_components=n_components,
        n_clusters=n_clusters,
        label_name="cluster",
    )

    avgs_dict = out.groupby("cluster").average()
    avgs = ip.asarray(
        np.stack(list(avgs_dict.values()), axis=0), axes=["cluster", "z", "y", "x"]
    ).set_scale(zyx=loader.scale, unit="nm")
    t0.toc()

    @thread_worker.callback
    def _on_return():
        layer.molecules = out.molecules  # update features
        pca_viewer = PcaViewer(pca)
        pca_viewer.native.setParent(self.native, pca_viewer.native.windowFlags())
        pca_viewer.show()
        self._show_rec(avgs, name=f"[PCA]{layer.name}", store=False)
        ACTIVE_WIDGETS.add(pca_viewer)

    return _on_return

construct_landscape(layer, template_path, mask_params=None, max_shifts=(0.8, 0.8, 0.8), rotations=((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)), cutoff=0.5, interpolation=3, bin_size=1, upsample_factor=5, method='zncc', norm=True)

Construct a cross-correlation landscape for subtomogram alignment.

Parameters:

Name	Type	Description	Default
layer	MoleculesLayer	Points layer of molecules to be used.	required
template_path	Path or str	Path to template image.	required
mask_params	str or (float, float)	Mask image path or dilation/Gaussian blur parameters. If a path is given, image must in the same shape as the template.	None
max_shifts	int or tuple of int	Maximum shift between subtomograms and template in nm. ZYX order.	(0.8, 0.8, 0.8)
rotations	((float, float), (float, float), (float, float))	Rotation in external degree around each axis.	((0.0, 0.0), (0.0, 0.0), (0.0, 0.0))
cutoff	float	Cutoff frequency of low-pass filter applied in each subtomogram.	0.5
interpolation	int	Interpolation order.	3
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1
upsample_factor	int	Upsampling factor of ZNCC landscape. Be careful not to set this parameter too large. Calculation will take much longer for larger upsample_factor.	5
method	str	Correlation metrics for alignment.	'zncc'
norm	bool	If true, each landscape will be normalized by its mean and standard deviation.	True

Source code in cylindra/widgets/sta.py

@set_design(text=capitalize, location=LandscapeMenu)
@dask_worker.with_progress(descs=_pdesc.construct_landscape_fmt)
def construct_landscape(
    self,
    layer: MoleculesLayerType,
    template_path: Annotated[_PathOrPathsOrNone, {"bind": _template_params}],
    mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
    max_shifts: _MaxShifts = (0.8, 0.8, 0.8),
    rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
    cutoff: _CutoffFreq = 0.5,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
    bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
    upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 5,
    method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
    norm: bool = True,
):
    """
    Construct a cross-correlation landscape for subtomogram alignment.

    Parameters
    ----------
    {layer}{template_path}{mask_params}{max_shifts}{rotations}{cutoff}
    {interpolation}{bin_size}{upsample_factor}{method}
    norm: bool, default True
        If true, each landscape will be normalized by its mean and standard
        deviation.
    """
    layer = assert_layer(layer, self.parent_viewer)
    lnd = self._construct_landscape(
        molecules=layer.molecules,
        template_path=template_path,
        mask_params=mask_params,
        max_shifts=max_shifts,
        rotations=rotations,
        cutoff=cutoff,
        order=interpolation,
        bin_size=bin_size,
        upsample_factor=upsample_factor,
        norm=norm,
        method=method,
    )
    surf = LandscapeSurface(lnd, name=f"{LANDSCAPE_PREFIX}{layer.name}")
    surf.source_component = layer.source_component

    @thread_worker.callback
    def _on_return():
        self.parent_viewer.add_layer(surf)
        self._get_main()._reserved_layers.to_be_removed.add(surf)
        layer.visible = False

    return _on_return

get_mask(mask_params, scale=None, template_path=None)

A non-GUI method to get the mask.

Source code in cylindra/widgets/sta.py

@nogui
@do_not_record
def get_mask(
    self,
    mask_params: Any,
    scale: float | None = None,
    template_path: str | Path | None = None,
) -> ip.ImgArray:
    """A non-GUI method to get the mask."""
    if scale is None:
        scale = self._get_main().tomogram.scale
    if isinstance(mask_params, tuple):
        if template_path is None:
            raise ValueError("Template path is required when using soft-Otsu mask.")
        template = self.params._norm_template_param(template_path).provide(scale)
        radius, sigma = mask_params
        mask = pipe.soft_otsu(radius=radius, sigma=sigma).convert(template, scale)
    elif isinstance(mask_params, (str, Path)):
        mask = pipe.from_file(mask_params).provide(scale)
    else:
        raise TypeError(
            f"Cannot create mask image using parameter: {mask_params!r}"
        )
    return ip.asarray(mask, axes="zyx").set_scale(zyx=scale)

get_subtomograms(layers, shape, bin_size=1, order=3)

A non-GUI method to get all the subtomograms as a dask array.

Parameters:

Name	Type	Description	Default
layers	str, MoleculesLayer or list of them	All the layers that will be used to construct the subtomogram array.	required
shape	(nm, nm, nm)	Shape of output subtomograms.	required
bin_size	(int, default)	Bin size of the subtomograms.	1
order	int	Interpolation order.	3

Returns:

Type	Description
Array	4D Dask array.

Source code in cylindra/widgets/sta.py

@nogui
@do_not_record
def get_subtomograms(
    self,
    layers: str | MoleculesLayer | list[str | MoleculesLayer],
    shape: tuple[nm, nm, nm],
    bin_size: int = 1,
    order: int = 3,
) -> "Array":
    """
    A non-GUI method to get all the subtomograms as a dask array.

    Parameters
    ----------
    layers : str, MoleculesLayer or list of them
        All the layers that will be used to construct the subtomogram array.
    shape : (nm, nm, nm)
        Shape of output subtomograms.
    bin_size : int, default
        Bin size of the subtomograms.
    order : int, default 3
        Interpolation order.

    Returns
    -------
    Array
        4D Dask array.
    """
    layers = assert_list_of_layers(layers, self.parent_viewer)
    parent = self._get_main()
    tomo = parent.tomogram
    loader = tomo.get_subtomogram_loader(
        _concat_molecules(layers), shape, binsize=bin_size, order=order
    )
    return loader.construct_dask()

get_template(template_path, scale=None)

A non-GUI method to get the template

Source code in cylindra/widgets/sta.py

@nogui
@do_not_record
def get_template(
    self, template_path: str | Path, scale: float | None = None
) -> ip.ImgArray:
    """A non-GUI method to get the template"""
    if scale is None:
        scale = self._get_main().tomogram.scale
    img = self.params._norm_template_param(template_path).provide(scale)
    return ip.asarray(img, axes="zyx").set_scale(zyx=scale)

normalize_landscape(landscape_layer, norm_sd=True)

Normalize the landscape.

Parameters:

Name	Type	Description	Default
landscape_layer	LandscapeSurface	Landscape layer to be used in this algorithm.	required
norm_sd	bool	If true, each landscape will also be normalized by its standard deviation.	True

Source code in cylindra/widgets/sta.py

@set_design(text=capitalize, location=LandscapeMenu)
def normalize_landscape(
    self,
    landscape_layer: _LandscapeLayer,
    norm_sd: bool = True,
):
    """
    Normalize the landscape.

    Parameters
    ----------
    {landscape_layer}
    norm_sd : bool, default True
        If true, each landscape will also be normalized by its standard deviation.
    """
    landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
    new = landscape_layer.landscape.normed(sd=norm_sd)
    surf = LandscapeSurface(new, name=f"{landscape_layer}-Norm")
    return self._add_new_landscape_layer(landscape_layer, surf)

remove_landscape_outliers(landscape_layer, lower=None, upper=None)

Remove outliers from the landscape.

This method will replace energy (inverse score) outliers with the thresholds. This method is useful for lattice with such as defects or strong artifacts.

Parameters:

Name	Type	Description	Default
landscape_layer	LandscapeSurface	Landscape layer to be used in this algorithm.	required
lower	float	Lower limit of the energy.	None
upper	float	Upper limit of the energy.	None

Source code in cylindra/widgets/sta.py

@set_design(text=capitalize, location=LandscapeMenu)
def remove_landscape_outliers(
    self,
    landscape_layer: _LandscapeLayer,
    lower: Annotated[Optional[float], {"text": "Do not process lower outliers"}] = None,
    upper: Annotated[Optional[float], {"text": "Do not process upper outliers"}] = None,
):  # fmt: skip
    """
    Remove outliers from the landscape.

    This method will replace energy (inverse score) outliers with the thresholds.
    This method is useful for lattice with such as defects or strong artifacts.

    Parameters
    ----------
    {landscape_layer}
    lower : float, optional
        Lower limit of the energy.
    upper : float, optional
        Upper limit of the energy.
    """
    landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
    new = landscape_layer.landscape.clip_energies(lower, upper)
    surf = LandscapeSurface(new, name=f"{landscape_layer}-Clip")
    return self._add_new_landscape_layer(landscape_layer, surf)

run_align_on_landscape(landscape_layer)

Find the optimal displacement for each molecule on the landscape.

Source code in cylindra/widgets/sta.py

@set_design(text="Run alignment on landscape", location=LandscapeMenu)
@dask_worker.with_progress(desc="Peak detection on landscape")
def run_align_on_landscape(self, landscape_layer: _LandscapeLayer):
    """Find the optimal displacement for each molecule on the landscape."""
    landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
    landscape = landscape_layer.landscape
    spl = landscape_layer.source_spline
    mole_opt, _ = landscape.run_min_energy(spl)
    return self._align_on_landscape_on_return.with_args(
        mole_opt, landscape_layer.name, spl
    )

run_annealing_on_landscape(landscape_layer, range_long=(4.0, 4.28), range_lat=(5.1, 5.3), angle_max=5.0, temperature_time_const=1.0, random_seeds=(0, 1, 2, 3, 4))

Run simulated annealing on the landscape.

Parameters:

Name	Type	Description	Default
landscape_layer	LandscapeSurface	Landscape layer to be used in this algorithm.	required
range_long	(float, float)	Minimum and maximum allowed distances between longitudinally consecutive monomers	(4.0, 4.28)
range_lat	(float, float)	Minimum and maximum allowed distances between laterally consecutive monomers	(5.1, 5.3)
angle_max	float	Maximum allowed angle between longitudinally consecutive monomers and the Y axis.	5.0
temperature_time_const	float	Time constant of the temperature decay during annealing. Larger value results in slower annealing. 1.0 is a moderate value.	1.0
random_seeds	iterable of int	Random seed integers. Number of integers will be the number of trials.	(0, 1, 2, 3, 4)

Source code in cylindra/widgets/sta.py

@set_design(text="Run annealing on landscape", location=LandscapeMenu)
@dask_worker.with_progress(desc="Running simulated annealing")
def run_annealing_on_landscape(
    self,
    landscape_layer: _LandscapeLayer,
    range_long: _DistRangeLon = (4.0, 4.28),
    range_lat: _DistRangeLat = (5.1, 5.3),
    angle_max: _AngleMaxLon = 5.0,
    temperature_time_const: Annotated[float, {"min": 0.01, "max": 10.0}] = 1.0,
    random_seeds: _RandomSeeds = (0, 1, 2, 3, 4),
):
    """
    Run simulated annealing on the landscape.

    Parameters
    ----------
    {landscape_layer}{range_long}{range_lat}{angle_max}{temperature_time_const}
    {random_seeds}
    """
    t0 = timer()
    landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
    spl = landscape_layer.source_spline
    if spl is None:
        raise ValueError("RMA requires a spline.")
    mole, results = landscape_layer.landscape.run_annealing_along_spline(
        spl=spl,
        range_long=range_long,
        range_lat=range_lat,
        angle_max=angle_max,
        temperature_time_const=temperature_time_const,
        random_seeds=random_seeds,
    )
    t0.toc()

    @thread_worker.callback
    def _plot_result():
        with _Logger.set_plt():
            _annealing.plot_annealing_result(results)

    yield _plot_result
    return self._align_on_landscape_on_return.with_args(
        mole,
        landscape_layer.name,
        source=spl,
        metadata={ANNEALING_RESULT: results[0]},
    )

run_viterbi_on_landscape(landscape_layer, range_long=(4.0, 4.28), angle_max=5.0)

Run Viterbi alignment on the landscape.

Parameters:

Name	Type	Description	Default
landscape_layer	LandscapeSurface	Landscape layer to be used in this algorithm.	required
range_long	(float, float)	Minimum and maximum allowed distances between longitudinally consecutive monomers	(4.0, 4.28)
angle_max	float	Maximum allowed angle between longitudinally consecutive monomers and the Y axis.	5.0

Source code in cylindra/widgets/sta.py

@set_design(text="Run Viterbi alignment on landscape", location=LandscapeMenu)
@dask_worker.with_progress(desc="Running Viterbi alignment")
def run_viterbi_on_landscape(
    self,
    landscape_layer: _LandscapeLayer,
    range_long: _DistRangeLon = (4.0, 4.28),
    angle_max: _AngleMaxLon = 5.0,
):
    """
    Run Viterbi alignment on the landscape.

    Parameters
    ----------
    {landscape_layer}{range_long}{angle_max}
    """
    t0 = timer()
    landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
    spl = landscape_layer.source_spline
    mole = landscape_layer.landscape.run_viterbi_along_spline(
        spl=spl,
        range_long=range_long,
        angle_max=angle_max,
    )
    t0.toc()
    return self._align_on_landscape_on_return.with_args(
        mole, landscape_layer.name, spl
    )

save_annealing_scores(layer, path)

Save RMA scores to a CSV file.

Source code in cylindra/widgets/sta.py

@set_design(text=capitalize, location=Alignment)
@do_not_record
def save_annealing_scores(
    self,
    layer: Annotated[MoleculesLayer, {"choices": _get_layers_with_annealing_result}],
    path: Path.Save[FileFilter.CSV],
):  # fmt: skip
    """Save RMA scores to a CSV file."""
    layer = assert_layer(layer, self.parent_viewer)
    try:
        result: AnnealingResult = layer.metadata[ANNEALING_RESULT]
    except KeyError:
        raise ValueError(
            f"Layer {layer!r} does not have annealing result."
        ) from None
    x = result.batch_size * np.arange(result.energies.size)
    df = pl.DataFrame({"iteration": x, "score": -result.energies})
    return df.write_csv(path, include_header=False)

save_last_average(path)

Save the lastly generated average image.

Source code in cylindra/widgets/sta.py

@set_design(text="Save last average", location=STAnalysis)
def save_last_average(self, path: Path.Save[FileFilter.IMAGE]):
    """Save the lastly generated average image."""
    path = Path(path)
    img = self.last_average
    if img is None:
        raise ValueError(
            "No average image is available. You have to average subtomograms first."
        )
    return img.imsave(path)

seam_search(layer, template_path, mask_params, anti_template_path=None, interpolation=3, npf=None, show_average='Filtered', cutoff=0.25)

Search for the best seam position.

Try all patterns of seam positions and compare cross correlation values. If molecule assembly has 13 protofilaments, this method will try 26 patterns.

Parameters:

Name	Type	Description	Default
layer	MoleculesLayer	Points layer of molecules to be used.	required
template_path	Path or str	Path to template image.	required
mask_params	str or (float, float)	Mask image path or dilation/Gaussian blur parameters. If a path is given, image must in the same shape as the template.	required
anti_template_path	Path	The anti-template used for seam search. For instance, if the template is beta-tubulin, the anti-template is alpha-tubulin.	None
interpolation	int	Interpolation order.	3
npf	int	Number of protofilaments. By default the global properties stored in the corresponding spline will be used.	None
show_average	bool	If true, all the subtomogram averages will be shown.	True
cutoff	float	Cutoff frequency of low-pass filter applied in each subtomogram.	0.25

Source code in cylindra/widgets/sta.py

@set_design(text="Seam search by correlation", location=STAnalysis.SeamSearch)
@dask_worker.with_progress(desc=_pdesc.fmt_layer("Seam search of {!r}"))
def seam_search(
    self,
    layer: MoleculesLayerType,
    template_path: Annotated[_PathOrNone, {"bind": _template_param}],
    mask_params: Annotated[Any, {"bind": _get_mask_params}],
    anti_template_path: Annotated[Optional[Path.Read[FileFilter.IMAGE]], {"text": "Do not use anti-template", "label": "anti-template path"}] = None,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
    npf: Annotated[Optional[int], {"text": "use global properties"}] = None,
    show_average: Annotated[str, {"label": "show averages as", "choices": [None, "Raw", "Filtered"]}] = "Filtered",
    cutoff: _CutoffFreq = 0.25,
):  # fmt: skip
    """
    Search for the best seam position.

    Try all patterns of seam positions and compare cross correlation values. If
    molecule assembly has 13 protofilaments, this method will try 26 patterns.

    Parameters
    ----------
    {layer}{template_path}{mask_params}
    anti_template_path : Path, optional
        The anti-template used for seam search. For instance, if the template is
        beta-tubulin, the anti-template is alpha-tubulin.
    {interpolation}
    npf : int, optional
        Number of protofilaments. By default the global properties stored in the
        corresponding spline will be used.
    show_average : bool, default True
        If true, all the subtomogram averages will be shown.
    {cutoff}
    """
    t0 = timer()
    layer = assert_layer(layer, self.parent_viewer)
    loader, npf = self._seam_search_input(layer, npf, interpolation)
    template, mask = loader.normalize_input(
        template=self.params._norm_template_param(template_path),
        mask=self.params._get_mask(params=mask_params),
    )
    if anti_template_path is not None:
        anti_template = ip.asarray(
            pipe.from_file(anti_template_path).provide(loader.scale),
            axes="zyx",
        )
    else:
        anti_template = None

    seam_searcher = CorrelationSeamSearcher(npf)
    result = seam_searcher.search(
        loader=loader,
        template=ip.asarray(template, axes="zyx"),
        anti_template=anti_template,
        mask=mask,
        cutoff=cutoff,
    )

    t0.toc()

    @thread_worker.callback
    def _seam_search_on_return():
        new_feat = result.as_series(loader.molecules.count())
        layer.features = layer.molecules.features.with_columns(new_feat)
        layer.metadata[SEAM_SEARCH_RESULT] = result
        if show_average is not None:
            if show_average == "Filtered":
                sigma = 0.25 / loader.scale
                result.averages.gaussian_filter(sigma=sigma, update=True)
            _imlayer: "Image" = self._show_rec(
                result.averages, layer.name, store=False
            )
            _imlayer.metadata[SEAM_SEARCH_RESULT] = result

        # plot all the correlation
        _Logger.print_html("<code>seam_search</code>")
        with _Logger.set_plt():
            _Logger.print(f"layer = {layer.name!r}")
            _Logger.print(f"template = {Path(template_path).as_posix()!r}")
            if anti_template_path is not None:
                _Logger.print(
                    f"anti_template = {Path(anti_template_path).as_posix()!r}"
                )
            widget_utils.plot_seam_search_result(result.scores, npf)

    return _seam_search_on_return

seam_search_by_feature(layer, by)

Search for seams by a feature.

Parameters:

Name	Type	Description	Default
layer	MoleculesLayer	Points layer of molecules to be used.	required
by	str	Name of the feature that will be used for seam search.	required

Source code in cylindra/widgets/sta.py

@set_design(text=capitalize, location=STAnalysis.SeamSearch)
def seam_search_by_feature(
    self,
    layer: MoleculesLayerType,
    by: Annotated[str, {"choices": _choice_getter("seam_search_by_feature")}],
):
    """
    Search for seams by a feature.

    Parameters
    ----------
    {layer}
    by : str
        Name of the feature that will be used for seam search.
    """
    layer = assert_layer(layer, self.parent_viewer)
    feat = layer.molecules.features
    if by not in feat.columns:
        raise ValueError(f"Column {by} does not exist.")
    npf = utils.roundint(layer.molecules.features[Mole.pf].max() + 1)
    seam_searcher = BooleanSeamSearcher(npf)
    result = seam_searcher.search(feat[by])
    new_feat = result.as_series(feat.shape[0])
    layer.features = layer.molecules.features.with_columns(new_feat)
    return undo_callback(layer.feature_setter(feat, layer.colormap_info))

seam_search_manually(layer, location=0)

Search for seams manually.

Seam location is represented by a number in the range [0, 2 * npf - 1].

Parameters:

Name	Type	Description	Default
layer	MoleculesLayer	Points layer of molecules to be used.	required
location	int	Seam location.	0

Source code in cylindra/widgets/sta.py

@set_design(text=capitalize, location=STAnalysis.SeamSearch)
def seam_search_manually(
    self,
    layer: MoleculesLayerType,
    location: int = 0,
):
    """
    Search for seams manually.

    Seam location is represented by a number in the range [0, 2 * npf - 1].

    Parameters
    ----------
    {layer}
    location : int
        Seam location.
    """
    layer = assert_layer(layer, self.parent_viewer)
    feat = layer.molecules.features
    npf = utils.roundint(layer.molecules.features[Mole.pf].max() + 1)
    seam_searcher = ManualSeamSearcher(npf)
    result = seam_searcher.search(location)
    new_feat = result.as_series(feat.shape[0])
    layer.features = layer.molecules.features.with_columns(new_feat)
    return undo_callback(layer.feature_setter(feat, layer.colormap_info))

show_mask()

Load and show mask image in the scale of the tomogram.

Source code in cylindra/widgets/sta.py

@set_design(icon="fluent:shape-organic-20-filled", location=STATools)
@do_not_record
def show_mask(self):
    """Load and show mask image in the scale of the tomogram."""
    mask = self._get_mask_image(self._template_params())
    self._show_rec(mask, name="Mask image", store=False, threshold=0.5)

show_template()

Load and show template image in the scale of the tomogram.

Source code in cylindra/widgets/sta.py

@set_design(icon="ic:baseline-view-in-ar", location=STATools)
@do_not_record
def show_template(self):
    """Load and show template image in the scale of the tomogram."""
    template = self._get_template_image()
    self._show_rec(template, name="Template image", store=False)

show_template_original()

Load and show template image in the original scale.

Source code in cylindra/widgets/sta.py

@set_design(icon="material-symbols:view-in-ar", location=STATools)
@do_not_record
def show_template_original(self):
    """Load and show template image in the original scale."""
    _input = self._template_params()
    if _input is None:
        raise ValueError("No template path provided.")
    elif isinstance(_input, Path):
        self._show_rec(ip.imread(_input), name="Template image", store=False)
    else:
        for i, fp in enumerate(_input):
            img = ip.imread(fp)
            self._show_rec(img, name=f"Template image [{i}]", store=False)

split_and_average(layers, n_pairs=1, size=None, interpolation=1, bin_size=1)

Split molecules into two groups and average separately.

Parameters:

Name	Type	Description	Default
layers	list of MoleculesLayer	All the points layers of molecules to be used.	required
n_pairs	int	How many pairs of average will be calculated.	1
size	nm	Size of the template in nm. Use the size of template image by default.	None
interpolation	int	Interpolation order.	1
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1

Source code in cylindra/widgets/sta.py

@set_design(text="Split and average molecules", location=Averaging)
@dask_worker.with_progress(desc=_pdesc.fmt_layers("Split-and-averaging of {!r}"))  # fmt: skip
def split_and_average(
    self,
    layers: MoleculesLayersType,
    n_pairs: Annotated[int, {"min": 1, "label": "number of image pairs"}] = 1,
    size: _SubVolumeSize = None,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 1,
    bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
):
    """
    Split molecules into two groups and average separately.

    Parameters
    ----------
    {layers}
    n_pairs : int, default 1
        How many pairs of average will be calculated.
    {size}{interpolation}{bin_size}
    """
    t0 = timer()
    layers = assert_list_of_layers(layers, self.parent_viewer)
    parent = self._get_main()
    molecules = _concat_molecules(layers)
    shape = self._get_shape_in_nm(size)
    loader = parent.tomogram.get_subtomogram_loader(
        molecules, shape, binsize=bin_size, order=interpolation
    )
    axes = "ipzyx" if n_pairs > 1 else "pzyx"
    img = ip.asarray(loader.average_split(n_set=n_pairs), axes=axes)
    img.set_scale(zyx=loader.scale)
    t0.toc()
    return self._show_rec.with_args(img, f"[Split]{_avg_name(layers)}", store=False)