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="pzyx")
        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[tuple[str, int]]:
        return self._get_main()._get_available_binsize()

    @set_design(text="Extract subtomograms", location=Averaging)
    @dask_worker.with_progress(desc=_pdesc.fmt_layers("Extract subtomograms of {!r}"))
    def extract_subtomograms(
        self,
        layers: MoleculesLayersType,
        size: _SubVolumeSize = None,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 1,
        bin_size: BinSizeType = 1,
        compute: bool = False,
    ):
        """Extract subtomograms around molecules in the selected layer(s).

        Extracted 3D images will be sent to the sub-viewer. Note that this extraction
        operation is NOT needed for subsequent averaging/alignment operations, as they
        internally extract subtomograms by themselves.

        Parameters
        ----------
        {layers}{size}{interpolation}{bin_size}
        compute : bool, default: False
            If True, the subtomograms will be computed immediately, and in-memory data
            will be sent to the sub-viewer.
        """
        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
        )
        subtomos = ip.lazy.asarray(loader.construct_dask(), axes="pzyx")
        if compute:
            subtomos = subtomos.compute()
        return self._show_rec.with_args(
            subtomos, f"[Subtomos]{_avg_name(layers)}", store=False
        )

    @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: BinSizeType = 1,
        wedge_weighting: bool = True,
    ):
        """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}{wedge_weighting}
        """
        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.order_optimize().average(wedge_weighting=wedge_weighting), axes="zyx"
        )
        img.set_scale(zyx=loader.scale, unit="nm")
        t0.toc()
        _Logger.print_html(
            f"{loader.molecules.count()} molecules. Image size: {shape[0]:.2f} nm ({img.shape[0]} pixel)"
        )
        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: BinSizeType = 1,
        wedge_weighting: bool = True,
    ):
        """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}{wedge_weighting}
        """
        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
        ).order_optimize()
        img = ip.asarray(
            loader.average(wedge_weighting=wedge_weighting), 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: BinSizeType = 1,
        wedge_weighting: bool = True,
    ):
        """Group-wise subtomogram averaging.

        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}{wedge_weighting}
        """
        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_polars_expr(by)
        groups = loader.order_optimize().groupby(expr)
        avg_dict = groups.average(wedge_weighting=wedge_weighting)
        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()
        mole_counts = [sub.molecules.count() for _, sub in groups]
        _Logger.print_html(
            f"Averages of {len(avg_dict)} groups, {mole_counts} molecules "
            f"respectively.\nImage size: {shape[0]:.2f} nm ({img.shape[-1]} pixel)"
        )
        return self._show_rec.with_args(img, f"[AVG]{_avg_name(layers)}")

    @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: BinSizeType = 1,
        wedge_weighting: bool = True,
    ):
        """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}{wedge_weighting}
        """
        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
        ).filter(widget_utils.norm_polars_expr(predicate))
        avg = loader.order_optimize().average(wedge_weighting=wedge_weighting)
        img = ip.asarray(avg, axes="zyx")
        img.set_scale(zyx=loader.scale, unit="nm")
        t0.toc()
        _Logger.print_html(
            f"Average of {loader.molecules.count()} molecules. Image size: {shape[0]:.2f} nm ({img.shape[0]} pixel)"
        )
        return self._show_rec.with_args(img, f"[AVG]{_avg_name(layers)}")

    @set_design(text="Align average to template", location=Alignment)
    @dask_worker.with_progress()
    def align_averaged(
        self,
        layers: MoleculesLayersType,
        template_path: Annotated[_PathOrNone, {"bind": _template_param}],
        mask_params: Annotated[Any, {"bind": None}] = None,  # deprecated
        max_shifts: Optional[_MaxShifts] = None,
        rotations: _Rotations = ((0.0, 0.0), (15.0, 1.0), (3.0, 1.0)),
        bin_size: BinSizeType = 1,
        method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
        cylindrical_mask: Annotated[bool, {"text": "Use cylindrical mask"}] = False,
    ):  # fmt: skip
        """Align the averaged image at current molecules to the template image.

        This function creates a new layer with transformed molecules, which should
        align well with the template image. Users usually run this function after using
        `map_monomers` on microtubules etc, because the initial coordinates have
        correct periodicity but the centers may not align well with the actual
        structure.

        Parameters
        ----------
        {layers}{template_path}{max_shifts}{rotations}{bin_size}{method}
        cylindrical_mask : bool, default: False
            If True, apply a cylindrical mask to the averaged image before alignment.
            This option is useful if the averaged image contains a lot of background
            density that may interfere with the alignment (e.g., MAPs decoration, in
            situ data, etc).
        """
        t0 = timer()
        layers = assert_list_of_layers(layers, self.parent_viewer)
        parent = self._get_main()

        new_layers = list[MoleculesLayer]()
        total = 2 * len(layers) + 1
        yield thread_worker.description(_pdesc.align_averaged_0(total))

        @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 = loader.normalize_template(
            self.params._norm_template_param(template_path, allow_multiple=False)
        )
        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,
            rotations=rotations,
            tilt=None,  # NOTE: because input is an average
        )
        _spl_globs = list[
            tuple[weakref.ReferenceType["CylSpline"], pl.DataFrame, pl.DataFrame]
        ]()
        for i, layer in enumerate(layers):
            mole = layer.molecules
            loader = self._get_loader(bin_size, mole, order=1)
            yield thread_worker.description(_pdesc.align_averaged_1(i, total, layer))
            avg = loader.average(template.shape)
            yield thread_worker.description(_pdesc.align_averaged_2(i, total, layer))
            max_shifts_px = [_s / _scale for _s in max_shifts]

            # apply cylindrical mask
            if cylindrical_mask:
                if spl := layer.source_spline:
                    _cyl_mask = _image_gen.make_cylindrical_mask(
                        spl, mole, avg.shape, loader.scale
                    )
                    avg *= _cyl_mask
                else:
                    _Logger.print(
                        f"Layer {layer.name} has no source spline, cannot apply "
                        "cylindrical mask."
                    )

            _img_trans, result = model.fit(avg, max_shifts=max_shifts_px)

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

            # write offsets to spline globalprops if available
            if spl := layer.source_spline:
                _Logger.print(
                    f"Layer {layer.name} has a source spline, update its globalprops."
                )
                _mole_trans = widget_utils.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, 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
            _Logger.print_table(
                _align_averages_table(svec, rotator.as_rotvec()),
                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.with_desc("Finished")

    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: BinSizeType = 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 = widget_utils.MoleculesCombiner()

        loader = self._get_loader(
            binsize=bin_size,
            molecules=combiner.concat(layer.molecules for layer in layers),
            order=interpolation,
        )
        _Logger.print(f"Aligning {loader.molecules.count()} molecules ...")
        aligned_loader = (
            loader.order_optimize()
            .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,
            )
            .order_restore()
        )
        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()
    def align_all_template_free(
        self,
        layers: MoleculesLayersType,
        mask_params: Annotated[Any, {"bind": _get_mask_params}],
        size: _SubVolumeSize = 12.0,
        max_shifts: _MaxShifts = (0.8, 0.8, 0.8),
        max_rotations: _MaxRotations = (3.0, 3.0, 3.0),
        min_rotation_step: Annotated[float, {"min": 0.1, "step": 0.1}] = 0.4,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
        method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
        bin_size: BinSizeType = 1,
        max_num_iters: Annotated[int, {"min": 3, "max": 100}] = 20,
    ):  # fmt: skip
        """Iteratively align molecules and validate using FSC without template.

        Parameters
        ----------
        {layers}{mask_params}{size}{max_shifts}{max_rotations}{min_rotation_step}
        {interpolation}{method}{bin_size}{max_num_iters}
        """
        t0 = timer()
        rng = np.random.default_rng(5926)
        layers = assert_list_of_layers(layers, self.parent_viewer)
        main = self._get_main()
        combiner = widget_utils.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)
            )
        loader = (
            self._get_loader(binsize=bin_size, molecules=molecules, order=interpolation)
            .reshape(shape=shape)
            .order_optimize()
        )
        _alignment_state = template_free.AlignmentState(
            rng=rng,
            mask=mask,
            alignment_model=_get_alignment(method),
            min_rotation_step=min_rotation_step,
        )
        _Logger.print(f"Starting template-free alignment ({loader.molecules.count()} molecules in total) ...")  # fmt: skip
        yield thread_worker.description(_pdesc.align_tf_0(_alignment_state))
        result = _alignment_state.fsc_step_init(loader, max_shifts, max_rotations)
        _Logger.print(f"Iteration {int(_alignment_state.num_iter)}")
        yield _plot_current_fsc.with_args(result.fsc, result.avg).with_desc(_pdesc.align_tf_1(_alignment_state))  # fmt: skip
        while True:
            _Logger.print(_alignment_state.next_params(loader.scale).format())
            _exceeded = max_num_iters <= _alignment_state.num_iter
            if _alignment_state.is_converged() or _exceeded:
                _Logger.print(
                    "Maximum iteration exceeded" if _exceeded else "FSC converged."
                )
                break
            loader = _alignment_state.align_step(loader)
            yield thread_worker.description(_pdesc.align_tf_0(_alignment_state))
            result = _alignment_state.fsc_step(loader)
            _Logger.print(f"Iteration {int(_alignment_state.num_iter)}")
            yield _plot_current_fsc.with_args(result.fsc, result.avg).with_desc(_pdesc.align_tf_1(_alignment_state))  # fmt: skip

        yield self._show_rec.with_args(result.avg, f"[Aligned]{_avg_name(layers)}")
        loader = loader.order_restore()
        molecules = combiner.split(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: BinSizeType = 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 molecules 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="RMA alignment", location=Alignment)
    @dask_worker.with_progress(descs=_pdesc.align_annealing_fmt)
    def align_all_rma(
        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 = ("d.mean() - 0.1", "d.mean() + 0.1"),
        angle_max: _AngleMaxLon = 5.0,
        bin_size: BinSizeType = 1,
        temperature_time_const: Annotated[float, {"min": 0.01, "max": 10.0}] = 1.0,
        lj_const: Annotated[float, {"min": 0.0, "max": 1000.0, "step": 0.1, "label": "LJ const"}] = 0.0,
        method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
        upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 5,
        num_trials: Annotated[int, {"min": 1, "max": 100}] = 5,
        seed: _SeedType = 0,
    ):  # fmt: skip
        """2D-constrained subtomogram alignment using Restricted Mesh 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}{lj_const}{method}{upsample_factor}{num_trials}{seed}
        """
        t0 = timer()
        layer = assert_layer(layer, self.parent_viewer)
        if (spl := layer.source_spline) is None:
            spl = CylSpline._dummy_spline()
            _Logger.print(
                "No source spline found, cylindrical boundary will not be considered."
            )
        main = self._get_main()
        _Logger.print(
            f"Constructing correlation landscape on {layer.name} ({layer.molecules.count()} molecules) for RMA ..."
        )
        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, method=method,
        )  # fmt: skip
        yield
        mole, results = landscape.run_annealing_along_spline(
            spl,
            range_long=range_long,
            range_lat=range_lat,
            angle_max=angle_max,
            temperature_time_const=temperature_time_const,
            lj_nstd=None if lj_const < 1e-3 else 1 / lj_const,
            random_seeds=utils.create_random_seeds(num_trials, seed),
        )
        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="MT-RMA alignment", location=Alignment)
    @dask_worker.with_progress(descs=_pdesc.align_annealing_fmt)
    def align_all_mt_rma(
        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 = ("d.mean() - 0.1", "d.mean() + 0.1"),
        range_angle: _AngleRange = (-1.5, 5.0),
        bin_size: BinSizeType = 1,
        temperature_time_const: Annotated[float, {"min": 0.01, "max": 10.0}] = 1.0,
        lj_const: Annotated[float, {"min": 0.0, "max": 1000.0, "step": 0.1, "label": "LJ const"}] = 0.0,
        method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
        upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 5,
        num_trials: Annotated[int, {"min": 1, "max": 100}] = 5,
        seed: _SeedType = 0,
    ):  # fmt: skip
        """RMA with microtubule-specialized energy potential.

        Parameters
        ----------
        {layer}{template_path}{mask_params}{max_shifts}{rotations}{cutoff}
        {interpolation}{range_long}{range_lat}
        range_angle : tuple[float, float], default: (-1.5, 5.0)
            The allowed range of protofilament bending angles, in degrees.
        {bin_size}{temperature_time_const}{lj_const}{method}{upsample_factor}
        {num_trials}{seed}
        """
        t0 = timer()
        layer = assert_layer(layer, self.parent_viewer)
        if (spl := layer.source_spline) is None:
            spl = CylSpline._dummy_spline()
            _Logger.print(
                "No source spline found, cylindrical boundary will not be considered."
            )
        main = self._get_main()
        _Logger.print(
            f"Constructing correlation landscape on {layer.name} ({layer.molecules.count()} molecules) for MT-RMA ..."
        )
        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, method=method,
            upsample_factor=upsample_factor,
        )  # fmt: skip
        yield
        mole, results = landscape.run_microtubule_annealing_along_spline(
            spl,
            range_long=range_long,
            range_lat=range_lat,
            range_angle=range_angle,
            temperature_time_const=temperature_time_const,
            lj_nstd=None if lj_const < 1e-3 else 1 / lj_const,
            random_seeds=utils.create_random_seeds(num_trials, seed),
        )
        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="RFA alignment", location=Alignment)
    @dask_worker.with_progress(descs=_pdesc.align_annealing_fmt)
    def align_all_rfa(
        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: BinSizeType = 1,
        temperature_time_const: Annotated[float, {"min": 0.01, "max": 10.0}] = 1.0,
        upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 5,
        num_trials: Annotated[int, {"min": 1, "max": 100}] = 5,
        seed: _SeedType = 0,
    ):
        """1D-constrained subtomogram alignment using Restricted Filament Annealing.

        This alignment method considers the distance between every adjacent monomers on
        the filament.

        Parameters
        ----------
        {layer}{template_path}{mask_params}{max_shifts}{rotations}{cutoff}
        {interpolation}{range_long}{angle_max}{bin_size}{temperature_time_const}
        {upsample_factor}{num_trials}{seed}
        """
        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()
        _Logger.print(
            f"Constructing correlation landscape on {layer.name} ({layer.molecules.count()} molecules) for RFA ..."
        )
        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,
        )  # fmt: skip
        yield
        mole, results = landscape.run_filamentous_annealing(
            range=range_long,
            angle_max=angle_max,
            temperature_time_const=temperature_time_const,
            random_seeds=utils.create_random_seeds(num_trials, seed),
        )
        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="Fit spline by RFA", location=Alignment)
    @dask_worker.with_progress(descs=_pdesc.fit_spline_rfa_fmt)
    def fit_spline_rfa(
        self,
        spline: SplineType,
        template_path: Annotated[_PathOrPathsOrNone, {"bind": _template_params}],
        forward_is: PlusOrMinusType = "MinusToPlus",
        interval: PolarsExprStrOrScalar = "4.1",
        err_max: Annotated[nm, {"label": "Max fit error (nm)", "step": 0.1}] = 0.5,
        mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
        max_shifts: _MaxShifts = (2.0, 2.0, 2.0),
        rotations: _Rotations = ((0.0, 0.0), (15.0, 5.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: BinSizeType = 1,
        temperature_time_const: Annotated[float, {"min": 0.01, "max": 100.0}] = 10.0,
        upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 5,
        num_trials: Annotated[int, {"min": 1, "max": 100}] = 5,
        seed: _SeedType = 0,
    ):
        """Fit spline by RFA.

        This algorithm uses a template image to precisely determine the center line of
        filaments. By comparing the score, the orientation of the filament will also be
        determined.

        Parameters
        ----------
        {spline}{template_path}
        forward_is : "PlusToMinus" or "MinusToPlus", default "MinusToPlus"
            Which orientation is the forward direction. Set "PlusToMinus" if the
            template image is oriented from the plus end to the minus end in the y
            direction.
        interval : float or str expression, default 4.1
            Interval of the sampling points along the spline.
        {err_max}{mask_params}{max_shifts}{rotations}{cutoff}{interpolation}{range_long}
        {angle_max}{bin_size}{temperature_time_const}{upsample_factor}{num_trials}{seed}
        """
        t0 = timer()
        main = self._get_main()
        tomo = main.tomogram
        spl_old = tomo.splines[spline]
        interv_expr = widget_utils.norm_scalar_expr(interval)
        mole_fw = tomo.map_centers(
            i=spline,
            interval=spl_old.props.get_glob(interv_expr),
            rotate_molecules=False,
        )
        mole_rv = mole_fw.rotate_by_rotvec_internal([np.pi, 0, 0])  # invert

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

        def _construct_landscape(mole_):
            landscape_, _ = self._construct_landscape(
                molecules=mole_, template_path=template_path, mask_params=mask_params,
                max_shifts=max_shifts, rotations=rotations, cutoff=cutoff, norm=False,
                order=interpolation, bin_size=bin_size, upsample_factor=upsample_factor,
            )  # fmt: skip
            yield
            mole, results = landscape_.run_filamentous_annealing(
                range=range_long,
                angle_max=angle_max,
                temperature_time_const=temperature_time_const,
                random_seeds=utils.create_random_seeds(num_trials, seed),
            )
            yield _on_yield.with_args(results)
            return mole, results

        mole_opt_fw, results_fw = yield from _construct_landscape(mole_fw)
        mole_opt_rv, results_rv = yield from _construct_landscape(mole_rv)

        if results_fw[0].energies[-1] < results_rv[0].energies[-1]:  # forward is better
            mole_opt = mole_opt_fw
            ori = Ori(forward_is)
        else:  # reverse is better
            mole_opt = mole_opt_rv
            ori = Ori.invert(Ori(forward_is))

        # calculate distances for logging
        _ds = np.diff(mole_opt.pos, axis=0)
        _dist: NDArray[np.float32] = np.sqrt(np.sum(_ds**2, axis=1))

        spl_new = CylSpline(
            order=spl_old.order,
            config=spl_old.config,
            extrapolate=spl_old.extrapolate,
        ).fit(mole_opt.pos, err_max=err_max)
        spl_new.props.loc = spl_old.props.loc
        spl_new.props.glob = spl_old.props.glob
        spl_new.orientation = ori
        tomo.splines[spline] = spl_new
        t0.toc()

        @thread_worker.callback
        def _on_return():
            main._update_splines_in_images()
            _Logger.print(f"Orientation: {ori.value}")
            _Logger.print(
                f"Distance: mean = {_dist.mean():.3f} nm (ranging from "
                f"{_dist.min():.3f} to {_dist.max():.3f}) nm "
            )
            main.reset_choices()

            @undo_callback
            def _out():
                tomo.splines[spline] = spl_old
                main._update_splines_in_images()

            @_out.with_redo
            def _out():
                tomo.splines[spline] = spl_new
                main._update_splines_in_images()

            return _out

        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: BinSizeType = 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, norm=norm,
            upsample_factor=upsample_factor, method=method,
        )  # fmt: skip
        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", location=LandscapeMenu)
    @dask_worker.with_progress(desc="Peak detection on landscape")
    def run_align_on_landscape(self, landscape_layer: LandscapeLayerType):
        """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", location=LandscapeMenu)
    @dask_worker.with_progress(desc="Running Viterbi alignment")
    def run_viterbi_on_landscape(
        self,
        landscape_layer: LandscapeLayerType,
        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 (RMA)", location=LandscapeMenu)
    @dask_worker.with_progress(desc="Running simulated annealing")
    def run_rma_on_landscape(
        self,
        landscape_layer: LandscapeLayerType,
        range_long: _DistRangeLon = (4.0, 4.28),
        range_lat: _DistRangeLat = ("d.mean() - 0.1", "d.mean() + 0.1"),
        angle_max: _AngleMaxLon = 5.0,
        temperature_time_const: Annotated[float, {"min": 0.01, "max": 10.0}] = 1.0,
        lj_const: Annotated[float, {"min": 0.00, "max": 1000.0, "step": 0.1, "label": "LJ const"}] = 0.0,
        num_trials: Annotated[int, {"min": 1, "max": 100}] = 5,
        seed: _SeedType = 0,
    ):  # fmt: skip
        """Run simulated annealing on the landscape, supposing a cylindric structure.

        Parameters
        ----------
        {landscape_layer}{range_long}{range_lat}{angle_max}{temperature_time_const}
        {lj_const}{num_trials}{seed}
        """
        t0 = timer()
        landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
        if (spl := landscape_layer.source_spline) is None:
            spl = CylSpline._dummy_spline()
            _Logger.print(
                "No source spline found, cylindrical boundary will not be considered."
            )

        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,
            lj_nstd=None if lj_const < 1e-3 else 1 / lj_const,
            random_seeds=utils.create_random_seeds(num_trials, seed),
        )
        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="Run annealing (RFA)", location=LandscapeMenu)
    @dask_worker.with_progress(desc="Running simulated annealing")
    def run_rfa_on_landscape(
        self,
        landscape_layer: LandscapeLayerType,
        range_long: _DistRangeLon = (4.0, 4.28),
        angle_max: _AngleMaxLon = 5.0,
        temperature_time_const: Annotated[float, {"min": 0.01, "max": 10.0}] = 1.0,
        num_trials: Annotated[int, {"min": 1, "max": 100}] = 5,
        seed: _SeedType = 0,
    ):
        """Run simulated annealing on the landscape, supposing a filamentous structure.

        Parameters
        ----------
        {landscape_layer}{range_long}{angle_max}{temperature_time_const}{num_trials}
        {seed}
        """
        t0 = timer()
        landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
        mole, results = landscape_layer.landscape.run_filamentous_annealing(
            range=range_long,
            angle_max=angle_max,
            temperature_time_const=temperature_time_const,
            random_seeds=utils.create_random_seeds(num_trials, seed),
        )
        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=landscape_layer.source_spline,
            metadata={ANNEALING_RESULT: results[0]},
        )

    @set_design(text=capitalize, location=LandscapeMenu)
    def remove_landscape_outliers(
        self,
        landscape_layer: LandscapeLayerType,
        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: LandscapeLayerType,
        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

    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.epoch_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: BinSizeType = 1,
        method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
        column_prefix: str = "score",
        use_tilt_series: bool = False,
    ):
        """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 = widget_utils.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)

        score_kwargs = {
            "templates": [pipe.from_file(t) for t in template_path],
            "mask": mask,
            "cutoff": cutoff,
        }
        if use_tilt_series:
            out = (
                main.tomogram.tilt_series.loader.replace(
                    molecules=all_mole, order=interpolation, binsize=bin_size
                )
                .score_on_tilt(**score_kwargs, method=method)
                .mean(axis=2)
            )
        else:
            out = main.tomogram.get_subtomogram_loader(
                all_mole,
                order=interpolation,
                binsize=bin_size,
            ).score(
                **score_kwargs,
                tilt=main.tomogram.tilt_model,
                alignment_model=_get_alignment(method),
            )
        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: _SeedType = 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,
        wedge_weighting: bool = True,
    ):
        """Calculate Fourier Shell Correlation using the selected 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.
        {wedge_weighting}
        """
        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,
            )
            .order_optimize()
            .fsc_with_halfmaps(
                mask,
                seed=seed,
                n_set=n_pairs,
                dfreq=dfreq,
                squeeze=False,
                wedge_weighting=wedge_weighting,
            )
        )

        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)
                plt.tight_layout()
                plt.show()
            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="3D classification", location=STAnalysis)
    @dask_worker.with_progress(desc="3D Classification")
    def classify_em(
        self,
        layers: MoleculesLayersType,
        templates: Annotated[list[Path], {"filter": FileFilter.IMAGE}],
        max_num_iters: Annotated[int, {"min": 1, "max": 100, "label": "iterations"}] = 10,
        min_num_molecules_per_class: _MinNumMoleculesPerClass = None,
        mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
        inverse_temperature: Annotated[float, {"min": 1, "max": 10000}] = 80,
        cutoff: _CutoffFreq = 0.5,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
        bin_size: BinSizeType = 1,
        seed: _SeedType = 0,
    ):  # fmt: skip
        """Classify molecules to the user given templates using EM algorithm.

        Parameters
        ----------
        {layers}{max_num_iters}{mask_params}{cutoff}{interpolation}
        {bin_size}
        templates : list of Path
            Paths to the template images used for classification. Must be a list of
            3D images or a 4D image containing multiple 3D images.
        inverse_temperature : float, default 100
            Temperature parameter for the softmax function used in the expectation step.
            Higher value (lower temperature) makes the classification harder, while
            lower value (higher temperature) makes it "softer" (i.e., more uncertain).
        """
        t0 = timer()
        layers = assert_list_of_layers(layers, self.parent_viewer)
        mask = self.params._get_mask(params=mask_params)
        tomo = self._get_main().tomogram
        if isinstance(templates, (Path, str)):
            templates = [Path(templates)]
        template_images: list[np.ndarray] = []
        for _path in templates:
            arr = ip.imread(_path)
            ratio = arr.scale.x / (tomo.scale * bin_size)
            if abs(ratio - 1) > 0.01:
                arr = arr.zoom(ratio, order=3, mode="reflect")
            if arr.ndim == 3:
                template_images.append(arr.value)
            elif arr.ndim == 4:
                template_images.extend(list(arr.value))
            else:
                raise ValueError(
                    f"Template image {_path!r} must be a 3D or 4D image, "
                    f"but got {arr.ndim}D image."
                )

        combiner = widget_utils.MoleculesCombiner()
        all_moles = combiner.concat(layer.molecules for layer in layers)

        shape = template_images[0].shape
        loader = (
            self._get_loader(
                binsize=bin_size,
                molecules=all_moles,
                order=interpolation,
            )
            .replace(output_shape=shape)
            .order_optimize()
        )
        _Logger.print(
            f"Starting EM classification into user supplied {len(template_images)} templates."
            f"\n{loader.molecules.count()} subtomograms with size {loader.output_shape}"
        )
        if min_num_molecules_per_class is None:
            min_num_molecules_per_class = all_moles.count()
        picker = ClassfySubsetPicker(
            max_num_iters,
            min_num_molecules_per_class * len(template_images),
            np.random.default_rng(seed),
        )
        yield thread_worker.description("Extracting subtomograms")

        with tempfile.TemporaryDirectory() as tmpdir:
            path = Path(tmpdir) / f"subtomograms_{id(t0)}"
            temp_loader = loader.extract_subtomograms(path)
            result = yield from _classify_em_impl(
                template_images, mask, temp_loader, max_num_iters=max_num_iters,
                inverse_temperature=inverse_temperature, cutoff=cutoff, picker=picker,
            )  # fmt: skip
        t0.toc()
        yield thread_worker.description("Classification finished")
        avgs, all_moles = _post_classify_em(result, loader, all_moles)
        moles = combiner.split(all_moles, layers)

        @thread_worker.callback
        def _on_return():
            for layer, each_mole in zip(layers, moles, strict=True):
                layer.set_molecules_with_new_features(each_mole)
            self._show_rec(avgs, name="Classification results", store=False)

        return _on_return

    @set_design(text="3D classification (template-free)", location=STAnalysis)
    @dask_worker.with_progress(desc="3D Classification (template-free)")
    def classify_em_template_free(
        self,
        layers: MoleculesLayersType,
        size: Annotated[float, {"value": 12.0, "max": 100.0, "label": "size (nm)"}] = None,
        num_classes: Annotated[int, {"min": 1, "max": 100}] = 2,
        min_num_molecules_per_class: _MinNumMoleculesPerClass = None,
        max_num_iters: Annotated[int, {"min": 1, "max": 100, "label": "iterations"}] = 10,
        mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
        inverse_temperature: Annotated[float, {"min": 1, "max": 10000}] = 80,
        cutoff: _CutoffFreq = 0.5,
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
        bin_size: BinSizeType = 1,
        seed: _SeedType = 0,
    ):  # fmt: skip
        """Classify molecules using EM algorithm.

        Random split averages will be used as initial templates.

        Parameters
        ----------
        {layers}{size}{max_num_iters}{mask_params}{cutoff}{interpolation}{bin_size}
        num_classes : int, default 2
            The number of classes to classify subtomograms into.
        inverse_temperature : float, default 100
            Temperature parameter for the softmax function used in the expectation step.
            Higher value (lower temperature) makes the classification harder, while
            lower value (higher temperature) makes it "softer" (i.e., more uncertain).
        seed : int, optional
            Random seed.
        """
        t0 = timer()
        layers = assert_list_of_layers(layers, self.parent_viewer)
        combiner = widget_utils.MoleculesCombiner()
        mask = self.params._get_mask(params=mask_params)
        all_moles = combiner.concat(layer.molecules for layer in layers)

        shape = (size,) * 3
        loader = self._get_loader(
            binsize=bin_size, molecules=all_moles, order=interpolation, shape=shape
        ).order_optimize()
        _Logger.print(
            f"Starting EM classification into {num_classes} classes.\n"
            f"{loader.molecules.count()} subtomograms with size {loader.output_shape}"
        )
        rng = np.random.default_rng(seed)
        if min_num_molecules_per_class is None:
            min_num_molecules_per_class = all_moles.count()
        picker = ClassfySubsetPicker(
            max_num_iters,
            min_num_molecules_per_class * num_classes,
            rng,
        )
        yield thread_worker.description("Extracting subtomograms")

        with tempfile.TemporaryDirectory() as tmpdir:
            path = Path(tmpdir) / f"subtomograms_{id(t0)}"
            temp_loader = loader.extract_subtomograms(path)
            yield thread_worker.description("Creating initial class templates")
            templates = temp_loader.average_split(n_split=num_classes, seed=rng)
            result = yield from _classify_em_impl(
                templates, mask, temp_loader, max_num_iters=max_num_iters,
                inverse_temperature=inverse_temperature, cutoff=cutoff, picker=picker,
            )  # fmt: skip

        t0.toc()
        yield thread_worker.description("Classification finished")
        avgs, all_moles = _post_classify_em(result, loader, all_moles)
        moles = combiner.split(all_moles, layers)

        @thread_worker.callback
        def _on_return():
            for layer, each_mole in zip(layers, moles, strict=True):
                layer.set_molecules_with_new_features(each_mole)
            self._show_rec(avgs, name="Classification results", store=False)

        return _on_return

    @set_design(text="Calculate S/N profiles", location=STAnalysis)
    @dask_worker.with_progress(desc="Preparing averaged image")
    def calculate_signal_noise_profiles(
        self,
        layer: MoleculesLayerType,
        diameter: Annotated[nm, {"min": 0.1, "max": 100.0}] = 6.0,
        max_shifts: tuple[nm, nm, nm] = (0.5, 0.5, 0.5),
        rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
        bin_size: BinSizeType = 1,
        upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 2,
        method: Annotated[str, {"choices": METHOD_CHOICES}] = "ncc",
        shift_to_background: Annotated[Optional[tuple[nm, nm, nm]], {"text": "Use half tomograms"}] = None,
    ):  # fmt: skip
        """Calculate the mean landscape energy of signal and noise.

        Parameters
        ----------
        {layer}
        diameter : float, default 6.0
            Diameter of the monomer.
        {max_shifts}{rotations}{interpolation}{bin_size}{upsample_factor}{method}
        shift_to_background : tuple of float, optional
            If given, the molecules will be internally shifted to get the background
            region. If not given, the background will be calculated using half1 - half2.
        """
        main = self._get_main()
        _layer = assert_layer(layer, main.parent_viewer)
        mole = _layer.molecules
        val_s, val_n = yield from self._calculate_signal_noise(
            _layer, diameter, max_shifts, rotations, cutoff=1.0,
            interpolation=interpolation, bin_size=bin_size,
            upsample_factor=upsample_factor, method=method,
            shift_to_background=shift_to_background,
        )  # fmt: skip

        @thread_worker.callback
        def _on_return():
            _layer.molecules = mole.with_features(
                pl.Series("profile_signal", val_s, dtype=pl.Float32),
                pl.Series("profile_noise", val_n, dtype=pl.Float32),
            )
            _Logger.print(f"Profile of {_layer.name!r} has been calculated.")
            with _Logger.set_plt():
                widget_utils.plot_signal_noise_profile(_layer.molecules.features)

        return _on_return

    @set_design(text="Detect defects", location=STAnalysis)
    @dask_worker.with_progress(desc="Preparing averaged image")
    def detect_defects(
        self,
        layer: MoleculesLayerType,
        diameter: Annotated[nm, {"min": 0.1, "max": 100.0}] = 5.8,
        prior_defect: Annotated[float, {"min": 0.001, "max": 0.999}] = 0.1,
        max_shifts: tuple[nm, nm, nm] = (0.5, 0.5, 0.5),
        rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
        interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
        bin_size: BinSizeType = 1,
        upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 2,
        method: Annotated[str, {"choices": METHOD_CHOICES}] = "ncc",
        shift_to_background: Annotated[Optional[tuple[nm, nm, nm]], {"text": "Use half tomograms"}] = None,
    ):  # fmt: skip
        """Detect lattice defect by Bayesian metrics.

        This method calculates the signal and noise distribution, and evaluate the
        probability of being a defect site for individual molecules.

        Parameters
        ----------
        {layer}
        diameter : float
            Diameter of monomer in nanometers.
        prior_defect : float
            Prior probability of a site being a defect. This value does not necessarily
            reflect the actual defect rate. Changing this value only rescales the output
            probability.
        {max_shifts}{rotations}{interpolation}{bin_size}{upsample_factor}{method}
        shift_to_background : tuple of float, optional
            Shift vector to move molecules to background. For example, use (7, 0, 0) for
            tubulin subunits.
        """
        main = self._get_main()
        _layer = assert_layer(layer, main.parent_viewer)
        mole = _layer.molecules
        val_s, val_n = yield from self._calculate_signal_noise(
            _layer, diameter, max_shifts, rotations, cutoff=1.0,
            interpolation=interpolation, bin_size=bin_size,
            upsample_factor=upsample_factor, method=method,
            shift_to_background=shift_to_background,
        )  # fmt: skip

        mean_s, std_s = val_s.mean(), val_s.std(ddof=0)
        mean_n, std_n = val_n.mean(), val_n.std(ddof=0)
        _log_mean_sd_for_detect_defect(mean_s, std_s, mean_n, std_n)

        post_prob = utils.calc_post_probability(
            val_s, 1 - prior_defect, mean_s, std_s, mean_n, std_n
        )

        prob_defect = 1 - post_prob

        @thread_worker.callback
        def _on_return():
            _layer.molecules = mole.with_features(
                pl.Series("profile_signal", val_s, dtype=pl.Float32),
                pl.Series("profile_noise", val_n, dtype=pl.Float32),
                pl.Series("prob_defect", prob_defect, dtype=pl.Float32),
            )
            if not isinstance(base_color := _layer.colormap_info, str):
                base_color = "lime"
            _layer.set_colormap("prob_defect", (0, 1), cmap={0: base_color, 1: "black"})

        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: Path.Read[FileFilter.IMAGE],
        mask_params: Annotated[Any, {"bind": _get_mask_params}],
        anti_template_path: Annotated[Path.Read[FileFilter.IMAGE] | None, {"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 between
        template(s) and averages. If molecule assembly has 13 protofilaments, this
        method will try 26 patterns of averages.

        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 : str, default "Filtered"
            If "Raw", all the raw subtomogram averages will be shown.
            If "Filtered", all the filtered subtomogram averages will be shown.
            If None, no subtomogram averages will be shown.
        {cutoff}
        """
        t0 = timer()
        layer = assert_layer(layer, self.parent_viewer)
        loader, npf = prep_seam_search_input(
            self._get_main().tomogram,
            layer.molecules,
            npf=npf,
            order=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 and str(anti_template_path) not in ("", "."):
            anti_template, anti_mask = loader.normalize_input(
                template=self.params._norm_template_param(anti_template_path),
                mask=self.params._get_mask(params=mask_params),
            )
            anti_template = ip.asarray(anti_template, axes="zyx")
            if mask is not None and anti_mask is not None:
                mask = np.maximum(mask, anti_mask)
        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.lower() == "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="Seam search on tilt series", location=STAnalysis.SeamSearch)
    @dask_worker.with_progress(desc=_pdesc.fmt_layer("Seam search of {!r} on tilt series"))  # fmt: skip
    def seam_search_on_tilt(
        self,
        layer: MoleculesLayerType,
        template_path: Path.Read[FileFilter.IMAGE],
        mask_params: Annotated[Any, {"bind": _get_mask_params}],
        anti_template_path: Annotated[Path.Read[FileFilter.IMAGE] | None, {"label": "anti-template path"}] = None,
        npf: Annotated[Optional[int], {"text": "use global properties"}] = None,
        bin_size: BinSizeType = 1,
        cutoff: _CutoffFreq = 0.25,
        ctf: Annotated[bool, {"label": "apply CTF"}] = True,
    ):  # fmt: skip
        """Search for the best seam position using the scores in the tilt images.

        The template (and the anti-template) is projected to each tilt image at the
        position and orientation of each molecule, and compared with the tilt images
        by ZNCC. The score of each molecule is the mean over the tilts. The class means
        of the scores are compared for all the seam positions. The tomogram must be
        opened from a tilt series.

        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.
        npf : int, optional
            Number of protofilaments. By default the global properties stored in the
            corresponding spline will be used.
        {bin_size}{cutoff}
        ctf : bool, default True
            If true, the CTF is applied to the projections of the templates, with the
            defocus of each molecule in each tilt, considering its depth. The tilt
            series must have the CTF parameters.
        """
        t0 = timer()
        layer = assert_layer(layer, self.parent_viewer)
        if (src := self._get_main().tomogram.tilt_series) is None:
            raise ValueError("The tomogram is not opened from a tilt series.")
        mole = layer.molecules
        if npf is None:
            npf = mole.features[Mole.pf].unique().len()
        loader = src.get_loader(mole, binsize=bin_size)
        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 and str(anti_template_path) not in ("", "."):
            anti_template, anti_mask = loader.normalize_input(
                template=self.params._norm_template_param(anti_template_path),
                mask=self.params._get_mask(params=mask_params),
            )
            # the same mask is used for both templates
            if mask is not None and anti_mask is not None:
                mask = np.maximum(mask, anti_mask)
        else:
            anti_template = None

        seam_searcher = TiltScoreSeamSearcher(npf)
        result = seam_searcher.search(
            loader=loader,
            template=template,
            anti_template=anti_template,
            mask=mask,
            cutoff=cutoff,
            ctf=ctf,
        )

        t0.toc()

        @thread_worker.callback
        def _seam_search_on_return():
            new_feat = result.as_series(mole.count())
            layer.features = layer.molecules.features.with_columns(new_feat)
            layer.metadata[SEAM_SEARCH_RESULT] = result

            with _Logger.set_plt():
                _Logger.print(f"layer = {layer.name!r}")
                _Logger.print(f"template = {Path(template_path).as_posix()!r}")
                if anti_template is not None:
                    _Logger.print(
                        f"anti_template = {Path(anti_template_path).as_posix()!r}"
                    )
                _Logger.print(f"seam position = {result.seam_pos}")
                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.
        """
        ui = self._get_main()
        layer = assert_layer(layer, ui.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)
        ui.reset_choices()  # choices regarding of features need update
        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.
        """
        ui = self._get_main()
        layer = assert_layer(layer, ui.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)
        ui.reset_choices()  # choices regarding of features need update
        return undo_callback(layer.feature_setter(feat, layer.colormap_info))

    @set_design(text=capitalize, location=Alignment)
    @dask_worker.with_progress(desc="Refining the tilt series")
    def refine_tilt_series(
        self,
        layers: MoleculesLayersType,
        size: _SubVolumeSize = None,
        max_shift: Annotated[nm, {"label": "max shift (nm)", "min": 0.1, "max": 50.0, "step": 0.1}] = 2.0,
        rotation: Annotated[tuple[float, float], {"label": "rotation max/step (deg)"}] = (0, 0),
        grid_shape: Annotated[tuple[int, int], {"label": "grid shape (y, x)", "options": {"min": 1, "max": 20}}] = (1, 1),
        smoothness: Annotated[float, {"min": 0.0, "max": 100.0, "step": 0.001}] = 0.01,
        ctf: Annotated[bool, {"label": "apply CTF"}] = True,
        bin_size: BinSizeType = 1,
        update_coarse_tomogram: bool = True,
    ):  # fmt: skip
        """Refine the alignment of the tilt series using the molecules.

        The subtomograms of each layer are averaged as the reference of the molecules
        of the layer. Projections of the averages placed at the molecules are compared
        with the tilt images, and the shifts of each tilt image are refined.

        Parameters
        ----------
        {layers}{size}
        max_shift : nm
            Maximum shift of the tilt images in nm.
        rotation : (float, float)
            Maximum rotation and the step of the tilt images in degrees.
        grid_shape : (int, int)
            Number of control points of the local shifts along y and x of the tilt
            images. (1, 1) means a uniform shift of each tilt image. Otherwise, the
            shifts are refined as a smooth field over each tilt image.
        smoothness : float
            Weight of the smoothness prior of the local shifts. Larger value is more
            robust against noise and fewer molecules.
        ctf : bool, default True
            If true, the CTF is applied to the projections of the averages, with the
            defocus of each molecule in each tilt, considering its depth. The tilt
            series must have the CTF parameters.
        {bin_size}
        update_coarse_tomogram : bool, default True
            If true, the coarse tomogram for visualization is also reconstructed from
            the refined tilt series.
        """
        t0 = timer()
        layers = assert_list_of_layers(layers, self.parent_viewer)
        main = self._get_main()
        tomo = main.tomogram
        if (src := tomo.tilt_series) is None:
            raise ValueError("The tomogram is not opened from a tilt series.")
        shape = self._get_shape_in_nm(size)
        refs = list[tuple[Molecules, "NDArray[np.float32]"]]()
        for layer in layers:
            # the average of the current tilt series is the reference of each layer
            mole = layer.molecules
            loader = tomo.get_subtomogram_loader(mole, shape, binsize=bin_size)
            yield thread_worker.description(f"Creating template for layer {layer.name}")
            refs.append((mole, loader.average(wedge_weighting=True)))
        nmole = sum(mole.count() for mole, _ in refs)
        _Logger.print(f"Refining the tilt series using {nmole} molecules ...")
        old = src.tilt_series
        loader = src.get_loader(Molecules.empty(), binsize=bin_size)
        yield thread_worker.description("Running tilt series refinement")
        new = loader.refine_tilt_series(
            refs,
            max_shift=max_shift,
            grid_shape=grid_shape,
            smoothness=smoothness,
            rotation=rotation,
            ctf=ctf,
        ).tilt_series
        if update_coarse_tomogram:
            yield thread_worker.description("Reconstructing coarse tomogram")
            new_ref = src.reconstruct_reference(new)
        else:
            new_ref = None
        t0.toc()

        old_ref = src.reference

        def _apply(tilt_series: "TiltSeries", reference: "ip.ImgArray | None"):
            removed = tomo.update_tilt_series(tilt_series)
            if new_ref is not None and reference is not None:
                src.set_reference(reference)
                self._show_coarse_tomogram(reference)
            if main.SplineControl.need_resample:
                main.sample_subtomograms()
            return removed

        @thread_worker.callback
        def _on_return():
            removed = _apply(new, new_ref)
            _log_tilt_series_shifts(old, new)
            if removed:
                _Logger.print(
                    f"Multi-scale images of bin size {removed} are removed. They are "
                    "reconstructed from the refined tilt series on demand."
                )
            # the tilt series and the coarse tomogram are restored by undo
            return (
                undo_callback(_apply)
                .with_args(old, old_ref)
                .with_redo(lambda: _apply(new, new_ref))
            )

        return _on_return

    @set_design(text=capitalize, location=Alignment)
    @thread_worker.with_progress(desc="Aligning the tilt series")
    def align_tilt_series(
        self,
        splines: SplinesType = None,
        interval: Annotated[nm, {"label": "segment interval (nm)", "min": 2.0, "max": 200.0, "step": 1.0}] = 25.0,
        length: Annotated[nm, {"label": "segment length (nm)", "min": 4.0, "max": 200.0, "step": 1.0}] = 32.0,
        max_shift: Annotated[nm, {"label": "max shift (nm)", "min": 1.0, "max": 200.0, "step": 1.0}] = 20.0,
        resolution: Annotated[nm, {"label": "low-pass resolution (nm)", "min": 0.5, "max": 100.0, "step": 0.5}] = 6.0,
        position_sigma: Annotated[Optional[nm], {"label": "position sigma (nm)", "text": "do not move splines", "options": {"min": 0.1, "max": 100.0, "step": 0.5, "value": 10.0}}] = 10.0,
        offset_per: Literal["spline", "segment"] = "segment",
        fit_tilt_axis: bool = True,
        max_tilt_angle_change: Annotated[float, {"label": "max tilt angle change (deg)", "min": 0.0, "max": 10.0, "step": 0.1}] = 0.5,
        progressive: bool = True,
        niter: Annotated[int, {"min": 1, "max": 20}] = 3,
        ctf: Annotated[bool, {"label": "apply CTF"}] = False,
        bin_size: Annotated[int, {"min": 1, "max": 16}] = 4,
        update_splines: bool = True,
        update_coarse_tomogram: bool = True,
    ):  # fmt: skip
        """Coarsely align the tilt series of large errors using the splines.

        Segments of tubes are placed along the splines, and their projections are
        compared with the tilt images. Call `refine_tilt_series` for further refinement.

        Parameters
        ----------
        {splines}
        interval : nm
            Interval between the centers of the segments along the splines.
        length : nm
            Length of each segment (tube) along the spline. The ends of the tubes are
            smoothly tapered.
        max_shift : nm
            Maximum shift of the segments in the tilt images from the predicted
            positions.
        resolution : nm
            Resolution of the low-pass filter applied to the projections and the tilt
            images. It should be coarser than the lattice spacing (~4 nm), because the
            ends of the segments correlate with the lattice.
        position_sigma : nm, optional
            Standard deviation of the prior of the 3D offsets of the splines. The
            splines picked in the tomogram of a wrong alignment are misplaced,
            especially along the z axis. If not given, the splines are not moved.
        offset_per : "spline" or "segment"
            The segments of the same spline share the 3D offset if "spline", otherwise
            each segment has its own offset.
        fit_tilt_axis : bool
            If true, the tilt axis angle common to all the tilts is also optimized.
        max_tilt_angle_change : float
            If positive, the tilt angles are also optimized within this range (deg)
            around the current ones. The mean change is fixed to zero, as it is a
            rotation of the tomogram.
        progressive : bool
            If true, the tilts are added one by one from the zero tilt, so that the
            errors accumulating over the tilts can be larger than `max_shift`.
        niter : int
            Maximum number of iterations of computing the correlation and optimizing
            the parameters.
        ctf : bool, default True
            If true, the CTF is applied to the projections of the segments, with the
            defocus of each segment in each tilt, considering its depth. The tilt
            series must have the CTF parameters.
        bin_size : int
            Bin size of the tilt images used for the alignment. It does not have to be
            a bin size of the multi-scale images.
        update_splines : bool, default True
            If true, the splines are moved by their offsets, so that they are at the
            filaments in the tomogram of the aligned tilt series. If the offsets are
            of each segment, the splines are fitted to the moved segments again. The
            components of the offsets along the splines are ignored.
        update_coarse_tomogram : bool, default True
            If true, the coarse tomogram for visualization is also reconstructed from
            the aligned tilt series.
        """
        t0 = timer()
        main = self._get_main()
        tomo = main.tomogram
        if (src := tomo.tilt_series) is None:
            raise ValueError("The tomogram is not opened from a tilt series.")
        indices = main._norm_splines(splines)
        if len(indices) == 0:
            raise ValueError("No spline is given.")
        loader = src.get_loader(Molecules.empty(), binsize=bin_size)
        refs = list[tuple[Molecules, "NDArray[np.float32]"]]()
        for i in indices:
            spl = tomo.splines[i]
            if spl.radius is None:
                raise ValueError(f"The radius of spline {i} is not measured yet.")
            mole = _spline_segments(spl, interval, length)
            refs.append(
                (mole, _image_gen.make_tube_template(spl, length, loader.scale))
            )
        nseg = sum(mole.count() for mole, _ in refs)
        _Logger.print(
            f"Aligning the tilt series using {nseg} segments of {len(indices)} "
            "splines ..."
        )
        old = src.tilt_series
        result = loader.align_tilt_series(
            refs,
            max_shift=max_shift,
            cutoff=loader.scale / resolution,  # cycles per pixel
            group_by=_SEGMENT_GROUP if offset_per == "spline" else None,
            position_sigma=position_sigma,
            tilt_axis="global" if fit_tilt_axis else "fixed",
            max_tilt_angle_change=max_tilt_angle_change,
            progressive=progressive,
            niter=niter,
            ctf=ctf,
        )
        new = result.loader.tilt_series
        offsets = [
            moved.pos - mole.pos
            for (mole, _), moved in zip(refs, result.molecules, strict=True)
        ]
        new_splines = dict[int, CylSpline]()
        if update_splines and position_sigma is not None:
            for i, (mole, _), offset in zip(indices, refs, offsets, strict=True):
                spl_old = tomo.splines[i]
                spl_new = _move_spline(spl_old, mole, offset)
                _u = np.linspace(0, 1, max(int(spl_old.length() / 4), 1))
                if spl_new.curvature(_u).max() > spl_old.curvature(_u).max() * 3:
                    _Logger.print(
                        f"Warning: Curvature of spline {i} increased significantly."
                    )
                new_splines[i] = spl_new
        if update_coarse_tomogram:
            yield thread_worker.description("Reconstructing coarse tomogram")
            new_ref = src.reconstruct_reference(new)
        else:
            new_ref = None
        t0.toc()

        old_ref = src.reference
        old_splines = {i: tomo.splines[i] for i in new_splines}

        def _apply(
            tilt_series: "TiltSeries",
            reference: "ip.ImgArray | None",
            spls: dict[int, CylSpline],
        ):
            removed = tomo.update_tilt_series(tilt_series)
            if new_ref is not None and reference is not None:
                src.set_reference(reference)
                self._show_coarse_tomogram(reference)
            for i, spl in spls.items():
                tomo.splines[i] = spl
            if spls:
                main._update_splines_in_images()
            if main.SplineControl.need_resample:
                main.sample_subtomograms()
            return removed

        @thread_worker.callback
        def _on_return():
            removed = _apply(new, new_ref, new_splines)
            _log_tilt_series_alignment(old, new, offsets, result.scores)
            if removed:
                _Logger.print(
                    f"Multi-scale images of bin size {removed} are removed. They are "
                    "reconstructed from the aligned tilt series on demand."
                )
            # the tilt series, the coarse tomogram and the splines are restored by undo
            return (
                undo_callback(_apply)
                .with_args(old, old_ref, old_splines)
                .with_redo(lambda: _apply(new, new_ref, new_splines))
            )

        return _on_return

    def _show_coarse_tomogram(self, img: ip.ImgArray):
        """Show the coarse tomogram with the current filter and inversion."""
        main = self._get_main()
        filt = main._reserved_layers.ref_filter
        inverted = main._reserved_layers.ref_inverted
        main._update_reference_image(-img if inverted else img)
        main._reserved_layers.ref_inverted = inverted
        if filt is not None:
            main.filter_reference_image(*filt)

    @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,
    ) -> tuple[Landscape, SubtomogramLoader]:
        parent = self._get_main()
        loader = parent.tomogram.get_subtomogram_loader(
            molecules, binsize=bin_size, order=order
        )
        model = _get_alignment(method)
        tmp_prov = self.params._norm_template_param(template_path, allow_multiple=True)
        landscape = Landscape.from_loader(
            loader=loader,
            template=tmp_prov.provide(loader.scale),
            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, loader

    def _calculate_signal_noise(
        self,
        layer: MoleculesLayer,
        diameter: nm = 6.0,
        max_shifts: tuple[nm, nm, nm] = (0.5, 0.5, 0.5),
        rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
        cutoff: _CutoffFreq = 0.5,
        interpolation: int = 3,
        bin_size: BinSizeType = 1,
        upsample_factor: int = 2,
        method: str = "ncc",
        shift_to_background: tuple[nm, nm, nm] | None = None,
    ):
        main = self._get_main()
        model = _get_alignment(method)
        mole = layer.molecules
        if shift_to_background is None and main.tomogram.image_halves is None:
            raise ValueError(
                "Tomogram does not have half tomograms. "
                "Please provide shift_to_background parameter."
            )

        loader = main.tomogram.get_subtomogram_loader(
            mole, binsize=bin_size, order=interpolation
        )
        _sigma = 0.8
        shape = (main.tomogram.nm2pixel(diameter + _sigma * 3, bin_size),) * 3
        _zeros = np.zeros(shape, dtype=np.float32)
        mask = spherical_mask(diameter / 2, _sigma)(_zeros, loader.scale * bin_size)
        fsc, (img_0, img_1), img_mask = (
            loader.reshape(shape=shape)
            .order_optimize()
            .fsc_with_halfmaps(mask, seed=0, n_set=1, dfreq=None, squeeze=True)
        )
        avg = (img_0 + img_1) / 2
        result = widget_utils.FscResult.from_dataframe(fsc, loader.scale)
        yield _plot_current_fsc(result, avg)

        kwargs = {
            "template": avg,
            "mask": mask,
            "max_shifts": max_shifts,
            "upsample_factor": upsample_factor,
            "alignment_model": model.with_params(
                rotations=rotations,
                cutoff=cutoff,
                tilt=main.tomogram.tilt_model,
            ),
        }
        if shift_to_background is not None:
            mole_bg = mole.translate_internal(shift_to_background)
            loader_s_n = main.tomogram.get_subtomogram_loader(
                mole.concat_with(mole_bg), binsize=bin_size, order=interpolation
            )
            # it is faster to load subtomograms at the same time
            yield thread_worker.description("Constructing signal and noise landscapes")
            landscape_s_n = Landscape.from_loader(loader_s_n, **kwargs)
            val_s_n = landscape_s_n.energies.mean(axis=(1, 2, 3))
            val_s = val_s_n[: mole.count()]
            val_n = val_s_n[mole.count() :]
        else:
            loaders = main.tomogram.get_subtomogram_and_noise_loader(
                mole, binsize=bin_size, order=interpolation
            )
            yield thread_worker.description("Constructing signal and noise landscapes")
            landscape_s, landscape_n = Landscape.diff_landscapes(*loaders, **kwargs)

            val_s = landscape_s.energies.mean(axis=(1, 2, 3))
            val_n = landscape_n.energies.mean(axis=(1, 2, 3))
        return val_s, val_n

last_average property

Last averaged image if exists.

sub_viewer 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: BinSizeType = 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 = widget_utils.MoleculesCombiner()

    loader = self._get_loader(
        binsize=bin_size,
        molecules=combiner.concat(layer.molecules for layer in layers),
        order=interpolation,
    )
    _Logger.print(f"Aligning {loader.molecules.count()} molecules ...")
    aligned_loader = (
        loader.order_optimize()
        .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,
        )
        .order_restore()
    )
    molecules = combiner.split(aligned_loader.molecules, layers)
    t0.toc()
    return self._align_all_on_return.with_args(molecules, layers)

align_all_mt_rma(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=('d.mean() - 0.1', 'd.mean() + 0.1'), range_angle=(-1.5, 5.0), bin_size=1, temperature_time_const=1.0, lj_const=0.0, method='zncc', upsample_factor=5, num_trials=5, seed=0)

RMA with microtubule-specialized energy potential.

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	('d.mean() - 0.1', 'd.mean() + 0.1')
range_angle	tuple[float, float]	The allowed range of protofilament bending angles, in degrees.	(-1.5, 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
lj_const	float	If greater than 0, the annealing will use a Lennard-Jones-like potential to allow molecules to separate further apart than the cutoff distance. Larger value will prefer longer distances between molecules. 0.25 - 0.5 is usually a good choice.	0.0
method	str	Correlation metrics for alignment.	'zncc'
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
num_trials	int	Number of independent trials to perform. The best result among all trials will be returned.	5
seed	int	Random seed for reproducibility.	0

Source code in cylindra/widgets/sta.py

@set_design(text="MT-RMA alignment", location=Alignment)
@dask_worker.with_progress(descs=_pdesc.align_annealing_fmt)
def align_all_mt_rma(
    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 = ("d.mean() - 0.1", "d.mean() + 0.1"),
    range_angle: _AngleRange = (-1.5, 5.0),
    bin_size: BinSizeType = 1,
    temperature_time_const: Annotated[float, {"min": 0.01, "max": 10.0}] = 1.0,
    lj_const: Annotated[float, {"min": 0.0, "max": 1000.0, "step": 0.1, "label": "LJ const"}] = 0.0,
    method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
    upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 5,
    num_trials: Annotated[int, {"min": 1, "max": 100}] = 5,
    seed: _SeedType = 0,
):  # fmt: skip
    """RMA with microtubule-specialized energy potential.

    Parameters
    ----------
    {layer}{template_path}{mask_params}{max_shifts}{rotations}{cutoff}
    {interpolation}{range_long}{range_lat}
    range_angle : tuple[float, float], default: (-1.5, 5.0)
        The allowed range of protofilament bending angles, in degrees.
    {bin_size}{temperature_time_const}{lj_const}{method}{upsample_factor}
    {num_trials}{seed}
    """
    t0 = timer()
    layer = assert_layer(layer, self.parent_viewer)
    if (spl := layer.source_spline) is None:
        spl = CylSpline._dummy_spline()
        _Logger.print(
            "No source spline found, cylindrical boundary will not be considered."
        )
    main = self._get_main()
    _Logger.print(
        f"Constructing correlation landscape on {layer.name} ({layer.molecules.count()} molecules) for MT-RMA ..."
    )
    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, method=method,
        upsample_factor=upsample_factor,
    )  # fmt: skip
    yield
    mole, results = landscape.run_microtubule_annealing_along_spline(
        spl,
        range_long=range_long,
        range_lat=range_lat,
        range_angle=range_angle,
        temperature_time_const=temperature_time_const,
        lj_nstd=None if lj_const < 1e-3 else 1 / lj_const,
        random_seeds=utils.create_random_seeds(num_trials, seed),
    )
    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_rfa(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, temperature_time_const=1.0, upsample_factor=5, num_trials=5, seed=0)

1D-constrained subtomogram alignment using Restricted Filament Annealing.

This alignment method considers the distance between every adjacent monomers on the filament.

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
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
num_trials	int	Number of independent trials to perform. The best result among all trials will be returned.	5
seed	int	Random seed for reproducibility.	0

Source code in cylindra/widgets/sta.py

@set_design(text="RFA alignment", location=Alignment)
@dask_worker.with_progress(descs=_pdesc.align_annealing_fmt)
def align_all_rfa(
    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: BinSizeType = 1,
    temperature_time_const: Annotated[float, {"min": 0.01, "max": 10.0}] = 1.0,
    upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 5,
    num_trials: Annotated[int, {"min": 1, "max": 100}] = 5,
    seed: _SeedType = 0,
):
    """1D-constrained subtomogram alignment using Restricted Filament Annealing.

    This alignment method considers the distance between every adjacent monomers on
    the filament.

    Parameters
    ----------
    {layer}{template_path}{mask_params}{max_shifts}{rotations}{cutoff}
    {interpolation}{range_long}{angle_max}{bin_size}{temperature_time_const}
    {upsample_factor}{num_trials}{seed}
    """
    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()
    _Logger.print(
        f"Constructing correlation landscape on {layer.name} ({layer.molecules.count()} molecules) for RFA ..."
    )
    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,
    )  # fmt: skip
    yield
    mole, results = landscape.run_filamentous_annealing(
        range=range_long,
        angle_max=angle_max,
        temperature_time_const=temperature_time_const,
        random_seeds=utils.create_random_seeds(num_trials, seed),
    )
    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_rma(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=('d.mean() - 0.1', 'd.mean() + 0.1'), angle_max=5.0, bin_size=1, temperature_time_const=1.0, lj_const=0.0, method='zncc', upsample_factor=5, num_trials=5, seed=0)

2D-constrained subtomogram alignment using Restricted Mesh 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	('d.mean() - 0.1', 'd.mean() + 0.1')
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
lj_const	float	If greater than 0, the annealing will use a Lennard-Jones-like potential to allow molecules to separate further apart than the cutoff distance. Larger value will prefer longer distances between molecules. 0.25 - 0.5 is usually a good choice.	0.0
method	str	Correlation metrics for alignment.	'zncc'
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
num_trials	int	Number of independent trials to perform. The best result among all trials will be returned.	5
seed	int	Random seed for reproducibility.	0

Source code in cylindra/widgets/sta.py

@set_design(text="RMA alignment", location=Alignment)
@dask_worker.with_progress(descs=_pdesc.align_annealing_fmt)
def align_all_rma(
    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 = ("d.mean() - 0.1", "d.mean() + 0.1"),
    angle_max: _AngleMaxLon = 5.0,
    bin_size: BinSizeType = 1,
    temperature_time_const: Annotated[float, {"min": 0.01, "max": 10.0}] = 1.0,
    lj_const: Annotated[float, {"min": 0.0, "max": 1000.0, "step": 0.1, "label": "LJ const"}] = 0.0,
    method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
    upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 5,
    num_trials: Annotated[int, {"min": 1, "max": 100}] = 5,
    seed: _SeedType = 0,
):  # fmt: skip
    """2D-constrained subtomogram alignment using Restricted Mesh 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}{lj_const}{method}{upsample_factor}{num_trials}{seed}
    """
    t0 = timer()
    layer = assert_layer(layer, self.parent_viewer)
    if (spl := layer.source_spline) is None:
        spl = CylSpline._dummy_spline()
        _Logger.print(
            "No source spline found, cylindrical boundary will not be considered."
        )
    main = self._get_main()
    _Logger.print(
        f"Constructing correlation landscape on {layer.name} ({layer.molecules.count()} molecules) for RMA ..."
    )
    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, method=method,
    )  # fmt: skip
    yield
    mole, results = landscape.run_annealing_along_spline(
        spl,
        range_long=range_long,
        range_lat=range_lat,
        angle_max=angle_max,
        temperature_time_const=temperature_time_const,
        lj_nstd=None if lj_const < 1e-3 else 1 / lj_const,
        random_seeds=utils.create_random_seeds(num_trials, seed),
    )
    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=(0.8, 0.8, 0.8), max_rotations=(3.0, 3.0, 3.0), min_rotation_step=0.4, interpolation=3, method='zncc', bin_size=1, max_num_iters=20)

Iteratively align molecules and validate using FSC without template.

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.	(0.8, 0.8, 0.8)
max_rotations	(float, float, float)	Maximum rotation in degree around Z, Y, X axis. Rotations are searched with steps of at most 30 degrees.	(3.0, 3.0, 3.0)
min_rotation_step	float	Minimum rotation search precision in degree.	0.4
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
max_num_iters	int	Maximum number of iterations to perform.	20

Source code in cylindra/widgets/sta.py

@set_design(text="Align all (template-free)", location=Alignment)
@dask_worker.with_progress()
def align_all_template_free(
    self,
    layers: MoleculesLayersType,
    mask_params: Annotated[Any, {"bind": _get_mask_params}],
    size: _SubVolumeSize = 12.0,
    max_shifts: _MaxShifts = (0.8, 0.8, 0.8),
    max_rotations: _MaxRotations = (3.0, 3.0, 3.0),
    min_rotation_step: Annotated[float, {"min": 0.1, "step": 0.1}] = 0.4,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
    method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
    bin_size: BinSizeType = 1,
    max_num_iters: Annotated[int, {"min": 3, "max": 100}] = 20,
):  # fmt: skip
    """Iteratively align molecules and validate using FSC without template.

    Parameters
    ----------
    {layers}{mask_params}{size}{max_shifts}{max_rotations}{min_rotation_step}
    {interpolation}{method}{bin_size}{max_num_iters}
    """
    t0 = timer()
    rng = np.random.default_rng(5926)
    layers = assert_list_of_layers(layers, self.parent_viewer)
    main = self._get_main()
    combiner = widget_utils.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)
        )
    loader = (
        self._get_loader(binsize=bin_size, molecules=molecules, order=interpolation)
        .reshape(shape=shape)
        .order_optimize()
    )
    _alignment_state = template_free.AlignmentState(
        rng=rng,
        mask=mask,
        alignment_model=_get_alignment(method),
        min_rotation_step=min_rotation_step,
    )
    _Logger.print(f"Starting template-free alignment ({loader.molecules.count()} molecules in total) ...")  # fmt: skip
    yield thread_worker.description(_pdesc.align_tf_0(_alignment_state))
    result = _alignment_state.fsc_step_init(loader, max_shifts, max_rotations)
    _Logger.print(f"Iteration {int(_alignment_state.num_iter)}")
    yield _plot_current_fsc.with_args(result.fsc, result.avg).with_desc(_pdesc.align_tf_1(_alignment_state))  # fmt: skip
    while True:
        _Logger.print(_alignment_state.next_params(loader.scale).format())
        _exceeded = max_num_iters <= _alignment_state.num_iter
        if _alignment_state.is_converged() or _exceeded:
            _Logger.print(
                "Maximum iteration exceeded" if _exceeded else "FSC converged."
            )
            break
        loader = _alignment_state.align_step(loader)
        yield thread_worker.description(_pdesc.align_tf_0(_alignment_state))
        result = _alignment_state.fsc_step(loader)
        _Logger.print(f"Iteration {int(_alignment_state.num_iter)}")
        yield _plot_current_fsc.with_args(result.fsc, result.avg).with_desc(_pdesc.align_tf_1(_alignment_state))  # fmt: skip

    yield self._show_rec.with_args(result.avg, f"[Aligned]{_avg_name(layers)}")
    loader = loader.order_restore()
    molecules = combiner.split(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 molecules 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: BinSizeType = 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 molecules 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=None, max_shifts=None, rotations=((0.0, 0.0), (15.0, 1.0), (3.0, 1.0)), bin_size=1, method='zncc', cylindrical_mask=False)

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

This function creates a new layer with transformed molecules, which should align well with the template image. Users usually run this function after using map_monomers on microtubules etc, because the initial coordinates have correct periodicity but the centers may not align well with the actual structure.

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
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'
cylindrical_mask	bool	If True, apply a cylindrical mask to the averaged image before alignment. This option is useful if the averaged image contains a lot of background density that may interfere with the alignment (e.g., MAPs decoration, in situ data, etc).	False

Source code in cylindra/widgets/sta.py

@set_design(text="Align average to template", location=Alignment)
@dask_worker.with_progress()
def align_averaged(
    self,
    layers: MoleculesLayersType,
    template_path: Annotated[_PathOrNone, {"bind": _template_param}],
    mask_params: Annotated[Any, {"bind": None}] = None,  # deprecated
    max_shifts: Optional[_MaxShifts] = None,
    rotations: _Rotations = ((0.0, 0.0), (15.0, 1.0), (3.0, 1.0)),
    bin_size: BinSizeType = 1,
    method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
    cylindrical_mask: Annotated[bool, {"text": "Use cylindrical mask"}] = False,
):  # fmt: skip
    """Align the averaged image at current molecules to the template image.

    This function creates a new layer with transformed molecules, which should
    align well with the template image. Users usually run this function after using
    `map_monomers` on microtubules etc, because the initial coordinates have
    correct periodicity but the centers may not align well with the actual
    structure.

    Parameters
    ----------
    {layers}{template_path}{max_shifts}{rotations}{bin_size}{method}
    cylindrical_mask : bool, default: False
        If True, apply a cylindrical mask to the averaged image before alignment.
        This option is useful if the averaged image contains a lot of background
        density that may interfere with the alignment (e.g., MAPs decoration, in
        situ data, etc).
    """
    t0 = timer()
    layers = assert_list_of_layers(layers, self.parent_viewer)
    parent = self._get_main()

    new_layers = list[MoleculesLayer]()
    total = 2 * len(layers) + 1
    yield thread_worker.description(_pdesc.align_averaged_0(total))

    @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 = loader.normalize_template(
        self.params._norm_template_param(template_path, allow_multiple=False)
    )
    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,
        rotations=rotations,
        tilt=None,  # NOTE: because input is an average
    )
    _spl_globs = list[
        tuple[weakref.ReferenceType["CylSpline"], pl.DataFrame, pl.DataFrame]
    ]()
    for i, layer in enumerate(layers):
        mole = layer.molecules
        loader = self._get_loader(bin_size, mole, order=1)
        yield thread_worker.description(_pdesc.align_averaged_1(i, total, layer))
        avg = loader.average(template.shape)
        yield thread_worker.description(_pdesc.align_averaged_2(i, total, layer))
        max_shifts_px = [_s / _scale for _s in max_shifts]

        # apply cylindrical mask
        if cylindrical_mask:
            if spl := layer.source_spline:
                _cyl_mask = _image_gen.make_cylindrical_mask(
                    spl, mole, avg.shape, loader.scale
                )
                avg *= _cyl_mask
            else:
                _Logger.print(
                    f"Layer {layer.name} has no source spline, cannot apply "
                    "cylindrical mask."
                )

        _img_trans, result = model.fit(avg, max_shifts=max_shifts_px)

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

        # write offsets to spline globalprops if available
        if spl := layer.source_spline:
            _Logger.print(
                f"Layer {layer.name} has a source spline, update its globalprops."
            )
            _mole_trans = widget_utils.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, 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
        _Logger.print_table(
            _align_averages_table(svec, rotator.as_rotvec()),
            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.with_desc("Finished")

align_tilt_series(splines=None, interval=25.0, length=32.0, max_shift=20.0, resolution=6.0, position_sigma=10.0, offset_per='segment', fit_tilt_axis=True, max_tilt_angle_change=0.5, progressive=True, niter=3, ctf=False, bin_size=4, update_splines=True, update_coarse_tomogram=True)

Coarsely align the tilt series of large errors using the splines.

Segments of tubes are placed along the splines, and their projections are compared with the tilt images. Call refine_tilt_series for further refinement.

Parameters:

Name	Type	Description	Default
splines	list of int	Indices of splines to be used.	None
interval	nm	Interval between the centers of the segments along the splines.	25.0
length	nm	Length of each segment (tube) along the spline. The ends of the tubes are smoothly tapered.	32.0
max_shift	nm	Maximum shift of the segments in the tilt images from the predicted positions.	20.0
resolution	nm	Resolution of the low-pass filter applied to the projections and the tilt images. It should be coarser than the lattice spacing (~4 nm), because the ends of the segments correlate with the lattice.	6.0
position_sigma	nm	Standard deviation of the prior of the 3D offsets of the splines. The splines picked in the tomogram of a wrong alignment are misplaced, especially along the z axis. If not given, the splines are not moved.	10.0
offset_per	spline or segment	The segments of the same spline share the 3D offset if "spline", otherwise each segment has its own offset.	'segment'
fit_tilt_axis	bool	If true, the tilt axis angle common to all the tilts is also optimized.	True
max_tilt_angle_change	float	If positive, the tilt angles are also optimized within this range (deg) around the current ones. The mean change is fixed to zero, as it is a rotation of the tomogram.	0.5
progressive	bool	If true, the tilts are added one by one from the zero tilt, so that the errors accumulating over the tilts can be larger than max_shift.	True
niter	int	Maximum number of iterations of computing the correlation and optimizing the parameters.	3
ctf	bool	If true, the CTF is applied to the projections of the segments, with the defocus of each segment in each tilt, considering its depth. The tilt series must have the CTF parameters.	True
bin_size	int	Bin size of the tilt images used for the alignment. It does not have to be a bin size of the multi-scale images.	4
update_splines	bool	If true, the splines are moved by their offsets, so that they are at the filaments in the tomogram of the aligned tilt series. If the offsets are of each segment, the splines are fitted to the moved segments again. The components of the offsets along the splines are ignored.	True
update_coarse_tomogram	bool	If true, the coarse tomogram for visualization is also reconstructed from the aligned tilt series.	True

Source code in cylindra/widgets/sta.py

@set_design(text=capitalize, location=Alignment)
@thread_worker.with_progress(desc="Aligning the tilt series")
def align_tilt_series(
    self,
    splines: SplinesType = None,
    interval: Annotated[nm, {"label": "segment interval (nm)", "min": 2.0, "max": 200.0, "step": 1.0}] = 25.0,
    length: Annotated[nm, {"label": "segment length (nm)", "min": 4.0, "max": 200.0, "step": 1.0}] = 32.0,
    max_shift: Annotated[nm, {"label": "max shift (nm)", "min": 1.0, "max": 200.0, "step": 1.0}] = 20.0,
    resolution: Annotated[nm, {"label": "low-pass resolution (nm)", "min": 0.5, "max": 100.0, "step": 0.5}] = 6.0,
    position_sigma: Annotated[Optional[nm], {"label": "position sigma (nm)", "text": "do not move splines", "options": {"min": 0.1, "max": 100.0, "step": 0.5, "value": 10.0}}] = 10.0,
    offset_per: Literal["spline", "segment"] = "segment",
    fit_tilt_axis: bool = True,
    max_tilt_angle_change: Annotated[float, {"label": "max tilt angle change (deg)", "min": 0.0, "max": 10.0, "step": 0.1}] = 0.5,
    progressive: bool = True,
    niter: Annotated[int, {"min": 1, "max": 20}] = 3,
    ctf: Annotated[bool, {"label": "apply CTF"}] = False,
    bin_size: Annotated[int, {"min": 1, "max": 16}] = 4,
    update_splines: bool = True,
    update_coarse_tomogram: bool = True,
):  # fmt: skip
    """Coarsely align the tilt series of large errors using the splines.

    Segments of tubes are placed along the splines, and their projections are
    compared with the tilt images. Call `refine_tilt_series` for further refinement.

    Parameters
    ----------
    {splines}
    interval : nm
        Interval between the centers of the segments along the splines.
    length : nm
        Length of each segment (tube) along the spline. The ends of the tubes are
        smoothly tapered.
    max_shift : nm
        Maximum shift of the segments in the tilt images from the predicted
        positions.
    resolution : nm
        Resolution of the low-pass filter applied to the projections and the tilt
        images. It should be coarser than the lattice spacing (~4 nm), because the
        ends of the segments correlate with the lattice.
    position_sigma : nm, optional
        Standard deviation of the prior of the 3D offsets of the splines. The
        splines picked in the tomogram of a wrong alignment are misplaced,
        especially along the z axis. If not given, the splines are not moved.
    offset_per : "spline" or "segment"
        The segments of the same spline share the 3D offset if "spline", otherwise
        each segment has its own offset.
    fit_tilt_axis : bool
        If true, the tilt axis angle common to all the tilts is also optimized.
    max_tilt_angle_change : float
        If positive, the tilt angles are also optimized within this range (deg)
        around the current ones. The mean change is fixed to zero, as it is a
        rotation of the tomogram.
    progressive : bool
        If true, the tilts are added one by one from the zero tilt, so that the
        errors accumulating over the tilts can be larger than `max_shift`.
    niter : int
        Maximum number of iterations of computing the correlation and optimizing
        the parameters.
    ctf : bool, default True
        If true, the CTF is applied to the projections of the segments, with the
        defocus of each segment in each tilt, considering its depth. The tilt
        series must have the CTF parameters.
    bin_size : int
        Bin size of the tilt images used for the alignment. It does not have to be
        a bin size of the multi-scale images.
    update_splines : bool, default True
        If true, the splines are moved by their offsets, so that they are at the
        filaments in the tomogram of the aligned tilt series. If the offsets are
        of each segment, the splines are fitted to the moved segments again. The
        components of the offsets along the splines are ignored.
    update_coarse_tomogram : bool, default True
        If true, the coarse tomogram for visualization is also reconstructed from
        the aligned tilt series.
    """
    t0 = timer()
    main = self._get_main()
    tomo = main.tomogram
    if (src := tomo.tilt_series) is None:
        raise ValueError("The tomogram is not opened from a tilt series.")
    indices = main._norm_splines(splines)
    if len(indices) == 0:
        raise ValueError("No spline is given.")
    loader = src.get_loader(Molecules.empty(), binsize=bin_size)
    refs = list[tuple[Molecules, "NDArray[np.float32]"]]()
    for i in indices:
        spl = tomo.splines[i]
        if spl.radius is None:
            raise ValueError(f"The radius of spline {i} is not measured yet.")
        mole = _spline_segments(spl, interval, length)
        refs.append(
            (mole, _image_gen.make_tube_template(spl, length, loader.scale))
        )
    nseg = sum(mole.count() for mole, _ in refs)
    _Logger.print(
        f"Aligning the tilt series using {nseg} segments of {len(indices)} "
        "splines ..."
    )
    old = src.tilt_series
    result = loader.align_tilt_series(
        refs,
        max_shift=max_shift,
        cutoff=loader.scale / resolution,  # cycles per pixel
        group_by=_SEGMENT_GROUP if offset_per == "spline" else None,
        position_sigma=position_sigma,
        tilt_axis="global" if fit_tilt_axis else "fixed",
        max_tilt_angle_change=max_tilt_angle_change,
        progressive=progressive,
        niter=niter,
        ctf=ctf,
    )
    new = result.loader.tilt_series
    offsets = [
        moved.pos - mole.pos
        for (mole, _), moved in zip(refs, result.molecules, strict=True)
    ]
    new_splines = dict[int, CylSpline]()
    if update_splines and position_sigma is not None:
        for i, (mole, _), offset in zip(indices, refs, offsets, strict=True):
            spl_old = tomo.splines[i]
            spl_new = _move_spline(spl_old, mole, offset)
            _u = np.linspace(0, 1, max(int(spl_old.length() / 4), 1))
            if spl_new.curvature(_u).max() > spl_old.curvature(_u).max() * 3:
                _Logger.print(
                    f"Warning: Curvature of spline {i} increased significantly."
                )
            new_splines[i] = spl_new
    if update_coarse_tomogram:
        yield thread_worker.description("Reconstructing coarse tomogram")
        new_ref = src.reconstruct_reference(new)
    else:
        new_ref = None
    t0.toc()

    old_ref = src.reference
    old_splines = {i: tomo.splines[i] for i in new_splines}

    def _apply(
        tilt_series: "TiltSeries",
        reference: "ip.ImgArray | None",
        spls: dict[int, CylSpline],
    ):
        removed = tomo.update_tilt_series(tilt_series)
        if new_ref is not None and reference is not None:
            src.set_reference(reference)
            self._show_coarse_tomogram(reference)
        for i, spl in spls.items():
            tomo.splines[i] = spl
        if spls:
            main._update_splines_in_images()
        if main.SplineControl.need_resample:
            main.sample_subtomograms()
        return removed

    @thread_worker.callback
    def _on_return():
        removed = _apply(new, new_ref, new_splines)
        _log_tilt_series_alignment(old, new, offsets, result.scores)
        if removed:
            _Logger.print(
                f"Multi-scale images of bin size {removed} are removed. They are "
                "reconstructed from the aligned tilt series on demand."
            )
        # the tilt series, the coarse tomogram and the splines are restored by undo
        return (
            undo_callback(_apply)
            .with_args(old, old_ref, old_splines)
            .with_redo(lambda: _apply(new, new_ref, new_splines))
        )

    return _on_return

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

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
wedge_weighting	bool	Whether to apply wedge weighting by Wiener filtering to the average(s).	True

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: BinSizeType = 1,
    wedge_weighting: bool = True,
):
    """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}{wedge_weighting}
    """
    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.order_optimize().average(wedge_weighting=wedge_weighting), axes="zyx"
    )
    img.set_scale(zyx=loader.scale, unit="nm")
    t0.toc()
    _Logger.print_html(
        f"{loader.molecules.count()} molecules. Image size: {shape[0]:.2f} nm ({img.shape[0]} pixel)"
    )
    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, wedge_weighting=True)

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
wedge_weighting	bool	Whether to apply wedge weighting by Wiener filtering to the average(s).	True

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: BinSizeType = 1,
    wedge_weighting: bool = True,
):
    """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}{wedge_weighting}
    """
    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
    ).filter(widget_utils.norm_polars_expr(predicate))
    avg = loader.order_optimize().average(wedge_weighting=wedge_weighting)
    img = ip.asarray(avg, axes="zyx")
    img.set_scale(zyx=loader.scale, unit="nm")
    t0.toc()
    _Logger.print_html(
        f"Average of {loader.molecules.count()} molecules. Image size: {shape[0]:.2f} nm ({img.shape[0]} pixel)"
    )
    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, wedge_weighting=True)

Group-wise subtomogram averaging.

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
wedge_weighting	bool	Whether to apply wedge weighting by Wiener filtering to the average(s).	True

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: BinSizeType = 1,
    wedge_weighting: bool = True,
):
    """Group-wise subtomogram averaging.

    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}{wedge_weighting}
    """
    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_polars_expr(by)
    groups = loader.order_optimize().groupby(expr)
    avg_dict = groups.average(wedge_weighting=wedge_weighting)
    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()
    mole_counts = [sub.molecules.count() for _, sub in groups]
    _Logger.print_html(
        f"Averages of {len(avg_dict)} groups, {mole_counts} molecules "
        f"respectively.\nImage size: {shape[0]:.2f} nm ({img.shape[-1]} pixel)"
    )
    return self._show_rec.with_args(img, f"[AVG]{_avg_name(layers)}")

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

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
wedge_weighting	bool	Whether to apply wedge weighting by Wiener filtering to the average(s).	True

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: BinSizeType = 1,
    wedge_weighting: bool = True,
):
    """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}{wedge_weighting}
    """
    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
    ).order_optimize()
    img = ip.asarray(
        loader.average(wedge_weighting=wedge_weighting), 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', use_tilt_series=False)

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: BinSizeType = 1,
    method: Annotated[str, {"choices": METHOD_CHOICES}] = "zncc",
    column_prefix: str = "score",
    use_tilt_series: bool = False,
):
    """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 = widget_utils.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)

    score_kwargs = {
        "templates": [pipe.from_file(t) for t in template_path],
        "mask": mask,
        "cutoff": cutoff,
    }
    if use_tilt_series:
        out = (
            main.tomogram.tilt_series.loader.replace(
                molecules=all_mole, order=interpolation, binsize=bin_size
            )
            .score_on_tilt(**score_kwargs, method=method)
            .mean(axis=2)
        )
    else:
        out = main.tomogram.get_subtomogram_loader(
            all_mole,
            order=interpolation,
            binsize=bin_size,
        ).score(
            **score_kwargs,
            tilt=main.tomogram.tilt_model,
            alignment_model=_get_alignment(method),
        )
    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, wedge_weighting=True)

