remove paint_cylinder

cylindra/components/cylindric.py

@@ -232,7 +232,12 @@ def locate_molecules(self, spl: Spline, coords: NDArray[np.int32]) -> Molecules:
         return mole
 
     def to_mesh(self, spl: Spline):
-        """Create a mesh data for cylinder visualization."""
+        """
+        Create a mesh data for cylinder visualization.
+
+        Returned mesh is a tuple of (nodes, vertices). Nodes is a (N, 3) array, where
+        N is the y-shape of this object.
+        """
         nodes = spl.cylindrical_to_world(
             self.replace(tilts=(0, 0))._get_regular_mesh().reshape(-1, 3)
         )

cylindra/cylmeasure.py

@@ -209,22 +209,27 @@ def transform_vector(
         start_pos = start[:, 3]
 
         # the displacement in the spline coordinate
-        dpos = self._integrate_dpos(start_pos, vec)
+        dy = self._integrate_yvec(start_pos, vec)
+        end_pos = start_pos + dy
 
         # surface normal vectors at the start and end points
         er0 = self._surface_vec(_concat(start_zyx, start_pos))
-        er1 = self._surface_vec(_concat(start_zyx + vec, start_pos + dpos))
+        er1 = self._surface_vec(_concat(start_zyx + vec, end_pos))
+        ey0 = self._spline_vec_norm(start_pos)
+        ey1 = self._spline_vec_norm(end_pos)
 
+        # the "mean" unit vectors
+        ey = _norm(ey0 + ey1)
         er = _norm(er0 + er1)
-        ey = self._spline_vec_norm(start_pos + dpos / 2)
-        ea = -np.cross(er, ey, axis=1)
+        ea = np.cross(er, ey, axis=1)
 
         er0y = _norm(_cancel_component(er0, ey))
         er1y = _norm(_cancel_component(er1, ey))
         v_ang_len = 2 * self._radius * _half_sin(er0y, er1y)
-        _factor = _dot(_norm(er0y - er1y), ea)
-        v_ang = v_ang_len * _factor
-        return np.stack([_dot(vec, er), _dot(vec, ey), v_ang], axis=1)
+        _sign = np.sign(_dot(_norm(er1y - er0y), ea))
+        v_ang = v_ang_len * _sign
+        v_ang[np.isnan(v_ang)] = 0.0
+        return np.stack([_dot(vec, er), dy, v_ang], axis=1)
 
     def _surface_vec(
         self,
@@ -253,12 +258,19 @@ def _spline_vec_norm(self, pos: NDArray[np.float32]) -> NDArray[np.float32]:
         _spl_vec = self._spl.map(u, der=1)
         return _norm(_spl_vec)
 
-    def _integrate_dpos(self, start_pos: NDArray[np.float32], vec: NDArray[np.float32]):
-        ds0 = _dot(self._spline_vec_norm(start_pos) / 4, vec)
-        ds1 = _dot(self._spline_vec_norm(start_pos + ds0) / 4, vec)
-        ds2 = _dot(self._spline_vec_norm(start_pos + ds0 + ds1) / 4, vec)
-        ds3 = _dot(self._spline_vec_norm(start_pos + ds0 + ds1 + ds2) / 4, vec)
-        return ds0 + ds1 + ds2 + ds3
+    def _integrate_yvec(
+        self,
+        start_pos: NDArray[np.float32],
+        vec: NDArray[np.float32],
+        num: int = 4,
+    ):
+        """Nonlinear dot product of vec and spline tangent vector."""
+        cur_pos = start_pos.copy()
+        for _ in range(num):
+            norm = self._spline_vec_norm(cur_pos) / num
+            ds0 = _dot(norm, vec)
+            cur_pos += ds0
+        return cur_pos - start_pos
 
 
 class RegionProfiler:
@@ -470,7 +482,8 @@ def _concat_groups(subsets: list[Molecules]) -> Molecules:
 
 
 def _half_sin(er0y, er1y):
+    """Calculate sin(a/2) from two vectors, where a is the angle between them."""
     # sin^2(a) = (1 - cos(2a)) / 2, when 0 < a < pi
     cos2a = _dot(er0y, er1y) / np.sqrt(_dot(er0y, er0y) * _dot(er1y, er1y))
-    sina = np.sqrt((1 - np.clip(cos2a, 0, 1)) / 2)
+    sina = np.sqrt((1 - np.clip(cos2a, -1, 1)) / 2)
     return sina

cylindra/widget_utils.py

@@ -352,114 +352,6 @@ def coordinates_with_extensions(
     return np.concatenate(coords, axis=0)
 
 
-class PaintDevice:
-    """
-    Device used for painting 3D images.
-
-    Parameters
-    ----------
-    shape : tuple[int, int, int]
-        Shape of the target image.
-    scale : nm
-        Scale of the target image.
-    """
-
-    def __init__(self, shape: tuple[int, int, int], scale: nm):
-        self._shape = shape
-        self._scale = scale
-
-    @property
-    def scale(self) -> nm:
-        return self._scale
-
-    def paint_cylinders(self, tomo: CylTomogram, prop: str):
-        from scipy import ndimage as ndi
-
-        lbl = np.zeros(self._shape, dtype=np.uint8)
-        all_df = tomo.splines.collect_localprops()
-        if all_df is None:
-            raise ValueError("No local property found.")
-        bin_scale = self._scale  # scale of binned reference image
-        binsize = utils.roundint(bin_scale / tomo.scale)
-
-        cylinders = list[list[NDArray[np.bool_]]]()
-        matrices = list[list[NDArray[np.float32]]]()
-        for i, spl in enumerate(tomo.splines):
-            depth = spl.props.window_size.get(prop, spl.config.fit_depth)
-            lz, ly, lx = (
-                utils.roundint(r / bin_scale * 1.73) * 2 + 1
-                for r in [15, depth / 2, 15]
-            )
-            center = np.array([lz, ly, lx]) / 2 + 0.5
-            z, _, x = np.indices((lz, ly, lx))
-            # Prepare template hollow image
-            _sq = (z - lz / 2 - 0.5) ** 2 + (x - lx / 2 - 0.5) ** 2
-            domains = list[NDArray[np.bool_]]()
-            dist = [-np.inf] + list(spl.distances()) + [np.inf]
-            if len(spl.props.get_loc(H.radius, [])) == spl.anchors.size:
-                radii = spl.props.get_loc(H.radius)
-            elif spl.props.has_glob(H.radius):
-                radii = [spl.props.get_glob(H.radius)] * spl.anchors.size
-            else:
-                raise RuntimeError(f"Radius not found in spline-{i}.")
-            for j, rc in enumerate(radii):
-                r0 = max(rc - spl.config.thickness_inner, 0.0) / tomo.scale / binsize
-                r1 = (rc + spl.config.thickness_outer) / tomo.scale / binsize
-                domain = (r0**2 < _sq) & (_sq < r1**2)
-                ry = (
-                    min(
-                        abs(dist[j + 1] - dist[j]) / 2,
-                        abs(dist[j + 2] - dist[j + 1]) / 2,
-                        depth / 2,
-                    )
-                    / bin_scale
-                    + 0.5
-                )
-
-                ry = max(utils.ceilint(ry), 1)
-                domain[:, : ly // 2 - ry] = 0
-                domain[:, ly // 2 + ry + 1 :] = 0
-                domain = domain.astype(np.float32)
-                domains.append(domain)
-
-            cylinders.append(domains)
-            matrices.append(spl.affine_matrix(center=center, inverse=True))
-            yield
-
-        cylinders = np.concatenate(cylinders, axis=0)
-        matrices = np.concatenate(matrices, axis=0)
-
-        out = np.empty_like(cylinders)
-        for i, (img, matrix) in enumerate(zip(cylinders, matrices, strict=True)):
-            out[i] = ndi.affine_transform(img, matrix, order=1, cval=0, prefilter=False)
-        out = out > 0.3
-
-        # paint roughly
-        for i, crd in enumerate(tomo.splines.iter_anchor_coords()):
-            center = crd / bin_scale
-            sl = list[slice]()
-            outsl = list[slice]()
-            # We should deal with the borders of image.
-            for c, l, size in zip(center, [lz, ly, lx], lbl.shape, strict=True):
-                _left = int(c - l / 2 - 0.5)
-                _right = _left + l
-                _sl, _pad = utils.make_slice_and_pad(_left, _right, size)
-                sl.append(_sl)
-                outsl.append(
-                    slice(
-                        _pad[0] if _pad[0] > 0 else None,
-                        -_pad[1] if _pad[1] > 0 else None,
-                    )
-                )
-
-            sl = tuple(sl)
-            outsl = tuple(outsl)
-            lbl[sl][out[i][outsl]] = i + 1
-            yield
-
-        return lbl
-
-
 def get_code_theme(self: MagicTemplate) -> str:
     """Get the theme for CodeEdit using the napari theme."""
     from napari.utils.theme import get_theme

cylindra/widgets/main.py

@@ -2329,65 +2329,6 @@ def regionprops_features(
         dock.setFloating(True)
         return undo_callback(dock.close).with_redo(dock.show)
 
-    @set_design(text=capitalize, location=_sw.ImageMenu)
-    @thread_worker.with_progress(desc="Paint cylinders ...")
-    def paint_cylinders(
-        self,
-        color_by: Annotated[str, {"choices": [H.spacing, H.twist, H.rise, H.npf]}] = H.spacing,
-        cmap: _CmapType = DEFAULT_COLORMAP,
-        limits: Optional[tuple[float, float]] = (3.95, 4.28),
-    ):  # fmt: skip
-        """
-        Paint cylinder fragments by its local properties.
-
-        Parameters
-        ----------
-        {color_by}{cmap}{limits}
-        """
-        tomo = self.tomogram
-        all_df = tomo.splines.collect_localprops()
-        if color_by not in all_df.columns:
-            raise ValueError(f"Column {color_by} does not exist.")
-
-        paint_device = widget_utils.PaintDevice(
-            self._reserved_layers.image.data.shape,
-            self._reserved_layers.image.scale[-1],
-        )
-        lbl = yield from paint_device.paint_cylinders(self.tomogram, color_by)
-
-        # Labels layer properties
-        _id = "ID"
-        _str = "structure"
-        columns = [_id, H.rise, H.spacing, H.twist, _str]
-        df = (
-            all_df.select([H.spline_id, H.pos_id, H.rise, H.spacing, H.twist, H.npf, H.start])
-            .with_columns(
-                pl.format("{}-{}", pl.col(H.spline_id), pl.col(H.pos_id)).alias(_id),
-                pl.format("{}_{}", pl.col(H.npf), pl.col(H.start)).alias(_str),
-                pl.col(H.rise),
-                pl.col(H.spacing),
-                pl.col(H.twist),
-            )
-        )  # fmt: skip
-        back = pl.DataFrame([pl.Series(_id, [None], dtype=pl.Utf8)])
-        props = pl.concat([back, df[columns]], how="diagonal")
-        if limits is None:
-            limits = float(all_df[color_by].min()), float(all_df[color_by].max())
-
-        @thread_worker.callback
-        def _on_return():
-            # Add labels layer
-            self._reserved_layers.add_paint(lbl, props)
-            if self._reserved_layers.paint not in self.parent_viewer.layers:
-                self.parent_viewer.add_layer(self._reserved_layers.paint)
-            self._reserved_layers.paint.set_colormap(color_by, limits, cmap)
-            # TODO: undo paint
-            return undo_callback(
-                lambda: _Logger.print("undoing paint_cylinders do nothing.")
-            )
-
-        return _on_return
-
     # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
     #   Non-GUI methods
     # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
@@ -2708,9 +2649,6 @@ def _refer_spline_config(self, cfg: SplineConfig):
         fgui.npf.min, fgui.npf.max = cfg.npf_range.astuple()
         fgui.npf.value = int(cfg.npf_range.center)
 
-        fgui = get_function_gui(self.paint_cylinders)
-        fgui.limits.value = cfg.spacing_range.astuple()
-
         for method in [self.map_monomers, self.map_monomers_with_extensions, self.map_along_pf, self.map_centers]:  # fmt: skip
             get_function_gui(method)["orientation"].value = cfg.clockwise
 

cylindra/widgets/sta.py

@@ -127,7 +127,7 @@ def _get_landscape_layers(self: "SubtomogramAveraging", *_) -> list[LandscapeSur
 _SubVolumeSize = Annotated[
     Optional[nm],
     {
-        "text": "Use template shape",
+        "text": "use template or mask shape",
         "options": {"value": 12.0, "max": 100.0},
         "label": "size (nm)",
         "validator": _get_template_shape,
@@ -985,7 +985,8 @@ def align_all_template_free(
 
         Parameters
         ----------
-        {layers}{mask_params}{size}{max_shifts}{rotations}{cutoff}{interpolation}{method}{bin_size}
+        {layers}{mask_params}{size}{max_shifts}{rotations}{cutoff}{interpolation}
+        {method}{bin_size}
         """
         t0 = timer()
         layers = assert_list_of_layers(layers, self.parent_viewer)
@@ -994,7 +995,6 @@ def align_all_template_free(
         molecules = combiner.concat(layer.molecules for layer in layers)
         mask = self.params._get_mask(params=mask_params)
         if size is None:
-            shape = None
             raise NotImplementedError("'size' must be given.")
         else:
             shape = tuple(
@@ -1098,8 +1098,9 @@ def align_all_annealing(
 
         Parameters
         ----------
-        {layer}{template_path}{mask_params}{max_shifts}{rotations}{cutoff}{interpolation}
-        {range_long}{range_lat}{angle_max}{temperature_time_const}{upsample_factor}{random_seeds}
+        {layer}{template_path}{mask_params}{max_shifts}{rotations}{cutoff}
+        {interpolation}{range_long}{range_lat}{angle_max}{temperature_time_const}
+        {upsample_factor}{random_seeds}
         """
         t0 = timer()
         layer = assert_layer(layer, self.parent_viewer)
@@ -1247,7 +1248,8 @@ def run_annealing_on_landscape(
 
         Parameters
         ----------
-        {landscape_layer}{range_long}{range_lat}{angle_max}{temperature_time_const}{random_seeds}
+        {landscape_layer}{range_long}{range_lat}{angle_max}{temperature_time_const}
+        {random_seeds}
         """
         t0 = timer()
         landscape_layer = _assert_landscape_layer(landscape_layer, self.parent_viewer)
@@ -1405,8 +1407,8 @@ def calculate_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.
+            Precision of frequency to calculate FSC. "0.02" means that FSC will be
+            calculated at frequency 0.01, 0.03, 0.05, ..., 0.45.
         """
         t0 = timer()
         layers = assert_list_of_layers(layers, self.parent_viewer)
@@ -1459,7 +1461,7 @@ def classify_pca(
         layer: MoleculesLayerType,
         template_path: Annotated[_PathOrNone, {"bind": _template_param}] = None,
         mask_params: Annotated[Any, {"bind": _get_mask_params}] = None,
-        size: Annotated[Optional[nm], {"text": "Use mask shape", "options": {"value": 12.0, "max": 100.0}, "label": "size (nm)"}] = None,
+        size: _SubVolumeSize = None,
         cutoff: _CutoffFreq = 0.5,
         interpolation: Annotated[int, {"choices": INTERPOLATION_CHOICES}] = 3,
         bin_size: Annotated[int, {"choices": _get_available_binsize}] = 1,
@@ -1538,14 +1540,14 @@ def seam_search(
         anti_template_path: Annotated[Optional[Path.Read[FileFilter.IMAGE]], {"text": "Do not use anti-template"}] = 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",
+        show_average: Annotated[str, {"label": "show averages as", "choices": [None, "Raw", "Filtered"]}] = "Filtered",
         cutoff: _CutoffFreq = 0.25,
     ):  # fmt: skip
         """
         Search for the best seam position.
 
-        Try all patterns of seam positions and compare cross correlation values. If molecule
-        assembly has 13 protofilaments, this method will try 26 patterns.
+        Try all patterns of seam positions and compare cross correlation values. If
+        molecule assembly has 13 protofilaments, this method will try 26 patterns.
 
         Parameters
         ----------

cylindra/widgets/subwidgets/menus.py

@@ -230,7 +230,6 @@ def open_simulator(self):
 
     sep1 = field(Separator)
     sample_subtomograms = abstractapi()
-    paint_cylinders = abstractapi()
 
     @set_design(text=capitalize)
     @do_not_record
@@ -294,15 +293,17 @@ def show_splines(self):
         @set_design(text=capitalize)
         def show_splines_as_meshes(self):
             """Show 3D spline cylinder as a surface layer."""
+            # TODO: after napari supports features in surface layer, add spline
+            # properties
             main = self._get_main()
             nodes = []
             vertices = []
             n_nodes = 0
             for i, spl in enumerate(main.tomogram.splines):
-                n, v = spl.cylinder_model().to_mesh(spl)
-                nodes.append(n)
-                vertices.append(v + i * n_nodes)
-                n_nodes += n.shape[0]
+                node, vert = spl.cylinder_model().to_mesh(spl)
+                nodes.append(node)
+                vertices.append(vert + i * n_nodes)
+                n_nodes += node.shape[0]
             nodes = np.concatenate(nodes, axis=0)
             vertices = np.concatenate(vertices, axis=0)
             return main.parent_viewer.add_surface([nodes, vertices], shading="smooth")