Calculate Fourier Shell Correlation using the selected 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
wedge_weighting	bool	Whether to apply wedge weighting by Wiener filtering to the average(s).	True

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: _SeedType = 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,
    wedge_weighting: bool = True,
):
    """Calculate Fourier Shell Correlation using the selected 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.
    {wedge_weighting}
    """
    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,
        )
        .order_optimize()
        .fsc_with_halfmaps(
            mask,
            seed=seed,
            n_set=n_pairs,
            dfreq=dfreq,
            squeeze=False,
            wedge_weighting=wedge_weighting,
        )
    )

    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)
            plt.tight_layout()
            plt.show()
        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

calculate_signal_noise_profiles(layer, diameter=6.0, max_shifts=(0.5, 0.5, 0.5), rotations=((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)), interpolation=3, bin_size=1, upsample_factor=2, method='ncc', shift_to_background=None)

Calculate the mean landscape energy of signal and noise.

Parameters:

Name	Type	Description	Default
layer	MoleculesLayer	Points layer of molecules to be used.	required
diameter	float	Diameter of the monomer.	6.0
max_shifts	int or tuple of int	Maximum shift between subtomograms and template in nm. ZYX order.	(0.5, 0.5, 0.5)
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))
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.	2
method	str	Correlation metrics for alignment.	'ncc'
shift_to_background	tuple of float	If given, the molecules will be internally shifted to get the background region. If not given, the background will be calculated using half1 - half2.	None

Source code in cylindra/widgets/sta.py

@set_design(text="Calculate S/N profiles", location=STAnalysis)
@dask_worker.with_progress(desc="Preparing averaged image")
def calculate_signal_noise_profiles(
    self,
    layer: MoleculesLayerType,
    diameter: Annotated[nm, {"min": 0.1, "max": 100.0}] = 6.0,
    max_shifts: tuple[nm, nm, nm] = (0.5, 0.5, 0.5),
    rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
    bin_size: BinSizeType = 1,
    upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 2,
    method: Annotated[str, {"choices": METHOD_CHOICES}] = "ncc",
    shift_to_background: Annotated[Optional[tuple[nm, nm, nm]], {"text": "Use half tomograms"}] = None,
):  # fmt: skip
    """Calculate the mean landscape energy of signal and noise.

    Parameters
    ----------
    {layer}
    diameter : float, default 6.0
        Diameter of the monomer.
    {max_shifts}{rotations}{interpolation}{bin_size}{upsample_factor}{method}
    shift_to_background : tuple of float, optional
        If given, the molecules will be internally shifted to get the background
        region. If not given, the background will be calculated using half1 - half2.
    """
    main = self._get_main()
    _layer = assert_layer(layer, main.parent_viewer)
    mole = _layer.molecules
    val_s, val_n = yield from self._calculate_signal_noise(
        _layer, diameter, max_shifts, rotations, cutoff=1.0,
        interpolation=interpolation, bin_size=bin_size,
        upsample_factor=upsample_factor, method=method,
        shift_to_background=shift_to_background,
    )  # fmt: skip

    @thread_worker.callback
    def _on_return():
        _layer.molecules = mole.with_features(
            pl.Series("profile_signal", val_s, dtype=pl.Float32),
            pl.Series("profile_noise", val_n, dtype=pl.Float32),
        )
        _Logger.print(f"Profile of {_layer.name!r} has been calculated.")
        with _Logger.set_plt():
            widget_utils.plot_signal_noise_profile(_layer.molecules.features)

    return _on_return

classify_em(layers, templates, max_num_iters=10, min_num_molecules_per_class=None, mask_params=None, inverse_temperature=80, cutoff=0.5, interpolation=3, bin_size=1, seed=0)

Classify molecules to the user given templates using EM algorithm.

Parameters:

Name	Type	Description	Default
layers	list of MoleculesLayer	All the points layers of molecules to be used.	required
max_num_iters	int	Maximum number of iterations to perform.	10
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
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
templates	list of Path	Paths to the template images used for classification. Must be a list of 3D images or a 4D image containing multiple 3D images.	required
inverse_temperature	float	Temperature parameter for the softmax function used in the expectation step. Higher value (lower temperature) makes the classification harder, while lower value (higher temperature) makes it "softer" (i.e., more uncertain).	100

Source code in cylindra/widgets/sta.py

@set_design(text="3D classification", location=STAnalysis)
@dask_worker.with_progress(desc="3D Classification")
def classify_em(
    self,
    layers: MoleculesLayersType,
    templates: Annotated[list[Path], {"filter": FileFilter.IMAGE}],
    max_num_iters: Annotated[int, {"min": 1, "max": 100, "label": "iterations"}] = 10,
    min_num_molecules_per_class: _MinNumMoleculesPerClass = None,
    mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
    inverse_temperature: Annotated[float, {"min": 1, "max": 10000}] = 80,
    cutoff: _CutoffFreq = 0.5,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
    bin_size: BinSizeType = 1,
    seed: _SeedType = 0,
):  # fmt: skip
    """Classify molecules to the user given templates using EM algorithm.

    Parameters
    ----------
    {layers}{max_num_iters}{mask_params}{cutoff}{interpolation}
    {bin_size}
    templates : list of Path
        Paths to the template images used for classification. Must be a list of
        3D images or a 4D image containing multiple 3D images.
    inverse_temperature : float, default 100
        Temperature parameter for the softmax function used in the expectation step.
        Higher value (lower temperature) makes the classification harder, while
        lower value (higher temperature) makes it "softer" (i.e., more uncertain).
    """
    t0 = timer()
    layers = assert_list_of_layers(layers, self.parent_viewer)
    mask = self.params._get_mask(params=mask_params)
    tomo = self._get_main().tomogram
    if isinstance(templates, (Path, str)):
        templates = [Path(templates)]
    template_images: list[np.ndarray] = []
    for _path in templates:
        arr = ip.imread(_path)
        ratio = arr.scale.x / (tomo.scale * bin_size)
        if abs(ratio - 1) > 0.01:
            arr = arr.zoom(ratio, order=3, mode="reflect")
        if arr.ndim == 3:
            template_images.append(arr.value)
        elif arr.ndim == 4:
            template_images.extend(list(arr.value))
        else:
            raise ValueError(
                f"Template image {_path!r} must be a 3D or 4D image, "
                f"but got {arr.ndim}D image."
            )

    combiner = widget_utils.MoleculesCombiner()
    all_moles = combiner.concat(layer.molecules for layer in layers)

    shape = template_images[0].shape
    loader = (
        self._get_loader(
            binsize=bin_size,
            molecules=all_moles,
            order=interpolation,
        )
        .replace(output_shape=shape)
        .order_optimize()
    )
    _Logger.print(
        f"Starting EM classification into user supplied {len(template_images)} templates."
        f"\n{loader.molecules.count()} subtomograms with size {loader.output_shape}"
    )
    if min_num_molecules_per_class is None:
        min_num_molecules_per_class = all_moles.count()
    picker = ClassfySubsetPicker(
        max_num_iters,
        min_num_molecules_per_class * len(template_images),
        np.random.default_rng(seed),
    )
    yield thread_worker.description("Extracting subtomograms")

    with tempfile.TemporaryDirectory() as tmpdir:
        path = Path(tmpdir) / f"subtomograms_{id(t0)}"
        temp_loader = loader.extract_subtomograms(path)
        result = yield from _classify_em_impl(
            template_images, mask, temp_loader, max_num_iters=max_num_iters,
            inverse_temperature=inverse_temperature, cutoff=cutoff, picker=picker,
        )  # fmt: skip
    t0.toc()
    yield thread_worker.description("Classification finished")
    avgs, all_moles = _post_classify_em(result, loader, all_moles)
    moles = combiner.split(all_moles, layers)

    @thread_worker.callback
    def _on_return():
        for layer, each_mole in zip(layers, moles, strict=True):
            layer.set_molecules_with_new_features(each_mole)
        self._show_rec(avgs, name="Classification results", store=False)

    return _on_return

classify_em_template_free(layers, size=None, num_classes=2, min_num_molecules_per_class=None, max_num_iters=10, mask_params=None, inverse_temperature=80, cutoff=0.5, interpolation=3, bin_size=1, seed=0)

Classify molecules using EM algorithm.

Random split averages will be used as initial templates.

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
max_num_iters	int	Maximum number of iterations to perform.	10
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
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
num_classes	int	The number of classes to classify subtomograms into.	2
inverse_temperature	float	Temperature parameter for the softmax function used in the expectation step. Higher value (lower temperature) makes the classification harder, while lower value (higher temperature) makes it "softer" (i.e., more uncertain).	100
seed	int	Random seed.	0

Source code in cylindra/widgets/sta.py

@set_design(text="3D classification (template-free)", location=STAnalysis)
@dask_worker.with_progress(desc="3D Classification (template-free)")
def classify_em_template_free(
    self,
    layers: MoleculesLayersType,
    size: Annotated[float, {"value": 12.0, "max": 100.0, "label": "size (nm)"}] = None,
    num_classes: Annotated[int, {"min": 1, "max": 100}] = 2,
    min_num_molecules_per_class: _MinNumMoleculesPerClass = None,
    max_num_iters: Annotated[int, {"min": 1, "max": 100, "label": "iterations"}] = 10,
    mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
    inverse_temperature: Annotated[float, {"min": 1, "max": 10000}] = 80,
    cutoff: _CutoffFreq = 0.5,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
    bin_size: BinSizeType = 1,
    seed: _SeedType = 0,
):  # fmt: skip
    """Classify molecules using EM algorithm.

    Random split averages will be used as initial templates.

    Parameters
    ----------
    {layers}{size}{max_num_iters}{mask_params}{cutoff}{interpolation}{bin_size}
    num_classes : int, default 2
        The number of classes to classify subtomograms into.
    inverse_temperature : float, default 100
        Temperature parameter for the softmax function used in the expectation step.
        Higher value (lower temperature) makes the classification harder, while
        lower value (higher temperature) makes it "softer" (i.e., more uncertain).
    seed : int, optional
        Random seed.
    """
    t0 = timer()
    layers = assert_list_of_layers(layers, self.parent_viewer)
    combiner = widget_utils.MoleculesCombiner()
    mask = self.params._get_mask(params=mask_params)
    all_moles = combiner.concat(layer.molecules for layer in layers)

    shape = (size,) * 3
    loader = self._get_loader(
        binsize=bin_size, molecules=all_moles, order=interpolation, shape=shape
    ).order_optimize()
    _Logger.print(
        f"Starting EM classification into {num_classes} classes.\n"
        f"{loader.molecules.count()} subtomograms with size {loader.output_shape}"
    )
    rng = np.random.default_rng(seed)
    if min_num_molecules_per_class is None:
        min_num_molecules_per_class = all_moles.count()
    picker = ClassfySubsetPicker(
        max_num_iters,
        min_num_molecules_per_class * num_classes,
        rng,
    )
    yield thread_worker.description("Extracting subtomograms")

    with tempfile.TemporaryDirectory() as tmpdir:
        path = Path(tmpdir) / f"subtomograms_{id(t0)}"
        temp_loader = loader.extract_subtomograms(path)
        yield thread_worker.description("Creating initial class templates")
        templates = temp_loader.average_split(n_split=num_classes, seed=rng)
        result = yield from _classify_em_impl(
            templates, mask, temp_loader, max_num_iters=max_num_iters,
            inverse_temperature=inverse_temperature, cutoff=cutoff, picker=picker,
        )  # fmt: skip

    t0.toc()
    yield thread_worker.description("Classification finished")
    avgs, all_moles = _post_classify_em(result, loader, all_moles)
    moles = combiner.split(all_moles, layers)

    @thread_worker.callback
    def _on_return():
        for layer, each_mole in zip(layers, moles, strict=True):
            layer.set_molecules_with_new_features(each_mole)
        self._show_rec(avgs, name="Classification results", store=False)

    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: BinSizeType = 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, norm=norm,
        upsample_factor=upsample_factor, method=method,
    )  # fmt: skip
    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

detect_defects(layer, diameter=5.8, prior_defect=0.1, max_shifts=(0.5, 0.5, 0.5), rotations=((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)), interpolation=3, bin_size=1, upsample_factor=2, method='ncc', shift_to_background=None)

Detect lattice defect by Bayesian metrics.

This method calculates the signal and noise distribution, and evaluate the probability of being a defect site for individual molecules.

Parameters:

Name	Type	Description	Default
layer	MoleculesLayer	Points layer of molecules to be used.	required
diameter	float	Diameter of monomer in nanometers.	5.8
prior_defect	float	Prior probability of a site being a defect. This value does not necessarily reflect the actual defect rate. Changing this value only rescales the output probability.	0.1
max_shifts	int or tuple of int	Maximum shift between subtomograms and template in nm. ZYX order.	(0.5, 0.5, 0.5)
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))
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.	2
method	str	Correlation metrics for alignment.	'ncc'
shift_to_background	tuple of float	Shift vector to move molecules to background. For example, use (7, 0, 0) for tubulin subunits.	None

Source code in cylindra/widgets/sta.py

@set_design(text="Detect defects", location=STAnalysis)
@dask_worker.with_progress(desc="Preparing averaged image")
def detect_defects(
    self,
    layer: MoleculesLayerType,
    diameter: Annotated[nm, {"min": 0.1, "max": 100.0}] = 5.8,
    prior_defect: Annotated[float, {"min": 0.001, "max": 0.999}] = 0.1,
    max_shifts: tuple[nm, nm, nm] = (0.5, 0.5, 0.5),
    rotations: _Rotations = ((0.0, 0.0), (0.0, 0.0), (0.0, 0.0)),
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
    bin_size: BinSizeType = 1,
    upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 2,
    method: Annotated[str, {"choices": METHOD_CHOICES}] = "ncc",
    shift_to_background: Annotated[Optional[tuple[nm, nm, nm]], {"text": "Use half tomograms"}] = None,
):  # fmt: skip
    """Detect lattice defect by Bayesian metrics.

    This method calculates the signal and noise distribution, and evaluate the
    probability of being a defect site for individual molecules.

    Parameters
    ----------
    {layer}
    diameter : float
        Diameter of monomer in nanometers.
    prior_defect : float
        Prior probability of a site being a defect. This value does not necessarily
        reflect the actual defect rate. Changing this value only rescales the output
        probability.
    {max_shifts}{rotations}{interpolation}{bin_size}{upsample_factor}{method}
    shift_to_background : tuple of float, optional
        Shift vector to move molecules to background. For example, use (7, 0, 0) for
        tubulin subunits.
    """
    main = self._get_main()
    _layer = assert_layer(layer, main.parent_viewer)
    mole = _layer.molecules
    val_s, val_n = yield from self._calculate_signal_noise(
        _layer, diameter, max_shifts, rotations, cutoff=1.0,
        interpolation=interpolation, bin_size=bin_size,
        upsample_factor=upsample_factor, method=method,
        shift_to_background=shift_to_background,
    )  # fmt: skip

    mean_s, std_s = val_s.mean(), val_s.std(ddof=0)
    mean_n, std_n = val_n.mean(), val_n.std(ddof=0)
    _log_mean_sd_for_detect_defect(mean_s, std_s, mean_n, std_n)

    post_prob = utils.calc_post_probability(
        val_s, 1 - prior_defect, mean_s, std_s, mean_n, std_n
    )

    prob_defect = 1 - post_prob

    @thread_worker.callback
    def _on_return():
        _layer.molecules = mole.with_features(
            pl.Series("profile_signal", val_s, dtype=pl.Float32),
            pl.Series("profile_noise", val_n, dtype=pl.Float32),
            pl.Series("prob_defect", prob_defect, dtype=pl.Float32),
        )
        if not isinstance(base_color := _layer.colormap_info, str):
            base_color = "lime"
        _layer.set_colormap("prob_defect", (0, 1), cmap={0: base_color, 1: "black"})

    return _on_return

extract_subtomograms(layers, size=None, interpolation=1, bin_size=1, compute=False)

Extract subtomograms around molecules in the selected layer(s).

Extracted 3D images will be sent to the sub-viewer. Note that this extraction operation is NOT needed for subsequent averaging/alignment operations, as they internally extract subtomograms by themselves.

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
compute	bool	If True, the subtomograms will be computed immediately, and in-memory data will be sent to the sub-viewer.	False

Source code in cylindra/widgets/sta.py

@set_design(text="Extract subtomograms", location=Averaging)
@dask_worker.with_progress(desc=_pdesc.fmt_layers("Extract subtomograms of {!r}"))
def extract_subtomograms(
    self,
    layers: MoleculesLayersType,
    size: _SubVolumeSize = None,
    interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 1,
    bin_size: BinSizeType = 1,
    compute: bool = False,
):
    """Extract subtomograms around molecules in the selected layer(s).

    Extracted 3D images will be sent to the sub-viewer. Note that this extraction
    operation is NOT needed for subsequent averaging/alignment operations, as they
    internally extract subtomograms by themselves.

    Parameters
    ----------
    {layers}{size}{interpolation}{bin_size}
    compute : bool, default: False
        If True, the subtomograms will be computed immediately, and in-memory data
        will be sent to the sub-viewer.
    """
    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
    )
    subtomos = ip.lazy.asarray(loader.construct_dask(), axes="pzyx")
    if compute:
        subtomos = subtomos.compute()
    return self._show_rec.with_args(
        subtomos, f"[Subtomos]{_avg_name(layers)}", store=False
    )

fit_spline_rfa(spline, template_path, forward_is='MinusToPlus', interval='4.1', err_max=0.5, mask_params=None, max_shifts=(2.0, 2.0, 2.0), rotations=((0.0, 0.0), (15.0, 5.0), (0.0, 0.0)), cutoff=0.5, interpolation=3, range_long=(4.0, 4.28), angle_max=5.0, bin_size=1, temperature_time_const=10.0, upsample_factor=5, num_trials=5, seed=0)

Fit spline by RFA.

This algorithm uses a template image to precisely determine the center line of filaments. By comparing the score, the orientation of the filament will also be determined.

Parameters:

Name	Type	Description	Default
spline	int	Index of splines to be used.	required
template_path	Path or str	Path to template image.	required
forward_is	PlusToMinus or MinusToPlus	Which orientation is the forward direction. Set "PlusToMinus" if the template image is oriented from the plus end to the minus end in the y direction.	"MinusToPlus"
interval	float or str expression	Interval of the sampling points along the spline.	4.1
err_max	float	S.D. allowed for spline fitting. Larger value will result in smoother spline, i.e. fewer spline knots.	0.5
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.	(2.0, 2.0, 2.0)
rotations	((float, float), (float, float), (float, float))	Rotation in external degree around each axis.	((0.0, 0.0), (15.0, 5.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
temperature_time_const	float	Time constant of the temperature decay during annealing. Larger value results in slower annealing. 1.0 is a moderate value.	10.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
num_trials	int	Number of independent trials to perform. The best result among all trials will be returned.	5
seed	int	Random seed for reproducibility.	0

Source code in cylindra/widgets/sta.py

@set_design(text="Fit spline by RFA", location=Alignment)
@dask_worker.with_progress(descs=_pdesc.fit_spline_rfa_fmt)
def fit_spline_rfa(
    self,
    spline: SplineType,
    template_path: Annotated[_PathOrPathsOrNone, {"bind": _template_params}],
    forward_is: PlusOrMinusType = "MinusToPlus",
    interval: PolarsExprStrOrScalar = "4.1",
    err_max: Annotated[nm, {"label": "Max fit error (nm)", "step": 0.1}] = 0.5,
    mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
    max_shifts: _MaxShifts = (2.0, 2.0, 2.0),
    rotations: _Rotations = ((0.0, 0.0), (15.0, 5.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: BinSizeType = 1,
    temperature_time_const: Annotated[float, {"min": 0.01, "max": 100.0}] = 10.0,
    upsample_factor: Annotated[int, {"min": 1, "max": 20}] = 5,
    num_trials: Annotated[int, {"min": 1, "max": 100}] = 5,
    seed: _SeedType = 0,
):
    """Fit spline by RFA.

    This algorithm uses a template image to precisely determine the center line of
    filaments. By comparing the score, the orientation of the filament will also be
    determined.

    Parameters
    ----------
    {spline}{template_path}
    forward_is : "PlusToMinus" or "MinusToPlus", default "MinusToPlus"
        Which orientation is the forward direction. Set "PlusToMinus" if the
        template image is oriented from the plus end to the minus end in the y
        direction.
    interval : float or str expression, default 4.1
        Interval of the sampling points along the spline.
    {err_max}{mask_params}{max_shifts}{rotations}{cutoff}{interpolation}{range_long}
    {angle_max}{bin_size}{temperature_time_const}{upsample_factor}{num_trials}{seed}
    """
    t0 = timer()
    main = self._get_main()
    tomo = main.tomogram
    spl_old = tomo.splines[spline]
    interv_expr = widget_utils.norm_scalar_expr(interval)
    mole_fw = tomo.map_centers(
        i=spline,
        interval=spl_old.props.get_glob(interv_expr),
        rotate_molecules=False,
    )
    mole_rv = mole_fw.rotate_by_rotvec_internal([np.pi, 0, 0])  # invert

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

    def _construct_landscape(mole_):
        landscape_, _ = self._construct_landscape(
            molecules=mole_, template_path=template_path, mask_params=mask_params,
            max_shifts=max_shifts, rotations=rotations, cutoff=cutoff, norm=False,
            order=interpolation, bin_size=bin_size, upsample_factor=upsample_factor,
        )  # fmt: skip
        yield
        mole, results = landscape_.run_filamentous_annealing(
            range=range_long,
            angle_max=angle_max,
            temperature_time_const=temperature_time_const,
            random_seeds=utils.create_random_seeds(num_trials, seed),
        )
        yield _on_yield.with_args(results)
        return mole, results

    mole_opt_fw, results_fw = yield from _construct_landscape(mole_fw)
    mole_opt_rv, results_rv = yield from _construct_landscape(mole_rv)

    if results_fw[0].energies[-1] < results_rv[0].energies[-1]:  # forward is better
        mole_opt = mole_opt_fw
        ori = Ori(forward_is)
    else:  # reverse is better
        mole_opt = mole_opt_rv
        ori = Ori.invert(Ori(forward_is))

    # calculate distances for logging
    _ds = np.diff(mole_opt.pos, axis=0)
    _dist: NDArray[np.float32] = np.sqrt(np.sum(_ds**2, axis=1))

    spl_new = CylSpline(
        order=spl_old.order,
        config=spl_old.config,
        extrapolate=spl_old.extrapolate,
    ).fit(mole_opt.pos, err_max=err_max)
    spl_new.props.loc = spl_old.props.loc
    spl_new.props.glob = spl_old.props.glob
    spl_new.orientation = ori
    tomo.splines[spline] = spl_new
    t0.toc()

    @thread_worker.callback
    def _on_return():
        main._update_splines_in_images()
        _Logger.print(f"Orientation: {ori.value}")
        _Logger.print(
            f"Distance: mean = {_dist.mean():.3f} nm (ranging from "
            f"{_dist.min():.3f} to {_dist.max():.3f}) nm "
        )
        main.reset_choices()

        @undo_callback
        def _out():
            tomo.splines[spline] = spl_old
            main._update_splines_in_images()

        @_out.with_redo
        def _out():
            tomo.splines[spline] = spl_new
            main._update_splines_in_images()

        return _out

    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: LandscapeLayerType,
    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)

refine_tilt_series(layers, size=None, max_shift=2.0, rotation=(0, 0), grid_shape=(1, 1), smoothness=0.01, ctf=True, bin_size=1, update_coarse_tomogram=True)

Refine the alignment of the tilt series using the molecules.

The subtomograms of each layer are averaged as the reference of the molecules of the layer. Projections of the averages placed at the molecules are compared with the tilt images, and the shifts of each tilt image are refined.

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
max_shift	nm	Maximum shift of the tilt images in nm.	2.0
rotation	(float, float)	Maximum rotation and the step of the tilt images in degrees.	(0, 0)
grid_shape	(int, int)	Number of control points of the local shifts along y and x of the tilt images. (1, 1) means a uniform shift of each tilt image. Otherwise, the shifts are refined as a smooth field over each tilt image.	(1, 1)
smoothness	float	Weight of the smoothness prior of the local shifts. Larger value is more robust against noise and fewer molecules.	0.01
ctf	bool	If true, the CTF is applied to the projections of the averages, with the defocus of each molecule in each tilt, considering its depth. The tilt series must have the CTF parameters.	True
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1
update_coarse_tomogram	bool	If true, the coarse tomogram for visualization is also reconstructed from the refined tilt series.	True

Source code in cylindra/widgets/sta.py

@set_design(text=capitalize, location=Alignment)
@dask_worker.with_progress(desc="Refining the tilt series")
def refine_tilt_series(
    self,
    layers: MoleculesLayersType,
    size: _SubVolumeSize = None,
    max_shift: Annotated[nm, {"label": "max shift (nm)", "min": 0.1, "max": 50.0, "step": 0.1}] = 2.0,
    rotation: Annotated[tuple[float, float], {"label": "rotation max/step (deg)"}] = (0, 0),
    grid_shape: Annotated[tuple[int, int], {"label": "grid shape (y, x)", "options": {"min": 1, "max": 20}}] = (1, 1),
    smoothness: Annotated[float, {"min": 0.0, "max": 100.0, "step": 0.001}] = 0.01,
    ctf: Annotated[bool, {"label": "apply CTF"}] = True,
    bin_size: BinSizeType = 1,
    update_coarse_tomogram: bool = True,
):  # fmt: skip
    """Refine the alignment of the tilt series using the molecules.

    The subtomograms of each layer are averaged as the reference of the molecules
    of the layer. Projections of the averages placed at the molecules are compared
    with the tilt images, and the shifts of each tilt image are refined.

    Parameters
    ----------
    {layers}{size}
    max_shift : nm
        Maximum shift of the tilt images in nm.
    rotation : (float, float)
        Maximum rotation and the step of the tilt images in degrees.
    grid_shape : (int, int)
        Number of control points of the local shifts along y and x of the tilt
        images. (1, 1) means a uniform shift of each tilt image. Otherwise, the
        shifts are refined as a smooth field over each tilt image.
    smoothness : float
        Weight of the smoothness prior of the local shifts. Larger value is more
        robust against noise and fewer molecules.
    ctf : bool, default True
        If true, the CTF is applied to the projections of the averages, with the
        defocus of each molecule in each tilt, considering its depth. The tilt
        series must have the CTF parameters.
    {bin_size}
    update_coarse_tomogram : bool, default True
        If true, the coarse tomogram for visualization is also reconstructed from
        the refined tilt series.
    """
    t0 = timer()
    layers = assert_list_of_layers(layers, self.parent_viewer)
    main = self._get_main()
    tomo = main.tomogram
    if (src := tomo.tilt_series) is None:
        raise ValueError("The tomogram is not opened from a tilt series.")
    shape = self._get_shape_in_nm(size)
    refs = list[tuple[Molecules, "NDArray[np.float32]"]]()
    for layer in layers:
        # the average of the current tilt series is the reference of each layer
        mole = layer.molecules
        loader = tomo.get_subtomogram_loader(mole, shape, binsize=bin_size)
        yield thread_worker.description(f"Creating template for layer {layer.name}")
        refs.append((mole, loader.average(wedge_weighting=True)))
    nmole = sum(mole.count() for mole, _ in refs)
    _Logger.print(f"Refining the tilt series using {nmole} molecules ...")
    old = src.tilt_series
    loader = src.get_loader(Molecules.empty(), binsize=bin_size)
    yield thread_worker.description("Running tilt series refinement")
    new = loader.refine_tilt_series(
        refs,
        max_shift=max_shift,
        grid_shape=grid_shape,
        smoothness=smoothness,
        rotation=rotation,
        ctf=ctf,
    ).tilt_series
    if update_coarse_tomogram:
        yield thread_worker.description("Reconstructing coarse tomogram")
        new_ref = src.reconstruct_reference(new)
    else:
        new_ref = None
    t0.toc()

    old_ref = src.reference

    def _apply(tilt_series: "TiltSeries", reference: "ip.ImgArray | None"):
        removed = tomo.update_tilt_series(tilt_series)
        if new_ref is not None and reference is not None:
            src.set_reference(reference)
            self._show_coarse_tomogram(reference)
        if main.SplineControl.need_resample:
            main.sample_subtomograms()
        return removed

    @thread_worker.callback
    def _on_return():
        removed = _apply(new, new_ref)
        _log_tilt_series_shifts(old, new)
        if removed:
            _Logger.print(
                f"Multi-scale images of bin size {removed} are removed. They are "
                "reconstructed from the refined tilt series on demand."
            )
        # the tilt series and the coarse tomogram are restored by undo
        return (
            undo_callback(_apply)
            .with_args(old, old_ref)
            .with_redo(lambda: _apply(new, new_ref))
        )

    return _on_return

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: LandscapeLayerType,
    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", location=LandscapeMenu)
@dask_worker.with_progress(desc="Peak detection on landscape")
def run_align_on_landscape(self, landscape_layer: LandscapeLayerType):
    """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_rfa_on_landscape(landscape_layer, range_long=(4.0, 4.28), angle_max=5.0, temperature_time_const=1.0, num_trials=5, seed=0)

Run simulated annealing on the landscape, supposing a filamentous structure.

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
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
num_trials	int	Number of independent trials to perform. The best result among all trials will be returned.	5
seed	int	Random seed for reproducibility.	0

Source code in cylindra/widgets/sta.py

@set_design(text="Run annealing (RFA)", location=LandscapeMenu)
@dask_worker.with_progress(desc="Running simulated annealing")
def run_rfa_on_landscape(
    self,
    landscape_layer: LandscapeLayerType,
    range_long: _DistRangeLon = (4.0, 4.28),
    angle_max: _AngleMaxLon = 5.0,
    temperature_time_const: Annotated[float, {"min": 0.01, "max": 10.0}] = 1.0,
    num_trials: Annotated[int, {"min": 1, "max": 100}] = 5,
    seed: _SeedType = 0,
):
    """Run simulated annealing on the landscape, supposing a filamentous structure.

    Parameters
    ----------
    {landscape_layer}{range_long}{angle_max}{temperature_time_const}{num_trials}
    {seed}
    """
    t0 = timer()
    landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
    mole, results = landscape_layer.landscape.run_filamentous_annealing(
        range=range_long,
        angle_max=angle_max,
        temperature_time_const=temperature_time_const,
        random_seeds=utils.create_random_seeds(num_trials, seed),
    )
    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=landscape_layer.source_spline,
        metadata={ANNEALING_RESULT: results[0]},
    )

run_rma_on_landscape(landscape_layer, range_long=(4.0, 4.28), range_lat=('d.mean() - 0.1', 'd.mean() + 0.1'), angle_max=5.0, temperature_time_const=1.0, lj_const=0.0, num_trials=5, seed=0)

Run simulated annealing on the landscape, supposing a cylindric structure.

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	('d.mean() - 0.1', 'd.mean() + 0.1')
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
lj_const	float	If greater than 0, the annealing will use a Lennard-Jones-like potential to allow molecules to separate further apart than the cutoff distance. Larger value will prefer longer distances between molecules. 0.25 - 0.5 is usually a good choice.	0.0
num_trials	int	Number of independent trials to perform. The best result among all trials will be returned.	5
seed	int	Random seed for reproducibility.	0

Source code in cylindra/widgets/sta.py

@set_design(text="Run annealing (RMA)", location=LandscapeMenu)
@dask_worker.with_progress(desc="Running simulated annealing")
def run_rma_on_landscape(
    self,
    landscape_layer: LandscapeLayerType,
    range_long: _DistRangeLon = (4.0, 4.28),
    range_lat: _DistRangeLat = ("d.mean() - 0.1", "d.mean() + 0.1"),
    angle_max: _AngleMaxLon = 5.0,
    temperature_time_const: Annotated[float, {"min": 0.01, "max": 10.0}] = 1.0,
    lj_const: Annotated[float, {"min": 0.00, "max": 1000.0, "step": 0.1, "label": "LJ const"}] = 0.0,
    num_trials: Annotated[int, {"min": 1, "max": 100}] = 5,
    seed: _SeedType = 0,
):  # fmt: skip
    """Run simulated annealing on the landscape, supposing a cylindric structure.

    Parameters
    ----------
    {landscape_layer}{range_long}{range_lat}{angle_max}{temperature_time_const}
    {lj_const}{num_trials}{seed}
    """
    t0 = timer()
    landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
    if (spl := landscape_layer.source_spline) is None:
        spl = CylSpline._dummy_spline()
        _Logger.print(
            "No source spline found, cylindrical boundary will not be considered."
        )

    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,
        lj_nstd=None if lj_const < 1e-3 else 1 / lj_const,
        random_seeds=utils.create_random_seeds(num_trials, seed),
    )
    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", location=LandscapeMenu)
@dask_worker.with_progress(desc="Running Viterbi alignment")
def run_viterbi_on_landscape(
    self,
    landscape_layer: LandscapeLayerType,
    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.epoch_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 between template(s) and averages. If molecule assembly has 13 protofilaments, this method will try 26 patterns of averages.

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	str	If "Raw", all the raw subtomogram averages will be shown. If "Filtered", all the filtered subtomogram averages will be shown. If None, no subtomogram averages will be shown.	"Filtered"
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: Path.Read[FileFilter.IMAGE],
    mask_params: Annotated[Any, {"bind": _get_mask_params}],
    anti_template_path: Annotated[Path.Read[FileFilter.IMAGE] | None, {"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 between
    template(s) and averages. If molecule assembly has 13 protofilaments, this
    method will try 26 patterns of averages.

    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 : str, default "Filtered"
        If "Raw", all the raw subtomogram averages will be shown.
        If "Filtered", all the filtered subtomogram averages will be shown.
        If None, no subtomogram averages will be shown.
    {cutoff}
    """
    t0 = timer()
    layer = assert_layer(layer, self.parent_viewer)
    loader, npf = prep_seam_search_input(
        self._get_main().tomogram,
        layer.molecules,
        npf=npf,
        order=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 and str(anti_template_path) not in ("", "."):
        anti_template, anti_mask = loader.normalize_input(
            template=self.params._norm_template_param(anti_template_path),
            mask=self.params._get_mask(params=mask_params),
        )
        anti_template = ip.asarray(anti_template, axes="zyx")
        if mask is not None and anti_mask is not None:
            mask = np.maximum(mask, anti_mask)
    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.lower() == "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.
    """
    ui = self._get_main()
    layer = assert_layer(layer, ui.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)
    ui.reset_choices()  # choices regarding of features need update
    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.
    """
    ui = self._get_main()
    layer = assert_layer(layer, ui.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)
    ui.reset_choices()  # choices regarding of features need update
    return undo_callback(layer.feature_setter(feat, layer.colormap_info))

seam_search_on_tilt(layer, template_path, mask_params, anti_template_path=None, npf=None, bin_size=1, cutoff=0.25, ctf=True)

Search for the best seam position using the scores in the tilt images.

The template (and the anti-template) is projected to each tilt image at the position and orientation of each molecule, and compared with the tilt images by ZNCC. The score of each molecule is the mean over the tilts. The class means of the scores are compared for all the seam positions. The tomogram must be opened from a tilt series.

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
npf	int	Number of protofilaments. By default the global properties stored in the corresponding spline will be used.	None
bin_size	int	Bin size of multiscale image to be used. Set to >1 to boost performance.	1
cutoff	float	Cutoff frequency of low-pass filter applied in each subtomogram.	0.25
ctf	bool	If true, the CTF is applied to the projections of the templates, with the defocus of each molecule in each tilt, considering its depth. The tilt series must have the CTF parameters.	True

Source code in cylindra/widgets/sta.py

@set_design(text="Seam search on tilt series", location=STAnalysis.SeamSearch)
@dask_worker.with_progress(desc=_pdesc.fmt_layer("Seam search of {!r} on tilt series"))  # fmt: skip
def seam_search_on_tilt(
    self,
    layer: MoleculesLayerType,
    template_path: Path.Read[FileFilter.IMAGE],
    mask_params: Annotated[Any, {"bind": _get_mask_params}],
    anti_template_path: Annotated[Path.Read[FileFilter.IMAGE] | None, {"label": "anti-template path"}] = None,
    npf: Annotated[Optional[int], {"text": "use global properties"}] = None,
    bin_size: BinSizeType = 1,
    cutoff: _CutoffFreq = 0.25,
    ctf: Annotated[bool, {"label": "apply CTF"}] = True,
):  # fmt: skip
    """Search for the best seam position using the scores in the tilt images.

    The template (and the anti-template) is projected to each tilt image at the
    position and orientation of each molecule, and compared with the tilt images
    by ZNCC. The score of each molecule is the mean over the tilts. The class means
    of the scores are compared for all the seam positions. The tomogram must be
    opened from a tilt series.

    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.
    npf : int, optional
        Number of protofilaments. By default the global properties stored in the
        corresponding spline will be used.
    {bin_size}{cutoff}
    ctf : bool, default True
        If true, the CTF is applied to the projections of the templates, with the
        defocus of each molecule in each tilt, considering its depth. The tilt
        series must have the CTF parameters.
    """
    t0 = timer()
    layer = assert_layer(layer, self.parent_viewer)
    if (src := self._get_main().tomogram.tilt_series) is None:
        raise ValueError("The tomogram is not opened from a tilt series.")
    mole = layer.molecules
    if npf is None:
        npf = mole.features[Mole.pf].unique().len()
    loader = src.get_loader(mole, binsize=bin_size)
    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 and str(anti_template_path) not in ("", "."):
        anti_template, anti_mask = loader.normalize_input(
            template=self.params._norm_template_param(anti_template_path),
            mask=self.params._get_mask(params=mask_params),
        )
        # the same mask is used for both templates
        if mask is not None and anti_mask is not None:
            mask = np.maximum(mask, anti_mask)
    else:
        anti_template = None

    seam_searcher = TiltScoreSeamSearcher(npf)
    result = seam_searcher.search(
        loader=loader,
        template=template,
        anti_template=anti_template,
        mask=mask,
        cutoff=cutoff,
        ctf=ctf,
    )

    t0.toc()

    @thread_worker.callback
    def _seam_search_on_return():
        new_feat = result.as_series(mole.count())
        layer.features = layer.molecules.features.with_columns(new_feat)
        layer.metadata[SEAM_SEARCH_RESULT] = result

        with _Logger.set_plt():
            _Logger.print(f"layer = {layer.name!r}")
            _Logger.print(f"template = {Path(template_path).as_posix()!r}")
            if anti_template is not None:
                _Logger.print(
                    f"anti_template = {Path(anti_template_path).as_posix()!r}"
                )
            _Logger.print(f"seam position = {result.seam_pos}")
            widget_utils.plot_seam_search_result(result.scores, npf)

    return _seam_search_on_return

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)